...
  View open merge request
Commits (79)
......@@ -140,6 +140,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 -DCOVERAGE_REPORT=On
-DCMAKE_BUILD_TYPE=None"
......
......@@ -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)
......
......@@ -371,4 +371,20 @@ A full list is found at https://www.w3.org/TR/REC-html40/sgml/entities.html
</parameter>
</category>
<category name="KnoFu">
<parameter name="saturationMin">
<definition> Minimum saturation value when transferring water content to
matric head </definition>
<values> float </values>
<suggestion> 1E-4 </suggestion>
</parameter>
<parameter name="saturationMax">
<definition> Maximum saturation value when transferring water content to
matric head </definition>
<values> float </values>
<suggestion> 1.0 </suggestion>
</parameter>
</category>
</dorie>
......@@ -154,7 +154,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
^^^^^^^^^
......@@ -165,8 +180,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)
......
#ifndef DUNE_DORIE_CALCULATION_CONTROLLER_HH
#define DUNE_DORIE_CALCULATION_CONTROLLER_HH
#include <iostream>
#include <cmath>
#include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>
#include <dune/dorie/common/logging.hh>
namespace Dune{
namespace Dorie{
/// Class to control time and maximum number of iterations for the Newton solver
template<typename R, typename TSC>
class CalculationController
{
private:
//! Typedef of logger for brevity
using Logger = std::shared_ptr<spdlog::logger>;
R time; //!< current time
R dt; //!< current time step
R suggested_dt; //!< suggestion for the next dt
bool use_dt_suggestion; //!< flag for suggested dt
int it; //!< current number of allowed iterations for Newton solver
//! logger of this instance
const Logger _log;
const R dtmax, dtmin; //!< time step limits
const R eps; //!< error margin
const R dtinc, dtdec; //!< time step mutliplicators
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()
public:
/// Read relevant time values, check for consistency and adapt time step to BC change, if necessary
/** \param tsc_ Class with public function getNextTimeStamp(time)
* \param log The logger to use for this instance,
* defaults to retrieving the base logger.
*/
CalculationController(const ParameterTree& config,
const TSC& tsc_,
const Logger log=get_logger(log_base))
: time(config.get<R>("time.start"))
, dt(config.get<R>("time.startTimestep"))
, suggested_dt(dt)
, use_dt_suggestion(false)
, _log(log)
, dtmax(config.get<R>("time.maxTimestep"))
, dtmin(config.get<R>("time.minTimestep"))
, eps(dtmin/10.0)
, dtinc(config.get<R>("time.timestepIncreaseFactor"))
, dtdec(config.get<R>("time.timestepDecreaseFactor"))
, tBegin(config.get<R>("time.start"))
, 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()) {
DUNE_THROW(IOError, "Invalid input for time controller");
}
// Check whether first time step needs adjustment
const R changeTime = tsc.getNextTimeStamp(time);
if(changeTime>0.0 && time+dt >= changeTime){
dt = changeTime-time-eps;
_log->trace("Adapting time step to change in boundary condition at "
"{:.2e}", changeTime);
}
it=calcIterations(dt);
}
/// suggest timestep
inline void suggest_timestep (R suggestion) {
if (FloatCmp::lt(suggestion,dt)) {
_log->trace("Adapting time step due to a suggestion: "
"{:.2e} -> {:.2e}", dt, suggestion);
dt = suggestion;
}
}
/// Return current time
inline R getTime () const { return time; }
/// Return timestep for next calculation
inline R getDT () const { return dt; }
/// Return begin time
inline R getBeginTime () const { return tBegin; }
/// Return end time
inline R getEndTime () const { return tEnd; }
/// Return maximum number of Newton iterations allowed for next calculation
inline int getIterations () const { return it; }
void set_time (const R _time)
{
time = _time;
}
/// Set a new finish time for the computation
void set_t_end (const R _t_end)
{
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
* function and the timestep will be reduced. Returning false causes the Problem
* solver to recalculate the timestep with new parameters provided by the controller.
* The next time step is also adjusted to the next time the BC changes
* \param exc Boolean if Newton Class threw an exception
* \return Boolean if solution is valid
* \throw Dune::Exception Newton solver does not converge at dt>=dtmin and it<=itmax
*/
bool validate (const bool exc)
{
R changeTime = 0.0;
if (!exc) // solution was successfully computed
{
// time has advanced after calculation
time += dt;
// 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);
changeTime = tsc.getNextTimeStamp(time);
if(time >= changeTime-eps && time <= changeTime){
// small step to ensure that next time step starts with new BC
time = changeTime;
changeTime = tsc.getNextTimeStamp(time);
}
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
{
if (dt <= dtmin) {
_log->error("Solver threw an exception at minimal time step "
"{:.2e}. Time step cannot be reduced", dt);
DUNE_THROW(Exception,"Solver fails at minimal time step");
}
auto failed_dt = dt;
dt = std::max(dt*dtdec,dtmin);
_log->trace("Adapting time step to improve solver: "
"{:.2e} -> {:.2e}", failed_dt, dt);
}
// Calculate new iterations based on adapted time step
it=calcIterations(dt);
return !exc;
}
private:
/// Calculate the allowed number of Newton iterations by logarithmic interpolation between max and min timestep
/** \param dt Timestep size
* \return Number of allowed Newton steps for inserted dt
*/
inline int calcIterations (const R dt)
{
return std::round( (itmin-itmax) * std::log10(dt-dtmin+1.0) / std::log10(dtmax-dtmin+1.0) + itmax );
}
/// Some checks if user input makes sense
/**
* \note This method does not throw exceptions by itself!
* \return True if the input variables are consistent
*/
bool input_valid () const
{
if (dtinc<=1.0) {
_log->error("timestepIncreaseFactor must be >1.0. Value is: {}",
dtinc);
return false;
}
if (dtdec>=1.0 || dtdec<=0.0) {
_log->error("timestepDecreaseFactor exceeds (0.0, 1.0). "
"Value is: {}", dtdec);
return false;
}
if (itmax<itmin) {
_log->error("maxIterations must be >=minIterations. "
"Values are: {}, {}", itmin, itmax);
return false;
}
if (dtmax<dtmin) {
_log->error("maxTimestep must be >=minTimestep. "
"Values are: {}, {}", dtmin, dtmax);
return false;
}
if (dt<dtmin || dt>dtmax) {
_log->error("startTimestep exceeds [minTimestep, maxTimestep]. "
"Value is: {}", dt);
return false;
}
return true;
}
};
} // namespace Dorie
} // namespace Dune
#endif
#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,22 +7,63 @@
#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
......@@ -30,16 +71,17 @@ public:
* \note In case the user wants the initial condition to be computed from
* the stationary problem, an empty pointer is returned.
*/
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)
)
{
if (config.get<std::string>("initial.type") == "stationary") {
// NOTE: Return empty pointer to communicate that the IC is
// computed internally.
return std::unique_ptr<InitialCondition<T>>();
return std::shared_ptr<InitialCondition<Traits>>();
}
const auto quantity = config.get<std::string>("initial.quantity");
......@@ -47,6 +89,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>
......
......@@ -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);
}
};
......
......@@ -58,7 +58,7 @@ ModelRichards<Traits>::ModelRichards (
const auto& boundary_map = _grid_creator.boundary_index_map();
fboundary = std::make_shared<FlowBoundary>(inifile, boundary_map, this->_log);
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;
......@@ -156,7 +156,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"));
......
......@@ -100,7 +100,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,
......@@ -325,7 +327,7 @@ protected:
std::shared_ptr<FlowParameters> fparam;
std::shared_ptr<FlowBoundary> fboundary;
std::shared_ptr<const FlowSource> fsource;
std::unique_ptr<FlowInitial> finitial;
std::shared_ptr<FlowInitial> finitial;
std::unique_ptr<SLOP> slop;
MBE mbe_slop;
......@@ -364,7 +366,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)
......
......@@ -364,7 +364,7 @@ protected:
std::shared_ptr<SoluteParameters> sparam;
std::shared_ptr<SoluteBoundary> sboundary;
std::unique_ptr<SoluteInitial> sinitial;
std::shared_ptr<SoluteInitial> sinitial;
std::shared_ptr<LSGFS> lsgfs;
std::shared_ptr<LSCC> lscc;
......
......@@ -53,6 +53,17 @@ dorie_add_metaini_test(
)
target_link_libraries(${tsc_target} 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(
UNIT_TEST
......@@ -73,12 +84,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
)
......@@ -90,7 +108,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);