...
  View open merge request
Commits (78)
......@@ -109,6 +109,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"
......
......@@ -39,6 +39,15 @@ set_target_properties(hdf5 PROPERTIES
# yaml-cpp
find_package (yaml-cpp 0.5.2 REQUIRED)
# 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)
......
......@@ -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 model 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>
......@@ -101,7 +102,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
^^^^^^^^^
......@@ -112,8 +128,12 @@ Initial condition tranformations for the Transport solver.
* ``quantity = soluteConcentration``
The input data is directly interpreted as solute concentration,
<<<<<<< HEAD
:math:`f = c_w [\text{kg}/\text{m}^3]`.
=======
:math:`f = c_w [\mathrm{kg}/\mathrm{m}^d]`, where :math:`d` indicates the
spatial dimensions of the grid.
>>>>>>> master
.. _H5: https://www.h5py.org/
.. _muparser: http://beltoforion.de/article.php?a=muparser&p=features
......@@ -646,7 +646,7 @@ namespace Dorie{
else if(buffer == "outflow"){
bcData[i][j].type.push_back(BoundaryCondition::Outflow);
instream >> buffer;
bcData[i][j].headValue.push_back(0.0);
bcData[i][j].headValue.push_back(stod(buffer));
bcData[i][j].fluxValue.push_back(stod(buffer));
}
else
......
......@@ -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)
......
......@@ -29,9 +29,10 @@ private:
const R dtmax, dtmin; //!< time step limits
const R eps; //!< error margin
const R dtinc, dtdec; //!< time step mutliplicators
const R tBegin; //!< model time limits
R tEnd; //!< model time limits
R tBegin; //!< simulation time limits
R tEnd; //!< simulation time limits
const int itmax, itmin; //!< newton iteration limits
const R dt_begin; //!< beginning time step
const TSC& tsc; //!< Class with Function getNextTimeStamp()
......@@ -60,6 +61,7 @@ public:
, tEnd(config.get<R>("time.end"))
, itmax(config.get<int>("time.maxIterations"))
, itmin(config.get<int>("time.minIterations"))
, dt_begin(dt)
, tsc(tsc_)
{
if (not input_valid()) {
......@@ -107,6 +109,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
......@@ -126,8 +134,17 @@ public:
// time has advanced after calculation
time += dt;
if(time >= tEnd-eps) // Model has ended. Do nothing
// Simulation has ended. Do nothing
if (time > tEnd) {
dt = dt_begin; // in case this was a data assimilation step
return true;
}
// NOTE: Simulation needs to know that we are actually at the end
else if (time <= tEnd && time >= tEnd-eps) {
time = tEnd;
dt = dt_begin; // in case this was a data assimilation step
return true;
}
dt = std::max(dt,dtmin);
dt = std::min(dt*dtinc,dtmax);
......@@ -139,18 +156,19 @@ public:
changeTime = tsc.getNextTimeStamp(time);
}
if(changeTime < 0.0 && time+dt > tEnd){
// time step adjustment to tEnd
dt = tEnd-time;
_log->trace("Adjusting time step to model end at: {:.2e}",
tEnd);
}
else if(changeTime > 0.0 && time+dt >= changeTime){
if(changeTime > 0.0 && time+dt >= changeTime){
// time step adjustment to BC Change
dt = changeTime-time-eps;
_log->trace("Adapting time step to change in boundary condition at "
"{:.2e}", changeTime);
}
if(time+dt > tEnd){
// time step adjustment to tEnd
dt = tEnd-time;
_log->trace("Adjusting time step to simulation end at: {:.2e}",
tEnd);
}
}
else // Newton solver threw error. Time is not advanced
......
#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.
......@@ -7,6 +7,7 @@
#include <dune/common/parametertree.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/float_cmp.hh>
#include <dune/geometry/referenceelements.hh>
......@@ -678,6 +679,11 @@ public:
bcType = BCType::Dirichlet;
}
else if (BoundaryCondition::isOutflow(bc))
{
bcType = BCType::Dirichlet;
}
else if (BoundaryCondition::isOther(bc))
bcType = BCType::none;
......@@ -741,6 +747,10 @@ public:
relCond = 2 * relCond_s * relCond_n / (relCond_s + relCond_n);
}
// No contribution if flux flows inward
if (BoundaryCondition::isOutflow(bc)
and Dune::FloatCmp::lt(numFlux, 0.0))
continue;
if constexpr (lopcase != LOPCase::RTVolume)
{
......
#ifndef DUNE_DORIE_RICHARDS_OPERATOR_FV_HH
#define DUNE_DORIE_RICHARDS_OPERATOR_FV_HH
#include <dune/common/float_cmp.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/pdelab/common/quadraturerules.hh>
......@@ -262,7 +264,7 @@ public:
const auto saturation_f_i = _param->saturation_f();
const auto conductivity_f_i = _param->conductivity_f();
if (BoundaryCondition::isDirichlet(bc))
if (BoundaryCondition::isDirichlet(bc) || BoundaryCondition::isOutflow(bc))
{
const auto geo_i = entity_i.geometry();
const auto center_position_i_g = geo_i.center();
......@@ -292,6 +294,11 @@ public:
// Water flux in normal direction w.r.t the intersection
const auto water_flux_n = - cond_i*dudn;
// No contribution for inward flux
if (BoundaryCondition::isOutflow(bc)
and Dune::FloatCmp::lt(double(water_flux_n), 0.0))
return;
// Contribution to residual from Dirichlet boundary
r_i.accumulate(lfsv_i, 0, water_flux_n*volume_f);
}
......
......@@ -5,6 +5,8 @@
#include <string>
#include <map>
#include <dune/common/float_cmp.hh>
namespace Dune {
namespace Dorie {
namespace Parameterization {
......@@ -118,47 +120,83 @@ 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
*/
virtual std::function<Conductivity(const Saturation)> conductivity_f ()
const = 0;
/// Return a multimap referecing all parameters by their names.
/// Return a map referecing all parameters by their names.
/** \return Map. Key: Name of parameter (string).
* Value(s): Value of parameter (double&)
*/
virtual std::multimap<std::string, double&> parameters () = 0;
virtual std::map<std::string, double&> parameters () = 0;
/// Return a map referecing all parameters by their names.
/** \return Map. Key: Name of parameter (string).
* Value: Value of parameter (const double&)
*/
virtual std::multimap<std::string, const double&> parameters () const = 0;
virtual std::map<std::string, const double&> parameters () const = 0;
/// Return a clone of this object
/** \return a unique pointer with a copy of this object.
/** \return Shared pointer with a copy of this object.
*/
virtual std::unique_ptr<Richards<Traits>> clone () const = 0;
virtual std::shared_ptr<Richards<Traits>> clone () const = 0;
};
......@@ -167,4 +205,4 @@ public:
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_PARAM_RICHARDS_INTERFACE_HH
\ No newline at end of file
#endif // DUNE_DORIE_PARAM_RICHARDS_INTERFACE_HH
......@@ -128,10 +128,26 @@ 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
*/
std::multimap<std::string, double&> parameters () override
std::map<std::string, double&> parameters () override
{
return {
{ResidualWaterContent::name, this->_theta_r.value},
......@@ -146,7 +162,7 @@ public:
/// Return a map of parameter names and values for manipulation
/** \return Map: Parameter name, parameter value in this object
*/
std::multimap<std::string, const double&> parameters () const override
std::map<std::string, const double&> parameters () const override
{
return {
{ResidualWaterContent::name, this->_theta_r.value},
......@@ -159,12 +175,12 @@ public:
}
/// Clone the plymorphic class
/** \return unique pointer to the cloned object
/** \return Shared pointer to the cloned object
*/
std::unique_ptr<Richards<Traits>> clone () const override
std::shared_ptr<Richards<Traits>> clone () const override
{
using ThisType = MualemVanGenuchten<Traits>;
return std::make_unique<ThisType>(*this);
return std::make_shared<ThisType>(*this);
}
};
......
......@@ -53,7 +53,7 @@ ModelRichards<Traits>::ModelRichards (
// --- 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);
// Read local operator settings
namespace OP = Dune::Dorie::Operator;
......@@ -139,7 +139,7 @@ void ModelRichards<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"));
......
......@@ -92,7 +92,9 @@ struct ModelRichardsTraits : 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 =
std::conditional_t<order==0,
......@@ -310,7 +312,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;
......@@ -351,7 +353,10 @@ public:
Dune::MPIHelper& _helper
);
~ModelRichards() = default;
/*------------------------------------------------------------------------*//**
* @brief Destroy this object
*/
virtual ~ModelRichards () 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)
......
......@@ -290,7 +290,7 @@ protected:
std::shared_ptr<SoluteParameters> sparam;
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;
......
......@@ -14,6 +14,17 @@ dorie_add_unit_test(SOURCES test-interpolators.cc NAME test-interpolators)
dorie_add_unit_test(SOURCES test-scaling-adapters.cc NAME test-scaling-adapters)
dorie_add_unit_test(SOURCES test-param-factory.cc NAME test-param-factory)
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(
UNIT_TEST
......@@ -34,12 +45,19 @@ 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
)
......@@ -51,7 +69,6 @@ dorie_add_metaini_test(
METAINI test-param-richards.mini.in
CREATED_TARGETS param_target
)
target_link_libraries(${param_target} spdlog)
dorie_add_metaini_test(
UNIT_TEST
......
#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-param-richards-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>
#include <dune/common/float_cmp.hh>
#include <dune/dorie/model/richards/initial_condition.hh>
#include <dune/dorie/model/transport/initial_condition.hh>
/// Traits stub used for this test
template<int d>
struct Traits
{
static constexpr int dim = d;
using Grid = Dune::YaspGrid<dim>;
using GridView = typename Grid::LeafGridView;
using GV = GridView;
using DF = typename Grid::ctype;
using Domain = Dune::FieldVector<DF, dim>;
using RF = double;
using Scalar = Dune::FieldVector<RF, 1>;
};
/// Test an initial condition against a lambda function on the grid
template<typename Grid, class IC, class Function>
void test_ic_against_function (
const std::shared_ptr<Grid> grid,
const std::unique_ptr<IC>& ic,
const Function func
)
{
using Range = typename IC::Traits::RangeType;
using namespace Dune::FloatCmp;
// iterate over the grid vertices
for (auto&& e : Dune::elements(grid->leafGridView()))
{
Range y;
const auto x = e.geometry().local(e.geometry().center());
ic->evaluate(e, x, y);
assert(eq(func(e, x), y));
}
}
#include "test-initial-condition.hh"
/// Test the basic initial conditions (without transformation)
/** \tparam T Traits
......@@ -97,53 +56,26 @@ int main (int argc, char** argv)
auto grid = gc.grid();
const auto index_map = gc.element_index_map();
{
auto config = Dune::Dorie::Setup::prep_ini_for_richards(inifile);
auto config = Dune::Dorie::Setup::prep_ini_for_transport(inifile);
// initial condition factory
using T = Traits<2>;
using ICF = Dune::Dorie::RichardsInitialConditionFactory<T>;
// initial condition factory
using T = Traits<2>;
using ICF = Dune::Dorie::TransportInitialConditionFactory<T>;
// 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 the transport logger
const auto log_level = config.get<std::string>("output.logLevel"); // FIXME
auto _log = Dune::Dorie::create_logger(Dune::Dorie::log_transport,
helper/*,spdlog::level::from_str(log_level)*/);
const auto quantity = config.get<std::string>("initial.quantity");
if (quantity == "matricHead") {
test_ic_base<T,ICF>(grid, config);
}
else {
log->error("Unsupported quantity for initial condition test: {}",
quantity);
DUNE_THROW(Dune::NotImplemented, "Unsupported initial condition "
"quantity");
}
const auto quantity = config.get<std::string>("initial.quantity");
if (quantity == "soluteConcentration") {
test_ic_base<T, ICF>(grid, config);
}
{
// auto config = Dune::Dorie::Setup::prep_ini_for_transport(inifile); // FIXME
auto config = inifile.sub("transport");
// initial condition factory
using T = Traits<2>;
using ICF = Dune::Dorie::TransportInitialConditionFactory<T>;
// create the transport logger
// const auto log_level = config.get<std::string>("output.logLevel"); // FIXME
auto _log = Dune::Dorie::create_logger(Dune::Dorie::log_transport,
helper/*,spdlog::level::from_str(log_level)*/);
const auto quantity = config.get<std::string>("initial.quantity");
if (quantity == "soluteConcentration") {
test_ic_base<T,ICF>(grid, config);
}
else {
log->error("Unsupported quantity for initial condition test: {}",
quantity);
DUNE_THROW(Dune::NotImplemented, "Unsupported initial condition "
"quantity");
}
else {
log->error("Unsupported quantity for initial condition test: {}",
quantity);
DUNE_THROW(Dune::NotImplemented, "Unsupported initial condition "
"quantity");
}
return 0;
......
include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini
_initial_type = analytic, data | expand
__name = {_initial_type}
__name = ic-transport-{_initial_type}
_asset_path = "${PROJECT_SOURCE_DIR}/test"
[grid]
......@@ -14,16 +14,6 @@ extensions = 1 1
[grid.mapping]
file = none
[richards.initial]
type = {_initial_type}
quantity = matricHead
file = ${CMAKE_CURRENT_LIST_DIR}/scaling.h5
dataset = twos
interpolation = nearest
equation = -1+2*y
[transport.initial]
type = {_initial_type}
quantity = soluteConcentration
......
#include <memory>
#include <cassert>
#include <dune/common/fvector.hh>
#include <dune/common/float_cmp.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/common/gridview.hh>
/// Traits stub used for this test
template<int d>
struct Traits
{
static constexpr int dim = d;
using Grid = Dune::YaspGrid<dim>;
using GridView = typename Grid::LeafGridView;
using GV = GridView;
using DF = typename Grid::ctype;
using Domain = Dune::FieldVector<DF, dim>;
using RF = double;
using Scalar = Dune::FieldVector<RF, 1>;
using Vector = Domain;
};
/// Test an initial condition against a lambda function on the grid
template<typename Grid, class IC, class Function>
void test_ic_against_function (
const std::shared_ptr<Grid> grid,
const std::shared_ptr<IC>& ic,
const Function func
)
{
using Range = typename IC::Traits::RangeType;
using namespace Dune::FloatCmp;
// iterate over the grid vertices
for (auto&& e : elements(grid->leafGridView()))
{
Range y;
const auto x = e.geometry().local(e.geometry().center());
ic->evaluate(e, x, y);
assert(eq(func(e, x), y));
}
}
spatial_resolution_north 0
spatial_resolution_south 0
spatial_resolution_west -1
spatial_resolution_east -1
number_BC_change_times 2
0 neumann 0 dirichlet 0
1000 neumann 0 dirichlet 0
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <random>
#include <iostream>
#include <iomanip>
#include <dune/common/float_cmp.hh>
#include <dune/common/exceptions.hh>
#include <dune/geometry/type.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/dorie/common/setup.hh>
#include <dune/dorie/common/util.hh>
#include <dune/dorie/common/utility.hh>
#include <dune/dorie/common/knofu/utility.hh>
#include <dune/dorie/model/richards/richards.hh>
#include <dune/dorie/model/richards/knofu/interface.hh>
#include <dune/dorie/model/richards/richards.cc>
#include <Eigen/Core>
using Grid = Dune::YaspGrid<2>;
using BaseTraits = Dune::Dorie::BaseTraits<
Grid, Dune::GeometryType::BasicType::cube>;
template<int order>
using Traits = Dune::Dorie::ModelRichardsTraits<BaseTraits, order>;
/// Assert that inserting and retrieving does not alter the vector
/** Draw a random state vector within given limits, store it into the model
* and retrieve it immediately. Compare both vectors. Only throw if deviations
* are larger than floating-point precision errors.
*
* \param sim The simulation object
* \param v_min Minimal value to insert
* \param v_max Maximal value to insert
* \param eps Margin from extreme values
*/
template<typename Sim>
void check_set_get_solution (Sim& sim,
const double v_min,
const double v_max,
const double eps=1E-6)
{
// retrieve the current state
using BackendVector = typename Sim::BackendVector;
auto in = Dune::Dorie::backend_to_eigen(sim.get_solution());
// insert random values (within limits)
std::uniform_real_distribution dist(v_min + eps, v_max - eps);
std::mt19937 ran;
std::vector<double> data(in.size());
std::generate(begin(data), end(data), std::bind(dist, std::ref(ran)));
in << Eigen::Map<Eigen::VectorXd>(data.data(), data.size());
// set/get and check
sim.set_solution(Dune::Dorie::eigen_to_backend<BackendVector>(in));
const auto out = Dune::Dorie::backend_to_eigen(sim.get_solution());
assert(in.size() == out.size());
if (in != out) {
const Eigen::ArrayXXd in_arr(in);
const Eigen::ArrayXXd out_arr(out);
std::cout << "Found " << (in_arr != out_arr).count()
<< " non-equal values" << std::endl;
std::cout << std::scientific
<< std::setprecision(std::numeric_limits<double>::digits10+1);
// compare again, neglecting floating-point precision discrepancies
size_t count = 0;
for (size_t i = 0; i < in_arr.size(); ++i) {
const auto val_in = in_arr(i);
const auto val_out = out_arr(i);
using namespace Dune::FloatCmp;
if (ne(val_in, val_out)) {
std::cout << "Dim: " << i
<< " In: " << val_in
<< " Out: " << val_out
<< std::endl;
count++;
}
}
std::cout << "FloatCmp returned " << count << " different values"
<< std::endl;
if (count > 0)
assert(false);
}
}
/// Verify the contents of the position vector
/** The following checks are performed:
* * All coordinates are within the extensions of the grid
* (assuming a rectangular grid)
*/
template<typename Sim, typename Grid>
void verify_position_vector (const Sim& sim,
const std::shared_ptr<Grid> grid)