Commit 0c8435e7 authored by Lukas Riedel's avatar Lukas Riedel Committed by Santiago Ospina De Los Ríos

Add option to insert water content as initial condition

* Implement matric head grid function adapter transforming water content
  into matric head.
* Add parameterization functions necessary to compute matric head from
  water content.
* Use shared_ptr when using instances of initial condition
  objects in models.
parent 2525bd1a
......@@ -38,6 +38,7 @@
* DG solver for solute transport model !112
* Cookbook tutorial on using the random field generator !184
* Outflow boundary condition for Richards model !191
* Specify water content as initial condition in Richards model !187
* Documentation about model solver loop and time step adaptation schemes !190
### Changed
......
......@@ -104,11 +104,11 @@ A full list is found at https://www.w3.org/TR/REC-html40/sgml/entities.html
<parameter name="quantity">
<definition> The physical quantity represented by the initial condition
data.
data or equation.
</definition>
<values> matricHead </values>
<values> matricHead, waterContent </values>
<suggestion> matricHead </suggestion>
<comment> Choose quantity represented: matricHead </comment>
<comment> Choose quantity represented: matricHead, waterContent </comment>
<versionadded version="unreleased"> </versionadded>
</parameter>
......
......@@ -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
^^^^^^^^^
......
......@@ -28,7 +28,7 @@ class InitialConditionFactory
protected:
/// Create a user defined initial condition
static auto create(
static std::unique_ptr<InitialCondition<T>> create(
const Dune::ParameterTree& ini_file,
const typename T::GV& grid_view,
std::shared_ptr<spdlog::logger> log=get_logger(log_base)
......
......@@ -374,6 +374,17 @@
- [Exception assertions](https://github.com/google/googletest/blob/master/googletest/docs/advanced.md#exception-assertions):
Check if a statement throws an exception and optionally check for its type.
### Notes for Unit-testing with DUNE structures
- The DUNE `GridView` objects and their associated structures like ID or index
sets are not necessarily thread-safe. As Google Test starts separate threads
for building fixtures and executing test cases, objects containing grid
views may lead to errors when built in and passed through fixtures. When in
doubt, create a factory-like function you call in the beginning of every
test case instead of using a fixture. You can also check the
[threading capabilities](https://www.dune-project.org/doxygen/2.6.0/structDune_1_1Capabilities_1_1viewThreadSafe.html)
of the respective grid.
@}
**/
#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 WaterContentToMatricHeadAdapter :
public Dune::PDELab::GridFunctionBase<
Dune::PDELab::GridFunctionTraits<typename T::GridView,
typename T::RF,
1,
typename T::Scalar>,
WaterContentToMatricHeadAdapter<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;
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
*/
WaterContentToMatricHeadAdapter (const GV& gv,
const std::shared_ptr<const GF> gf,
const std::shared_ptr<const P> param)
:
_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(sat);
}
const GV& getGridView () const
{
return _gv;
}
};
#endif // DUNE_DORIE_RICHARDS_ADAPTERS_HEAD_HH
......@@ -250,6 +250,50 @@ public:
};
}
/// Return a scaled version of the inverse Water Content function
/** Uses the function of the underlying parameterization and applies
* scaling as specified by ScalingType. Cast to operator-internal RF.
* \return Function: Water content -> Saturation
* \see water_content_f() for the inverse function
*/
std::function<RF(RF)> water_content_f_inv () const
{
verify_cache();
const auto& par =
std::get<std::shared_ptr<ParameterizationType>>(_cache);
// get water content function and apply the scaling
const auto& scale_por = std::get<Scaling>(_cache).scale_por;
const auto wc_f_inv = par->water_content_f_inv(scale_por);
using WaterContent = typename ParameterizationType::WaterContent;
return [wc_f_inv](const RF water_content) {
return wc_f_inv(WaterContent{water_content}).value;
};
}
/// Return a scaled version of the Saturation function
/** Uses the function of the underlying parameterization and applies
* scaling as specified by ScalingType. Cast to operator-internal RF.
* \return Function: Saturation -> Matric head
* \see saturation_f() for the inverse function
*/
std::function<RF(RF)> matric_head_f () const
{
verify_cache();
const auto& par =
std::get<std::shared_ptr<ParameterizationType>>(_cache);
const auto head_f = par->matric_head_f();
using Saturation = typename ParameterizationType::Saturation;
const auto& scale_head = std::get<Scaling>(_cache).scale_head;
return [head_f, &scale_head](const RF saturation){
return head_f(Saturation{saturation}).value / scale_head;
};
}
private:
/// Build the parameterization from an element mapping and the input file.
......
......@@ -7,22 +7,62 @@
#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 = WaterContentToMatricHeadAdapter<Traits, Param, Base>;
public:
/// Construct the transformation from a grid view and an initial condition
WaterContentInitialCondition (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>(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 +70,19 @@ 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::unique_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)
)
{
using IC = InitialCondition<Traits>;
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::unique_ptr<IC>();
}
const auto quantity = config.get<std::string>("initial.quantity");
......@@ -47,6 +90,14 @@ public:
if (quantity == "matricHead") {
return Base::create(config, grid_view, log);
}
else if (quantity == "waterContent") {
// create initial condition grid function from input data
std::shared_ptr<IC> ic_gf = Base::create(config, grid_view, log);
// plug into adapter
using WCIC = WaterContentInitialCondition<Traits, Param>;
return std::make_unique<WCIC>(grid_view, ic_gf, param);
}
else {
log->error("Unsupported quantity for initial condition: {}",
quantity);
......
......@@ -131,12 +131,41 @@ public:
};
}
/// 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
{
return [this, &scale_por](const WaterContent theta)
{
// Small value to avoid numerical issues
constexpr double eps_min = 1E-6;
const double sat = (theta.value - _theta_r.value)
/ (_theta_s.value + scale_por - _theta_r.value);
return Saturation{std::clamp(sat, eps_min, 1.0)};
};
}
/// Return a bound version of the saturation function
/** \return Function: Matric Head -> Water content
* \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
*/
......@@ -167,4 +196,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,6 +128,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
*/
......
......@@ -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;
......
......@@ -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 +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;
......
This diff is collapsed.
include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini
_initial_type = analytic, data | expand
__name = {_initial_type}
_asset_path = "${PROJECT_SOURCE_DIR}/test"
_quantity = matricHead, waterContent | expand quantity_type
{_quantity} == waterContent and {_initial_type} == analytic | exclude
__name = ic-{_initial_type}-{_quantity}
[grid]
dimensions = 2
......@@ -16,7 +17,7 @@ file = none
[richards.initial]
type = {_initial_type}
quantity = matricHead
quantity = {_quantity}
file = ${CMAKE_CURRENT_LIST_DIR}/scaling.h5
dataset = twos
......@@ -24,6 +25,9 @@ interpolation = nearest
equation = -1+2*y
[richards.parameters]
file = ${CMAKE_CURRENT_LIST_DIR}/test-param-richards-unscaled.yml
[transport.initial]
type = {_initial_type}
quantity = soluteConcentration
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment