Merge branch '137-add-a-finite-volume-local-operator-for-the-richards-equation' into 'master'

Resolve "Add a Finite Volume local operator for the Richards equation"

Closes #137

See merge request !132

CHANGELOG.md

@@ -148,6 +148,7 @@
 * Structure and setup of Sphinx user docs !126
 * Switch to stable `dune-randomfield` release branch !151, !153
 * System tests for executing `dorie pfg` module !153
+* Finite volume solver for the Richards equation !132
 
 ### Fixed
 * Solver in `RichardsSimulation` was using the wrong time variable.

doc/default_files/richards-parameters.xml

@@ -52,9 +52,8 @@ adding an empty line, make text **bold** or ``monospaced``.
 
     <parameter name="vertexData">
       <definition> Plot vertex based (``true``) or cell-centered (``false``)
-       data into VTK files. Vertex based data might render sharp
-       parameterization boundaries inappropriately.
-       System tests and plotting functions (``dorie plot``) require
+       data into VTK files. Prefer vertex over cell data for full-precision
+       output. System tests and plotting functions (``dorie plot``) require
        cell-centered data.
       </definition>
       <values> true, false </values>
@@ -240,10 +239,12 @@ 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> Polynomial order of the DG method used. Setting this value
+        to 0 (zero) selects the finite volume (FV) solver (only compatible to
+        structured rectangular grids).
+      </definition>
       <suggestion> 1 </suggestion>
-      <values> 1, 2, 3 </values>
+      <values> 0, 1, 2, 3 </values>
     </parameter>
   </category>
 

doc/introduction/features.rst

@@ -39,6 +39,17 @@ singularities.
 Use the value ``richards`` in the ``simulation.mode`` keyword to run the
 standalone Richards solver.
 
+.. _richards_solver_options:
+
+Solver Options
+--------------
+
+DORiE includes a finite volume (FV) and a discontinuous Galerkin (DG) solver
+for the Richards equation. The FV solver can be selected for structured
+rectangular grids by setting ``richards.numerics.FEorder = 0`` when adaptivity
+is disabled. The DG solver is used for orders :math:`k \geq 1` and is available
+for all grid options.
+
 Solute Transport Solver
 =======================
 

doc/manual/fluxes.rst

@@ -10,34 +10,44 @@ Firstly, we have to recall that DORiE solves a Discontinuous Galerking finite el
 
 From the description above one can infer that one has to distinguish between * fluxes* at the interior of each element and at the intersections of all elements (we call these intersections skeleton of the grid). Unfortunately, there is no a standard form to write the skeleton fluxes on formats like VTK and that's the main reason why DORiE only provides the interior fluxes. However, assuming one can write both fluxes into some output format, they are still discontinuous (notice that direct use of discontinuous fluxes are useless for conservative computations since the transported quantities are very likely to get stagnated or over-transported in the nearby of intersections between elements). It means that it is needed some sort of post-processing that ensures that the *mass* is still locally and globally conserved.
 
-.. DORiE allow finite volume computations under certain specific conditions. In such case, if generated, the raw flux output generated by DORiE has no meaning. The reason is that finite volumes are computed with finite elements of order 0 where gradients are 0.
+When employing a finite volume solver, the regular flux output is omitted
+because the basis function gradients are always zero. In this case, the only
+option for evaluating flux data is flux reconstruction which has to be enabled
+in the :doc:`config file <config-file>`.
 
 Flux reconstruction
 -------------------
 The flux reconstruction is a projection of the fluxes used in the Discontinuous Galerkin method into a vector field function. Using correct elements, this procedure can ensure that fluxes in normal direction to the element are *equivalent* to those computed by the Discontinuous Galerkin method, and most importantly, it can also ensure the continuity of them. Hence, the resulting vector field is useful to compute other problems that rely on the fluxes of the water (i.e. solute transport).
 
