...
 
Commits (56)
......@@ -107,6 +107,11 @@ build:system-tests: &build-tests
build:unit-tests:
<<: *build-tests
# for KnoFu util testing
# TODO: Remove before_script, better handling of Eigen3(?)
before_script:
- cd /opt/dune
- apt-get update && apt-get install -y libeigen3-dev && apt-get clean
script:
- CMAKE_FLAGS="$CMAKE_FLAGS
-DCMAKE_BUILD_TYPE=Debug"
......
cmake_minimum_required(VERSION 2.8.12)
project(dorie CXX)
cmake_minimum_required(VERSION 3.10)
project(dorie C CXX)
# set build type
if(NOT CMAKE_BUILD_TYPE)
......
......@@ -103,6 +103,9 @@ by CI tests.
#### Step-by-step Instructions
These instructions are suitable for a clean **Ubuntu** or **macOS** setup. The main difference between the two systems is the package manager. Debian-based systems have the APT manager already built in. On Mac, we recommend installing [Homebrew](http://brew.sh/). If you prefer to use [MacPorts](https://www.macports.org/), you need to check if some of the packages require different installation options than displayed here.
Manual installations on macOS require installing HDF5 from source. This can
be tricky, but the following instructions should work on a clean system.
If you installed [Anaconda](https://conda.io/docs/user-guide/install/download.html) on your machine, you don't need to install Python or Pip. Simply skip these packages when using the package managers for installing the software. However, notice the warnings when compiling DORiE below!
1. **macOS** users need to start by installing the Apple Command Line Tools by executing
......@@ -125,10 +128,33 @@ If you installed [Anaconda](https://conda.io/docs/user-guide/install/download.ht
**macOS:**
brew update
brew install cmake doxygen gcc libpng open-mpi muparser \
brew install cmake doxygen fftw gcc libpng open-mpi muparser \
pkg-config python3 superlu yaml-cpp
brew install hdf5 --with-mpi
brew install fftw --with-mpi
2. **macOS only:** Install HDF5 with MPI support from source.
1. Download an archive of the
[HDF5 source code](https://www.hdfgroup.org/downloads/hdf5/source-code/),
and extract it.
2. Enter the extracted folder. In there, create a `build` directory, and
enter it:
mkdir build && cd build
3. Configure your build. If you followed the instructions above, the
OpenMPI C compiler is reachable via the command `mpicc`. If not, you have
to specify a full path to it. Use the option `prefix` to specify where
you want the package to be installed. This should *not* be a
system-reserved path like `/usr/local`, and *not* be located in a
sub-directory of the source code. Execute the configuration script:
./../configure CC=mpicc --prefix=<hdf5-path> --enable-parallel
4. Build and install the library:
make && make install
3. The parallel linear solver of DORiE can make use of the ParMETIS package. If you want to run DORiE in parallel on multiple processes, additionally install METIS and ParMETIS:
......@@ -158,11 +184,18 @@ If you installed [Anaconda](https://conda.io/docs/user-guide/install/download.ht
If you installed software into paths not appended to your `PATH` variable, you will have to add `CMAKE_FLAGS` to the call to make sure that CMake finds all packages. Alternatively, you can add a custom options file. See the [DUNE Installation Instructions](https://dune-project.org/doc/installation/) for details. CMake will throw an error if required packages are not found.
**Warning:** Anacoda supplies its own version of HDF5 which is typically found first by CMake. If you have Anaconda installed on your machine, add
If you installed HDF5 from source (all macOS users) or use Anaconda,
specify the path to your HDF5 installation by using the `HDF5_ROOT`
variable. On Ubuntu, add the path to the APT package,
-DHDF5_ROOT=/usr/
and on macOS, add
HDF5_ROOT=/usr/local
-DHDF5_ROOT=<hdf5-path>
to the `CMAKE_FLAGS` in the call to `dunecontrol` above.
to the `CMAKE_FLAGS` in the above command, replacing `<hdf5-path>` with the
path chosen as installation prefix when configuring HDF5.
#### Experimental Features
The local operator implementing Richards equation's discretization supports
......
# find all required packages
FIND_PACKAGE (HDF5 REQUIRED)
# These macros check for the following packages, yielding the respective
# targets
#
# HDF5 ... hdf5
# yaml-cpp ... yaml-cpp
# muparser ... muparser::muparser
#
# All of these are required an will be autmatically linked to executables when
# using the DorieTesting CMake macros.
# Add CMake policy stack
cmake_policy(PUSH)
# Allow using `<Package>_ROOT` on CMake >3.12
if(POLICY CMP0074)
cmake_policy(SET CMP0074 NEW)
endif()
# HDF5
set(HDF5_PREFER_PARALLEL TRUE)
find_package (HDF5 REQUIRED COMPONENTS C)
if(NOT HDF5_IS_PARALLEL)
message(SEND_ERROR "Parallel HDF5 must be installed!")
endif()
add_definitions(-DHDF5_PARALLEL)
FIND_PACKAGE (FFTW REQUIRED)
FIND_PACKAGE (SuperLU REQUIRED)
FIND_PACKAGE (MPI REQUIRED)
# Define a target because the CMake module does not
add_library(hdf5 INTERFACE IMPORTED GLOBAL)
# TODO: Use `target_XXX` functions after raising CMake requirements
# Sanitize the definitions because we cannot use the proper CMake function...
string(REPLACE -D "" HDF5_DEFINITIONS "${HDF5_DEFINITIONS}")
# Set properties directly because of a bug in CMake 3.10
set_target_properties(hdf5 PROPERTIES
INTERFACE_LINK_LIBRARIES "${HDF5_LIBRARIES}"
INTERFACE_INCLUDE_DIRECTORIES "${HDF5_INCLUDE_DIRS}"
INTERFACE_COMPILE_DEFINITIONS "${HDF5_DEFINITIONS};H5_HAVE_PARALLEL")
# yaml-cpp
find_package (yaml-cpp 0.5.2 REQUIRED)
FIND_PACKAGE (muparser REQUIRED)
FIND_PACKAGE (METIS)
FIND_PACKAGE (ParMETIS)
# Eigen for KnoFu interface
find_package(Eigen3 3.3.4)
# add library target (header only-library)
if(EIGEN3_FOUND AND NOT TARGET Eigen3::Eigen)
add_library(Eigen3::Eigen INTERFACE IMPORTED)
set_target_properties(Eigen3::Eigen PROPERTIES
INTERFACE_INCLUDE_DIRECTORIES "${EIGEN3_INCLUDE_DIR}")
endif()
# muparser
find_package (muparser REQUIRED)
include_directories(${FFTW_INCLUDES}
${HDF5_INCLUDE_DIRS}
${YAML_CPP_INCLUDE_DIR})
# Report the DUNE libs
message (STATUS "DUNE Libraries: ${DUNE_LIBS}")
list (APPEND DUNE_LIBS
${FFTW_LIBRARIES}
${HDF5_LIBRARIES}
${YAML_CPP_LIBRARIES}
muparser::muparser)
MESSAGE(STATUS "DUNE Libraries: ${DUNE_LIBS}")
# Remove CMake policy stack
cmake_policy(POP)
# Add DORiE testing functions
include(DorieTesting)
......@@ -53,6 +53,9 @@ endfunction()
# tests within DORiE and adds flags for coverage reports. Notice that
# `dune_add_test` requires more parameters than this function alone.
#
# The registered executable is automatically linked to the libraries DORiE
# DORiE depends on.
#
# Use this function exactly like `dune_add_test`.
#
function(dorie_add_unit_test)
......@@ -75,6 +78,8 @@ function(dorie_add_unit_test)
endif()
# add to build target and employ compile options
target_link_libraries(${UNIT_TEST_TARGET}
muparser::muparser hdf5 yaml-cpp spdlog)
add_coverage_links(${UNIT_TEST_TARGET})
add_dependencies(build_unit_tests ${UNIT_TEST_TARGET})
endfunction()
......@@ -129,6 +134,9 @@ endfunction()
# Use this function like `dune_add_system_test`, considering the options
# given above.
#
# When registering a new executable (`TARGET` is not specified), this
# executable will be linked to the libraries DORiE depends on.
#
function(dorie_add_metaini_test)
set(OPTIONS UNIT_TEST)
set(SINGLE TARGET METAINI SCRIPT BASENAME CREATED_TARGETS)
......@@ -177,9 +185,16 @@ function(dorie_add_metaini_test)
SCRIPT ${SYSTEM_TEST_SCRIPT}
CREATED_TARGETS created_targets
)
# report created targets to parent scope
set(${SYSTEM_TEST_CREATED_TARGETS} ${created_targets} PARENT_SCOPE)
# Link to dependencies
if(NOT SYSTEM_TEST_TARGET)
target_link_libraries(${created_targets}
muparser::muparser hdf5 yaml-cpp spdlog)
endif()
# add dependencies and flags
if(SYSTEM_TEST_UNIT_TEST)
add_coverage_links(${created_targets})
......
# - Find the FFTW library
#
# Usage:
# find_package(FFTW [REQUIRED] [QUIET] )
#
# It sets the following variables:
# FFTW_FOUND ... true if fftw is found on the system
# FFTW_LIBRARIES ... full path to fftw library
# FFTW_INCLUDES ... fftw include directory
#
# The following variables will be checked by the function
# FFTW_USE_STATIC_LIBS ... if true, only static libraries are found
# FFTW_ROOT ... if set, the libraries are exclusively searched
# under this path
# FFTW_LIBRARY ... fftw library to use
# FFTW_INCLUDE_DIR ... fftw include directory
#
#If environment variable FFTWDIR is specified, it has same effect as FFTW_ROOT
if( NOT FFTW_ROOT AND ENV{FFTWDIR} )
set( FFTW_ROOT $ENV{FFTWDIR} )
endif()
# Check if we can use PkgConfig
find_package(PkgConfig)
#Determine from PKG
if( PKG_CONFIG_FOUND AND NOT FFTW_ROOT )
pkg_check_modules( PKG_FFTW QUIET "fftw3" )
endif()
#Check whether to search static or dynamic libs
set( CMAKE_FIND_LIBRARY_SUFFIXES_SAV ${CMAKE_FIND_LIBRARY_SUFFIXES} )
if( FFTW_ROOT )
#find libs
find_library(
FFTW_LIB
NAMES "fftw3"
PATHS ${FFTW_ROOT}
PATH_SUFFIXES "lib" "lib64"
NO_DEFAULT_PATH
)
find_library(
FFTWF_LIB
NAMES "fftw3_mpi"
PATHS ${FFTW_ROOT}
PATH_SUFFIXES "lib" "lib64"
NO_DEFAULT_PATH
)
find_library(
FFTWL_LIB
NAMES "fftw3l"
PATHS ${FFTW_ROOT}
PATH_SUFFIXES "lib" "lib64"
NO_DEFAULT_PATH
)
#find includes
find_path(
FFTW_INCLUDES
NAMES "fftw3.h"
PATHS ${FFTW_ROOT}
PATH_SUFFIXES "include"
NO_DEFAULT_PATH
)
else()
find_library(
FFTW_LIB
NAMES "fftw3"
PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR}
)
find_library(
FFTWF_LIB
NAMES "fftw3_mpi"
PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR}
)
find_library(
FFTWL_LIB
NAMES "fftw3l"
PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR}
)
find_path(
FFTW_INCLUDES
NAMES "fftw3.h"
PATHS ${PKG_FFTW_INCLUDE_DIRS} ${INCLUDE_INSTALL_DIR}
)
endif( FFTW_ROOT )
set(FFTW_LIBRARIES ${FFTWF_LIB} ${FFTW_LIB})
if(FFTWL_LIB)
set(FFTW_LIBRARIES ${FFTWL_LIB} ${FFTW_LIBRARIES})
endif()
set( CMAKE_FIND_LIBRARY_SUFFIXES ${CMAKE_FIND_LIBRARY_SUFFIXES_SAV} )
MESSAGE(STATUS "Found FFTW library: ${FFTW_LIBRARIES}")
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(FFTW DEFAULT_MSG
FFTW_INCLUDES FFTW_LIBRARIES)
mark_as_advanced(FFTW_INCLUDES FFTW_LIBRARIES)
# - Find the muParser library
# Find the muparser library
#
# Usage:
# find_package(muParser [REQUIRED] [QUIET] )
# find_package(muparser [REQUIRED] [QUIET] )
#
# Hints may be given by the (environment) variables
# MUPARSER_ROOT ... Path to the installation library
# muparser_ROOT ... Path to the installation library
#
# It sets the following variables:
# MUPARSER_FOUND ... true if muParser is found on the system
# MUPARSER_LIBRARIES ... full path to muParser library
# MUPARSER_INCLUDES ... muParser include directory
# muparser_FOUND ... true if muparser is found on the system
# muparser_LIBRARIES ... full path to muparser library
# muparser_INCLUDES ... muparser include directory
#
# It defines the following targets:
# muparser::muparser ... muparser library to link against
#
find_library(MUPARSER_LIBRARY
NAMES muparser muparserd
HINTS ${MUPARSER_ROOT} ENV MUPARSER_ROOT
find_library(muparser_LIBRARY
NAMES muparser muparserd
HINTS ${muparser_ROOT} ENV muparser_ROOT
)
find_path(MUPARSER_INCLUDE_DIR
muParserDef.h
HINTS ${MUPARSER_ROOT} ENV MUPARSER_ROOT
find_path(muparser_INCLUDE_DIR
muParserDef.h
HINTS ${muparser_ROOT} ENV muparser_ROOT
)
mark_as_advanced(MUPARSER_INCLUDE_DIR MUPARSER_LIBRARY)
mark_as_advanced(muparser_INCLUDE_DIR muparser_LIBRARY)
find_package_handle_standard_args(MUPARSER
DEFAULT_MSG
MUPARSER_LIBRARY MUPARSER_INCLUDE_DIR)
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(muparser
DEFAULT_MSG
muparser_LIBRARY muparser_INCLUDE_DIR)
if (MUPARSER_FOUND)
set(MUPARSER_LIBRARIES ${MUPARSER_LIBRARY} )
set(MUPARSER_INCLUDES ${MUPARSER_INCLUDE_DIR} )
if (muparser_FOUND)
set(muparser_LIBRARIES ${muparser_LIBRARY} )
set(muparser_INCLUDES ${muparser_INCLUDE_DIR} )
# add the target
add_library(muparser::muparser SHARED IMPORTED)
set_target_properties(muparser::muparser
PROPERTIES IMPORTED_LOCATION ${MUPARSER_LIBRARIES}
INTERACE_INCLUDE_DIRECTORIES ${MUPARSER_INCLUDES}
PROPERTIES IMPORTED_LOCATION ${muparser_LIBRARIES}
INTERACE_INCLUDE_DIRECTORIES ${muparser_INCLUDES}
)
endif()
......@@ -363,4 +363,24 @@ adding an empty line, make text **bold** or ``monospaced``.
</parameter>
</category>
<category name="KnoFu" hidden="true">
<parameter name="saturationMin">
<definition> Minimum saturation when transforming water content to
matric head.
</definition>
<values> float </values>
<suggestion> 1E-4 </suggestion>
<comment> Settings for the KnoFu data assimilation interface
</comment>
</parameter>
<parameter name="saturationMax">
<definition> Maximum saturation when transforming water content to
matric head.
</definition>
<values> float </values>
<suggestion> 1.0 </suggestion>
</parameter>
</category>
</dorie>
......@@ -31,16 +31,16 @@ They are controlled by the ``initial.type`` key and available for every model.
* ``type = analytic``
An analytic function :math:`f(\vec{p})` which depends on the physical
position :math:`\vec{p}`. The function must be defined via the key
An analytic function :math:`f(\mathbf{p})` which depends on the physical
position :math:`\mathbf{p}`. The function must be defined via the key
``initial.equation``. For parsing the input expression, we use muparser_
which supports a set of common mathematical functions. Additionally, the
following variables can be used:
Available variables:
* ``x``: X-coordinate :math:`p_1 \, [\mathrm{m}]`.
* ``y``: Y-coordinate :math:`p_2 \, [\mathrm{m}]`.
* ``z``: Z-coordinate :math:`p_3 \, [\mathrm{m}]` (only in 3D).
* ``x``: X-coordinate :math:`p_1 \, [\text{m}]`.
* ``y``: Y-coordinate :math:`p_2 \, [\text{m}]`.
* ``z``: Z-coordinate :math:`p_3 \, [\text{m}]` (only in 3D).
* ``h``: Height above origin. Synonymous to ``y`` in 2D and ``z`` in 3D.
* ``pi``: Mathematical constant :math:`\pi`.
* ``dim``: Number of physical dimensions.
......@@ -50,8 +50,9 @@ They are controlled by the ``initial.type`` key and available for every model.
:ref:`initial-transformation`), typical initial conditions for a
Richards simulation are
* Hydrostatic equilibrium: A vertical gradient of :math:`-1` and a
fixed value ``<v>`` at height :math:`h = 0 \, \mathrm{m}`::
* Hydrostatic equilibrium: A vertical gradient of
:math:`\partial_h h_m = -1` and a fixed value ``<v>`` at height
:math:`h = 0 \, \text{m}`:
initial.equation = -h + <v>
......@@ -59,7 +60,7 @@ They are controlled by the ``initial.type`` key and available for every model.
.. tip::
The expression for a gaussian pulse of solute concentration centered at
:math:`\vec{p} = [0.5, 0.5]^T \, \mathrm{m}` is::
:math:`\mathbf{p} = [0.5, 0.5]^T \, \text{m}` is::
initial.equation = exp(-sqrt((x-0.5)^2+(y-0.5)^2)/(4.*0.002))/(4*pi*0.002)^(2/dim).
......@@ -95,7 +96,22 @@ Initial condition tranformations for the Richards solver.
* ``quantity = matricHead``
The input data is directly interpreted as matric head
:math:`f = h_m [\mathrm{m}]`.
:math:`f = h_m \, [\text{m}]`.
.. object:: Water Content to Matric Head
* ``quantity = waterContent``
The input data is interpreted as water content,
:math:`f = \theta_w \, [\text{-}]`, and transformed into matric head via
the :doc:`parameterization <parameter-file>` of the medium.
Values greater than the porosity :math:`\phi` and less than the residual
water content :math:`\theta_r` are automatically clamped to fit the allowed
range. Additionally, any input value :math:`f(x_0)` at some position
:math:`x_0` on the grid will result in a saturation greater zero,
:math:`\Theta (x_0) > 0`, to avoid divergence of the matric head towards
negative infinity.
Transport
^^^^^^^^^
......@@ -106,7 +122,7 @@ Initial condition tranformations for the Transport solver.
* ``quantity = soluteConcentration``
The input data is directly interpreted as solute concentration,
:math:`f = c_w [\mathrm{kg}/\mathrm{m}^3]`.
:math:`f = c_w [\text{kg}/\text{m}^3]`.
.. _H5: https://www.h5py.org/
.. _muparser: http://beltoforion.de/article.php?a=muparser&p=features
add_executable("dorie-rfg" dorie-rfg.cc)
dune_target_link_libraries(dorie-rfg ${DUNE_LIBS})
\ No newline at end of file
add_executable(dorie-rfg dorie-rfg.cc)
dune_target_link_libraries(dorie-rfg ${DUNE_LIBS})
......@@ -56,11 +56,6 @@ int main(int argc, char** argv)
DUNE_THROW(Dune::IOError,"No parameter file specified!");
const std::string inifilename = argv[1];
#if !(HDF5_PARALLEL)
if (helper.size() > 1)
DUNE_THROW(Dune::Exception,"Parallel HDF5 library is needed to run in parallel");
#endif
// Read ini file
Dune::ParameterTree inifile;
Dune::ParameterTreeParser ptreeparser;
......
......@@ -4,16 +4,19 @@ if(dune-testtools_FOUND)
add_subdirectory("test")
endif()
add_executable("richards" richards.cc)
dune_target_link_libraries(richards dorie-richards ${DUNE_LIBS} spdlog)
add_executable(richards richards.cc)
dune_target_link_libraries(richards ${DUNE_LIBS})
target_link_libraries(richards
dorie-richards spdlog muparser::muparser hdf5 yaml-cpp)
add_executable("transport" transport.cc)
dune_target_link_libraries(transport dorie-richards ${DUNE_LIBS} spdlog)
dune_target_link_libraries(transport dorie-transport ${DUNE_LIBS} spdlog)
add_executable(transport transport.cc)
dune_target_link_libraries(transport ${DUNE_LIBS})
target_link_libraries(transport
dorie-richards dorie-transport spdlog muparser::muparser hdf5 yaml-cpp)
add_custom_target("dorie" DEPENDS richards transport)
# enable setting operator scheme from config file
if(EXPERIMENTAL_DG_FEATURES)
target_compile_definitions("dorie" PUBLIC -DEXPERIMENTAL_DG_FEATURES)
endif()
\ No newline at end of file
endif()
......@@ -34,7 +34,7 @@ protected:
std::shared_ptr<spdlog::logger> log=get_logger(log_base)
)
{
std::unique_ptr<InitialCondition<T>> ic;
std::shared_ptr<InitialCondition<T>> ic;
auto ic_type = ini_file.get<std::string>("initial.type");
log->debug("Creating initial condition of type: {}", ic_type);
......@@ -56,7 +56,7 @@ protected:
const auto file_type = ic_datafile.substr(ext_start + 1);
if (file_type == "h5") {
using ICH5 = InitialConditionH5<T>;
ic = std::make_unique<ICH5>(grid_view, ini_file, log);
ic = std::make_shared<ICH5>(grid_view, ini_file, log);
} else if (file_type == "vtu") {
DUNE_THROW(NotImplemented,
......@@ -73,7 +73,7 @@ protected:
auto equation_str = ini_file.get<std::string>("initial.equation");
using ICAnalytic = InitialConditionAnalytic<T>;
ic = std::make_unique<ICAnalytic>(log,grid_view,equation_str);
ic = std::make_shared<ICAnalytic>(log,grid_view,equation_str);
}
else {
log->error("Unsupported initial condition type: {}", ic_type);
......
This diff is collapsed.
......@@ -39,9 +39,9 @@ namespace Setup
* \param greetings Initial message of the logger (log-level `info`,
* will *always* be displayed)
*/
auto init (int& argc,
char**& argv,
const std::string greetings="Initializing program")
inline auto init (int& argc,
char**& argv,
const std::string greetings="Initializing program")
-> std::tuple<Dune::ParameterTree,
std::shared_ptr<spdlog::logger>,
Dune::MPIHelper&>
......@@ -92,7 +92,7 @@ namespace Setup
*
* The respective process ID is displayed in the log messages.
*/
void debug_hook ()
inline void debug_hook ()
{
// cache log level
const auto log = get_logger(log_base);
......@@ -122,7 +122,7 @@ namespace Setup
*
* @param ini The inifile
*/
Dune::ParameterTree prep_ini_for_richards(const Dune::ParameterTree& ini)
inline Dune::ParameterTree prep_ini_for_richards(const Dune::ParameterTree& ini)
{
const auto log = get_logger(log_base);
Dune::ParameterTree richards_ini(ini.sub("richards"));
......@@ -147,7 +147,7 @@ namespace Setup
*
* @param ini The inifile
*/
Dune::ParameterTree prep_ini_for_transport(const Dune::ParameterTree& ini)
inline Dune::ParameterTree prep_ini_for_transport(const Dune::ParameterTree& ini)
{
Dune::ParameterTree transport_ini(ini.sub("transport"));
// move grid extensions into transport category (needed for reading parameters)
......
......@@ -48,6 +48,7 @@ public:
{ }
virtual ~SimulationBase () = default;
/*-----------------------------------------------------------------------*//**
* @brief Sets the output policy.
*
......
......@@ -29,7 +29,7 @@ private:
const R dtmax, dtmin; //!< time step limits
const R eps; //!< error margin
const R dtinc, dtdec; //!< time step mutliplicators
const R tBegin; //!< simulation time limits
R tBegin; //!< simulation time limits
R tEnd; //!< simulation time limits
const int itmax, itmin; //!< newton iteration limits
......@@ -107,6 +107,12 @@ public:
tEnd = _t_end;
}
/// Set a new begin time
void set_t_begin (const R _t_begin)
{
tBegin = _t_begin;
}
/// Validation of solution and time step calculation
/** If the Newton solver throws an exception (because the solution did not converge
* in the allowed steps or the timestep is too large), exc=true is passed to this
......
#ifndef DUNE_DORIE_RICHARDS_ADAPTERS_HEAD_HH
#define DUNE_DORIE_RICHARDS_ADAPTERS_HEAD_HH
#include <memory>
#include <algorithm>
#include <dune/pdelab/common/function.hh>
/// An adapter for translating water content to matric head
/** \tparam T BaseTraits
* \tparam GF Grid function representing water content
* \tparam P Type of the parameter interface
*/
template<typename T, typename P, typename GF>
class MatricHeadAdapter :
public Dune::PDELab::GridFunctionBase<
Dune::PDELab::GridFunctionTraits<typename T::GridView,
typename T::RF,
1,
typename T::Scalar>,
MatricHeadAdapter<T, GF, P>
>
{
public:
using Traits = Dune::PDELab::GridFunctionTraits<typename T::GridView,
typename T::RF,
1,
typename T::Scalar>;
private:
using GV = typename Traits::GridViewType;
using RF = typename Traits::RangeType;
/// Lower boundary to clamp the saturation to
const RF _sat_min;
/// Upper boundary to clamp the saturation to
const RF _sat_max;
const GV& _gv;
const std::shared_ptr<const GF> _gf;
const std::shared_ptr<const P> _param;
public:
/// Constructor
/*
* \param config Dune config file tree
* \param gv GridView
* \param gf Grid function to interpolate from
* \param param Parameter object
*/
MatricHeadAdapter (const Dune::ParameterTree& config,
const GV& gv,
const std::shared_ptr<const GF> gf,
const std::shared_ptr<const P> param)
:
_sat_min(config.get<RF>("KnoFu.saturationMin")),
_sat_max(config.get<RF>("KnoFu.saturationMax")),
_gv(gv),
_gf(gf),
_param(param)
{ }
/**
* \brief Evaluation of grid function at certain position
* \param e Element entity
* \param x Position in local entity coordinates
* \param y Return value
* \return Value of grid function
*/
void evaluate (const typename Traits::ElementType& e,
const typename Traits::DomainType& x,
typename Traits::RangeType& y) const
{
// evaluate water content
typename Traits::RangeType water_content;
_gf->evaluate(e, x, water_content);
// evaluate matric head
_param->bind(e);
const auto saturation_f = _param->water_content_f_inv();
const auto matric_head_f = _param->matric_head_f();
const auto sat = saturation_f(water_content);
y = matric_head_f(std::clamp<double>(sat, _sat_min, _sat_max));
}
const GV& getGridView () const
{
return _gv;
}
};
#endif // DUNE_DORIE_RICHARDS_ADAPTERS_HEAD_HH
......@@ -7,28 +7,70 @@
#include <dune/common/exceptions.hh>
#include <dune/dorie/common/initial_condition/factory.hh>
#include <dune/dorie/model/richards/adapters/matric_head.hh>
namespace Dune {
namespace Dorie {
/// Transforms water content into matric head, serving as initial condition
template<typename Traits, typename Param>
class WaterContentInitialCondition :
public InitialConditionTransformation<Traits>
{
public:
using ICTransf = InitialConditionTransformation<Traits>;
using GFTraits = typename ICTransf::Traits;
using Base = InitialCondition<Traits>;
protected:
using GV = typename ICTransf::GV;
using MHA = MatricHeadAdapter<Traits, Param, Base>;
public:
/// Construct the transformation from a grid view and an initial condition
WaterContentInitialCondition (const Dune::ParameterTree& config,
const GV& grid_view,
std::shared_ptr<Base> ic,
const std::shared_ptr<const Param> param ):
ICTransf(grid_view, ic),
_adapter(std::make_shared<MHA>(config, grid_view, ic, param))
{ }
void evaluate ( const typename GFTraits::ElementType& e,
const typename GFTraits::DomainType& x,
typename GFTraits::RangeType& y) const override
{
_adapter->evaluate(e, x, y);
};
private:
std::shared_ptr<MHA> _adapter;
};
/// An initial condition factory for the Richards solver
/**
* \tparam T Traits
* \tparam Traits
* \tparam Param
*
* \ingroup InitialConditions
* \author Lukas Riedel
* \date 2019
*/
template<typename T>
class RichardsInitialConditionFactory : public InitialConditionFactory<T>
template<typename Traits, typename Param>
class RichardsInitialConditionFactory : public InitialConditionFactory<Traits>
{
private:
using Base = InitialConditionFactory<T>;
using Base = InitialConditionFactory<Traits>;
public:
/// Create an initial condition for the Richards solver
static std::unique_ptr<InitialCondition<T>> create (
static std::shared_ptr<InitialCondition<Traits>> create (
const Dune::ParameterTree& config,
const typename T::GridView& grid_view,
const typename Traits::GridView& grid_view,
const std::shared_ptr<const Param> param,
const std::shared_ptr<spdlog::logger> log=get_logger(log_richards)
)
{
......@@ -37,6 +79,14 @@ public:
if (quantity == "matricHead") {
return Base::create(config, grid_view, log);
}
else if (quantity == "waterContent") {
// create initial condition grid function from input data
auto ic_gf = Base::create(config, grid_view, log);
// plug into adapter
using WCIC = WaterContentInitialCondition<Traits, Param>;
return std::make_shared<WCIC>(config, grid_view, ic_gf, param);
}
else {
log->error("Unsupported quantity for initial condition: {}",
quantity);
......
This diff is collapsed.
#ifndef DUNE_DORIE_RICHARDS_KNOFU_PLOTTING_HH
#define DUNE_DORIE_RICHARDS_KNOFU_PLOTTING_HH
#include <string>
#include <memory>
#include <vector>
#include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/common/vtkexport.hh>
namespace Dune {
namespace Dorie {
/// Simulation used for plotting a state returned by KnoFu
template<typename Traits>
class KnoFuPlotter : public KnoFuInterface<Traits>
{
private:
/// Type of the base class
using Base = KnoFuInterface<Traits>;
/// Type of the DOF vector
using U = typename Base::U;
/// Type of the discrete grid function based on the solution vector
using GFSolution = typename Base::GFSolution;
/// Type of the VTK function adapter based on the solution grid function
using VTKFunction = Dune::PDELab::VTKGridFunctionAdapter<GFSolution>;
/// VTK writer
Dune::SubsamplingVTKWriter<typename Traits::GV> vtkwriter;
/// State storage
std::vector<std::shared_ptr<U>> member_states;
/// Name storage associated with respective state
std::vector<std::string> member_names;
/// Output file encoding
const Dune::VTK::OutputType output_type;
public:
using GridCreator = typename Traits::GridCreator;
using BackendVector = typename Base::BackendVector;
/// Constructor. Set up Simulation base class and VTK writer
KnoFuPlotter (
Dune::ParameterTree& inifile,
const GridCreator& grid_creator,
Dune::MPIHelper& helper,
const FilterSolutionType solution_type)
:
KnoFuInterface<Traits>(inifile, grid_creator, helper, solution_type),
vtkwriter(this->gv, Dune::refinementLevels(
inifile.get<int>("output.subsamplingLevel"))),
output_type(Dune::VTK::OutputType::base64)
{ }
/// Default destructor
virtual ~KnoFuPlotter () override = default;
/// Add a solution to the VTK output.
/** Create a new backend vector based on GFS and the solution vector
* supplied.
* \param state State vector (water content)
* \param name Name of the state (i.e. member number)
*/
void add_solution_to_output (
const BackendVector& state,
const std::string& name)
{
auto solution = std::make_shared<U>(*this->gfs, 0.0);
Dune::PDELab::Backend::native(*solution) = state;
member_states.push_back(solution);
member_names.push_back(name);
}
/// Write all stored states into a .vtu file. Then clear all data
/** Create discrete grid functions and VTK functions from the saved
* states.
* \param path Path to write into
* \param name Filename to write
*/
void write_output (const std::string& path, const std::string& name)
{
// stack VTKFunctions into writer
for(size_t i = 0; i < member_states.size(); ++i) {
auto dgf = std::make_shared<GFSolution>(*this->gfs,
*member_states.at(i));
auto dgf_vtk = std::make_shared<VTKFunction>(dgf,
member_names.at(i));
if (this->inifile.template get<bool>("output.vertexData"))
vtkwriter.addVertexData(dgf_vtk);
else
vtkwriter.addCellData(dgf_vtk);
}
// write output
vtkwriter.pwrite(name, path, "", output_type);
// release data
vtkwriter.clear();
member_states.clear();
member_names.clear();
}
};
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_RICHARDS_KNOFU_PLOTTING_HH
......@@ -82,7 +82,7 @@ inline auto read_experimental_operator_settings(
RichardsDGUpwinding::Type upwinding;
RichardsDGWeights::Type weights;
const auto method_str = inifile.get<std::string>("dg.experimental.method");
const auto method_str = inifile.get<std::string>("numerics.experimental.method");
log->debug(" DG method: {}", method_str);
if (method_str == "SIPG")
method = RichardsDGMethod::SIPG;
......@@ -98,7 +98,7 @@ inline auto read_experimental_operator_settings(
}
const auto upwinding_str
= inifile.get<std::string>("dg.experimental.upwinding");
= inifile.get<std::string>("numerics.experimental.upwinding");
log->debug(" Upwinding scheme: {}", upwinding_str);
if (upwinding_str == "none")
upwinding = RichardsDGUpwinding::none;
......@@ -111,7 +111,7 @@ inline auto read_experimental_operator_settings(
DUNE_THROW(Dune::IOError, "Unknown upwinding scheme");
}
const auto weights_str = inifile.get<bool>("dg.experimental.weights");
const auto weights_str = inifile.get<bool>("numerics.experimental.weights");
log->debug(" Harmonic weights at interfaces: {}", weights_str);
if (weights_str)
weights = RichardsDGWeights::weightsOn;
......
......@@ -5,6 +5,8 @@
#include <string>
#include <map>
#include <dune/common/float_cmp.hh>
namespace Dune {
namespace Dorie {
......@@ -117,25 +119,61 @@ public:
/// Return a bound version of the water content function
/** \param scale_por Porosity scaling factor
* \return Function: Saturation -> Water content
* \see water_content_f_inv() for the inverse function
*/
std::function<WaterContent(const Saturation)> water_content_f (
const RF& scale_por
) const
{
return [this, &scale_por](const Saturation sat) {
const auto porosity = _theta_s.value + scale_por;
return [this, porosity](const Saturation sat) {
return WaterContent{
sat.value * (_theta_s.value - _theta_r.value + scale_por)
+ _theta_r.value
sat.value * (porosity - _theta_r.value) + _theta_r.value
};
};
}
/// Return a bound version of the inverse water content function
/** \param scale_por Porosity scaling factor
* \return Function: Water content -> Saturation
* \see water_content_f() for the inverse function
*/
std::function<Saturation(const WaterContent)> water_content_f_inv (
const RF& scale_por
) const
{
const auto porosity = _theta_s.value + scale_por;
return [this, porosity](const WaterContent theta) {
// sanitize result
using namespace FloatCmp;
if (le(theta.value, _theta_r.value)) {
return Saturation{1E-6};
}
else if (ge(theta.value, porosity)) {
return Saturation{1.0};
}
// actual computation
return Saturation{
(theta.value - _theta_r.value) / (porosity - _theta_r.value)
};
};
}
/// Return a bound version of the saturation function
/** \return Function: Matric Head -> Water content
/** \return Function: Matric Head -> Saturation
* \see matric_head_f for the inverse function
*/
virtual std::function<Saturation(const MatricHead)> saturation_f () const
= 0;
/// Return a bound version of the inverse saturation function
/** \return Function: Saturation -> Matric Head
* \see saturation_f for the inverse function
*/
virtual std::function<MatricHead(const Saturation)> matric_head_f () const
= 0;
/// Return a bound version of the conductivity function
/** \return Function: Saturation -> Hydraulic conductivity
*/
......
......@@ -127,6 +127,22 @@ public:
};
}
/// Return a scaled inverse saturation function
/** \return Function: Water content -> Matric Head
* \see saturation_f for the inverse function
*/
std::function<MatricHead(const Saturation)> matric_head_f () const
override
{
const auto m = 1.0 - 1.0 / _n.value;
return [this, m](const Saturation sat) {
return MatricHead {
std::pow(_alpha.value, -1)
* std::pow(std::pow(sat.value, -1.0/m) - 1.0, 1.0/_n.value)
};
};
}
/// Return a map of parameter names and values for manipulation
/** \return Map: Parameter name, parameter value in this object
*/
......
......@@ -53,7 +53,7 @@ RichardsSimulation<Traits>::RichardsSimulation (
// --- Operator Helper Classes ---
fboundary = std::make_shared<const FlowBoundary>(inifile);
fsource = std::make_shared<const FlowSource>(inifile);
finitial = FlowInitialFactory::create(inifile, gv, this->_log);
finitial = FlowInitialFactory::create(inifile, gv, fparam, this->_log);
// --- Local Operators ---
#ifdef EXPERIMENTAL_DG_FEATURES
......@@ -129,7 +129,7 @@ void RichardsSimulation<Traits>::operator_setup()
// --- Solvers ---
lsgfs = std::make_unique<LSGFS>(LSGFSHelper::create(gv));
lscc = std::make_unique<LSCC>();
ls = std::make_unique<LS>(*igo,*cc,*lsgfs,*lscc,1000,0,true,true);
ls = std::make_unique<LS>(*igo, *cc, *lsgfs, *lscc, 1000, 0, true, true);
pdesolver = std::make_unique<PDESOLVER>(*igo,*ls);
pdesolver->setParameters(inifile.sub("NewtonParameters"));
......
......@@ -86,7 +86,9 @@ struct RichardsSimulationTraits : public BaseTraits
/// Class defining the initial condition
using FlowInitial = Dune::Dorie::InitialCondition<BaseTraits>;
/// Class defining the initial condition factory
using FlowInitialFactory = Dune::Dorie::RichardsInitialConditionFactory<BaseTraits>;
using FlowInitialFactory
= Dune::Dorie::RichardsInitialConditionFactory<
BaseTraits, FlowParameters>;
/// Local spatial operator
using SLOP = Dune::Dorie::Operator::RichardsDGSpatialOperator<BaseTraits,FlowParameters,FlowBoundary,FlowSource,typename GFSHelper::FEM,false>;
/// Local temporal operator
......@@ -296,7 +298,7 @@ protected:
std::shared_ptr<FlowParameters> fparam;
std::shared_ptr<const FlowBoundary> fboundary;
std::shared_ptr<const FlowSource> fsource;
std::unique_ptr<FlowInitial> finitial;
std::shared_ptr<FlowInitial> finitial;
std::unique_ptr<CalculationController> controller;
std::unique_ptr<SLOP> slop;
......@@ -337,7 +339,10 @@ public:
Dune::MPIHelper& _helper
);
~RichardsSimulation() = default;
/*------------------------------------------------------------------------*//**
* @brief Destroy this object
*/
virtual ~RichardsSimulation () override = default;
/*------------------------------------------------------------------------*//**
* @brief Compute a time step. Models the time step requirement of
......
......@@ -26,7 +26,7 @@ private:
public:
/// Create an initial condition for the Richards solver
static std::unique_ptr<InitialCondition<T>> create (
static std::shared_ptr<InitialCondition<T>> create (
const Dune::ParameterTree& config,
const typename T::GridView& grid_view,
const std::shared_ptr<spdlog::logger> log=get_logger(log_transport)
......
......@@ -274,7 +274,7 @@ protected:
std::unique_ptr<CC> cc;
std::shared_ptr<SoluteBoundary> sboundary;
std::unique_ptr<SoluteInitial> sinitial;
std::shared_ptr<SoluteInitial> sinitial;
std::unique_ptr<CalculationController> controller;
std::unique_ptr<SLOP> slop;
......
......@@ -7,13 +7,20 @@ file(COPY scaling.h5
# scaling adapter test
dorie_add_unit_test(SOURCES test-simulation-base.cc NAME test-simulation-base)
target_link_libraries(ut-test-simulation-base spdlog)
dorie_add_unit_test(SOURCES test-grid-function-container.cc NAME test-grid-function-container)
target_link_libraries(ut-test-grid-function-container spdlog)
dorie_add_unit_test(SOURCES test-interpolators.cc NAME test-interpolators)
target_link_libraries(ut-test-interpolators spdlog)
dorie_add_unit_test(SOURCES test-scaling-adapters.cc NAME test-scaling-adapters)
target_link_libraries(ut-test-scaling-adapters spdlog)
if (EIGEN3_FOUND)
dorie_add_metaini_test(
UNIT_TEST
SOURCE test-knofu-interface-richards.cc
BASENAME test-knofu-interface-richards
METAINI test-knofu-interface-richards.mini.in
CREATED_TARGETS kir_target
)
dune_target_link_libraries(${kir_target} spdlog Eigen3::Eigen)
endif()
# grid creator test
dorie_add_metaini_test(
......@@ -23,7 +30,6 @@ dorie_add_metaini_test(
METAINI test-grid-creation.mini.in
CREATED_TARGETS gc_target
)
target_link_libraries(${gc_target} spdlog)
configure_file(test-grid-creation-global.mini.in test-grid-creation-global.mini)
dune_add_system_test(
......@@ -36,15 +42,21 @@ add_custom_target(test_grid_creator
)
add_dependencies(test_grid_creator prepare_testing)
# initial condition test
# initial condition tests
dorie_add_metaini_test(
UNIT_TEST
SOURCE test-initial-condition-richards.cc
BASENAME test-initial-condition-richards
METAINI test-initial-condition-richards.mini.in
CREATED_TARGETS ic_target
)
dorie_add_metaini_test(
UNIT_TEST
SOURCE test-initial-condition.cc
BASENAME test-initial-condition
METAINI test-initial-condition.mini.in
SOURCE test-initial-condition-transport.cc
BASENAME test-initial-condition-transport
METAINI test-initial-condition-transport.mini.in
CREATED_TARGETS ic_target
)
target_link_libraries(${ic_target} spdlog)
# parameterization and scaling test
dorie_add_metaini_test(
......@@ -54,7 +66,6 @@ dorie_add_metaini_test(
METAINI test-parameterization.mini.in
CREATED_TARGETS param_target
)
target_link_libraries(${param_target} spdlog)
dorie_add_metaini_test(
UNIT_TEST
......@@ -63,7 +74,6 @@ dorie_add_metaini_test(
METAINI test-parameterization-scaled.mini.in
CREATED_TARGETS param_target
)
target_link_libraries(${param_target} spdlog)
add_custom_target(test-parameterization
COMMAND ctest --output-on-failure --tests-regex ^.+test-parameterization.+$
......@@ -86,4 +96,4 @@ dune_add_test(
MPI_RANKS 2
TIMEOUT 30
CMD_ARGS ${CMAKE_CURRENT_BINARY_DIR}/gc_files_0003.ini
)
\ No newline at end of file
)
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/dorie/common/setup.hh>
#include <dune/dorie/common/grid_creator.hh>
#include <dune/dorie/model/richards/flow_parameters.hh>
#include <dune/dorie/model/richards/initial_condition.hh>
#include "test-initial-condition.hh"
/// Test the basic initial conditions (without transformation)
/** \tparam T Traits
* \tparam IFC InitialConditionFactory
*/
template<typename T, typename ICF, typename Param>
void test_ic_base (const std::shared_ptr<typename T::Grid> grid,
const Dune::ParameterTree config,
const std::shared_ptr<Param> param)
{
// create the initial condition
const auto ic = ICF::create(config, grid->leafGridView(), param);
const auto type = config.get<std::string>("initial.type");
if (type == "analytic") {
// define testing function
const auto equation = config.get<std::string>("initial.equation");
assert(equation == "-1+2*y");
auto test_f = [](const auto e, const auto x){
using Range = typename T::Scalar;
constexpr int dim = T::Grid::dimension;
const auto x_global = e.geometry().global(x);
Range y = -1.0 + x_global[dim-1] * 2.0;
return y;
};
test_ic_against_function(grid, ic, test_f);
}
else if (type == "file") {
const auto dataset = config.get<std::string>("initial.dataset");
// define testing function
auto test_f = [&dataset](const auto e, const auto x){
typename T::Scalar y;
if (dataset == "twos") // case: matric head input
y = 2.0;
else if (dataset == "ones") // case: water content input
y = 0.0;
return y;
};
test_ic_against_function(grid, ic, test_f);
}
}
int main (int argc, char** argv)
{
try{
// Initialize ALL the things!
auto [inifile, log, helper] = Dune::Dorie::Setup::init(argc, argv);
// create the grid
Dune::Dorie::GridCreator<Traits<2>::Grid> gc(inifile, helper);
auto grid = gc.grid();
const auto index_map = gc.element_index_map();
// fetch config for Richards
auto config = Dune::Dorie::Setup::prep_ini_for_richards(inifile);
// create the Richards logger
const auto log_level = config.get<std::string>("output.logLevel");
auto _log = Dune::Dorie::create_logger(
Dune::Dorie::log_richards,
helper,
spdlog::level::from_str(log_level));
// create Richards parameterization
using Param = Dune::Dorie::FlowParameters<Traits<2>>;
auto param = std::make_shared<Param>(config, grid, index_map);
// initial condition factory
using T = Traits<2>;
using ICF = Dune::Dorie::RichardsInitialConditionFactory<T, Param>;
const auto quantity = config.get<std::string>("initial.quantity");
if (quantity == "matricHead" || quantity == "waterContent") {
test_ic_base<T, ICF>(grid, config, param);
}
else {
log->error("Unsupported quantity for initial condition test: {}",
quantity);
DUNE_THROW(Dune::NotImplemented, "Unsupported initial condition "
"quantity");
}
return 0;
}
catch (Dune::Exception &e) {
std::cerr << "Dune reported error: " << e << std::endl;
return 1;
}
catch (std::exception &e) {
std::cerr << "Exception thrown: " << e.what() << std::endl;
return 1;
}
catch (...) {
std::cerr << "Unknown exception!" << std::endl;
return 1;
}
}
include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini
_initial_type = analytic, data | expand
_quantity = matricHead, waterContent | expand quantity_type
{_quantity} == waterContent and {_initial_type} == analytic | exclude
__name = ic-richards-{_initial_type}-{_quantity}
_asset_path = "${PROJECT_SOURCE_DIR}/test"
[grid]
dimensions = 2
initialLevel = 1
gridType = rectangular
cells = 50 50
extensions = 1 1
[grid.mapping]
file = none
[richards.initial]
type = {_initial_type}
quantity = {_quantity}
file = ${CMAKE_CURRENT_LIST_DIR}/scaling.h5
dataset = twos
interpolation = nearest
equation = -1+2*y
[richards.parameters]
file = test-parameterization-unscaled.yml
......@@ -2,53 +2,12 @@
#include "config.h"
#endif
#include <dune/common/fvector.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/common/gridview.hh>
#include <dune/dorie/common/setup.hh>
#include <dune/dorie/common/grid_creator.hh>