Merge branch 'water-content-as-initial-condition' into...

Merge branch 'water-content-as-initial-condition' into 127-add-data-assimilation-interface-to-richardssimulation

doc/manual/initial.rst

@@ -31,16 +31,16 @@ They are controlled by the ``initial.type`` key and available for every model.
 
     * ``type = analytic``
 
-    An analytic function :math:`f(\vec{p})` which depends on the physical
-    position :math:`\vec{p}`. The function must be defined via the key
+    An analytic function :math:`f(\mathbf{p})` which depends on the physical
+    position :math:`\mathbf{p}`. The function must be defined via the key
     ``initial.equation``. For parsing the input expression, we use muparser_
     which supports a set of common mathematical functions. Additionally, the
     following variables can be used:
 
     Available variables:
-      * ``x``: X-coordinate :math:`p_1 \, [\mathrm{m}]`.
-      * ``y``: Y-coordinate :math:`p_2 \, [\mathrm{m}]`.
-      * ``z``: Z-coordinate :math:`p_3 \, [\mathrm{m}]` (only in 3D).
+      * ``x``: X-coordinate :math:`p_1 \, [\text{m}]`.
+      * ``y``: Y-coordinate :math:`p_2 \, [\text{m}]`.
+      * ``z``: Z-coordinate :math:`p_3 \, [\text{m}]` (only in 3D).
       * ``h``: Height above origin. Synonymous to ``y`` in 2D and ``z`` in 3D.
       * ``pi``: Mathematical constant :math:`\pi`.
       * ``dim``: Number of physical dimensions.
@@ -50,8 +50,9 @@ They are controlled by the ``initial.type`` key and available for every model.
        :ref:`initial-transformation`), typical initial conditions for a
        Richards simulation are
 
-       * Hydrostatic equilibrium: A vertical gradient of :math:`-1` and a
-         fixed value ``<v>`` at height :math:`h = 0 \, \mathrm{m}`::
+       * Hydrostatic equilibrium: A vertical gradient of
+         :math:`\partial_h h_m = -1` and a fixed value ``<v>`` at height
+         :math:`h = 0 \, \text{m}`:
 
             initial.equation = -h + <v>
 
@@ -59,7 +60,7 @@ They are controlled by the ``initial.type`` key and available for every model.
 
     .. tip::
        The expression for a gaussian pulse of solute concentration centered at
-       :math:`\vec{p} = [0.5, 0.5]^T \, \mathrm{m}` is::
+       :math:`\mathbf{p} = [0.5, 0.5]^T \, \text{m}` is::
 
           initial.equation = exp(-sqrt((x-0.5)^2+(y-0.5)^2)/(4.*0.002))/(4*pi*0.002)^(2/dim).
 
@@ -97,7 +98,22 @@ Initial condition tranformations for the Richards solver.
     * ``quantity = matricHead``
 
     The input data is directly interpreted as matric head
-    :math:`f = h_m [\mathrm{m}]`.
+    :math:`f = h_m \, [\text{m}]`.
+
+.. object:: Water Content to Matric Head
+
+    * ``quantity = waterContent``
+
+    The input data is interpreted as water content,
+    :math:`f = \theta_w \, [\text{-}]`, and transformed into matric head via
+    the :doc:`parameterization <parameter-file>` of the medium.
+
+    Values greater than the porosity :math:`\phi` and less than the residual
+    water content :math:`\theta_r` are automatically clamped to fit the allowed
+    range. Additionally, any input value :math:`f(x_0)` at some position
+    :math:`x_0` on the grid will result in a saturation greater zero,
+    :math:`\Theta (x_0) > 0`, to avoid divergence of the matric head towards
+    negative infinity.
 
 Transport
 ^^^^^^^^^
@@ -108,7 +124,7 @@ Initial condition tranformations for the Transport solver.
     * ``quantity = soluteConcentration``
 
     The input data is directly interpreted as solute concentration,
-    :math:`f = c_w [\mathrm{kg}/\mathrm{m}^3]`.
+    :math:`f = c_w [\text{kg}/\text{m}^3]`.
 
 .. _H5: https://www.h5py.org/
 .. _muparser: http://beltoforion.de/article.php?a=muparser&p=features