-The flux reconstruction technique always use Raviar Thomas finite elements of one degree less than the one set for the Richards model. It can be identified in the vtk file by the name ``flux_RT{k-1}``, where ``k`` is the finite element order set for the Richards model. Flux reconstruction is not available for non-conforming grids (i.e. Cube-adaptive grids).
-
-+---------------------------+---+---+---+
-| Richards FEorder          | 1 | 2 | 3 |
-+============+====+=========+===+===+===+
-|            | 2D | Simplex | ✓ | ✓ | ✓ |
-|            |    +---------+---+---+---+
-|            |    | Cube    | ✓ | ✓ | ✓ |
-| Non-adapt. +----+---------+---+---+---+
-|            | 3D | Simplex | ✓ | ✓ | ✓ |
-|            |    +---------+---+---+---+
-|            |    | Cube    | ✓ | ✓ |   |
-+------------+----+---------+---+---+---+
-|            | 2D | Simplex | ✓ | ✓ | ✓ |
-|            |    +---------+---+---+---+
-|            |    | Cube    |   |   |   |
-| Adapt.     +----+---------+---+---+---+
-|            | 3D | Simplex | ✓ | ✓ | ✓ |
-|            |    +---------+---+---+---+
-|            |    | Cube    |   |   |   |
-+------------+----+---------+---+---+---+
+The flux reconstruction technique always use Raviar Thomas finite elements of one degree less than the one set for the Richards model. It can be identified in the vtk file by the name ``flux_RT{min(k-1,0)}``, where ``k`` is the finite element order set for the Richards model. Flux reconstruction is not available for non-conforming grids (i.e. Cube-adaptive grids).
+
++---------------------------+---+---+---+---+
+| Richards FEorder          | 0 | 1 | 2 | 3 |
++============+====+=========+===+===+===+===+
+|            | 2D | Simplex | ✗ | ✓ | ✓ | ✓ |
+|            |    +---------+---+---+---+---+
+|            |    | Cube    | ✓ | ✓ | ✓ | ✓ |
+| Non-adapt. +----+---------+---+---+---+---+
+|            | 3D | Simplex | ✗ | ✓ | ✓ | ✓ |
+|            |    +---------+---+---+---+---+
+|            |    | Cube    | ✓ | ✓ | ✓ |   |
++------------+----+---------+---+---+---+---+
+|            | 2D | Simplex | ✗ | ✓ | ✓ | ✓ |
+|            |    +---------+---+---+---+---+
+|            |    | Cube    | ✗ |   |   |   |
+| Adapt.     +----+---------+---+---+---+---+
+|            | 3D | Simplex | ✗ | ✓ | ✓ | ✓ |
+|            |    +---------+---+---+---+---+
+|            |    | Cube    | ✗ |   |   |   |
++------------+----+---------+---+---+---+---+
+
+Legend:
+
+* ✓: Flux reconstruction available.
+* ( ): Flux reconstruction unavailable.
+* ✗: Invalid setting. Finite volume solvers only works on regular grids,
+  see :ref:`Richards Solver Options <richards_solver_options>`.
 
 Usage
 -----
-To activate/deactivate flux reconstruction use the keyword ``richards.numeric.fluxReconstruction = true/false`` in the :doc:`config file<man-config-file>`.
+To activate/deactivate flux reconstruction use the keyword ``richards.numeric.fluxReconstruction = true/false`` in the :doc:`config file <config-file>`.

dune/dorie/model/richards/adaptivity/operator_error_indicator.hh

@@ -15,7 +15,7 @@
 #include <dune/pdelab/localoperator/defaultimp.hh>
 #include <dune/pdelab/finiteelement/localbasiscache.hh>
 
