Merge branch '153-add-upwinding-option-to-richards-fv-solver' into 'master'

Resolve "Add upwinding option to Richards FV solver"

Closes #153

See merge request !161

CHANGELOG.md

@@ -46,6 +46,7 @@
 * Define compile-time settings via CPP macros !144
 * [Google Test](https://github.com/google/googletest) unit test framework
     as Git Submodule !159
+* Upwinding options for Richards finite volume local operator !161
 
 ### Changed
 * `Simulation` is renamed `RichardsSimulation` and moved to 

CMakeLists.txt

@@ -14,12 +14,6 @@ endif()
 # add extra flags to debug compiler flags
 set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall")
 
-# option to change DG scheme via config file
-option(EXPERIMENTAL_DG_FEATURES
-  "Enable experimental DG settings through the config file"
-  OFF
-)
-
 # 
 if(NOT (dune-common_DIR OR dune-common_ROOT OR
       "${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))

README.md

@@ -199,16 +199,6 @@ If you installed [Anaconda](https://conda.io/docs/user-guide/install/download.ht
     path chosen as installation prefix when configuring HDF5.
 
 
-### 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
-`-DEXPERIMENTAL_DG_FEATURES=On`, and rebuilding it.
-
-The configuration settings in the section `[dg.experimental]` will then override
-the default settings.
-
-
 ### Recommended Third-Party Software
 
 The following software packages are cross-platform, so you should be able to find a release that fits your operating system:

doc/default_files/richards-parameters.xml

@@ -231,13 +231,6 @@ adding an empty line, make text **bold** or ``monospaced``.
   </category>
 
   <category name="numerics">
-    <parameter name="penaltyFactor">
-      <definition> Penalty factor to be used in the Discontinuous Galerkin scheme
-      </definition>
-      <values> float </values>
-      <suggestion> 10 </suggestion>
-    </parameter>
-
     <parameter name="FEorder">
       <definition> Polynomial order of the DG method used. Setting this value
         to 0 (zero) selects the finite volume (FV) solver (only compatible to
@@ -245,6 +238,52 @@ adding an empty line, make text **bold** or ``monospaced``.
       </definition>
       <suggestion> 1 </suggestion>
       <values> 0, 1, 2, 3 </values>
+      <comment> Select '0' for the finite volume (FV) solver </comment>
+    </parameter>
+
+    <parameter name="upwinding">
+      <definition> Upwinding method for skeleton terms. Upwinding typically
+        increases numeric stability while reducing accuracy.
+
+        **semiUpwind:** Apply upwinding to conductivity factor (only).
+
+        **fullUpwind:** Apply upwinding to conductivity.
+        **Not recommended for DG**.
+      </definition>
+      <values> none, semiUpwind, fullUpwind </values>
+      <suggestion> none </suggestion>
+      <comment> Choose upwinding type: 'none', 'semiUpwind', 'fullUpwind'
+      </comment>
+    </parameter>
+
+    <parameter name="DGMethod">
+      <definition> DG discretization method for skeleton terms.
+
+        **SIPG:** Symmetric Interior Penalty
+
+        **NIPG:** Non-Symmetric Interior Penalty
+
+        **OBB:** Oden, Babuska, Baumann: no penalty term
+
+        **IIP:** Incomplete Interior Penalty: no symmetry term
+      </definition>
+      <values> SIPG, NIPG, OBB, IIP </values>
+      <suggestion> SIPG </suggestion>
+    </parameter>
+
+    <parameter name="DGWeights">
+      <definition> Apply harmonic weighting to skeleton term contributions
+        in DG.
+      </definition>
+      <values> true, false </values>
+      <suggestion> true </suggestion>
+    </parameter>
+
+    <parameter name="penaltyFactor">
+      <definition> Penalty factor to be used in the Discontinuous Galerkin scheme
+      </definition>
+      <values> float </values>
+      <suggestion> 10 </suggestion>
     </parameter>
   </category>
 
@@ -328,41 +367,4 @@ adding an empty line, make text **bold** or ``monospaced``.
     </parameter>
   </category>
 
-   <category name="numerics.experimental">
-    <parameter name="method">
-      <definition> DG discretization method for skeleton terms.
-
-        **SIPG:** Symmetric Interior Penalty
-
-        **NIPG:** Non-Symmetric Interior Penalty
-
-        **OOB:** Oden, Babuska, Baumann: no penalty term
-
-        **IIP:** Incomplete Interior Penalty: no symmetry term
-      </definition>
-      <values> SIPG, NIPG, OOB, IIP </values>
-      <suggestion> SIPG </suggestion>
-      <comment> Experimental settings are enabled by the appropriate CMake flag.
-      </comment>
-    </parameter>
-
-    <parameter name="upwinding">
-      <definition> Upwinding method for skeleton terms.
-
-        **semiUpwind:** Apply upwinding to conductivity factor (only).
-
-        **fullUpwind:** Apply upwinding on numeric flux and conductivity.
-      </definition>
-      <values> none, semiUpwind, fullUpwind </values>
-      <suggestion> none </suggestion>
-    </parameter>
-
-    <parameter name="weights">
-      <definition> Apply harmonic weighting to skeleton term contributions.
-      </definition>
-      <values> true, false </values>
-      <suggestion> true </suggestion>
-    </parameter>
-  </category>
-
 </dorie>

dune/dorie/CMakeLists.txt

@@ -22,8 +22,3 @@ endforeach()
 
 # Global target
 add_custom_target("dorie" DEPENDS richards transport)
-
-# enable setting operator scheme from config file
-if(EXPERIMENTAL_DG_FEATURES)
-    target_compile_definitions("dorie" PUBLIC -DEXPERIMENTAL_DG_FEATURES)
-endif()

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

@@ -15,7 +15,8 @@
 #include <dune/pdelab/localoperator/defaultimp.hh>
 #include <dune/pdelab/finiteelement/localbasiscache.hh>
 
-#include <dune/dorie/model/richards/local_operator_DG.hh>
+#include <dune/dorie/common/typedefs.hh>
+#include <dune/dorie/model/richards/local_operator_cfg.hh>
 
 namespace Dune {
   namespace Dorie {
@@ -63,14 +64,14 @@ namespace Dune {
 
       FluxErrorEstimator (const Dune::ParameterTree& config,
         const Parameter& param_, const Boundary& boundary_,
-        RichardsDGMethod::Type method_ = RichardsDGMethod::SIPG,
+        RichardsDGMethod method_ = RichardsDGMethod::SIPG,
         int intorderadd_ = 2, int quadrature_factor_ = 2)
       : param(param_), boundary(boundary_),
         intorderadd(intorderadd_), quadrature_factor(quadrature_factor_),
         penalty_factor(config.get<RF>("numerics.penaltyFactor")),
         time(0.0)
       {
-        if (method_ == RichardsDGMethod::OOB)
+        if (method_ == RichardsDGMethod::OBB)
           penalty_factor = 0.0;
       }
 
@@ -310,4 +311,4 @@ namespace Dune {
   } // namespace Dorie
 } // namespace Dune
 
-#endif // DUNE_DORIE_ERROR_INDICATOR_HH
\ No newline at end of file
+#endif // DUNE_DORIE_ERROR_INDICATOR_HH

dune/dorie/model/richards/local_operator_FV.hh

@@ -10,6 +10,7 @@
 #include <dune/pdelab/localoperator/defaultimp.hh>
 
 #include <dune/dorie/common/typedefs.hh>
+#include <dune/dorie/model/richards/local_operator_cfg.hh>
 
 namespace Dune {
 namespace Dorie {
@@ -60,9 +61,11 @@ public:
    */
   RichardsFVSpatialOperator(
     const std::shared_ptr<const Parameter>& param,
-    const std::shared_ptr<const Boundary>& boundary
+    const std::shared_ptr<const Boundary>& boundary,
+    const RichardsUpwinding upwinding = RichardsUpwinding::none
   ) : _param(param)
     , _boundary(boundary)
+    , _upwinding(upwinding)
     , _time(0.)
   {}
 
@@ -126,10 +129,21 @@ public:
     const auto center_position_i_g = geo_i.center();
     const auto center_position_o_g = geo_o.center();
 
+    // 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();
+
     // Inside/outside solute value
     const RangeU u_i = x_i(lfsu_i,0);
     const RangeU u_o = x_o(lfsu_o,0);
 
+    // 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;
+
     // bind parameterization and retrieve functions
     _param->bind(entity_i);
     const auto saturation_f_i = _param->saturation_f();
@@ -143,28 +157,49 @@ public:
     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;
+    // Interface conductivity factor (Upwinding applies)
+    RangeU cond_f(0.0);
+    if (_upwinding == RichardsUpwinding::none)
+    {
+      // Harmonic average for interface conductivity
+      cond_f = 2.0/(   1.0/(cond_i + 1E-30)
+                     + 1.0/(cond_o + 1E-30) );
+    }
+    else if (_upwinding == RichardsUpwinding::semiUpwind)
+    {
+      // Harmonic average of saturated conductivities
+      const RangeU cond_sat_i = conductivity_f_i(1.0);
+      const RangeU cond_sat_o = conductivity_f_o(1.0);
+
+      // Upwinding conductivity factor
+      RangeU cond_upwind_factor(0.0);
+      if (dudn < 0.0) {
+        // inward flux
+        cond_upwind_factor = cond_o / cond_sat_o;
+      }
+      else {
+        // outward flux
+        cond_upwind_factor = cond_i / cond_sat_i;
+      }
+
+      cond_f = 2.0 * cond_upwind_factor / (   1.0/(cond_sat_i + 1E-30)
+                                            + 1.0/(cond_sat_o + 1E-30) );
+    }
+    else // RichardsFVUpwinding::fullUpwind
+    {
+      // No average: Use upwind conductivity
+      if (dudn < 0.0) {
+        // inward flux
+        cond_f = cond_o;
+      }
+      else {
+        // outward flux
+        cond_f = cond_i;
+      }
+    }
 
     // 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;
+    const auto water_flux_n = - cond_f * dudn;
 
     // Symmetric contribution to residual on inside element
     r_i.accumulate(lfsv_i, 0,  water_flux_n*volume_f );
@@ -290,6 +325,9 @@ private:
   const std::shared_ptr<const Parameter> _param;
   /// Class providing boundary conditions
   const std::shared_ptr<const Boundary> _boundary;
+  /// Upwinding setting
+  const RichardsUpwinding _upwinding;
+  /// Operator internal time
   double _time;
 };
 

dune/dorie/model/richards/local_operator_cfg.hh

+#ifndef DUNE_DORIE_RICHARDS_LOP_CFG_HH
+#define DUNE_DORIE_RICHARDS_LOP_CFG_HH
+
+#include <tuple>
+#include <string>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/dorie/common/logging.hh>
+
+namespace Dune {
+namespace Dorie {
+namespace Operator {
+
+/**
+ *  \addtogroup RichardsModel
+ *  \{
+ */
+
+//! Interior penalty type
+enum class RichardsDGMethod
+{
+    NIPG, //!< Non-symmetric interior penalty
+    SIPG, //!< Symmetric interior penalty (default)
+    OBB,  //!< Oden, Babuska, Baumann (no penalty term)
+    IIP   //!< Incomplete interior penalty (no symmetry term)
+};
+
+//! Gradient weighting
+enum class RichardsDGWeights
+{
+    weightsOn, //!< harmonic weighting of flux contributions
+    weightsOff //!< no weighting of flux contributions
+};
+
+//! Upwinding type
+enum class RichardsUpwinding
+{
+    fullUpwind, //!< Upwind conductivity
+    semiUpwind, //!< Upwind conductivity factor
+    none //!< No upwinding
+};
+
+
+/// Read operator settings from a parameter file
+/* \return Tuple of RichardsUpwinding, RichardsDGMethod, and
+ *    RichardsDGWeights to be inserted into the local operator constructor.
+ */
+inline auto read_operator_settings(
+    const Dune::ParameterTree& inifile)
+        -> std::tuple<Operator::RichardsUpwinding,
+                      Operator::RichardsDGMethod,
+                      Operator::RichardsDGWeights>
+{
+    const auto log = Dorie::get_logger(log_richards);
+    log->debug("Reading local operator settings:");
+
+    RichardsUpwinding upwinding;
+    RichardsDGMethod method = RichardsDGMethod::SIPG;
+    RichardsDGWeights weights = RichardsDGWeights::weightsOn;
+
+    // Upwinding
+    const auto upwinding_str
+        = inifile.get<std::string>("numerics.upwinding");
+    log->debug("  Upwinding scheme: {}", upwinding_str);
+    if (upwinding_str == "none")
+        upwinding = RichardsUpwinding::none;
+    else if (upwinding_str == "semiUpwind")
+        upwinding = RichardsUpwinding::semiUpwind;
+    else if (upwinding_str == "fullUpwind")
+        upwinding = RichardsUpwinding::fullUpwind;
+    else {
+        log->error("Unknown upwinding scheme: {}", upwinding_str);
+        DUNE_THROW(Dune::IOError, "Unknown upwinding scheme");
+    }
+
+    // Return here if the local operator is FV only
+    if (inifile.get<int>("numerics.FEorder") == 0)
+    {
+        log->debug("  Ignoring settings 'numerics.DGMethod' and  "
+                   "'numerics.DGWeights' for finite volume solver.");
+        return std::make_tuple(upwinding, method, weights);
+    }
+
+    // DG Method
+    const auto method_str = inifile.get<std::string>("numerics.DGMethod");
+    log->debug("  DG method: {}", method_str);
+    if (method_str == "SIPG")
+        method = RichardsDGMethod::SIPG;
+    else if (method_str == "NIPG")
+        method = RichardsDGMethod::NIPG;
+    else if (method_str == "OBB")
+        method = RichardsDGMethod::OBB;
+    else if (method_str == "IIP")
+        method = RichardsDGMethod::IIP;
+    else {
+        log->error("Unknown DG method: {}", method_str);
+        DUNE_THROW(Dune::IOError, "Unknown DG method");
+    }
+
+    // Harmonic weighting of fluxes
+    const auto weights_str = inifile.get<bool>("numerics.DGWeights");
+    log->debug("  Harmonic weights for interface fluxes (DG): {}", 
+               weights_str);
+    if (weights_str)
+        weights = RichardsDGWeights::weightsOn;
+    else
+        weights = RichardsDGWeights::weightsOff;
+
+    return std::make_tuple(upwinding, method, weights);
+}
+
+/**
+ *  \}
+ *  // group RichardsModel
+ */
+
+} // namespace Operator
+} // namespace Dorie
+} // namespace Dune
+
+#endif // DUNE_DORIE_RICHARDS_LOP_CFG_HH

dune/dorie/model/richards/richards.cc

@@ -55,30 +55,27 @@ ModelRichards<Traits>::ModelRichards (
 	fsource = std::make_shared<const FlowSource>(inifile);
 	finitial = FlowInitialFactory::create(inifile, gv, this->_log);
 
+	// Read local operator settings
+	namespace OP = Dune::Dorie::Operator;
+	const auto settings = OP::read_operator_settings(inifile);
+
 	// --- Local Operators ---
-  if constexpr (order>0)
-  {
+	if constexpr (order>0)
+	{
 		this->_log->debug("Setting up local grid operators: DG method");
 
-#ifdef EXPERIMENTAL_DG_FEATURES
-		// read experimental settings from inifile
-		namespace OP = Dune::Dorie::Operator;
-		const auto settings = OP::read_experimental_operator_settings(inifile);
-		const auto method = std::get<OP::RichardsDGMethod::Type>(settings);
-		const auto upwinding = std::get<OP::RichardsDGUpwinding::Type>(settings);
-		const auto weights = std::get<OP::RichardsDGWeights::Type>(settings);
+		const auto upwinding = std::get<OP::RichardsUpwinding>(settings);
+		const auto method = std::get<OP::RichardsDGMethod>(settings);
+		const auto weights = std::get<OP::RichardsDGWeights>(settings);
 
 		slop = std::make_unique<SLOP>(inifile, fparam, fboundary, fsource,
 									  method, upwinding, weights);
-#else
-		slop = std::make_unique<SLOP>(inifile, fparam, fboundary, fsource);
-#endif // EXPERIMENTAL_DG_FEATURES
-
 		tlop = std::make_unique<TLOP>(inifile, fparam);
 	}
 	else {
 		this->_log->debug("Setting up local grid operators: FV method");
-		slop = std::make_unique<SLOP>(fparam, fboundary);
+		const auto upwinding = std::get<OP::RichardsUpwinding>(settings);
+		slop = std::make_unique<SLOP>(fparam, fboundary, upwinding);
 		tlop = std::make_unique<TLOP>(fparam);
 	}
 

dune/dorie/model/richards/richards.hh

@@ -33,6 +33,7 @@
 #include <dune/dorie/model/richards/flow_source.hh>
 #include <dune/dorie/model/richards/local_operator_DG.hh>
 #include <dune/dorie/model/richards/local_operator_FV.hh>
+#include <dune/dorie/model/richards/local_operator_cfg.hh>
 
 namespace Dune{
 namespace Dorie{

test/ode_homogeneous_sand.mini.in

@@ -5,7 +5,10 @@ _asset_path = "${CMAKE_CURRENT_LIST_DIR}"
 _evaluation = ode
 
 _order = 0, 1, 2, 3 | expand prec
-__name = ode_homogeneous_sand_k{_order}
+_upwind = none, semiUpwind, fullUpwind | expand
+__name = ode_homogeneous_sand_k{_order}_{_upwind}
+
+{_order} > 0 and {_upwind} == fullUpwind | exclude
 
 [grid]
 gridType = rectangular
@@ -34,6 +37,7 @@ NewtonParameters.Reduction = 1E-10
 
 numerics.penaltyFactor = 10
 numerics.FEorder = 0, 1, 2, 3 | expand prec
+numerics.upwinding = {_upwind}
 
 [_ode]
 
@@ -41,4 +45,4 @@ headLower = 0.
 flux = -5.55e-6
 head_abstol = 3E-4, 3E-6, 2E-8, 6E-8 | expand prec
 flux_abstol = 1E100, 5E-11, 1E-9, 3E-12 | expand prec
-flux_rt_abstol = 4E-14
\ No newline at end of file
+flux_rt_abstol = 4E-14

test/ode_homogeneous_silt.mini.in

@@ -5,7 +5,10 @@ _asset_path = "${CMAKE_CURRENT_LIST_DIR}"
 _evaluation = ode
 
 _order = 0, 1, 2, 3 | expand prec
-__name = ode_homogeneous_silt_k{_order}
+_upwind = none, semiUpwind, fullUpwind | expand
+__name = ode_homogeneous_silt_k{_order}_{_upwind}
+
+{_order} > 0 and {_upwind} == fullUpwind | exclude
 
 [grid]
 gridType = rectangular
@@ -38,6 +41,7 @@ NewtonParameters.Reduction = 1E-10
 
 numerics.penaltyFactor = 10
 numerics.FEorder = 0, 1, 2, 3 | expand prec
+numerics.upwinding = {_upwind}
 
 [_ode]
 

test/ode_layered.mini.in

@@ -5,7 +5,10 @@ _asset_path = "${CMAKE_CURRENT_LIST_DIR}"
 _evaluation = ode
 
 _order = 0, 1, 2, 3 | expand prec
-__name = ode_layered_k{_order}
+_upwind = none, semiUpwind, fullUpwind | expand
+__name = ode_layered_k{_order}_{_upwind}
+
+{_order} > 0 and {_upwind} != none | exclude
 
 [grid]
 gridType = rectangular
@@ -38,6 +41,7 @@ NewtonParameters.Reduction = 1E-10
 
 numerics.penaltyFactor = 10
 numerics.FEorder = {_order}
+numerics.upwinding = {_upwind}
 
 [_ode]
 
@@ -45,4 +49,4 @@ headLower = 0.
 flux = -5.55e-6
 head_abstol = 1E-3, 3E-5, 4E-8, 4E-6 | expand prec
 flux_abstol = 1E100, 5E-10, 4E-9, 5E-11 | expand prec
-flux_rt_abstol = 4E-14
\ No newline at end of file
+flux_rt_abstol = 4E-14