Merge branch '110-add-scaling-field-input' into 'master'

Resolve "Add scaling field input"

Closes #18 and #110

See merge request !110

CHANGELOG.md

@@ -23,6 +23,7 @@
 * New class [`SimulationTransport`](dune/dorie/interface/transport_simulation.hh)
     which solves the transport equation with a finite volume scheme for give
     grid functions of water content and velocity.
+* Infrastructure for the input of Miller scaling fields. !110
 
 ### Changed
 * `Simulation` is renamed `RichardsSimulation` and moved to 

doc/default_files/CMakeLists.txt

@@ -116,9 +116,8 @@ function(create_default_config INPUT OUTPUT SOURCE_DIR CSS)
     endif()
 endfunction()
 
-# copy BC files
-file(COPY 2d_infiltr.bcdat DESTINATION .)
-file(COPY 3d_infiltr.bcdat DESTINATION .)
+# copy BC & parameter files
+file(COPY 2d_infiltr.bcdat 3d_infiltr.bcdat param.yml DESTINATION .)
 
 # Random field generator
 scrape_parameters(

doc/default_files/field-parameters.xml

@@ -50,10 +50,12 @@ adding an empty line, make text **bold** or ``monospaced``.
     <parameter name="converter">
       <definition> Identifier of the random field converter. The converter is
         applied after the random field is created and modifies it.
+
+        Use ``exponential`` for creating a suitable Miller scaling field.
       </definition>
-      <values> none, binary </values>
+      <values> none, binary, exponential </values>
       <suggestion> binary </suggestion>
-      <comment> none, binary </comment>
+      <comment> none, binary, exponential </comment>
     </parameter>
 
     <parameter name="tempDir">
@@ -136,4 +138,17 @@ adding an empty line, make text **bold** or ``monospaced``.
       <suggestion> 0 1 </suggestion>
     </parameter>
   </category>
+
+  <category name="converter.exponential">
+    <parameter name="varianceScaling">
+      <definition> Substract the variance from the random field values before
+        applying the exponential. For a random field with ``gaussian``
+        covariance, this will result in a conductivity field with the same
+        macroscopic properties that a homogeneous field, in the sense that
+        average water flux is the same (Roth 1995).
+      </definition>
+      <values> bool </values>
+      <suggestion> true </suggestion>
+    </parameter>
+  </category>
 </dorie>

doc/default_files/param.yml

+volumes:
+  sand:
+    index: 0
+    type: MvG
+    parameters:
+      alpha: -2.3
+      n: 4.17
+      k0: 2.2E-5
+      theta_r: 0.03
+      theta_s: 0.31
+      tau: -1.1
+
+  silt:
+    index: 1
+    type: MvG
+    parameters:
+      alpha: -0.7
+      n: 1.3
+      k0: 1.0E-5
+      theta_r: 0.01
+      theta_s: 0.41
+      tau: 0.0
+
+scaling:
+  type: None

doc/doxygen/Doxylocal

@@ -21,4 +21,4 @@ USE_MATHJAX = YES
 # EXCLUDE 	+= @top_srcdir@/...
 
 # Include documents with other extesions
-FILE_PATTERNS += *.dox
\ No newline at end of file
+FILE_PATTERNS += *.dox

doc/man-parameter-file.rst

@@ -14,6 +14,11 @@ If you insert an unstructured grid with GMSH, you have to incorporate the
 medium information. If you create a rectangular grid, you need to add a mapping
 dataset, specifying the medium index for each (initial) grid cell.
 
+Additionally, you can specify settings for the small-scale heterogeneities
+of the soil. They require an input via H5 datasets that contain appropriate
+scaling factors. You can conveniently create such datasets using the random
+field generator of the :doc:`command line interface <python-dorie-wrapper>`.
+
 YAML Parameter File
 -------------------
 This file is used to specify the parameterization for each medium inside the
@@ -24,24 +29,49 @@ name must be specified via the ``parameters.file`` key of the
 The top-level mapping must contain the key ``volumes``. This key contains a
 sequence of arbitrary parameterization names. These, in turn, contain the
 parameterization ``type``, the medium ``index``, and the actual ``parameters``.
+The medium index must be an **integer**.
+
+Heterogeneities throughout the entire domain are set via the top level key
+``scaling``. It must contain a supported scaling ``type``, which may default
+to ``none``. In case an actual scaling is to be applied, the section must
+contain the key ``data``, which specifies the datasets required for this type
+of scaling.
+
+Every ``scaling_section`` has the same keys: The H5 filepath is given by
+``file``, and the internal dataset path by ``dataset``. You can then choose an
+``interpolation`` method for the dataset, which requires you to insert values
+for the physical ``extensions`` of the dataset and the ``offset`` of its (front)
+lower left corner from the grid origin.
 
 .. code-block:: yaml
 
     volumes:
       <name-0>:
-        type: <type>
-        index: <index>
+        type: <string>
+        index: <int>
         parameters:
-          <param-name-0>: <value>
-          # ...
+          <param-name-0>: <float>
+          # more parameters ...
+
+      # more volumes ...
+
+    scaling:
+      type: <string>
+      data:
+        <scaling_section>:
+          file: <string>
+          dataset: <string>
+          interpolation: <string>
+          extensions: <sequence>
+          offset: <sequence>
 
-      # ...
+        # more scaling sections ...
 
 You find a documentation of the parameterizations implemented in DORiE
 along with the parameters they require below.
 The parameterization ``type`` must match the static ``type`` member of one of
 them. Likewise, the parameters must match the names of parameters
-these objects use. The medium index must be an **integer**.
+these objects use.
 
 Medium Mapping Dataset
 ----------------------
@@ -135,10 +165,120 @@ script. It is readily installed into the
     Physical Line(5) = {4}; // top
     Physical Line(6) = {5, 6}; // left
 
+Scaling Field Datasets
+----------------------
+One or more datasets in possibly multiple HDF5_ files are required for creating
+a medium with small-scale heterogeneities. The datasets must consist of
+floating-point values and must have the same dimensionality as the grid.
+Other than that, there is no restriction on their shape because they are
+inserted into interpolators. The interpolators supply scaling factor values
+based on positions on the grid and thus create a grid-independent
+representation. Scaling field input works the same for any supported grid type.
+
+During setup, the program reads the interpolator values at the **barycenter**
+of every grid cell on the **coarsest** grid configuration. These values are
+stored and referenced with the grid cell. If the cell is refined, the same
+scaling factors apply in all its child cells.
+
+Supported Parameter File Types
+------------------------------
+This section lists the available types for parameterizations, scalings, and
+interpolators in a understandable format. You can also dig through the code
+documentation below, or the code itself, to retrieve this information.
+
+Parameterizations
+^^^^^^^^^^^^^^^^^
+
+.. object:: Mualem–van Genuchten parameterization
+
+    Implements the following parameterization functions:
+
+    .. math::
+
+        \Theta (h_m) &= \left[ 1 + \left[ \alpha h_m \right]^n \right]^{-1+1/n}
+        \\
+        K (\Theta) &= K_0 \Theta^\tau
+                        \left[ 1 -
+                            \left[ 1 - \Theta^{n / \left[ n-1 \right]}
+                            \right]^{1-1/n}
+                        \right]^2
+
+    * ``type: MvG``
+
+    Parameters:
+        * ``theta_r``: Residual water content :math:`\theta_r`
+        * ``theta_s``: Saturated water content (porosity) :math:`\theta_s`
+        * ``k0``: Saturated conductivity :math:`K_0`
+        * ``alpha``: Air-entry value :math:`\alpha`
+        * ``n``: Retention curve shape factor :math:`n`
+        * ``tau``: Skew factor :math:`\tau`
+
+    YAML template:
+
+    .. code-block:: yaml
+
+        type: MvG
+        parameters:
+          theta_r:
+          theta_s:
+          k0:
+          alpha:
+          n:
+          tau:
+
+Scalings
+^^^^^^^^
+Every ``scaling_section`` has the following layout:
+
+.. code-block:: yaml
+
+    <scaling_section>:
+      file: <string>
+      dataset: <string>
+      interpolation: <string>
+      extensions: <sequence>
+      offset: <sequence>
+
+The values of ``extensions`` and ``offset`` are sequences containing the
+coordinates in the respective spatial dimensions.
+
+.. object:: Miller scaling
+
+    Assumes geometric similarity between scaled regions. Applies the following
+    scaling:
+
+    .. math ::
+
+        \Theta &= \Theta (h_m \cdot \xi_M)\\
+        K &= K (\Theta) \cdot \xi_M^2
+
+    * ``type: Miller``
+
+    Scaling sections:
+        * ``scale_miller``: Scaling factor :math:`\xi_M` to be applied onto
+          matric head and conductivity simultaneously.
+
+    YAML template:
+
+    .. code-block:: yaml
+
+        type: Miller
+        data:
+          scale_miller:
+            # ...
+
+Interpolators
+^^^^^^^^^^^^^
+
+.. object:: Nearest neighbor interpolator
+
+    Interprets dataset values as cell values.
+
+    * ``interpolation: nearest``
+
 
 Parameterization class documentation
 ------------------------------------
-
 The following is the doxygen documentation of the classes implementing the soil
 parameterization. Parameterizations specified in the YAML parameter file
 must match one of the instances derived from ``RichardsParameterization``.

dune/dorie/common/grid_creator.hh

@@ -315,6 +315,9 @@ private:
                           cells
                          );
 
