Commit 9b244ced authored by Lukas Riedel's avatar Lukas Riedel

Introduce new parameterization structures

* FlowParameters will serve as interface to the operators
    * It also manages storage and access of parameters
* RichardsParameterization is the base class for all parameterizations
* MualemVanGenuchtenParameterization is a derived class defining all
    necessary parameterization functions and stores the parameter values
parent f176bfc3
#ifndef DUNE_DORIE_PARAM_INTERFACE_HH
#define DUNE_DORIE_PARAM_INTERFACE_HH
#include <memory>
#include <dune/common/exceptions.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/dorie/solver/param_richards.hh>
#include <dune/dorie/solver/param_mvg.hh>
namespace Dune {
namespace Dorie {
template <typename Traits>
class FlowParameters
{
private:
using RF = typename Traits::RF;
using Grid = typename Traits::Grid;
using LevelGridView = typename Grid::LevelGridView;
using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<
LevelGridView
>;
using Index = typename Mapper::Index;
using ParameterStorage = std::map<
Index, std::shared_ptr<Dorie::RichardsParameterization<Traits>>
>;
/// Need this gruesome typedef because 'map::value_type' has const key
using Cache = std::pair<typename ParameterStorage::key_type,
typename ParameterStorage::mapped_type
>;
/// Configuration file tree
const Dune::ParameterTree& _config;
/// Grid view of the coarsest grid configuration (level 0)
const LevelGridView _gv;
/// Mapper for mapping grid elements to indices
Mapper _mapper;
/// Map for storing parameterization information for each grid entity
ParameterStorage _param;
/// Currently cached parameterization (bound to element)
mutable Cache _cache;
private:
/// Check if an entity has been cached
void verify_cache () const
{
if (not _cache.second){
DUNE_THROW(Dune::Exception, "No parameterization cached."
<< "Call 'bind' before accessing the parameterization.");
}
}
public:
FlowParameters (
const Dune::ParameterTree& config,
std::shared_ptr<Grid> grid
):
_config(config),
_gv(grid->levelGridView(0)),
_mapper(_gv, Dune::mcmgElementLayout())
{ }
/// Bind the parameterization of an Entity, caching its values
/** \param entity Entity to bind to
*/
template <class Entity>
void bind (const Entity& entity) const
{
const auto index = _mapper.index(entity);
const auto it = _param.find(index);
if (it == _param.end()) {
DUNE_THROW(Dune::Exception, "Could not retrieve parameterization"
<< "for entity " << index);
}
_cache = *it;
}
/// Return the conductivity function. Cast to RangeField
std::function<RF(RF)> conductivity_f () const
{
verify_cache();
auto& par = _cache.second;
return [&par](const RF saturation){
return par->conductivity_f()(Saturation<RF>{saturation}).value;
};
}
/// Return the saturaton function. Cast to RangeField
std::function<RF(RF)> saturation_f () const
{
verify_cache();
auto& par = _cache.second;
return [&par](const RF matric_head){
return par->saturation_f()(MatricHead<RF>{matric_head}).value;
};
}
template <class Param>
void copy_from_old_parameters (const Param& p_base)
{
// iterate over grid elements
for (auto&& entity: elements(_gv))
{
const auto pos = entity.geometry().center();
auto get_param_value = [&p_base, &pos](auto param) {
return p_base.interpolate(*param, pos);
};
std::vector<RF> p_values;
std::transform(p_base.parameters.begin(), p_base.parameters.end(),
std::back_inserter(p_values),
get_param_value
);
// read parameter values
using R = RichardsParameterization<Traits>;
using MvG = MualemVanGenuchtenParameterization<Traits>;
typename R::ResidualWaterContent theta_r{p_values.at(0)};
typename R::SaturatedWaterContent theta_s{p_values.at(1)};
typename MvG::Alpha alpha{p_values.at(2)};
typename MvG::N n{p_values.at(3)};
typename MvG::Tau tau{p_values.at(4)};
typename MvG::SaturatedConductivity k0{p_values.at(5)};
auto mvg =
std::make_shared<MualemVanGenuchtenParameterization<Traits>>
(
theta_r,
theta_s,
alpha,
tau,
n,
k0
);
// insert parameters into map
const auto index = _mapper.index(entity);
auto res = _param.emplace(index, mvg);
if (not res.second) {
DUNE_THROW(Dune::Exception, "Parameterization for entity"
<< index << " already exists."
);
}
}
}
};
} // namespace Dune
} // namespace Dorie
#endif // DUNE_DORIE_PARAM_INTERFACE_HH
\ No newline at end of file
#ifndef DUNE_DORIE_PARAM_MVG_HH
#define DUNE_DORIE_PARAM_MVG_HH
#include <cmath>
#include <dune/dorie/solver/param_richards.hh>
namespace Dune {
namespace Dorie {
/// Class specifying the Mualem–van Genuchten parameterization
template <typename Traits>
class MualemVanGenuchtenParameterization :
public RichardsParameterization<Traits>
{
private:
using RF = typename Traits::RF;
using Base = RichardsParameterization<Traits>;
public:
using Saturation = Saturation<RF>;
using Conductivity = Conductivity<RF>;
using MatricHead = MatricHead<RF>;
using ResidualWaterContent = typename Base::ResidualWaterContent;
using SaturatedWaterContent = typename Base::SaturatedWaterContent;
/// Parameter defining the saturated hydraulic conductivity
struct SaturatedConductivity
{
RF value;
inline static const std::string name = "k_0";
};
/// Parameter defining the value of alpha (air-entry value)
struct Alpha
{
RF value;
inline static const std::string name = "a";
};
/// Parameter defining the value of n (retention curve shape)
struct N
{
RF value;
inline static const std::string name = "n";
};
/// Parameter defining the value of tau (skew factor)
struct Tau
{
RF value;
inline static const std::string name = "tau";
};
private:
std::function<Conductivity(const Saturation,
const SaturatedConductivity,
const N,
const Tau)>
_conductivity_f = [](const Saturation sat,
const SaturatedConductivity k0,
const N n,
const Tau tau) {
const auto m = 1.0 - 1.0 / n.value;
return Conductivity {
k0.value
* std::pow(sat.value, tau.value)
* std::pow(1 - std::pow(1 - std::pow(sat.value, 1.0/m), m), 2)
};
};
std::function<Saturation(const MatricHead, const Alpha, const N)>
_saturation_f = [](const MatricHead head,
const Alpha alpha,
const N n) {
if (head.value >= 0.0) {
return Saturation{1.0};
}
const auto m = 1.0 - 1.0 / n.value;
return Saturation {
std::pow(1 + std::pow(alpha.value * head.value, n.value), -m)
};
};
Alpha _alpha;
Tau _tau;
N _n;
SaturatedConductivity _k0;
public:
MualemVanGenuchtenParameterization (
const ResidualWaterContent theta_r,
const SaturatedWaterContent theta_s,
const Alpha alpha,
const Tau tau,
const N n,
const SaturatedConductivity k0
):
RichardsParameterization<Traits>(theta_r, theta_s),
_alpha(alpha),
_tau(tau),
_n(n),
_k0(k0)
{ }
std::function<Conductivity(const Saturation)> conductivity_f () const
override
{
return std::bind(_conductivity_f,
std::placeholders::_1,
std::cref(_k0), std::cref(_n), std::cref(_tau)
);
}
std::function<Saturation(const MatricHead)> saturation_f () const
override
{
return std::bind(_saturation_f,
std::placeholders::_1,
std::cref(_alpha), std::cref(_n)
);
}
};
} // namespace Dune
} // namespace Dorie
#endif // DUNE_DORIE_PARAM_MVG_HH
\ No newline at end of file
#ifndef DUNE_DORIE_PARAM_RICHARDS_HH
#define DUNE_DORIE_PARAM_RICHARDS_HH
#include <functional>
namespace Dune {
namespace Dorie {
// --- Strong types for parameterization --- //
// These are standard types for all parameterizations of Richards Equation
/// Type of water saturation
template <typename T>
struct Saturation
{
T value;
};
/// Type of volumetric water content
template <typename T>
struct WaterContent
{
T value;
};
/// Type of hydraulic conductivity
template <typename T>
struct Conductivity
{
T value;
};
/// Type of matric head
template <typename T>
struct MatricHead
{
T value;
};
/// Class defining the interface for parameterizations of Richards equation
template <typename Traits>
class RichardsParameterization
{
private:
using RF = typename Traits::RF;
public:
using Saturation = Saturation<RF>;
using WaterContent = WaterContent<RF>;
using Conductivity = Conductivity<RF>;
using MatricHead = MatricHead<RF>;
/// Parameter defining the residual water content
struct ResidualWaterContent
{
RF value;
inline static const std::string name = "theta_r";
};
/// Parameter defining the saturated water content
struct SaturatedWaterContent
{
RF value;
inline static const std::string name = "theta_s";
};
private:
ResidualWaterContent _theta_r;
SaturatedWaterContent _theta_s;
/// Functional relation for Saturation depending on water content
/** \param WaterContent Volumetric water content
* \return Saturation
*/
std::function<
Saturation(const WaterContent, const ResidualWaterContent, const SaturatedWaterContent)>
_saturation_f = [](const WaterContent _wc,
const ResidualWaterContent _theta_r,
const SaturatedWaterContent _theta_s) {
return Saturation{
(_wc.value - _theta_r.value) / (_theta_s.value - _theta_r.value)
};
};
/// Functional relation for Saturation depending on water content
/** \param WaterContent Volumetric water content
* \return Saturation
*/
std::function<
WaterContent(const Saturation, const ResidualWaterContent, const SaturatedWaterContent)>
_water_content_f = [](const Saturation _sat,
const ResidualWaterContent _theta_r,
const SaturatedWaterContent _theta_s) {
return WaterContent{
_sat.value * (_theta_s.value - _theta_r.value) + _theta_r.value
};
};
public:
RichardsParameterization (
const ResidualWaterContent theta_r,
const SaturatedWaterContent theta_s
):
_theta_r(theta_r),
_theta_s(theta_s)
{ }
virtual ~RichardsParameterization () = default;
/// Return a bound version of the water content function
/** \return Function: Saturation -> Water content
*/
std::function<WaterContent(const Saturation)> water_content_f() const
{
return std::bind(_water_content_f,
std::placeholders::_1,
std::cref(_theta_r), std::cref(_theta_s)
);
}
/// Return a bound version of the saturation function
/** \return Function: Matric Head -> Water content
*/
virtual std::function<Saturation(const MatricHead)> saturation_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;
};
} // namespace Dune
} // namespace Dorie
#endif // DUNE_DORIE_PARAM_RICHARDS_HH
\ No newline at end of file
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