Merge branch '142-add-parameterization-for-the-effective-hydrodynamic-dispersion-tensor'

.gitignore

@@ -11,3 +11,6 @@ test/maps/cell_ids.h5
 # Ignore auxiliary Python files
 __pycache__/
 *.egg-info/
+
+# Ignore Desktop-Service-Store files from MacOS
+**/.DS_Store
\ No newline at end of file

README.md

@@ -164,7 +164,35 @@ If you installed [Anaconda](https://conda.io/docs/user-guide/install/download.ht
 
     to the `CMAKE_FLAGS` in the call to `dunecontrol` above.
 
-#### Experimental Features
+### Richards and Transport builds
+DORiE can solve the solute transport equation together with the Richards equation. However, for users that do not require this feature, it is possible to build a standalone solver for the Richards equation by using the make taget `richards`. This will half the compilation time. Notice that the other way around is not possible: transport model always require the solution of the richards equation.
+
+| Target     | Dependencies  | Detail |
+| ---------- | ------------- | ------ |
+| `richards` | `richards-impl` | Executable for Richards models
+| `transport` | `richards-impl`, `transport-impl` | Executable for solute transport models
+| `richards-impl` |  | Internal Richards library
+| `transport-impl` |  | Internal solute transport library
+| `dorie`    | `richards`, `transport` | Full DORiE compilation 
+
+For doing so, dune modules have to be built in more steps. Make sure you have completed the first 5 steps in the step-by-step instructions, and follow the instructions below:
+
+6.1. Building all DUNE dependencies for DORiE. If you have installed `dune-testtools, add it to the list of modules.
+
+        CMAKE_FLAGS="-DDUNE_PYTHON_VIRTUALENV_SETUP=True -DDUNE_PYTHON_ALLOW_GET_PIP=True" ./dune-common/bin/dunecontrol --module=dune-pdelab,dune-uggrid,dune-randomfield all
+
+6.2 Configure DORiE:
+
+        CMAKE_FLAGS="-DDUNE_PYTHON_VIRTUALENV_SETUP=True -DDUNE_PYTHON_ALLOW_GET_PIP=True" ./dune-common/bin/dunecontrol --only=dorie configure
+
+6.3 Build a DORiE target, `<dorie-target>`:
+
+        CMAKE_FLAGS="-DDUNE_PYTHON_VIRTUALENV_SETUP=True -DDUNE_PYTHON_ALLOW_GET_PIP=True" ./dune-common/bin/dunecontrol --only=dorie make <dorie-target>
+
+### Parallel compilation
+DUNE, and hence, DORiE support compilation of targets in parallel using the environmental variable `MAKE_FLAGS="-j X"`, where `X` is the number of processes. Some care has to be taken when setting this number. We estimate that each process will use about 2GB of memory RAM. If the memory RAM is overflown during compilation, build process may take much longer due to use of swap memory. Hence, make sure take this into account when assigning a higher number of processors.
+
+### 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

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: superposition
+      eff_diff:
+        type: milligton_quirk_1
+      eff_hydromechanic_disp:
+        type: isotropic
+        long_disp:
+          type: const
+        trans_disp:
+          type: const
+      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: superposition
+      eff_diff:
+        type: milligton_quirk_1
+      eff_hydromechanic_disp:
+        type: isotropic
+        long_disp:
+          type: const
+        trans_disp:
+          type: const
+      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/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
 

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,43 @@
 namespace Dune{
 namespace Dorie{
 
+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:
+
+  VTKGridFunctionAdapter(std::shared_ptr<const GF> gf)
+    : _gf(gf)
+  {}
+
+  void bind(const Entity& e)
+  {
+    _entity = &e;
+  }
+
+  void unbind()
+  {
+    _entity = nullptr;
+  }
+
+  Range operator()(const Domain& x) const
+  {
+    Range y;
+    _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 +76,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,8 +106,15 @@ 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);
   }
 
   /*-----------------------------------------------------------------------*//**

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  { description }
+ */
+template<class Param>
+struct ParameterizationFactory
+{
+  virtual ~ParameterizationFactory() = default;
+  /**
+   * @brief      Parameterization selector
+   *
+   * @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
+   *
+   * @param[in]  param_file  YAML node with information about the type and the
+   *                         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 name = 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 material '{}' of parameterization '{}'"
+                      "already exists.",
+                      index, name, 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, name);
+      auto parameter_values = parameterization->parameters();
+
+      // fill parameterization with parameter entries
+      log->trace("  Reading parameters. Volume: {}, Index: {}",
+                  name, index);
+
+      log->trace("    Required parameters: {}", parameter_values.size());
+      for (auto p : parameter_values)
+        log->trace("      {}", p.first);
+
+      auto yaml_parameters = yaml_sequence_to_map(type_node["parameters"]);
+      log->trace("    Available parameters: {}",
+                  yaml_parameters.size());
+      unsigned param_used = 0;
+      for (auto&& yaml_param : yaml_parameters)
+      {
+        const auto p_name = yaml_param.first.as<std::string>();
+        const auto p_value = yaml_param.second.as<double>();
+        
+        // try to insert the value into the parameter map
+        auto [begin,end] = parameter_values.equal_range(p_name);
+        if (begin==end) { // mark unused parameter
+          log->trace("      - {}: {}", p_name, p_value);
+        } else {          // mark and set parameter
+          for (auto p = begin; p != end; ++p) {
+            p->second = p_value;
+            log->trace("      + {}: {}", p->first, p->second);
+            param_used++;
+          }
+        }
+      }
+
+      if (param_used!=parameter_values.size()) {
+        log->error("    Parameters used: {} of {}. "
+                    "Volume: {}, Index: {}.",
+                    param_used,parameter_values.size(),name, 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;
 	using RangeField = RF;

dune/dorie/common/utility.hh

@@ -8,6 +8,9 @@
 #include <memory>
 #include <utility>
 
+#include <yaml-cpp/yaml.h>
+
+#include <dune/common/exceptions.hh>
 #include <dune/common/fvector.hh>
 
 #include <dune/grid/common/gridview.hh>
@@ -113,6 +116,61 @@ inline bool is_none (const std::string& in)
     return false;
 }
 
+/*-------------------------------------------------------------------------*//**
+ * @brief      Transform a Yaml node with sequences to map
+ * @details    This function check all the nodes within 'in' and transform every
+ *             sequence into a map. In particular it transform vectors and
+ *             tensors into single key-value map pair by adding a suffix
+ *             indicating the component.
+ *
+ * @param[in]  in    Map Yaml node.
+ *
+ * @return     Map Yaml node with only key-double pairs.
+ */
+inline YAML::Node yaml_sequence_to_map(const YAML::Node& in)
+{
+  YAML::Node out;
+  for (auto&& node : in)
+  {
+    const auto key = node.first.as<std::string>();
+    const YAML::Node& value = node.second;
+
+    if (value.IsScalar()) {
+      out[key] = value;
+    } else if (value.IsSequence()) {
+      if (value.size() == 1) { // scalar
+        out[key] = value[0].as<double>();
+      } else if (value.size() == 2) { // assume vector in 2D
+        out[key+"_x"] = value[0].as<double>();
+        out[key+"_y"] = value[1].as<double>();
+      } else if (value.size() == 3) { // assume vector in 3D
+        out[key+"_x"] = value[0].as<double>();
+        out[key+"_y"] = value[1].as<double>();
+        out[key+"_z"] = value[2].as<double>();
+      } else if (value.size() == 4) { // assume tensor in 2D
+        out[key+"_xx"] = value[0].as<double>();
+        out[key+"_xy"] = value[1].as<double>();
+        out[key+"_yx"] = value[2].as<double>();
+        out[key+"_yy"] = value[3].as<double>();
+      } else if (value.size() == 9) { // assume tensor in 3D
+        out[key+"_xx"] = value[0].as<double>();
+        out[key+"_xy"] = value[1].as<double>();
+        out[key+"_xz"] = value[2].as<double>();
+        out[key+"_yx"] = value[3].as<double>();
+        out[key+"_yy"] = value[4].as<double>();
+        out[key+"_yz"] = value[5].as<double>();
+        out[key+"_zx"] = value[6].as<double>();
+        out[key+"_zy"] = value[7].as<double>();
+        out[key+"_zz"] = value[8].as<double>();
+      }
+    } else {
+      DUNE_THROW(Dune::IOError,
+        "Yaml sequence '" << key << "' is not convertible to map pair!");
+    }
+  }
+  return out;
+}
+
 } // namespace Dorie
 } // namespace Dune
 

dune/dorie/model/richards/parameterization/interface.hh

-#ifndef DUNE_DORIE_PARAM_RICHARDS_HH
-#define DUNE_DORIE_PARAM_RICHARDS_HH
+#ifndef DUNE_DORIE_PARAM_RICHARDS_INTERFACE_HH
+#define DUNE_DORIE_PARAM_RICHARDS_INTERFACE_HH
 
 #include <functional>
 #include <string>
@@ -7,6 +7,7 @@
 
 namespace Dune {
 namespace Dorie {