Merge branch 'master' into 127-add-data-assimilation-interface-to-richardssimulation

Adjustments to comply to data assimilation interface:
* Switch back to `std::map` for Richards parameterization
  (introduced in !141).
* Rename typedef for `Parameterization` to `ParameterizationType` to
  avoid clash with new namespace.

.gitlab-ci.yml

@@ -10,7 +10,7 @@ variables:
   BASE_IMAGE: dorie/dune-env
   # Use semantic versioning (not the version of DUNE) and bump according to
   # to whether changes are backwards-compatible or not.
-  IMAGE_VERSION: "1.1"
+  IMAGE_VERSION: "1.2"
   DUNE_ENV_IMAGE: ${BASE_IMAGE}:img-v${IMAGE_VERSION}
 
   CMAKE_FLAGS:
@@ -69,6 +69,7 @@ setup:dune-env-clang:
 prep:update-dune: &update
   <<: *setup
   stage: prep
+  allow_failure: true
   only:
     - master
     - tags

CHANGELOG.md

@@ -40,6 +40,8 @@
 * Coupling between transient models for water flow and solute transport !96
 * Initial conditions generated from H5 input data !130
 * Generic Python VTK file reader !143
+* Parameterizations for hydrodynamic dispersion in solute transport !141
+* Generic Python VTK file reader !143, !150
 
 ### Changed
 * `Simulation` is renamed `RichardsSimulation` and moved to 
@@ -143,6 +145,7 @@
 * Specifying scaling field `extensions` and `offset` is now optional !133
 * Generalized initial condition specification in config file !129
 * Structure and setup of Sphinx user docs !126
+* Switch to stable `dune-randomfield` release branch !151
 
 ### Fixed
 * Solver in `RichardsSimulation` was using the wrong time variable.

README.md