dune/dorie/common/initial_condition/factory.hh

@@ -34,7 +34,7 @@ protected:
     std::shared_ptr<spdlog::logger> log=get_logger(log_base)
   )
   {
-    std::unique_ptr<InitialCondition<T>> ic;
+    std::shared_ptr<InitialCondition<T>> ic;
 
     auto ic_type = ini_file.get<std::string>("initial.type");
     log->debug("Creating initial condition of type: {}", ic_type);
@@ -56,7 +56,7 @@ protected:
       const auto file_type = ic_datafile.substr(ext_start + 1);
       if (file_type == "h5") {
         using ICH5 = InitialConditionH5<T>;
-        ic = std::make_unique<ICH5>(grid_view, ini_file, log);
+        ic = std::make_shared<ICH5>(grid_view, ini_file, log);
 
       } else if (file_type == "vtu") {
         DUNE_THROW(NotImplemented,
@@ -73,7 +73,7 @@ protected:
       auto equation_str = ini_file.get<std::string>("initial.equation");
 
       using ICAnalytic = InitialConditionAnalytic<T>;
-      ic = std::make_unique<ICAnalytic>(log,grid_view,equation_str);
+      ic = std::make_shared<ICAnalytic>(log,grid_view,equation_str);
     }
     else {
       log->error("Unsupported initial condition type: {}", ic_type);

dune/dorie/model/richards/adapters/matric_head.hh

@@ -14,7 +14,7 @@ template<typename T, typename P, typename GF>
 class MatricHeadAdapter :
     public Dune::PDELab::GridFunctionBase<
             Dune::PDELab::GridFunctionTraits<typename T::GridView,
-                                             typename T::RangeField,
+                                             typename T::RF,
                                              1,
                                              typename T::Scalar>,
             MatricHeadAdapter<T, GF, P>
@@ -22,7 +22,7 @@ class MatricHeadAdapter :
 {
 public:
     using Traits = Dune::PDELab::GridFunctionTraits<typename T::GridView,
-                                                    typename T::RangeField,
+                                                    typename T::RF,
                                                     1,
                                                     typename T::Scalar>;
 

dune/dorie/model/richards/initial_condition.hh

@@ -7,28 +7,70 @@
 #include <dune/common/exceptions.hh>
 
 #include <dune/dorie/common/initial_condition/factory.hh>
+#include <dune/dorie/model/richards/adapters/matric_head.hh>
 
 namespace Dune {
 namespace Dorie {
 
+/// Transforms water content into matric head, serving as initial condition
+template<typename Traits, typename Param>
+class WaterContentInitialCondition :
+    public InitialConditionTransformation<Traits>
+{
+public:
+  using ICTransf = InitialConditionTransformation<Traits>;
+  using GFTraits = typename ICTransf::Traits;
+
+  using Base = InitialCondition<Traits>;
+
+protected:
+  using GV = typename ICTransf::GV;
+  using MHA = MatricHeadAdapter<Traits, Param, Base>;
+
+public:
+
+    /// Construct the transformation from a grid view and an initial condition
+    WaterContentInitialCondition (const Dune::ParameterTree& config,
+                                  const GV& grid_view,
+                                  std::shared_ptr<Base> ic,
+                                  const std::shared_ptr<const Param> param ):
+        ICTransf(grid_view, ic),
+        _adapter(std::make_shared<MHA>(config, grid_view, ic, param))
+    { }
+
+    void evaluate ( const typename GFTraits::ElementType& e,
+                    const typename GFTraits::DomainType& x,
+                    typename GFTraits::RangeType& y) const override
+    {
+        _adapter->evaluate(e, x, y);
+    };
+
+private:
+    std::shared_ptr<MHA> _adapter;
+};
+
+
 /// An initial condition factory for the Richards solver
 /**
- *  \tparam T Traits
+ *  \tparam Traits
+ *  \tparam Param
+ *
  *  \ingroup InitialConditions
  *  \author Lukas Riedel
  *  \date 2019
  */
-template<typename T>
-class RichardsInitialConditionFactory : public InitialConditionFactory<T>
+template<typename Traits, typename Param>
+class RichardsInitialConditionFactory : public InitialConditionFactory<Traits>
 {
 private:
-    using Base = InitialConditionFactory<T>;
+    using Base = InitialConditionFactory<Traits>;
 
 public:
     /// Create an initial condition for the Richards solver
-    static std::unique_ptr<InitialCondition<T>> create (
+    static std::shared_ptr<InitialCondition<Traits>> create (
         const Dune::ParameterTree& config,
-        const typename T::GridView& grid_view,
+        const typename Traits::GridView& grid_view,
+        const std::shared_ptr<const Param> param,
         const std::shared_ptr<spdlog::logger> log=get_logger(log_richards)
     )
     {
@@ -37,6 +79,14 @@ public:
         if (quantity == "matricHead") {
             return Base::create(config, grid_view, log);
         }
+        else if (quantity == "waterContent") {
+            // create initial condition grid function from input data
+            auto ic_gf = Base::create(config, grid_view, log);
+
+            // plug into adapter
+            using WCIC = WaterContentInitialCondition<Traits, Param>;
+            return std::make_shared<WCIC>(config, grid_view, ic_gf, param);
+        }
         else {
             log->error("Unsupported quantity for initial condition: {}",
                        quantity);

dune/dorie/model/richards/richards.cc

@@ -53,7 +53,7 @@ RichardsSimulation<Traits>::RichardsSimulation (
 	// --- Operator Helper Classes ---
 	fboundary = std::make_shared<const FlowBoundary>(inifile);
 	fsource = std::make_shared<const FlowSource>(inifile);
-	finitial = FlowInitialFactory::create(inifile, gv, this->_log);
+	finitial = FlowInitialFactory::create(inifile, gv, fparam, this->_log);
 
 	// --- Local Operators ---
 #ifdef EXPERIMENTAL_DG_FEATURES

dune/dorie/model/richards/richards.hh

@@ -86,7 +86,9 @@ struct RichardsSimulationTraits : public BaseTraits
 	/// Class defining the initial condition
 	using FlowInitial = Dune::Dorie::InitialCondition<BaseTraits>;
 	/// Class defining the initial condition factory
-	using FlowInitialFactory = Dune::Dorie::RichardsInitialConditionFactory<BaseTraits>;
+	using FlowInitialFactory
+		= Dune::Dorie::RichardsInitialConditionFactory<
+			BaseTraits, FlowParameters>;
 	/// Local spatial operator
 	using SLOP = Dune::Dorie::Operator::RichardsDGSpatialOperator<BaseTraits,FlowParameters,FlowBoundary,FlowSource,typename GFSHelper::FEM,false>;
 	/// Local temporal operator
@@ -296,7 +298,7 @@ protected:
 	std::shared_ptr<FlowParameters> fparam;
 	std::shared_ptr<const FlowBoundary> fboundary;
 	std::shared_ptr<const FlowSource> fsource;
-	std::unique_ptr<FlowInitial> finitial;
+	std::shared_ptr<FlowInitial> finitial;
 	std::unique_ptr<CalculationController> controller;
 
 	std::unique_ptr<SLOP> slop;

dune/dorie/model/transport/initial_condition.hh

@@ -26,7 +26,7 @@ private:
 
 public:
     /// Create an initial condition for the Richards solver
-    static std::unique_ptr<InitialCondition<T>> create (
+    static std::shared_ptr<InitialCondition<T>> create (
         const Dune::ParameterTree& config,
         const typename T::GridView& grid_view,
         const std::shared_ptr<spdlog::logger> log=get_logger(log_transport)

dune/dorie/model/transport/transport.hh

@@ -290,7 +290,7 @@ protected:
 
   std::shared_ptr<SoluteParameters> sparam;
   std::shared_ptr<SoluteBoundary> sboundary;
-  std::unique_ptr<SoluteInitial> sinitial;
+  std::shared_ptr<SoluteInitial> sinitial;
   std::unique_ptr<CalculationController> controller;
 
   std::unique_ptr<SLOP> slop;

dune/dorie/test/CMakeLists.txt

@@ -45,12 +45,19 @@ add_custom_target(test_grid_creator
 )
 add_dependencies(test_grid_creator prepare_testing)
 
-# initial condition test
+# initial condition tests
 dorie_add_metaini_test(
 	UNIT_TEST
-	SOURCE test-initial-condition.cc
-	BASENAME test-initial-condition
-	METAINI test-initial-condition.mini.in
+	SOURCE test-initial-condition-richards.cc
+	BASENAME test-initial-condition-richards
+	METAINI test-initial-condition-richards.mini.in
+	CREATED_TARGETS ic_target
+)
+dorie_add_metaini_test(
+	UNIT_TEST
+	SOURCE test-initial-condition-transport.cc
+	BASENAME test-initial-condition-transport
+	METAINI test-initial-condition-transport.mini.in
 	CREATED_TARGETS ic_target
 )
 

dune/dorie/test/test-initial-condition-richards.cc

+#ifdef HAVE_CONFIG_H
+    #include "config.h"
+#endif
+
+#include <dune/dorie/common/setup.hh>
+#include <dune/dorie/common/grid_creator.hh>
+
+#include <dune/dorie/model/richards/flow_parameters.hh>
+#include <dune/dorie/model/richards/initial_condition.hh>
+
+#include "test-initial-condition.hh"
+
+/// Test the basic initial conditions (without transformation)
+/** \tparam T Traits
+ *  \tparam IFC InitialConditionFactory
+ */
+template<typename T, typename ICF, typename Param>
+void test_ic_base (const std::shared_ptr<typename T::Grid> grid,
+                   const Dune::ParameterTree config,
+                   const std::shared_ptr<Param> param)
+{
+    // create the initial condition
+    const auto ic = ICF::create(config, grid->leafGridView(), param);
+
+    const auto type = config.get<std::string>("initial.type");
+    if (type == "analytic") {
+        // define testing function
+        const auto equation = config.get<std::string>("initial.equation");
+        assert(equation == "-1+2*y");
+        auto test_f = [](const auto e, const auto x){
+            using Range = typename T::Scalar;
+            constexpr int dim = T::Grid::dimension;
+            const auto x_global = e.geometry().global(x);
+            Range y = -1.0 + x_global[dim-1] * 2.0;
+            return y;
+        };
+        test_ic_against_function(grid, ic, test_f);
+    }
+    else if (type == "file") {
+        const auto dataset = config.get<std::string>("initial.dataset");
+
+        // define testing function
+        auto test_f = [&dataset](const auto e, const auto x){
+            typename T::Scalar y;
+            if (dataset == "twos") // case: matric head input
+                y = 2.0;
+            else if (dataset == "ones") // case: water content input
+                y = 0.0;
+
+            return y;
+        };
+        test_ic_against_function(grid, ic, test_f);
+    }
+}
+
+int main (int argc, char** argv)
+{
+    try{
+        // Initialize ALL the things!
+		auto [inifile, log, helper] = Dune::Dorie::Setup::init(argc, argv);
+
+        // create the grid
+        Dune::Dorie::GridCreator<Traits<2>::Grid> gc(inifile, helper);
+        auto grid = gc.grid();
+        const auto index_map = gc.element_index_map();
+
+        // fetch config for Richards
+        auto config = Dune::Dorie::Setup::prep_ini_for_richards(inifile);
+
+        // create the Richards logger
+        const auto log_level = config.get<std::string>("output.logLevel");
+        auto _log = Dune::Dorie::create_logger(
+            Dune::Dorie::log_richards,
+            helper,
+            spdlog::level::from_str(log_level));
+
+        // create Richards parameterization
+        using Param = Dune::Dorie::FlowParameters<Traits<2>>;
+        auto param = std::make_shared<Param>(config, grid, index_map);
+
+        // initial condition factory
+        using T = Traits<2>;
+        using ICF = Dune::Dorie::RichardsInitialConditionFactory<T, Param>;
+
+        const auto quantity = config.get<std::string>("initial.quantity");
+        if (quantity == "matricHead" || quantity == "waterContent") {
+            test_ic_base<T, ICF>(grid, config, param);
+        }
+        else {
+            log->error("Unsupported quantity for initial condition test: {}",
+                        quantity);
+            DUNE_THROW(Dune::NotImplemented, "Unsupported initial condition "
+                                             "quantity");
+        }
+
+        return 0;
+    }
+
+    catch (Dune::Exception &e) {
+        std::cerr << "Dune reported error: " << e << std::endl;
+        return 1;
+    }
+    catch (std::exception &e) {
+        std::cerr << "Exception thrown: " << e.what() << std::endl;
+        return 1;
+    }
+    catch (...) {
+        std::cerr << "Unknown exception!" << std::endl;
+        return 1;
+    }
+}

dune/dorie/test/test-initial-condition-richards.mini.in

+include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini
+
+_initial_type = analytic, data | expand
+_quantity = matricHead, waterContent | expand quantity_type
+{_quantity} == waterContent and {_initial_type} == analytic | exclude
+
+__name = ic-richards-{_initial_type}-{_quantity}
+_asset_path = "${PROJECT_SOURCE_DIR}/test"
+
+[grid]
+dimensions = 2
+initialLevel = 1
+gridType = rectangular
+cells = 50 50
+extensions = 1 1
+
+[grid.mapping]
+file = none
+
+[richards.initial]
+type = {_initial_type}
+quantity = {_quantity}
+
+file = ${CMAKE_CURRENT_LIST_DIR}/scaling.h5
+dataset = twos
+interpolation = nearest
+
+equation = -1+2*y
+
+[richards.parameters]
+file = test-parameterization-unscaled.yml

dune/dorie/test/test-initial-condition.cc → dune/dorie/test/test-initial-condition-transport.cc

@@ -2,53 +2,12 @@
     #include "config.h"
 #endif
 
-#include <dune/common/fvector.hh>
-
-#include <dune/grid/yaspgrid.hh>
-#include <dune/grid/common/gridview.hh>
-
 #include <dune/dorie/common/setup.hh>
 #include <dune/dorie/common/grid_creator.hh>
-#include <dune/common/float_cmp.hh>
 
-#include <dune/dorie/model/richards/initial_condition.hh>
 #include <dune/dorie/model/transport/initial_condition.hh>
 
-/// Traits stub used for this test
-template<int d>
-struct Traits
-{
-    static constexpr int dim = d;
-    using Grid = Dune::YaspGrid<dim>;
-    using GridView = typename Grid::LeafGridView;
-    using GV = GridView;
-    using DF = typename Grid::ctype;
-    using Domain = Dune::FieldVector<DF, dim>;
-
-    using RF = double;
-    using Scalar = Dune::FieldVector<RF, 1>;
-};
-
-/// Test an initial condition against a lambda function on the grid
-template<typename Grid, class IC, class Function>
-void test_ic_against_function (
-    const std::shared_ptr<Grid> grid,
-    const std::unique_ptr<IC>& ic,
-    const Function func
-)
-{
-    using Range = typename IC::Traits::RangeType;
-    using namespace Dune::FloatCmp;
-
-    // iterate over the grid vertices
-    for (auto&& e : Dune::elements(grid->leafGridView()))
-    {
-        Range y;
-        const auto x = e.geometry().local(e.geometry().center());
-        ic->evaluate(e, x, y);
-        assert(eq(func(e, x), y));
-    }
-}
+#include "test-initial-condition.hh"
 
 /// Test the basic initial conditions (without transformation)
 /** \tparam T Traits
@@ -97,53 +56,26 @@ int main (int argc, char** argv)
         auto grid = gc.grid();
         const auto index_map = gc.element_index_map();
 
-        {
-            auto config = Dune::Dorie::Setup::prep_ini_for_richards(inifile);
+        auto config = Dune::Dorie::Setup::prep_ini_for_transport(inifile);
 
-            // initial condition factory
-            using T = Traits<2>;
-            using ICF = Dune::Dorie::RichardsInitialConditionFactory<T>;
+        // initial condition factory
+        using T = Traits<2>;
+        using ICF = Dune::Dorie::TransportInitialConditionFactory<T>;
 
-            // create the Richards logger
-            const auto log_level = config.get<std::string>("output.logLevel");
-            auto _log = Dune::Dorie::create_logger(Dune::Dorie::log_richards,
-                                    helper,spdlog::level::from_str(log_level));
+        // create the transport logger
+        const auto log_level = config.get<std::string>("output.logLevel"); // FIXME
+        auto _log = Dune::Dorie::create_logger(Dune::Dorie::log_transport,
+                                helper/*,spdlog::level::from_str(log_level)*/);
 
-            const auto quantity = config.get<std::string>("initial.quantity");
-            if (quantity == "matricHead") {
-                test_ic_base<T,ICF>(grid, config);
-            }
-            else {
-                log->error("Unsupported quantity for initial condition test: {}",
-                           quantity);
-                DUNE_THROW(Dune::NotImplemented, "Unsupported initial condition "
-                                                 "quantity");
-            }
+        const auto quantity = config.get<std::string>("initial.quantity");
+        if (quantity == "soluteConcentration") {
+            test_ic_base<T, ICF>(grid, config);
         }
-
-        {
-            // auto config = Dune::Dorie::Setup::prep_ini_for_transport(inifile); // FIXME
-            auto config = inifile.sub("transport");
-
-            // initial condition factory
-            using T = Traits<2>;
-            using ICF = Dune::Dorie::TransportInitialConditionFactory<T>;
-
-            // create the transport logger
-            // const auto log_level = config.get<std::string>("output.logLevel"); // FIXME
-            auto _log = Dune::Dorie::create_logger(Dune::Dorie::log_transport,
-                                    helper/*,spdlog::level::from_str(log_level)*/);
-
-            const auto quantity = config.get<std::string>("initial.quantity");
-            if (quantity == "soluteConcentration") {
-                test_ic_base<T,ICF>(grid, config);
-            }
-            else {
-                log->error("Unsupported quantity for initial condition test: {}",
-                           quantity);
-                DUNE_THROW(Dune::NotImplemented, "Unsupported initial condition "
-                                                 "quantity");
-            }
+        else {
+            log->error("Unsupported quantity for initial condition test: {}",
+                       quantity);
+            DUNE_THROW(Dune::NotImplemented, "Unsupported initial condition "
+                                             "quantity");
         }
 
         return 0;

dune/dorie/test/test-initial-condition.mini.in → dune/dorie/test/test-initial-condition-transport.mini.in

 include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini
 
 _initial_type = analytic, data | expand
-__name = {_initial_type}
+__name = ic-transport-{_initial_type}
 _asset_path = "${PROJECT_SOURCE_DIR}/test"
 
 [grid]
@@ -14,16 +14,6 @@ extensions = 1 1
 [grid.mapping]
 file = none
 
-[richards.initial]
-type = {_initial_type}
-quantity = matricHead
-
-file = ${CMAKE_CURRENT_LIST_DIR}/scaling.h5
-dataset = twos
-interpolation = nearest
-
-equation = -1+2*y
-
 [transport.initial]
 type = {_initial_type}
 quantity = soluteConcentration

dune/dorie/test/test-initial-condition.hh

+#include <memory>
+#include <cassert>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/float_cmp.hh>
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/common/gridview.hh>
+
+/// Traits stub used for this test
+template<int d>
+struct Traits
+{
+    static constexpr int dim = d;
+    using Grid = Dune::YaspGrid<dim>;
+    using GridView = typename Grid::LeafGridView;
+    using GV = GridView;
+    using DF = typename Grid::ctype;
+    using Domain = Dune::FieldVector<DF, dim>;
+
+    using RF = double;
+    using Scalar = Dune::FieldVector<RF, 1>;
+    using Vector = Domain;
+};
+
+/// Test an initial condition against a lambda function on the grid
+template<typename Grid, class IC, class Function>
+void test_ic_against_function (
+    const std::shared_ptr<Grid> grid,
+    const std::shared_ptr<IC>& ic,
+    const Function func
+)
+{
+    using Range = typename IC::Traits::RangeType;
+    using namespace Dune::FloatCmp;
+
+    // iterate over the grid vertices
+    for (auto&& e : elements(grid->leafGridView()))
+    {
+        Range y;
+        const auto x = e.geometry().local(e.geometry().center());
+        ic->evaluate(e, x, y);
+        assert(eq(func(e, x), y));