Commit 940254dd authored by Lukas Riedel's avatar Lukas Riedel 📝

Restructure parameterization source files and remove deprecated ones

* Delete old parameterization, factory, and parameterization base class
* Delete FlowParameters method for copying parameters from the old
    parameterization instances
* Rename 'param_interface.hh' to 'richards_param.hh', to comply to
    other source files directly related to the local operator.
* Rename 'param_richards.hh' to 'param_interface.hh'.
parent e3743fea
......@@ -17,7 +17,7 @@
#include "util.hh"
#include "adaptivity.hh"
#include "../solver/util_interpolator.hh"
#include "../solver/param_interface.hh"
#include "../solver/richards_param.hh"
#include "../solver/richards_initial.hh"
#include "../solver/richards_boundary.hh"
#include "../solver/richards_source.hh"
......
#ifndef DUNE_DORIE_CONDUCTIVITY_ADAPTER_HH
#define DUNE_DORIE_CONDUCTIVITY_ADAPTER_HH
#include "../param_richards.hh"
namespace Dune{
namespace Dorie{
......
// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_DORIE_PARAMETERS_BASE_HH
#define DUNE_DORIE_PARAMETERS_BASE_HH
#include <string>
#include <cmath>
#include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parallel/mpihelper.hh>
#include "util_h5tools.hh"
namespace Dune {
namespace Dorie {
/// Simple class providing a data structure for storing parameters
template<typename T>
class Parameter
{
private:
typedef typename T::RangeField RF;
typedef typename T::Array Array;
const RF min, max;
Array values;
const bool check_bounds;
bool initialized;
public:
const std::string name; //!< Name of the parameter
Parameter(const std::string name_, const RF min_, const RF max_)
: min(min_), max(max_), check_bounds(true), initialized(false), name(name_)
{}
Parameter(const std::string name_)
: min(0), max(0), check_bounds(false), initialized(false), name(name_)
{}
void set_values(const Array values_)
{
values = values_;
if (check_bounds) {
for (std::size_t i = 0; i < values.size(); i++){
const auto& v = values[i];
if ((v < min) || (v > max))
{
DUNE_THROW(Dune::IOError, "Invalid parameter value encontered for"
<< " parameter " << name << " at (" << i << "): " << v
<< " (must be between "<< min << " and "<< max <<")");
}
}
}
initialized = true;
}
const Array& get_values() const
{
if(!initialized)
DUNE_THROW(Exception,"Parameter must be initialized before accessing it");
return values;
}
};
/// Base class for all parameterizations
template<typename T, class Interpolator>
class ParameterizationBase
{
private:
enum {dim = T::dim};
typedef typename T::RangeField RF;
typedef typename T::Vector Vector;
typedef typename T::Domain Domain;
typedef typename T::Element Element;
typedef typename T::Array Array;
typedef typename T::Index Index;
typedef std::vector<Parameter<T>*> ParameterArray;
protected:
const ParameterTree& config; //!< ini file parser
const std::string parName; //!< Identifier of the parameterization (main group name in H5 file)
public:
ParameterArray& parameters; //!< Array containing all defined parameters
protected:
const Dune::MPIHelper& helper; //!< MPI helper
Interpolator& interpolator; //!< Interpolation handler instance
H5File& h5file; //!< HDF5 interface
const int verbose; //!< verbosity level
const Domain pf_extensions; //!< Extensions of the parameter field
public:
const bool millerSimilarity; //!< Switches miller similarity on or off
const RF variance; //!< Random field variance, needed for Miller Similarity
protected:
const Domain grid_extensions; //!< Spatial extensions of the model grid
const Domain scale; //!< Scaling factor of the random field
const Domain offset; //!< Position of the origin of the random field
public:
Parameter<T> millerField; //!< Random field used for Miller similarity
protected:
Vector gVector; //!< Unit length vector pointing towards gravitational force
Domain fieldCells; //!< Random field extensions in cells
bool array_open;
public:
explicit ParameterizationBase(const ParameterTree& config_, const std::string parName_,
ParameterArray& parameters_, const Dune::MPIHelper& helper_,
Interpolator& interpolator_, H5File& h5file_,
const int verbose_)
: config(config_), parName(parName_), parameters(parameters_), helper(helper_)
, interpolator(interpolator_), h5file(h5file_), verbose(verbose_)
, pf_extensions(h5file.read_attribute<Domain>("extensions","parameters"))
, millerSimilarity(h5file.read_attribute<bool>("millerSimilarity","parameters"))
, variance(h5file.read_attribute<RF>("variance","parameters"))
, grid_extensions(config.get<Domain>("grid.extensions"))
, scale(config.get<Domain>("parameters.scale"))
, offset(config.get<Domain>("parameters.offset"))
, millerField("raw_field")
{
gVector=0;
gVector[dim-1]=-1.;
for(Index i=0; i<dim; i++)
{
if(pf_extensions[i]*scale[i] - offset[i] < grid_extensions[i])
DUNE_THROW(IOError,"Random field extensions are smaller than model grid extensions");
}
}
virtual ~ParameterizationBase () = default;
/// Iterate through parameter array and read each one from the H5 file
void read_parameters()
{
open_array("parameters/" + parName);
read_parameter(millerField);
for (auto &p: parameters)
{
Domain fieldCells_old = fieldCells;
read_parameter(*p);
if (fieldCells_old != fieldCells) // consistency check
DUNE_THROW(Exception,"Varying array sizes encountered in random field array");
}
close_array();
Domain scaled_extensions;
for(typename Domain::size_type i=0; i<dim; i++)
scaled_extensions[i] = pf_extensions[i] * scale[i];
interpolator.initialize(fieldCells, scaled_extensions);
}
// Purely virtual functions queried by the DG operator
virtual RF saturation(const RF& h, const Domain& pos) const =0;
virtual RF waterContent(const RF& h, const Domain& pos) const =0;
virtual RF condFactor(const RF& h, const Domain& pos) const =0;
virtual RF K(const Domain& pos) const =0;
virtual RF matricHead(const RF waterContent, const Domain& pos) const = 0;
/******************************************
* MILLER SIMILARITY QUERIES
* According to Miller and Roth
******************************************/
/// Transformation: Head: Reference to Miller similar (if activated)
/** \f[ h^* = h \exp\left( - (v - \text{Var}(v))\right) \f]
* with the raw random field value \f$v\f$ and its variance.
* \param h Matric head
* \param pos Position in global coordinates
* \return Transformed value
*/
RF headRefToMiller (const RF& h, const Domain& pos) const
{
if (!millerSimilarity)
return h;
return h * exp(-(interpolate(millerField, pos)-variance));
}
/// Transformation: Head: Miller similar to Reference (if activated)
/** \see headRefToMiller
*/
RF headMillerToRef (const RF& h, const Domain& pos) const
{
if (!millerSimilarity)
return h;
return h * exp(interpolate(millerField, pos)-variance);
}
/// Transformation: Conductivity: Reference to Miller similar (if activated)
/** \f[ K^* = K \exp\left( 2 (v - \text{Var}(v))\right) \f]
* with the raw random field value \f$v\f$ and its variance.
* \param h Matric head
* \param pos Position in global coordinates
* \return Transformed value
*/
RF condRefToMiller (const RF& K, const Domain& pos) const
{
if (!millerSimilarity)
return K;
return K * exp (2.0*(interpolate(millerField, pos)-variance));
}
/// Transformation: Conductivity: Miller similar to Reference (if activated)
/** \see condRefToMiller
*/
RF condMillerToRef (const RF& k, const Domain& pos) const
{
if (!millerSimilarity)
return k;
return k * exp (-2.0*(interpolate(millerField, pos)-variance));
}
/// Return gravity vector
const Vector& gravity() const
{
return gVector;
}
public:
/******************************************
* UTILITY FUNCTIONS
******************************************/
/// Thin wrapper around the interpolator evaluation function
RF interpolate(const Parameter<T>& par, const Domain& pos) const
{
return interpolator.evaluate(par.get_values(),pos + offset);
}
protected:
/// Opens H5 array for reading
void open_array(std::string groupName)
{
h5file.open(groupName);
array_open = true;
}
/// Reads parameter array for Parameter p from H5 file
void read_parameter(Parameter<T>& p)
{
if (!array_open)
DUNE_THROW(Exception,"Array file must be opened before reading a parameter");
Array fieldArray;
h5file.read_dataset(fieldArray, fieldCells, p.name);
p.set_values(fieldArray);
}
/// Closes H5 file
void close_array()
{
if (!array_open)
DUNE_THROW(Exception,"Array file already closed");
h5file.close();
array_open = false;
}
};
}
}
#endif // DUNE_DORIE_PARAMETERS_BASE_HH
// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_DORIE_PARAMETERS_FACTORY_HH
#define DUNE_DORIE_PARAMETERS_FACTORY_HH
#include <string>
#include <memory>
#include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parallel/mpihelper.hh>
#include "param_van_genuchten.hh"
namespace Dune {
namespace Dorie {
/// Factory class for a generic parameterization (base class of ParameterizationBase)
template<typename T, class Interpolator>
class ParameterizationFactory
{
protected:
const ParameterTree& config;
const Dune::MPIHelper& helper;
Interpolator& interpolator;
const int verbose;
const std::string arrayPath;
const std::string interpMethod;
H5File h5array;
std::string parType;
public:
explicit ParameterizationFactory(const ParameterTree& config_,
const Dune::MPIHelper& helper_,
Interpolator& interpolator_,
const int verbose_)
: config(config_), helper(helper_), interpolator(interpolator_), verbose(verbose_)
, arrayPath(config.get<std::string>("parameters.arrayFile"))
, interpMethod(config.get<std::string>("parameters.interpolation"))
, h5array(arrayPath,verbose)
, parType(h5array.read_attribute<std::string>("parameterization","parameters"))
{}
std::unique_ptr<ParameterizationBase<T, Interpolator>> create()
{
std::unique_ptr<ParameterizationBase<T, Interpolator>> param;
if(parType == "vanGenuchten")
param = std::make_unique<MualemVanGenuchten<T, Interpolator>>(config, helper, interpolator, h5array, verbose);
else
DUNE_THROW(Dune::IOError,"Unrecognized parameterization type encountered in array file: " + parType);
return param;
}
};
}
}
#endif
#ifndef DUNE_DORIE_PARAM_INTERFACE_HH
#define DUNE_DORIE_PARAM_INTERFACE_HH
#ifndef DUNE_DORIE_PARAM_RICHARDS_HH
#define DUNE_DORIE_PARAM_RICHARDS_HH
#include <cmath>
#include <vector>
#include <map>
#include <utility>
#include <functional>
#include <memory>
#include <yaml-cpp/yaml.h>
#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>
#include <string>
#include <map>
namespace Dune {
namespace Dorie {
/// Top-level parameterization interface.
/** This class stores the parameters, maps them to grid entities, and provides
* functions to access the parameterization. For doing so, the
* FlowParameters::bind function has to be called first to cache the
* parameters of the given grid entity. Parameters are only stored for entities
* of codim 0.
*
* \detail The parameterization consists of three structures:
* - This top-level class serves as interface and storage. It also casts to
* data types used by the DUNE operators
* - The Dorie::RichardsParameterization interface that provides strong
* data types for the base parameters and purely virtual functions that
* need to be defined by derived parameterizations
* - The actual parameterization defining all parameters and functions
*/
/// Class defining the interface for parameterizations of Richards equation
template <typename Traits>
class FlowParameters
class RichardsParameterization
{
private:
using RF = typename Traits::RF;
using Grid = typename Traits::Grid;
using LevelGridView = typename Grid::LevelGridView;
using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<
LevelGridView
>;
using Vector = typename Traits::Vector;
// Define the storage and cache for parameterizations and scaling
/// Index type used by the element mapper
using Index = typename Mapper::Index;
/// Base class of the parameterization
using Parameterization = Dorie::RichardsParameterization<Traits>;
using ScalingType = typename Parameterization::Scaling::Type;
using ScalingFactor = typename Parameterization::ScalingFactor;
public:
/// Storage for parameters and skaling factors
using ParameterStorage = std::vector<
std::tuple<std::shared_ptr<Parameterization>,
ScalingType,
ScalingFactor
>
>;
// These are standard types for all parameterizations of Richards Equation
/// Type of water saturation
struct Saturation
{
RF value;
};
/// Need this gruesome typedef because 'map::value_type' has const key
using Cache = typename ParameterStorage::value_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;
/// Gravity vector
Vector _gravity;
/// Type of volumetric water content
struct WaterContent
{
RF value;
};
private:
/// Type of hydraulic conductivity
struct Conductivity
{
RF value;
};
/// Check if an entity has been cached
void verify_cache () const
/// Type of matric head
struct MatricHead
{
if (not std::get<std::shared_ptr<Parameterization>>(_cache)) {
DUNE_THROW(Dune::Exception, "No parameterization cached. "
<< "Call 'bind' before accessing the parameterization.");
}
}
RF value;
};
public:
/// Constructor.
/** Create a level grid view of level 0, create a mapper on it, and store
* the config tree.
* \param config Configuration file tree
* \param grid Shared pointer to the grid
* \param element_index_map The mapping from grid entity index to
* parameterization index.
*/
FlowParameters (
const Dune::ParameterTree& config,
const std::shared_ptr<Grid> grid,
const std::vector<int>& element_index_map
):
_config(config),
_gv(grid->levelGridView(0)),
_mapper(_gv, Dune::mcmgElementLayout()),
_gravity(0.0)
/// Generic scaling factor
struct ScalingFactor
{
_gravity[Traits::dim - 1] = -1.0;
build_parameterization(element_index_map);
}
RF value;
};
/// Return normalized gravity vector
const Vector& gravity() const { return _gravity; }
/// Define scaling types for parameterization
struct Scaling
{
enum Type
{
Miller, //!< Geometric similarity
Warrick,
None //!< No scaling applied
};
};
/// Return the current cache
const Cache& cache() const
/// Parameter defining the saturated hydraulic conductivity
struct SaturatedConductivity
{
verify_cache();
return _cache;
}
RF value;
inline static const std::string name = "k0";
};
/// Bind to a grid entity. Required to call parameterization functions!
/** \param entity Grid entity (codim 0) to bind to
*/
template <class Entity>
void bind (Entity entity) const
/// Parameter defining the residual water content
struct ResidualWaterContent
{
// retrieve index of top father element of chosen entity
while(entity.hasFather()) {
entity = entity.father();
}
const auto index = _mapper.index(entity);
_cache = _param.at(index);
}
RF value;
inline static const std::string name = "theta_r";
};
/// Return a scaled version of the Conductivity function
/** Uses the function of the underlying parameterization and applies
* scaling as specified by SkalingType. Cast to operator-internal RF.
* \return Function: Saturation -> Conductivity
*/
std::function<RF(RF)> conductivity_f () const
/// Parameter defining the saturated water content
struct SaturatedWaterContent
{
verify_cache();
const auto& par =
std::get<std::shared_ptr<Parameterization>>(_cache);
const auto& type = std::get<ScalingType>(_cache);
const auto cond_f = par->conductivity_f();
using Saturation = typename Parameterization::Saturation;
RF value;
inline static const std::string name = "theta_s";
};
//! Value of the residual water content.
ResidualWaterContent _theta_r;
//! Value of the saturated water content.
SaturatedWaterContent _theta_s;
//! Value of the saturated hydraulic conductivity.
SaturatedConductivity _k0;
//! The name of this parameterization instance, associated with the layer.
const std::string _name;
if (type == ScalingType::Miller) {
auto& scale = std::get<ScalingFactor>(_cache).value;
return [cond_f, &scale](const RF saturation){
return cond_f(Saturation{saturation}).value * scale * scale;
};
}
public:
return [cond_f](const RF saturation){
return cond_f(Saturation{saturation}).value;
};
}
/// Construct with default-initialized parameters
/** \param name The name associated with this soil layer
* \param parameters Tuple of parameters to use in this parameterization
*/
RichardsParameterization (const std::string name) :
_name(name)
{ }
/// Return a scaled version of the Saturation function
/** Uses the function of the underlying parameterization and applies
* scaling as specified by SkalingType. Cast to operator-internal RF.
* \return Function: Matric Head -> Saturation
/// Construct from a tuple of parameters
/** \param name The name associated with this soil layer
* \param parameters Tuple of parameters to use in this parameterization
*/
template<typename... Args>
RichardsParameterization (const std::string name,
const std::tuple<Args...> parameters) :
_theta_r(std::get<ResidualWaterContent>(parameters)),
_theta_s(std::get<SaturatedWaterContent>(parameters)),
_k0(std::get<SaturatedConductivity>(parameters)),
_name(name)
{ }
/// Default constructor (virtual).
virtual ~RichardsParameterization () = default;
/// Return the name of this parameterization instance.
const std::string& get_name() const { return _name; }
/// Return a bound version of the water content function
/** \return Function: Saturation -> Water content
*/
std::function<RF(RF)> saturation_f () const
std::function<WaterContent(const Saturation)> water_content_f () const
{
verify_cache();
const auto& par =
std::get<std::shared_ptr<Parameterization>>(_cache);
const auto& type = std::get<ScalingType>(_cache);
const auto sat_f = par->saturation_f();
using MatricHead = typename Parameterization::MatricHead;
if (type == ScalingType::Miller) {
const auto& scale = std::get<ScalingFactor>(_cache).value;
return [sat_f, &scale](const RF matric_head){
return sat_f(MatricHead{matric_head * scale}).value;
return [this](const Saturation sat) {
return WaterContent{
sat.value * (_theta_s.value - _theta_r.value) + _theta_r.value
};
}
return [sat_f](const RF matric_head){