@@ -88,7 +88,7 @@ by CI tests.
 | [dune-functions](https://gitlab.dune-project.org/staging/dune-functions) | releases/2.6
 | [dune-typetree](https://gitlab.dune-project.org/staging/dune-typetree) | releases/2.6
 | [dune-pdelab](https://gitlab.dune-project.org/pdelab/dune-pdelab) | releases/2.6
-| [dune-randomfield](https://gitlab.dune-project.org/oklein/dune-randomfield) | master
+| [dune-randomfield](https://gitlab.dune-project.org/oklein/dune-randomfield) | releases/2.6
 
 
 #### Optional Packages
@@ -197,7 +197,8 @@ If you installed [Anaconda](https://conda.io/docs/user-guide/install/download.ht
     to the `CMAKE_FLAGS` in the above command, replacing `<hdf5-path>` with the
     path chosen as installation prefix when configuring HDF5.
 
-#### Experimental Features
+
+### Experimental Features
 The local operator implementing Richards equation's discretization supports
 multiple scheme settings. Setting these via the config file is disabled by
 default. You can enable this feature by reconfiguring DORiE with the CMake flag

cmake/modules/DorieTesting.cmake

@@ -7,7 +7,8 @@ add_custom_target(unit_tests
 # Add the system test build target and test command
 add_custom_target(build_system_tests)
 add_custom_target(system_tests
-    COMMAND ctest --output-on-failure --exclude-regex ^ut.+$
+    COMMAND ctest
+            --output-on-failure --exclude-regex "^\\(ut.+\\|python_tests\\)$"
 )
 
 # create the mapping datafile as preparation for all tests

doc/default_files/CMakeLists.txt

@@ -119,7 +119,7 @@ endfunction()
 # copy BC & parameter files
 file(COPY   2d_infiltr.bcdat 3d_infiltr.bcdat 
             2d_solute.bcdat 3d_solute.bcdat 
-            param.yml 
+            richards_param.yml  transport_param.yml
     DESTINATION .)
 
 # Random field generator

doc/default_files/common-parameters.xml

@@ -32,10 +32,15 @@ adding an empty line, make text **bold** or ``monospaced``.
 <dorie>
   <category name="simulation">
     <parameter name="mode">
-      <definition> Mode of simulation in DORiE. </definition>
+      <definition> Sets the simulation mode of DORiE. 
+        (``richards``) mode calculates the matric head with a
+        DG scheme and produce water fluxes with Raviart Thomas reconstruction.
+        (``richards+transport``) mode calculates (``richards``) mode and use 
+        the generated water flux and saturation at each step to calculate the 
+        solute transport model for unsaturated media with a DG scheme.
+        </definition>
       <suggestion> richards </suggestion>
       <values> richards, richards+transport </values>
-      <comment> richards, richards+transport </comment>
     </parameter>
   </category>
 

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/default_files/richards-parameters.xml

@@ -227,6 +227,7 @@ adding an empty line, make text **bold** or ``monospaced``.
       <definition> YAML file containing the parameter definitions.
       </definition>
       <values> path </values>
+      <suggestion> richards_param.yml </suggestion>
     </parameter>
   </category>
 

doc/default_files/richards_param.yml

+volumes:
+  sand:
+    index: 0
+    richards: 
+      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
+    richards: 
+      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/default_files/transport-parameters.xml

@@ -67,6 +67,16 @@ adding an empty line, make text **bold** or ``monospaced``.
       <values> endOfTransportStep, endOfRichardsStep, none </values>
     </parameter>
 
+    <parameter name="writeDispersionTensor">
+      <definition> Defines whether VTK files should include the hydrodynamic 
+       dispersion tensor. Tensors are written in 3D and have 9 componentents 
+       independently of the world dimension. This can be easily be visualizated
+       in Paraview with the ``Tensor Glyph`` filter.
+      </definition>
+      <values> true, false </values>
+      <suggestion> false </suggestion>
+    </parameter>
+
     <parameter name="subsamplingLevel">
       <definition> Plot VTK files with virtually refined grids. VTK only
         supports bilinear triangulations and displays higher-order solutions
@@ -109,6 +119,13 @@ adding an empty line, make text **bold** or ``monospaced``.
       <values> true, false </values>
       <suggestion> false </suggestion>
     </parameter>
+
+    <parameter name="dirichletMode">
+      <definition> Type of the input value for the dirichlet condition. 
+      </definition>
+      <values> soluteConcentration, totalSolute </values>
+      <suggestion> soluteConcentration </suggestion>
+    </parameter>
   </category>
 
   <category name="initial">
@@ -222,11 +239,12 @@ adding an empty line, make text **bold** or ``monospaced``.
     </parameter>
   </category>
 
-  <category name="parameters">
-    <parameter name="effHydrDips">
-      <definition> Effective hydrodynamic dispersion </definition>
-      <values> float &gt; 0 </values>
-      <suggestion> 2E-9 </suggestion>
+ <category name="parameters">
+    <parameter name="file">
+      <definition> YAML file containing the parameter definitions.
+      </definition>
+      <values> path </values>
+      <suggestion> transport_param.yml </suggestion>
     </parameter>
   </category>
 
@@ -253,8 +271,7 @@ adding an empty line, make text **bold** or ``monospaced``.
     </parameter>
 
     <parameter name="FEorder">
-      <definition> Order of the finite element method used. Values &gt; 1 are not
-        thoroughly tested. </definition>
+      <definition> Order of the finite element method used. </definition>
       <suggestion> 0 </suggestion>
       <values> 0 </values>
     </parameter>

doc/default_files/transport_param.yml

+volumes:
+  sand:
+    index: 0
+    transport: 
+      type: Deff_MQ1 + Dhm_iso
+      parameters:
+        char_length: 1.5E-11
+        mol_diff: 2.E-9
+        phi: 0.34
+        lambda_t: 0.0025
+        lambda_l: 0.025
+      
+  silt:
+    index: 1
+    transport: 
+      type: Deff_MQ1 + Dhm_iso
+      parameters:
+        char_length: 1.5E-11
+        mol_diff: 2.E-9
+        phi: 0.34
+        lambda_t: 0.0025
+        lambda_l: 0.025
\ No newline at end of file

doc/introduction/data-io.rst

@@ -64,11 +64,16 @@ Output Files
   .. object:: Transport
 
     - ``solute``: Solute concentration in water phase
-      :math:`c_w \, [\mathrm{kg}/\mathrm{m}^3]`
+      :math:`c_w \, [\mathrm{kg}\mathrm{m}^{-3}]`
     - ``solute_total``: Total solute concentration
-      :math:`c_t \, [\mathrm{kg}/\mathrm{m}^3]`
+      :math:`c_t \, [\mathrm{kg}\mathrm{m}^{-3}]`
+    - ``micro_peclet``: Microscopic peclet number
+      :math:`\mathsf{pe_\mu} \, [-]`
     - ``flux_RT<k-1>``: Reconstructed solute flux
-      :math:`\mathbf{j}_{s, \mathrm{rc}} \, [\mathrm{kg}/\mathrm{m}^{-2}\mathrm{s}^{-1}]`
+      :math:`\mathbf{j}_{s, \mathrm{rc}} \, [\mathrm{kg}\,\mathrm{m}^{-2}\mathrm{s}^{-1}]`
+      (if enabled!)
+    - ``eff_hd_dispersion``: Hydrodynamic dispersion tensor
+      :math:`\mathsf{D}_\mathsf{hd} \, [\mathrm{m}^2\mathrm{s}^{-1}]`
       (if enabled!)
 
 * VTK Parallel Collection Files (``.pvtu``): Merging multiple VTK files in case

doc/man-cookbook.rst

@@ -7,7 +7,7 @@ Cook Book
 
 Introduction
 ============
-This part of the documentation is intended for first-time DORiE users. It explaines the basic usage of the program, how to execute a simulation and how to analyze its results. It will also showcase more complex features like Adaptive Grid Refinement.
+This part of the documentation is intended for first-time DORiE users. It explains the basic usage of the program, how to execute a simulation and how to analyze its results. It will also showcase more complex features like Adaptive Grid Refinement.
 
 Prerequisites
 -------------
@@ -57,7 +57,7 @@ Let's begin with the ``output``. We want the program to give us at least some ou
 .. code-block:: none
 
 	[output]
-	verbose = 1
+	logLevel = info
 	outputPath = ./sand
 	fileName = sand
 

doc/python/utilities.rst

@@ -53,12 +53,17 @@ Generic VTK Reader
 ------------------
 
 This is a flexible VTK XML file reader which can evaluate all data sets written
-by DORiE.
+by DORiE. It is derived from the abstract base class ``Mapping`` and thus
+supports accessing the data arrays stored inside like a Python ``dict``. Use
+the data array names as keys.
 
 .. autoclass:: dorie.utilities.vtktools.vtkreader.VTKReader
-   :members:
 
-   .. automethod:: __init__
+The ``VTKReader`` returns instances of ``VTKDataArray``, which is a wrapper
+around the underlying VTK data array.
+
+.. autoclass:: dorie.utilities.vtktools.vtkreader.VTKDataArray
+   :members:
 
 dorie.utilities.check_path
 --------------------------

doc/requirements.txt

-Sphinx<2.0
+Sphinx
 sphinx_rtd_theme
 sphinx-argparse
 sphinx-markdown-tables

docker/dune-env.dockerfile

@@ -43,7 +43,7 @@ ENV LANG en_US.utf8
 
 WORKDIR /opt/dune
 RUN git clone https://gitlab.dune-project.org/staging/dune-uggrid.git -b releases/2.6 \
-    && git clone https://gitlab.dune-project.org/oklein/dune-randomfield.git -b master \
+    && git clone https://gitlab.dune-project.org/oklein/dune-randomfield.git -b releases/2.6 \
     && git clone https://gitlab.dune-project.org/core/dune-common.git -b releases/2.6 \
     && git clone https://gitlab.dune-project.org/core/dune-geometry.git -b releases/2.6 \
     && git clone https://gitlab.dune-project.org/core/dune-grid.git -b releases/2.6 \

dune/dorie/common/boundary_condition.hh

@@ -816,6 +816,8 @@ namespace Dorie{
 								std::cout << "outflow : ";
 							else if(bcData[i][j].type[k]==BoundaryCondition::Dirichlet)
 								std::cout << "dirichlet : ";
+							else if(bcData[i][j].type[k]==BoundaryCondition::Outflow)
+								std::cout << "outflow : ";
 						}
 						index++;
 					}

dune/dorie/common/grid_function_writer.hh

@@ -16,6 +16,72 @@
 namespace Dune{
 namespace Dorie{
 
+/**
+ * @brief      Converts PDELab grid function to match a dune-function.
+ * @note       This adapter differs from the one provided in PDELab in that this
+ *             one has the current interface of dune-grid to write vtk files.
+ *             The main difference using this adapter is that it allows to write
+ *             more than 3 components.
+ *
+ * @tparam     GF    PDELab grid function type.
+ */
+template<class GF>
+class VTKGridFunctionAdapter 
+{
+  using GridFunction = typename std::decay<GF>::type;
+  using GridView = typename GF::Traits::GridViewType;
+  using Entity = typename GridView::template Codim<0>::Entity;
+  using Range = typename GF::Traits::RangeType;
+  using Domain = typename GF::Traits::DomainType;
+public:
+
+  /**
+   * @brief      Constructs the object
+   *
+   * @param[in]  gf    pointer to a constant grid function
+   */
+  VTKGridFunctionAdapter(std::shared_ptr<const GF> gf)
+    : _gf(gf)
+    , _entity(nullptr)
+  {}
+
+  /**
+   * @brief      Binds an entity to the grid function
+   *
+   * @param[in]  e     Entity
+   */
+  void bind(const Entity& e)
+  {
+    _entity = &e;
+  }
+
+  /**
+   * @brief      Unbinds the previously binded entity
+   */
+  void unbind()
+  {
+    _entity = nullptr;
+  }
+
+  /**
+   * @brief      Evaluates the grid function at a given coordinate
+   *
+   * @param[in]  x     Local coordinate with respect to the last binded entity
+   */
+  Range operator()(const Domain& x) const
+  {
+    Range y;
+    assert(_entity);
+    _gf->evaluate(*_entity,x,y);
+    return y;
+  }
+
+private:
+  const Entity* _entity;
+  std::shared_ptr<const GF> _gf;
+};
+
+
 /*-------------------------------------------------------------------------*//**
  * @brief      Class for grid fucntion vtk sequence writer. This class works
  *             exacly the same way as a VTKSequence writer but it receives
@@ -39,13 +105,20 @@ public:
    * @tparam     GF      Type of the grid function
    */
   template<class GF>
-  void addVertexData(const std::shared_ptr<const GF> &gf_ptr, std::string name)
+  void addVertexData(std::shared_ptr<const GF> gf_ptr, std::string name)
   {
     static_assert(std::is_same<typename GF::Traits::GridViewType,GV>::value, 
       "GridFunctionVTKSequenceWriter only works with only one GridView type.");
 
-    auto vtk_gf_ptr = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<GF> >(gf_ptr, name);
-    Dune::VTKSequenceWriter<GV>::addVertexData(vtk_gf_ptr);
+    Dune::Dorie::VTKGridFunctionAdapter<GF> vtk_gf(gf_ptr);
+
+    const int vtk_ncomps = GF::Traits::dimRange;
+    const int dim = GF::Traits::dimDomain;
+    auto vtk_type = (vtk_ncomps == dim) ? VTK::FieldInfo::Type::vector 
+                                        : VTK::FieldInfo::Type::scalar;
+    auto vtk_info = VTK::FieldInfo(name,vtk_type,vtk_ncomps);
+
+    this->vtkWriter()->addVertexData(vtk_gf,vtk_info);
   }
 
   /*-----------------------------------------------------------------------*//**
@@ -62,14 +135,21 @@ public:
     static_assert(std::is_same<typename GF::Traits::GridViewType,GV>::value, 
       "GridFunctionVTKSequenceWriter only works with only one GridView type.");
     
-    auto vtk_gf_ptr = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<GF> >(gf_ptr, name);
-    Dune::VTKSequenceWriter<GV>::addCellData(vtk_gf_ptr);
+    Dune::Dorie::VTKGridFunctionAdapter<GF> vtk_gf(gf_ptr);
+
+    const int vtk_ncomps = GF::Traits::dimRange;
+    const int dim = GF::Traits::dimDomain;
+    auto vtk_type = (vtk_ncomps == dim) ? VTK::FieldInfo::Type::vector 
+                                        : VTK::FieldInfo::Type::scalar;
+    auto vtk_info = VTK::FieldInfo(name,vtk_type,vtk_ncomps);
+
+    this->vtkWriter()->addCellData(vtk_gf,vtk_info);
   }
 
   /*-----------------------------------------------------------------------*//**
-   * @brief      Clear all the attached VTKFunctions. (Necessary for each
-   *             iteration if VTKFunctions are managed internaly with pointers
-   *             instead of references)
+   * @brief      Clear all the attached VTKFunctions.
+   * @details    This is necessary for each iteration if VTKFunctions are
+   *             managed internaly with pointers instead of references
    */
   void clear()
   {

dune/dorie/common/initial_condition/analytic.hh

@@ -20,7 +20,7 @@ namespace Dorie {
  * @brief      Class for analytic initial conditions.
  *
  * @ingroup    InitialConditions
- * @author     Santiago Ospina
+ * @author     Santiago Ospina De Los Ríos
  * @date       2019
  *
  * @tparam     T     BaseTraits

dune/dorie/common/initial_condition/factory.hh

@@ -17,7 +17,7 @@ namespace Dorie{
  * @brief      Abstract Base Class for initial conditions.
  *
  * @ingroup    InitialConditions
- * @author     Santiago Ospina
+ * @author     Santiago Ospina De Los Ríos
  * @date       2019
  *
  * @tparam     T     BaseTraits

dune/dorie/common/initial_condition/h5.hh

@@ -18,7 +18,7 @@ namespace Dorie{
  * @brief      Class for initial conditions from hdf5 data files.
  *
  * @ingroup    InitialConditions
- * @author     Santiago Ospina
+ * @author     Santiago Ospina De Los Ríos
  * @date       2019
  *
  * @tparam     T     BaseTraits

dune/dorie/common/initial_condition/initial_condition.hh

@@ -13,7 +13,7 @@ namespace Dorie{
  * @brief      Abstract Base Class for initial conditions.
  *
  * @ingroup    InitialConditions
- * @author     Santiago Ospina
+ * @author     Santiago Ospina De Los Ríos
  * @date       2019
  *
  * @tparam     T     BaseTraits

dune/dorie/common/parameterization_factory.hh

+#ifndef DUNE_DORIE_PARAMETERIZATION_READER_HH
+#define DUNE_DORIE_PARAMETERIZATION_READER_HH
+
+#include <unordered_map>
+#include <memory>
+#include <string>
+#include <set>
+
+#include <yaml-cpp/yaml.h>
+
+#include <dune/common/exceptions.hh>
+
+#include <dune/dorie/common/utility.hh>
+#include <dune/dorie/common/logging.hh>
+
+namespace Dune{
+namespace Dorie{
+
+/*-------------------------------------------------------------------------*//**
+ * @brief      Interface and default reader for parameterization factories
+ *
+ * @tparam     Param  The parameter type to build
+ */
+template<class Param>
+struct ParameterizationFactory
+{
+  virtual ~ParameterizationFactory() = default;
+
+  /**
+   * @brief      Parameterization selector
+   * @details    Given a type_node in yaml, it selects a -polymorphic-
+   *             parameter.
+   *
+   * @param[in]  type_node  YAML node describing the type of parameterization to
+   *                        choose from.
+   * @param[in]  name       Name that describes the parametrized material.
+   *
+   * @return     Shared pointer to the selected parameter.
+   */
+  virtual std::shared_ptr<Param> selector(
+              const YAML::Node& type_node,
+              const std::string& name) const = 0;
+
+
+  /**
+   * @brief      Parameterization reader
+   * @details    It reads parameterizations and its parameters. For each volume,
+   *             it assings a index to a parameterization (choosen by the
+   *             parameterization selector) and later filled with the paramaters
+   *             in the yaml file.
+   *
+   * @param[in]  param_file  YAML node with information about type and
+   *                         parameters
+   * @param[in]  param_name  Name of the parameterization.
+   * @param[in]  log         Logger
+   *
+   * @return     (Unsorted) Map with parameterizations for each material.
+   */
+  std::unordered_map<int, std::shared_ptr<Param>> 
+  reader( const YAML::Node& param_file, 
+          const std::string& param_name,
+          const std::shared_ptr<spdlog::logger> log) const
+  {
+    // mapping: index on grid to parameterization object (will be returned)
+    std::unordered_map<int, std::shared_ptr<Param>> ret;
+
+    // set to check for duplicate indices
+    std::set<int> param_index_set;
+
+    // iterate over all volume nodes
+    for (auto&& node : param_file["volumes"])
+    {
+      const auto volume = node.first.as<std::string>();
+      const YAML::Node param_info = node.second;
+
+      // read relevant data
+      const auto index = param_info["index"].as<int>();
+      if (param_index_set.count(index) != 0) {
+          log->error("Index '{}' in volume '{}' of parameterization '{}'"
+                      "already exists.",
+                      index, volume, param_name);
+          DUNE_THROW(IOError, "Parameter file indices must be unique");
+      }
+      param_index_set.emplace(index);
+
+      // get parameterization object
+      const YAML::Node& type_node = param_info[param_name];
+      auto parameterization = this->selector(type_node, volume);
+      auto parameter_values = parameterization->parameters();
+
+      // fill parameterization with parameter entries
+      log->trace("  Reading parameters. Volume: {}, Index: {}",
+                  volume, index);
+
+      auto yaml_parameters = yaml_sequence_to_map(type_node["parameters"]);
+
+      for (auto&& [name, value] : parameter_values)
+      {
+        try {
+            value = yaml_parameters[name].template as<double>();
+        }
+        catch (...) {
+          log->error("    Parameter '{}' is missing in '{}' parameterization,"
+                     " volume {}, index {}.",
+                      name, param_name, volume, index);
+          DUNE_THROW(IOError, "Missing parameters!");
+        }
+      }
+
+      ret.emplace(index, parameterization);
+    }
+
+    return ret;
+  }
+};
+
+}
+}
+
+#endif // DUNE_DORIE_PARAMETERIZATION_READER_HH
\ No newline at end of file

dune/dorie/common/setup.hh

@@ -63,7 +63,7 @@ namespace Setup
 		const std::string inifilename = argv[1];
 
         // Read ini file
-		log->info("Reading parameter file: {}", inifilename);
+		log->info("Reading configuration file: {}", inifilename);
 		Dune::ParameterTree inifile;
 		Dune::ParameterTreeParser ptreeparser;
 		ptreeparser.readINITree(inifilename, inifile);

dune/dorie/common/typedefs.hh

@@ -13,7 +13,7 @@ namespace Dune{
 			enum Type {
 				Neumann, //!< Fixed flux at boundary
 				Dirichlet, //!< Fixed matric head at boundary
-				Outflow, //!< Fixed flux at boundary
+				Outflow, //!< Fixed outflow boundary condition
 				Other //!< UNUSED
 			};
 

dune/dorie/common/util.hh

@@ -195,16 +195,6 @@ public:
 	}
 };
 
-/// A wrapper for begin and end time stamps of a time interval
-/** \tparam TF Arithmetic data type for time stamps
- */
-template<typename TF>
-struct TimeInterval {
-  TF begin;
-  TF end;
-};
-
-
 /// Traits struct defining basic types based on grid an geometry
 template<class GridType, Dune::GeometryType::BasicType GeometryType>
 struct BaseTraits
@@ -214,7 +204,6 @@ struct BaseTraits
 
 	using TF = double;
 	using TimeField = TF;
-	using TimeInterval = Dune::Dorie::TimeInterval<TimeField>;
 
 	using RF = double;