-#include <dune/dorie/model/richards/local_operator.hh>
+#include <dune/dorie/model/richards/local_operator_DG.hh>
 
 namespace Dune {
   namespace Dorie {

dune/dorie/model/richards/impl/CMakeLists.txt

 add_library(dorie-richards STATIC
+	sim_yasp_2_0.cc
 	sim_yasp_2_1.cc
-	sim_yasp_2_2.cc
 	sim_ug_2_1.cc
+	sim_yasp_2_2.cc
 	sim_yasp_2_3.cc
-	sim_yasp_3_1.cc
 	sim_ug_2_2.cc
+	sim_yasp_3_0.cc
+	sim_yasp_3_1.cc
+	sim_ug_2_3.cc
 	sim_yasp_3_2.cc
 	sim_yasp_3_3.cc
-	sim_ug_2_3.cc
 	sim_ug_3_1.cc
 	sim_ug_3_2.cc
 	sim_ug_3_3.cc)

dune/dorie/model/richards/impl/sim_yasp_2_0.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/dorie/model/richards/impl/impl.hh>
+#include <dune/dorie/model/richards/richards.cc>
+
+namespace Dune{
+namespace Dorie{
+
+template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,0>>;
+
+} // namespace Dorie
+} // namespace Dune
\ No newline at end of file

dune/dorie/model/richards/impl/sim_yasp_3_0.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/dorie/model/richards/impl/impl.hh>
+#include <dune/dorie/model/richards/richards.cc>
+
+namespace Dune{
+namespace Dorie{
+
+template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,0>>;
+
+} // namespace Dorie
+} // namespace Dune
\ No newline at end of file

dune/dorie/model/richards/local_operator.hh → dune/dorie/model/richards/local_operator_DG.hh

@@ -133,7 +133,8 @@ inline auto read_experimental_operator_settings(
  *  enough to support p-adaptivity.  (Another likely candidate would be
  *  differing geometry types, i.e. hybrid meshes.)
  *
- * @ingroup LocalOperators
+ * @ingroup    LocalOperators
+ * @ingroup    RichardsModel
  */
 template<typename Traits, typename Parameter, typename Boundary, typename SourceTerm, typename FEM, bool adjoint>
   class RichardsDGSpatialOperator
@@ -937,6 +938,7 @@ public:
  * @brief DG local operator for the temporal derivative of the Richards equation
  *
  * @ingroup LocalOperators
+ * @ingroup RichardsModel
  */
 template<typename Traits, typename Parameter, typename FEM, bool adjoint>
   class RichardsDGTemporalOperator

dune/dorie/model/richards/local_operator_FV.hh

+#ifndef DUNE_DORIE_RICHARDS_OPERATOR_FV_HH
+#define DUNE_DORIE_RICHARDS_OPERATOR_FV_HH
+
+#include <dune/geometry/referenceelements.hh>
+
+#include <dune/pdelab/common/quadraturerules.hh>
+#include <dune/pdelab/localoperator/pattern.hh>
+#include <dune/pdelab/localoperator/flags.hh>
+#include <dune/pdelab/localoperator/idefault.hh>
+#include <dune/pdelab/localoperator/defaultimp.hh>
+
+#include <dune/dorie/common/typedefs.hh>
+
+namespace Dune {
+namespace Dorie {
+namespace Operator {
+
+/*-------------------------------------------------------------------------*//**
+ * @brief      Spatial local operator for the Richards equation in a finite 
+ *             volume scheme.
+ * @details    It solves the spatial part of the Richards equation.
+ * 
+ * @author     Santiago Ospina De Los Ríos
+ * @date       2019
+ * @ingroup    LocalOperators
+ * @ingroup    RichardsModel
+ * 
+ * @tparam     Parameter                 Type of the class providing parameters
+ * @tparam     Boundary                  Type of the class providing boundary
+ *                                       conditions
+ */
+template<class Parameter, class Boundary>
+class RichardsFVSpatialOperator
+  : public Dune::PDELab::NumericalJacobianVolume<
+              RichardsFVSpatialOperator<Parameter,Boundary>>
+  , public Dune::PDELab::NumericalJacobianSkeleton<
+              RichardsFVSpatialOperator<Parameter,Boundary>>
+  , public Dune::PDELab::NumericalJacobianBoundary<
+              RichardsFVSpatialOperator<Parameter,Boundary>>
+  , public Dune::PDELab::FullSkeletonPattern
+  , public Dune::PDELab::FullVolumePattern
+  , public Dune::PDELab::LocalOperatorDefaultFlags
+  , public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<double>
+{
+public:
+  // Pattern assembly flags
+  enum { doPatternVolume = true };
+  enum { doPatternSkeleton = true };
+
+  // Residual assembly flags
+  enum { doAlphaBoundary = true };
+  enum { doAlphaSkeleton  = true };
+
+
+  /**
+   * @brief      Constructor of RichardsFVSpatialOperator
+   *
+   * @param[in]  param     The parameter class
+   * @param[in]  boundary  The boundary class
+   */
+  RichardsFVSpatialOperator(
+    const std::shared_ptr<const Parameter>& param,
+    const std::shared_ptr<const Boundary>& boundary
+  ) : _param(param)
+    , _boundary(boundary)
+    , _time(0.)
+  {}
+
+   /*----------------------------------------------------------------------*//**
+    * @brief      Skeleton integral depending on test and ansatz functions. Each
+    *             face is only visited once since this method is symmetric
+    *
+    * @param[in]  ig      The intersection entity of the grid (inside + outside
+    *                     entities)
+    * @param[in]  lfsu_i  The inside ansatz local function space
+    * @param[in]  x_i     The coefficients of the lfsu_i
+    * @param[in]  lfsv_i  The inside test local function space
+    * @param[in]  lfsu_o  The outside ansatz local function space
+    * @param[in]  x_o     The coefficients of the lfsu_o
+    * @param[in]  lfsv_o  The outside test local function space
+    * @param[out] r_i     The view of the residual vector w.r.t lfsu_i
+    * @param[out] r_o     The view of the residual vector w.r.t lfsu_o
+    *
+    * @tparam     IG      The type for ig
+    * @tparam     LFSU    The type for lfsu_i and lfsu_o
+    * @tparam     X       The type for x_i and x_o
+    * @tparam     LFSV    The type for lfsv_i and lfsv_o
+    * @tparam     R       The type for r_i and r_o
+    */
+  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
+  void alpha_skeleton(const IG& ig, 
+    const LFSU& lfsu_i, const X& x_i, const LFSV& lfsv_i,
+    const LFSU& lfsu_o, const X& x_o, const LFSV& lfsv_o, 
+    R& r_i, R& r_o) const
+  {
+    // Get local basis traits from local function space
+    using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType::
+                        Traits::LocalBasisType::Traits;
+
+    // Get range, range field and gradient for trial space
+    using RangeU = typename LocalBasisTraitsU::RangeType;
+
+    // Assert that polynomial degree is always 0
+    assert(lfsu_i.finiteElement().localBasis().order()==0);
+    assert(lfsu_o.finiteElement().localBasis().order()==0);
+    assert(lfsv_i.finiteElement().localBasis().order()==0);
+    assert(lfsv_o.finiteElement().localBasis().order()==0);
+
+    // Define entities for each frame of reference
+    const auto& entity_f = ig.intersection();
+    const auto& entity_i = ig.inside();
+    const auto& entity_o = ig.outside();
+
+    // Get geometries
+    const auto geo_f = entity_f.geometry();
+    const auto geo_i = entity_i.geometry();
+    const auto geo_o = entity_o.geometry();
+
+    // Get volume of entities
+    const auto volume_f = geo_f.volume();
+
+    // Get normal
+    const auto normal_f = entity_f.centerUnitOuterNormal();
+
+    // Entity centers in global coordinates
+    const auto center_position_i_g = geo_i.center();
+    const auto center_position_o_g = geo_o.center();
+
+    // Inside/outside solute value
+    const RangeU u_i = x_i(lfsu_i,0);
+    const RangeU u_o = x_o(lfsu_o,0);
+
+    // bind parameterization and retrieve functions
+    _param->bind(entity_i);
+    const auto saturation_f_i = _param->saturation_f();
+    const auto conductivity_f_i = _param->conductivity_f();
+    
+    _param->bind(entity_o);
+    const auto saturation_f_o = _param->saturation_f();
+    const auto conductivity_f_o = _param->conductivity_f();
+
+    // Compute the conductivity
+    const RangeU cond_i = conductivity_f_i(saturation_f_i(u_i));
+    const RangeU cond_o = conductivity_f_o(saturation_f_o(u_o));
+
+    // Harmonic average for interface conductivity
+    const RangeU cond_f = 2.0/(   1.0/(cond_i + 1E-30) 
+                                + 1.0/(cond_o + 1E-30) );
+
+    // Distance between the two entity centers
+    const auto center_position_diff = center_position_i_g - center_position_o_g;
+    const auto distance = center_position_diff.two_norm();
+
+    // Finite difference of u between the two entities
+    RangeU dudn = (u_o-u_i)/distance;
+
+    // Update gradient with gravity vector
+    dudn -= _param->gravity()*normal_f;
+
+    // Water flux in normal direction w.r.t the intersection
+    auto water_flux_n = - cond_f*dudn;
+
+    // // Upwinding
+    // if (water_flux_n > 0.)
+    //   water_flux_n = - cond_i*dudn;
+    // else
+    //   water_flux_n = - cond_o*dudn;
+
+    // Symmetric contribution to residual on inside element
+    r_i.accumulate(lfsv_i, 0,  water_flux_n*volume_f );
+    r_o.accumulate(lfsv_o, 0, -water_flux_n*volume_f );
+  }
+
+  /*-----------------------------------------------------------------------*//**
+   * @brief      Boundary integral depending on test and ansatz functions.
+   *
+   * @param[in]  ig      The intersection entity of the grid (inside + outside
+   *                     entities)
+   * @param[in]  lfsu_i  The inside ansatz local function space
+   * @param[in]  x_i     The coefficients of the lfsu_i
+   * @param[in]  lfsv_i  The inside test local function space
+   * @param[out] r_i     The view of the residual vector w.r.t lfsu_i
+   *
+   * @tparam     IG      The type for ig
+   * @tparam     LFSU    The type for lfsu_i
+   * @tparam     X       The type for x_i
+   * @tparam     LFSV    The type for lfsv_i
+   * @tparam     R       The type for r_i
+   */
+  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
+  void alpha_boundary (const IG& ig, const LFSU& lfsu_i, const X& x_i,
+                       const LFSV& lfsv_i, R& r_i) const
+  {
+    // Get local basis traits from local function space
+    using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType::
+                        Traits::LocalBasisType::Traits;
+
+    // Get range, range field and gradient for trial space
+    using RangeU = typename LocalBasisTraitsU::RangeType;
+    using RangeFieldU = typename LocalBasisTraitsU::RangeFieldType;
+    
+    // Assert that polynomial degree is always 0
+    assert(lfsu_i.finiteElement().localBasis().order()==0);
+    assert(lfsv_i.finiteElement().localBasis().order()==0);
+
+    // Get entities
+    const auto& entity_f = ig.intersection(); 
+    const auto& entity_i = ig.inside();
+
+    // Get geometries
+    const auto geo_f = entity_f.geometry();
+    const auto& ref_el_f = referenceElement(geo_f);
+    const auto& center_position_f = ref_el_f.position(0,0);
+
+    // Face volume for integration
+    const auto volume_f = geo_f.volume();
+
+    // Normal vector of intersection
+    const auto normal_f = ig.centerUnitOuterNormal();
+
+    // Retrieve boundary condition
+    auto bc = _boundary->bc(entity_f,center_position_f,_time);
+
+    // bind parameterization and retrieve functions
+    _param->bind(entity_i);
+    const auto saturation_f_i = _param->saturation_f();
+    const auto conductivity_f_i = _param->conductivity_f();
+
+    if (BoundaryCondition::isDirichlet(bc))
+    {
+      const auto geo_i = entity_i.geometry();
+      const auto center_position_i_g = geo_i.center();
+
+      // Face center in global coordinates
+      const auto center_position_f_g = geo_f.center();
+
+      // Compute distance of these two points
+      const auto position_diff = center_position_i_g - center_position_f_g;
+      const auto distance = position_diff.two_norm();
+
+      // Evaluate Dirichlet condition
+      const auto g = _boundary->g(entity_f,center_position_f,_time);
+
+      // Inside unknown value
+      const RangeU u_i = x_i(lfsu_i,0);
+
+      // Compute the conductivity
+      const RangeU cond_i = conductivity_f_i(saturation_f_i(u_i));
+
+      // Finite difference of u between the two entities
+      RangeU dudn = (g-u_i)/distance;
+
+      // Update gradient with gravity vector
+      dudn -= _param->gravity()*normal_f;
+
+      // Water flux in normal direction w.r.t the intersection
+      const auto water_flux_n = - cond_i*dudn;
+
+      // Contribution to residual from Dirichlet boundary
+      r_i.accumulate(lfsv_i, 0, water_flux_n*volume_f);
+    }
+    else if (BoundaryCondition::isNeumann(bc))
+    {
+      // Contribution to residual from Neumann boundary
+      RangeFieldU water_flux_f = _boundary->j(entity_f,center_position_f,_time);
+
+      // Convert water flux into units of matric head flux
+      water_flux_f *= std::abs( normal_f * _param->gravity() );
+      // water_flux_f *= _param->gravity().two_norm();
+
+      r_i.accumulate(lfsv_i, 0, water_flux_f*volume_f );
+    }
+  }
+
+  /*-----------------------------------------------------------------------*//**
+   * @brief      Sets the time.
+   *
+   * @param[in]  t           time of type RangeField
+   *
+   * @tparam     RangeField  type of the range field
+   */
+  template<class RangeField>
+  void setTime (RangeField t)
+  {
+    _time = t;
+  }
+
+private:
+  /// Class providing equation parameters
+  const std::shared_ptr<const Parameter> _param;
+  /// Class providing boundary conditions
+  const std::shared_ptr<const Boundary> _boundary;
+  double _time;
+};
+
+/**
+ * @brief      Temporal local operator Richards equation in a finite volume
+ *             scheme.
+ * @details    It solves the temporal part of the Richards equation.
+ * @author     Santiago Ospina De Los Ríos
+ * @date       2018
+ * @ingroup    LocalOperators
+ * @ingroup    RichardsModel
+ *
+ * @tparam     Parameter                 Type of the class providing parameters
+ */
+template<class Parameter>
+  class RichardsFVTemporalOperator
+  : public Dune::PDELab::NumericalJacobianVolume<
+              RichardsFVTemporalOperator<Parameter>>
+  , public Dune::PDELab::NumericalJacobianApplyVolume<
+              RichardsFVTemporalOperator<Parameter>>
+  , public Dune::PDELab::FullSkeletonPattern
+  , public Dune::PDELab::FullVolumePattern
+  , public Dune::PDELab::LocalOperatorDefaultFlags
+  , public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<double>
+{
+public:
+
+  // Pattern assembly flags
+  enum { doPatternVolume = true };
+
+  // Residual assembly flags
+  enum { doAlphaVolume = true };
+
+public:
+
+  /**
+   * @brief      Constructor of RichardsFVTemporalOperator
+   *
+   * @param[in]  param     The parameter class
+   */
+  RichardsFVTemporalOperator(const std::shared_ptr<const Parameter>& param)
+    : _param(param)
+    , _time(0.)
+  {}
+
+  /*-----------------------------------------------------------------------*//**
+   * @brief      Volume integral depending on test and ansatz functions
+   *
+   * @param[in]  eg    The entity of the grid
+   * @param[in]  lfsu  The ansatz local function space
+   * @param[in]  x     The coefficients of the lfsu
+   * @param[in]  lfsv  The test local function space
+   * @param[out] r     The view of the residual vector w.r.t lfsu
+   *
+   * @tparam     EG    The type for eg
+   * @tparam     LFSU  The type for lfsu
+   * @tparam     X     The type for x
+   * @tparam     LFSV  The type for lfsv
+   * @tparam     R     The type for r
+   */
+  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
+  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, 
+                     const LFSV& lfsv, R& r) const
+  {
+    // Get local basis traits from local function space
+    using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType::
+                        Traits::LocalBasisType::Traits;
+    // Get range for trial space
+    using RangeU = typename LocalBasisTraitsU::RangeType;
+    const RangeU u = x(lfsu,0);
+
+    // entity geometry
+    const auto geo = eg.geometry();
+
+    // bind parameterization and retrieve functions
+    _param->bind(eg.entity());
+    const auto saturation_f = _param->saturation_f();
+    const auto water_content_f = _param->water_content_f();
+
+    // Calculate water content from matric head
+    const RangeU water_content = water_content_f(saturation_f(u));
+
+    // update residual
+    r.accumulate(lfsv ,0 , water_content*geo.volume());
+  }
+
+  /**
+   * @brief      Sets the time.
+   *
+   * @param[in]  t     time of type RangeField
+   *
+   * @tparam     RF    Type of the range field
+   */
+  template<class RF>
+  void setTime (RF t)
+  {
+    _time = t;
+  }
+
+private:
+  /// class providing equation parameters
+  const std::shared_ptr<const Parameter> _param;
+  double _time;
+};
+
+
+} // namespace Operator
+} // namespace Dorie
+} // namespace Dune
+