diff --git a/dune/dorie/interface/richards_simulation.hh b/dune/dorie/interface/richards_simulation.hh index 674055e13044ffe8b7f66b764491a59bf02c5a13..d017a7e8b12b418257adb672df2a97d71d29296d 100644 --- a/dune/dorie/interface/richards_simulation.hh +++ b/dune/dorie/interface/richards_simulation.hh @@ -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" diff --git a/dune/dorie/solver/adapters/conductivity.hh b/dune/dorie/solver/adapters/conductivity.hh index 8537e9c185db44f5375f89aa8d98994357d2d6e8..bd1878597ab27a911eb3fd7f9d2fa2361d0c06c7 100644 --- a/dune/dorie/solver/adapters/conductivity.hh +++ b/dune/dorie/solver/adapters/conductivity.hh @@ -1,8 +1,6 @@ #ifndef DUNE_DORIE_CONDUCTIVITY_ADAPTER_HH #define DUNE_DORIE_CONDUCTIVITY_ADAPTER_HH -#include "../param_richards.hh" - namespace Dune{ namespace Dorie{ diff --git a/dune/dorie/solver/param_base.hh b/dune/dorie/solver/param_base.hh deleted file mode 100644 index b9ef5cbc6fe273e2fd3d66bde30292c1cc1a5eb6..0000000000000000000000000000000000000000 --- a/dune/dorie/solver/param_base.hh +++ /dev/null @@ -1,270 +0,0 @@ -// -*- tab-width: 4; indent-tabs-mode: nil -*- -#ifndef DUNE_DORIE_PARAMETERS_BASE_HH -#define DUNE_DORIE_PARAMETERS_BASE_HH - -#include -#include - -#include -#include -#include - -#include "util_h5tools.hh" - -namespace Dune { - namespace Dorie { - - /// Simple class providing a data structure for storing parameters - template - 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 - 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*> 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 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("extensions","parameters")) - , millerSimilarity(h5file.read_attribute("millerSimilarity","parameters")) - , variance(h5file.read_attribute("variance","parameters")) - , grid_extensions(config.get("grid.extensions")) - , scale(config.get("parameters.scale")) - , offset(config.get("parameters.offset")) - , millerField("raw_field") - { - gVector=0; - gVector[dim-1]=-1.; - - for(Index i=0; i& 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& 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 diff --git a/dune/dorie/solver/param_factory.hh b/dune/dorie/solver/param_factory.hh deleted file mode 100644 index 75d8f4d1d02503e821780d425101c387b7107e77..0000000000000000000000000000000000000000 --- a/dune/dorie/solver/param_factory.hh +++ /dev/null @@ -1,60 +0,0 @@ -// -*- tab-width: 4; indent-tabs-mode: nil -*- -#ifndef DUNE_DORIE_PARAMETERS_FACTORY_HH -#define DUNE_DORIE_PARAMETERS_FACTORY_HH - -#include -#include - -#include -#include -#include - -#include "param_van_genuchten.hh" - -namespace Dune { - namespace Dorie { - - /// Factory class for a generic parameterization (base class of ParameterizationBase) - template - 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("parameters.arrayFile")) - , interpMethod(config.get("parameters.interpolation")) - , h5array(arrayPath,verbose) - , parType(h5array.read_attribute("parameterization","parameters")) - {} - - std::unique_ptr> create() - { - std::unique_ptr> param; - if(parType == "vanGenuchten") - param = std::make_unique>(config, helper, interpolator, h5array, verbose); - else - DUNE_THROW(Dune::IOError,"Unrecognized parameterization type encountered in array file: " + parType); - - return param; - } - }; - } -} - -#endif diff --git a/dune/dorie/solver/param_interface.hh b/dune/dorie/solver/param_interface.hh index 6059a86b7eb2654e72c8998977223009a62b9c6b..ff3bd649697160df8aa30448e236600d950f9422 100644 --- a/dune/dorie/solver/param_interface.hh +++ b/dune/dorie/solver/param_interface.hh @@ -1,359 +1,156 @@ -#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 -#include -#include -#include #include -#include - -#include - -#include -#include - -#include -#include +#include +#include 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 -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; - - using ScalingType = typename Parameterization::Scaling::Type; - using ScalingFactor = typename Parameterization::ScalingFactor; +public: - /// Storage for parameters and skaling factors - using ParameterStorage = std::vector< - std::tuple, - 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>(_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, - const std::vector& 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 - 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 conductivity_f () const + /// Parameter defining the saturated water content + struct SaturatedWaterContent { - verify_cache(); - const auto& par = - std::get>(_cache); - const auto& type = std::get(_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(_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 + RichardsParameterization (const std::string name, + const std::tuple parameters) : + _theta_r(std::get(parameters)), + _theta_s(std::get(parameters)), + _k0(std::get(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 saturation_f () const + std::function water_content_f () const { - verify_cache(); - const auto& par = - std::get>(_cache); - const auto& type = std::get(_cache); - const auto sat_f = par->saturation_f(); - using MatricHead = typename Parameterization::MatricHead; - - if (type == ScalingType::Miller) { - const auto& scale = std::get(_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){ - return sat_f(MatricHead{matric_head}).value; }; } - /// Return a scaled version of the Water Content function - /** Uses the function of the underlying parameterization and applies - * scaling as specified by SkalingType. Cast to operator-internal RF. - * \return Function: Saturation -> Water Content + /// Return a bound version of the saturation function + /** \return Function: Matric Head -> Water content */ - std::function water_content_f () const - { - verify_cache(); - const auto& par = - std::get>(_cache); - const auto wc_f = par->water_content_f(); - using Saturation = typename Parameterization::Saturation; - - return [wc_f](const RF saturation) { - return wc_f(Saturation{saturation}).value; - }; - } + virtual std::function saturation_f () const + = 0; - /// Copy parameters from the 'old' parameterization implementation - /** \todo Remove once new interface is established - * \param p_base Old parameterization base class + /// Return a bound version of the conductivity function + /** \return Function: Saturation -> Hydraulic conductivity */ - template - 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 p_values; - std::transform(p_base.parameters.begin(), p_base.parameters.end(), - std::back_inserter(p_values), - get_param_value - ); - - // read parameter values - using RP = Parameterization; - using MvG = MualemVanGenuchtenParameterization; - typename RP::ResidualWaterContent theta_r{p_values.at(0)}; - typename RP::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 RP::SaturatedConductivity k0{p_values.at(5)}; - - // can create parameters out of order by using tuple - auto mvg = - std::make_shared> - ( - std::make_tuple("layer", theta_s, theta_r, k0, alpha, tau, n) - ); - - ScalingType type = RP::Scaling::None; - ScalingFactor scale{1.0}; - - if (p_base.millerSimilarity) { - type = ScalingType::Miller; - const auto val = get_param_value(&p_base.millerField); - scale = ScalingFactor{ std::exp(val - p_base.variance) }; - } - - // insert parameters into map - const auto index = _mapper.index(entity); - const auto res = _param.emplace(index, - std::make_tuple(mvg, type, scale) - ); - - if (not res.second) { - DUNE_THROW(Dune::Exception, "Parameterization for entity" - << index << " already exists." - ); - } - } - } - -private: - - /// Build the parameterization from an element mapping and the input file. - void build_parameterization (const std::vector& element_index_map) - { - const auto parameterization_map = read_parameter_file(); - _param.reserve(element_index_map.size()); - - // build dummy scaling - ScalingType s_type{Parameterization::Scaling::None}; - ScalingFactor s_factor{1.0}; - - // insert parameterization and dummy scaling for each element - for (auto&& index : element_index_map) { - auto p = parameterization_map.at(index); - _param.push_back(std::make_tuple(p, s_type, s_factor)); - } + virtual std::function conductivity_f () + const = 0; - // check that mapper can map to parameterization - assert(_param.size() == _mapper.size()); - } - - /// Read the YAML parameter file and insert information into a map. - std::map> read_parameter_file () - const - { - const auto param_file_name = _config.get("parameters.file"); - std::cout << "Reading file " << param_file_name << std::endl; - YAML::Node param_file = YAML::LoadFile(param_file_name); - - // mapping: index on grid to parameterization object (will be returned) - std::map> ret; - - // set to check for duplicate indices - std::set param_index_set; - - // iterate over all volume nodes - for (auto&& node : param_file["volumes"]) - { - const auto name = node.first.as(); - const YAML::Node param_info = node.second; - - // read relevant data - const auto type = param_info["type"].as(); - const auto index = param_info["index"].as(); - if (param_index_set.count(index) != 0) { - DUNE_THROW(IOError, "Parameter file contains index " - + std::to_string(index) + " at least twice!"); - } - param_index_set.emplace(index); - - // get parameterization object - auto parameterization = parameterization_selector(type, name); - auto parameter_values = parameterization->parameters(); - - // fill parameterization with parameter entries - const YAML::Node yaml_parameters = param_info["parameters"]; - for (auto&& yaml_param : yaml_parameters) - { - const auto p_name = yaml_param.first.as(); - const auto p_value = yaml_param.second.as(); - parameter_values.at(p_name) = p_value; - } - - ret.emplace(index, parameterization); - } - - return ret; - } - - /// Return a default-initialized parameterization object - /** \param type The type indicating the parameterization implementation - * \param name The name of the parameterization object (layer type) + /// Return a map referecing all parameters by their names. + /** \return Map. Key: Name of parameter (string). + * Value: Value of parameter (double&) */ - std::shared_ptr parameterization_selector ( - const std::string type, - const std::string name) const - { - if (type == MualemVanGenuchtenParameterization::type) { - return std::make_shared>(name); - } - else - DUNE_THROW(IOError, "Unknown parameterization type " + type); - } + virtual std::map parameters () = 0; + }; } // namespace Dune } // namespace Dorie -#endif // DUNE_DORIE_PARAM_INTERFACE_HH \ No newline at end of file +#endif // DUNE_DORIE_PARAM_RICHARDS_HH \ No newline at end of file diff --git a/dune/dorie/solver/param_mvg.hh b/dune/dorie/solver/param_mvg.hh index 3d8bf3d8693e866b7d86b61655646d728a67421c..d34dedbb616dee7c3e99049832c1eb1b86a894dd 100644 --- a/dune/dorie/solver/param_mvg.hh +++ b/dune/dorie/solver/param_mvg.hh @@ -5,7 +5,7 @@ #include #include -#include +#include namespace Dune { namespace Dorie { diff --git a/dune/dorie/solver/param_richards.hh b/dune/dorie/solver/param_richards.hh deleted file mode 100644 index ff3bd649697160df8aa30448e236600d950f9422..0000000000000000000000000000000000000000 --- a/dune/dorie/solver/param_richards.hh +++ /dev/null @@ -1,156 +0,0 @@ -#ifndef DUNE_DORIE_PARAM_RICHARDS_HH -#define DUNE_DORIE_PARAM_RICHARDS_HH - -#include -#include -#include - -namespace Dune { -namespace Dorie { - -/// Class defining the interface for parameterizations of Richards equation -template -class RichardsParameterization -{ -private: - using RF = typename Traits::RF; - -public: - - // These are standard types for all parameterizations of Richards Equation - /// Type of water saturation - struct Saturation - { - RF value; - }; - - /// Type of volumetric water content - struct WaterContent - { - RF value; - }; - - /// Type of hydraulic conductivity - struct Conductivity - { - RF value; - }; - - /// Type of matric head - struct MatricHead - { - RF value; - }; - - /// Generic scaling factor - struct ScalingFactor - { - RF value; - }; - - /// Define scaling types for parameterization - struct Scaling - { - enum Type - { - Miller, //!< Geometric similarity - Warrick, - None //!< No scaling applied - }; - }; - - /// Parameter defining the saturated hydraulic conductivity - struct SaturatedConductivity - { - RF value; - inline static const std::string name = "k0"; - }; - - /// 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"; - }; - - //! 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; - -public: - - /// 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) - { } - - /// 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 - RichardsParameterization (const std::string name, - const std::tuple parameters) : - _theta_r(std::get(parameters)), - _theta_s(std::get(parameters)), - _k0(std::get(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 water_content_f () const - { - return [this](const Saturation sat) { - return WaterContent{ - sat.value * (_theta_s.value - _theta_r.value) + _theta_r.value - }; - }; - } - - /// Return a bound version of the saturation function - /** \return Function: Matric Head -> Water content - */ - virtual std::function saturation_f () const - = 0; - - /// Return a bound version of the conductivity function - /** \return Function: Saturation -> Hydraulic conductivity - */ - virtual std::function conductivity_f () - const = 0; - - /// Return a map referecing all parameters by their names. - /** \return Map. Key: Name of parameter (string). - * Value: Value of parameter (double&) - */ - virtual std::map parameters () = 0; - -}; - -} // namespace Dune -} // namespace Dorie - -#endif // DUNE_DORIE_PARAM_RICHARDS_HH \ No newline at end of file diff --git a/dune/dorie/solver/richards_param.hh b/dune/dorie/solver/richards_param.hh new file mode 100644 index 0000000000000000000000000000000000000000..279d7de8907cfe0582dd6d87711649e99db5d10a --- /dev/null +++ b/dune/dorie/solver/richards_param.hh @@ -0,0 +1,298 @@ +#ifndef DUNE_DORIE_PARAM_INTERFACE_HH +#define DUNE_DORIE_PARAM_INTERFACE_HH + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include + +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 + */ +template +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 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; + + using ScalingType = typename Parameterization::Scaling::Type; + using ScalingFactor = typename Parameterization::ScalingFactor; + + /// Storage for parameters and skaling factors + using ParameterStorage = std::vector< + std::tuple, + ScalingType, + ScalingFactor + > + >; + + /// 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; + +private: + + /// Check if an entity has been cached + void verify_cache () const + { + if (not std::get>(_cache)) { + DUNE_THROW(Dune::Exception, "No parameterization cached. " + << "Call 'bind' before accessing the parameterization."); + } + } + +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, + const std::vector& element_index_map + ): + _config(config), + _gv(grid->levelGridView(0)), + _mapper(_gv, Dune::mcmgElementLayout()), + _gravity(0.0) + { + _gravity[Traits::dim - 1] = -1.0; + build_parameterization(element_index_map); + } + + /// Return normalized gravity vector + const Vector& gravity() const { return _gravity; } + + /// Return the current cache + const Cache& cache() const + { + verify_cache(); + return _cache; + } + + /// Bind to a grid entity. Required to call parameterization functions! + /** \param entity Grid entity (codim 0) to bind to + */ + template + void bind (Entity entity) const + { + // 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); + } + + /// 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 conductivity_f () const + { + verify_cache(); + const auto& par = + std::get>(_cache); + const auto& type = std::get(_cache); + const auto cond_f = par->conductivity_f(); + using Saturation = typename Parameterization::Saturation; + + if (type == ScalingType::Miller) { + auto& scale = std::get(_cache).value; + return [cond_f, &scale](const RF saturation){ + return cond_f(Saturation{saturation}).value * scale * scale; + }; + } + + return [cond_f](const RF saturation){ + return cond_f(Saturation{saturation}).value; + }; + } + + /// 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 + */ + std::function saturation_f () const + { + verify_cache(); + const auto& par = + std::get>(_cache); + const auto& type = std::get(_cache); + const auto sat_f = par->saturation_f(); + using MatricHead = typename Parameterization::MatricHead; + + if (type == ScalingType::Miller) { + const auto& scale = std::get(_cache).value; + return [sat_f, &scale](const RF matric_head){ + return sat_f(MatricHead{matric_head * scale}).value; + }; + } + + return [sat_f](const RF matric_head){ + return sat_f(MatricHead{matric_head}).value; + }; + } + + /// Return a scaled version of the Water Content function + /** Uses the function of the underlying parameterization and applies + * scaling as specified by SkalingType. Cast to operator-internal RF. + * \return Function: Saturation -> Water Content + */ + std::function water_content_f () const + { + verify_cache(); + const auto& par = + std::get>(_cache); + const auto wc_f = par->water_content_f(); + using Saturation = typename Parameterization::Saturation; + + return [wc_f](const RF saturation) { + return wc_f(Saturation{saturation}).value; + }; + } + +private: + + /// Build the parameterization from an element mapping and the input file. + void build_parameterization (const std::vector& element_index_map) + { + const auto parameterization_map = read_parameter_file(); + _param.reserve(element_index_map.size()); + + // build dummy scaling + ScalingType s_type{Parameterization::Scaling::None}; + ScalingFactor s_factor{1.0}; + + // insert parameterization and dummy scaling for each element + for (auto&& index : element_index_map) { + auto p = parameterization_map.at(index); + _param.push_back(std::make_tuple(p, s_type, s_factor)); + } + + // check that mapper can map to parameterization + assert(_param.size() == _mapper.size()); + } + + /// Read the YAML parameter file and insert information into a map. + std::map> read_parameter_file () + const + { + const auto param_file_name = _config.get("parameters.file"); + std::cout << "Reading file " << param_file_name << std::endl; + YAML::Node param_file = YAML::LoadFile(param_file_name); + + // mapping: index on grid to parameterization object (will be returned) + std::map> ret; + + // set to check for duplicate indices + std::set param_index_set; + + // iterate over all volume nodes + for (auto&& node : param_file["volumes"]) + { + const auto name = node.first.as(); + const YAML::Node param_info = node.second; + + // read relevant data + const auto type = param_info["type"].as(); + const auto index = param_info["index"].as(); + if (param_index_set.count(index) != 0) { + DUNE_THROW(IOError, "Parameter file contains index " + + std::to_string(index) + " at least twice!"); + } + param_index_set.emplace(index); + + // get parameterization object + auto parameterization = parameterization_selector(type, name); + auto parameter_values = parameterization->parameters(); + + // fill parameterization with parameter entries + const YAML::Node yaml_parameters = param_info["parameters"]; + for (auto&& yaml_param : yaml_parameters) + { + const auto p_name = yaml_param.first.as(); + const auto p_value = yaml_param.second.as(); + parameter_values.at(p_name) = p_value; + } + + ret.emplace(index, parameterization); + } + + return ret; + } + + /// Return a default-initialized parameterization object + /** \param type The type indicating the parameterization implementation + * \param name The name of the parameterization object (layer type) + */ + std::shared_ptr parameterization_selector ( + const std::string type, + const std::string name) const + { + if (type == MualemVanGenuchtenParameterization::type) { + return std::make_shared>(name); + } + else + DUNE_THROW(IOError, "Unknown parameterization type " + type); + } +}; + +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_INTERFACE_HH \ No newline at end of file diff --git a/dune/dorie/test/test-parameterization.hh b/dune/dorie/test/test-parameterization.hh index ed9467591f5144c34194aaf91c3130324ca6c37b..0656796f2110a5bca33a4d8ae7226c73839c78b1 100644 --- a/dune/dorie/test/test-parameterization.hh +++ b/dune/dorie/test/test-parameterization.hh @@ -11,7 +11,7 @@ #include "../interface/richards_simulation.cc" -#include "../solver/param_interface.hh" +#include "../solver/param_mvg.hh" namespace Dune { namespace Dorie {