+        // shape has to be reversed
+        std::reverse(begin(cells), end(cells));
+
         // check if cells has correct size and extensions
         if (cells.size() != Grid::dimension) {
             DUNE_THROW(IOError, "Mapping dataset has wrong dimensions!");

dune/dorie/common/groups.dox

@@ -12,9 +12,24 @@
   @todo document models
 @}
 
+@defgroup Interpolators Interpolators
+@{
+  @ingroup Common
+  @brief Interpolation of structured data for evaluation on any point
+      inside the grid.
+
+  ## Overview
+
+  Interpolators have the task to store input data and evaluate it for global
+  positions on the grid. The data is not shared or scattered across processes,
+  but exists independently on every process. Therefore, using these structures
+  is only efficient in a setup process. Interpolators should not be used to
+  retrieve data repeatedly, especially not from a PDE solver.
+@}
+
 @defgroup Models Models
 
-@defgroup LocalOperators Local operators
+@defgroup LocalOperators Local Operators
 @{
   Local operators are in many senses the heart of dorie; in them, we transform
   finite element basis into a vector of residuals and a mass matrix taking into

dune/dorie/common/h5file.hh

@@ -100,9 +100,11 @@ public:
      *  \param[in] data_type The H5 internal data type of the array to read.
      *  \param[out] data The vector for reading the data into. Is resized
      *      automatically.
-     *  \param[out] shape The vector containing the extensions of the dataset.
-     *      This shape is the reverse information of the H5 internal shape.
-     *      It indicated directions x, y[, z].
+     *  \param[out] shape The vector containing the shape of the dataset.
+     *      It indicates directions [z ,] y, x.
+     *
+     *  \todo Check if ``ext_t`` is convertible to ``hsize_t``. We currently
+     *      do a copy of values with implicit coversion.
      */
     template<typename data_t, typename ext_t>
     void read_dataset(const std::string& dataset_path,
@@ -147,9 +149,9 @@ public:
                                                 std::multiplies<hsize_t>());
         const std::vector<hsize_t> offset(arr_dim, 0);
 
-        // reverse order of dimensions!
-        shape.resize(arr_dim);
-        std::copy(dims.rbegin(), dims.rend(), shape.begin());
+        // FIXME: Converting unsigned to signed?
+        shape.resize(dims.size());
+        std::copy(begin(dims), end(dims), begin(shape));
 
         // create the memory space
         hid_t memspace_id = H5Screate_simple(arr_dim, dims.data(), NULL);

dune/dorie/dorie.cc

@@ -28,17 +28,6 @@
 // Main program with grid setup
 //===============================================================
 
-/**
-* \mainpage DORiE \n
-* __Dune Operated Richards Equation Solving Environment__ \n
-* _by Lukas Riedel, Felix Riexinger, Dion Häfner_ \n
-* \n
-* Easy access links for the most important documentation pages:
-* \see main Main Program Function
-* \see Dune::Dorie::FlowBoundary Class handling Boundary Condition query functions
-* \see Dune::Dorie::FlowSource Class handling Source Term query functions 
-*/
-
 template<typename Traits>
 using Sim = Dune::Dorie::RichardsSimulation<Traits>;
 
@@ -54,21 +43,7 @@ template<int dim, int order>
 using CubeAdaptive = Dune::Dorie::RichardsSimulationTraits<Dune::Dorie::BaseTraits<Dune::UGGrid<dim>,
 	Dune::GeometryType::BasicType::cube>,order>;
 
-/// Main Program Function: Initialize parameter file reader, build grid and call Richards Solver.
-/** As simplex and rectangular grids demand different FiniteElementMaps, the program calls different functions for these tasks.
- *  The objects and types are then passed to the generic solver function.\n
- *  The UG Grid and FEM templates need the dimension as _constexpr_ at compile time,
- *  so we need the tedious dim and FEorder queries.
- *  \see RichardsSolver Main functions for assembling operators and solving the problem
- *  \see RichardsSolverSimplex Helper for building simplex Grid Function Spaces
- *  \see RichardsSolverRectangular Helper for building rectangular Grid Function Spaces
- *  \see Dune::Dorie::Traits Variable type definitions
- *  \check Parameter file is specified
- *  \check Output directory is writable
- *  \check Finite Element Order is supported
- *  \check Grid type is supported
- *  \check Dimensions are supported
- */
+/// Main Program Function: Read config file, build grid and call Richards Solver
 int main(int argc, char** argv)
 {
   try{

dune/dorie/mainpage.dox

+/**
+
+@mainpage DORiE Doxygen Documentation
+
+Dune Operated Richards Equation Solving Environment
+
+For authors and license information, please refer to the license file in the
+top-level directory.
+
+@section Introduction
+
+DORiE implements a solver for Richards Equation based on DUNE and DUNE-PDELab.
+
+This is the developers' documentation of the source code. Refer to the user
+documentation built by `make doc` for detailed information on how to use the
+complete program.
+
+A good starting point for diving into the code is the Modules page.
+You find it in the navigation bar above.
+
+*/

dune/dorie/model/richards/parameterization/richards.hh → dune/dorie/model/richards/flow_parameters.hh

@@ -15,24 +15,37 @@
 
 #include <dune/dorie/model/richards/parameterization/interface.hh>
 #include <dune/dorie/model/richards/parameterization/mualem_van_genuchten.hh>
+#include <dune/dorie/model/richards/parameterization/scaling.hh>
 
 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
+/// Top-level parameterization interface for the local operator.
+/** This class loads and 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:
+ *
+ *  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 
+ *      - 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
+ *
+ *  \warning Load-balancing the grid *after* this object has been instantiated
+ *      is currently not supported.
+ *
+ *  \todo Make this object support load-balancing. This would require a global
+ *      ID mapper and communicating the data in the _param storage. The latter
+ *      is rather involved because of non-standard data types like maps and
+ *      Dorie::ScalingFactors
+ *
+ *  \author Lukas Riedel
+ *  \date 2018
+ *  \ingroup RichardsParam
  */
 template <typename Traits>
 class FlowParameters
@@ -51,16 +64,14 @@ private:
     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;
+    /// Struct for storing scaling factors
+    using Scaling = Dorie::ScalingFactors<RF>;
+    /// Base class for scaling adapters
+    using ScaleAdapter = Dorie::ScalingAdapter<Traits>;
 
     /// Storage for parameters and skaling factors
     using ParameterStorage = std::vector<
-        std::tuple<std::shared_ptr<Parameterization>,
-                   ScalingType,
-                   ScalingFactor
-        >
+        std::tuple<std::shared_ptr<Parameterization>, Scaling>
     >;
 
     /// Need this gruesome typedef because 'map::value_type' has const key
@@ -91,7 +102,7 @@ private:
     }
 
 public:
-    /// Constructor.
+    /// Create this object and load the data.
     /** Create a level grid view of level 0, create a mapper on it, and store
      *  the config tree.
      *  \param config Configuration file tree
@@ -129,6 +140,9 @@ public:
     template <class Entity>
     void bind (Entity entity) const
     {
+        static_assert(Entity::codimension == 0,
+                      "Parameters are mapped to codim 0 entities only!");
+
         // retrieve index of top father element of chosen entity
         while(entity.hasFather()) {
             entity = entity.father();
@@ -147,19 +161,13 @@ public:
         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;
 
-        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;
-            };
-        }
+        const auto& xi_cond = std::get<Scaling>(_cache).scale_cond;
 
-        return [cond_f](const RF saturation){
-            return cond_f(Saturation{saturation}).value;
+        return [cond_f, &xi_cond](const RF saturation){
+            return cond_f(Saturation{saturation}).value * xi_cond * xi_cond;
         };
     }
 
@@ -171,21 +179,15 @@ public:
     std::function<RF(RF)> saturation_f () const
     {
         verify_cache();
-        const auto& par = 
+        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;
-            };
-        }
+        const auto& scale_head = std::get<Scaling>(_cache).scale_head;
 
-        return [sat_f](const RF matric_head){
-            return sat_f(MatricHead{matric_head}).value;
+        return [sat_f, &scale_head](const RF matric_head){
+            return sat_f(MatricHead{matric_head * scale_head}).value;
         };
     }
 
@@ -199,9 +201,12 @@ public:
         verify_cache();
         const auto& par =
             std::get<std::shared_ptr<Parameterization>>(_cache);
-        const auto wc_f = par->water_content_f();
-        using Saturation = typename Parameterization::Saturation;
 
+        // get water content function and apply the scaling
+        const auto& scale_por = std::get<Scaling>(_cache).scale_por;
+        const auto wc_f = par->water_content_f(scale_por);
+
+        using Saturation = typename Parameterization::Saturation;
         return [wc_f](const RF saturation) {
             return wc_f(Saturation{saturation}).value;
         };
@@ -212,17 +217,32 @@ private:
     /// Build the parameterization from an element mapping and the input file.
     void build_parameterization (const std::vector<int>& element_index_map)
     {
-        const auto parameterization_map = read_parameter_file();
-        _param.reserve(element_index_map.size());
+        // Open the YAML file
+        const auto param_file_name = _config["parameters.file"];
+        std::cout << "Opening file " << param_file_name << std::endl;
+        YAML::Node param_file = YAML::LoadFile(param_file_name);
 
-        // build dummy scaling
-        ScalingType s_type{Parameterization::Scaling::None};
-        ScalingFactor s_factor{1.0};
+        // create the parameterization data
+        const auto parameterization_map = read_parameter_file(param_file);
+        _param.resize(element_index_map.size());
 
-        // 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));
+        // build global scaling
+        using SAF = Dorie::ScalingAdapterFactory<Traits>;
+        auto scale_cfg = param_file["scaling"];
+        auto scaling_adapter = SAF::create(scale_cfg["type"].as<std::string>(),
+                                           scale_cfg["data"]);
+
+        // insert parameterization and global scaling for each element
+        for (auto&& elem : elements(_gv))
+        {
+            // evaluate element index and position
+            const auto index = _mapper.index(elem);
+            const auto pos = elem.geometry().center();
+
+            // read values and insert
+            auto p = parameterization_map.at(element_index_map.at(index));
+            const auto scale = scaling_adapter->evaluate(pos);
+            _param.at(index) = std::make_tuple(p, scale);
         }
 
         // check that mapper can map to parameterization
@@ -230,13 +250,12 @@ private:
     }
 
     /// Read the YAML parameter file and insert information into a map.
-    std::map<int, std::shared_ptr<Parameterization>> read_parameter_file ()
-        const
+    /** \param param_file The top YAML node to read from (file root)
+     */
+    std::map<int, std::shared_ptr<Parameterization>> read_parameter_file (
+        const YAML::Node& param_file
+    )
     {
-        const auto param_file_name = _config.get<std::string>("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<int, std::shared_ptr<Parameterization>> ret;