Merge branch '43-add-options-for-altering-dg-formulation-in-local-operator' into 'master'

Resolve "add options for altering DG formulation in local operator"

Closes #43 and #41

See merge request !33

CMakeLists.txt

@@ -12,6 +12,12 @@ if(CMAKE_BUILD_TYPE_UPPER MATCHES DEBUG)
 endif()
 set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -Werror")
 
+# 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

@@ -143,6 +143,15 @@ If you installed [Anaconda](https://conda.io/docs/user-guide/install/download.ht
 
     **Installing DUNE and DORiE via `make install` into your usual work environment is strongly discouraged!**
 
+### 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
 

doc/default_files/parameters.xml

@@ -412,4 +412,41 @@ adding an empty line, make text **bold** or ``monospaced``.
     </parameter>
   </category>
 
+   <category name="dg.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

@@ -4,4 +4,9 @@ add_subdirectory(impl)
 add_subdirectory(test)
 
 add_executable("dorie" dorie.cc)
-dune_target_link_libraries(dorie dorie-impl ${DUNE_LIBS})
\ No newline at end of file
+dune_target_link_libraries(dorie dorie-impl ${DUNE_LIBS})
+
+# enable setting operator scheme from config file
+if(EXPERIMENTAL_DG_FEATURES)
+    target_compile_definitions("dorie" PUBLIC -DEXPERIMENTAL_DG_FEATURES)
+endif()
\ No newline at end of file

dune/dorie/interface/adaptivity.hh

@@ -42,7 +42,7 @@ private:
 	/// Adaptivity GFS
 	using AGFS = typename AGFSHelper::Type;
 	/// Error estimator local operator
-	using ESTLOP = Dune::Dorie::FluxErrorEstimator<Traits,Param,Boundary>;
+	using ESTLOP = Dune::Dorie::Operator::FluxErrorEstimator<Traits,Param,Boundary>;
 	/// Empty constraints container
 	using NoTrafo = Dune::PDELab::EmptyTransformation;
 	/// Solution vector type

dune/dorie/interface/simulation.cc

@@ -31,7 +31,20 @@ Simulation<Traits>::Simulation (Dune::MPIHelper& _helper, std::shared_ptr<Grid>
 	finitial = std::make_unique<FlowInitial>(inifile,*param,gv);
 
 	// --- Local Operators ---
+#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);
+
+	slop = std::make_unique<SLOP>(inifile, *param, gv, *fboundary, *fsource,
+																method, upwinding, weights);
+#else
 	slop = std::make_unique<SLOP>(inifile,*param,gv,*fboundary,*fsource);
+#endif // EXPERIMENTAL_DG_FEATURES
+
 	tlop = std::make_unique<TLOP>(inifile,*param,gv);
 	controller = std::make_unique<CalculationController>(inifile,*fboundary,helper);
 
@@ -128,6 +141,7 @@ bool Simulation<Traits>::compute_time_step ()
 	const RF t = controller->getTime();
 	const RF dt = controller->getDT();
 	bool exception = false;
+
 	try
 	{
 		// solve in parallel
@@ -140,7 +154,7 @@ bool Simulation<Traits>::compute_time_step ()
 		// solve sequentially
 		else {
 			pdesolver_seq->setMaxIterations(controller->getIterations());
-			osm_seq->apply(t,dt,*uold,*unew);
+			osm_seq->apply(t, dt, *uold, *unew);
 		}
 
 		*uold = *unew;

dune/dorie/solver/operator_DG.hh

@@ -19,16 +19,29 @@ namespace Operator {
 struct RichardsDGMethod
 {
   enum Type { NIPG, //!< Non-symmetric interior penalty
-              SIPG //!< Symmetric interior penalty (default)
+              SIPG, //!< Symmetric interior penalty (default)
+              OOB, //!< Oden, Babuska, Baumann (no penalty term)
+              IIP //!< Incomplete interior penalty (no symmetry)
             };
 };
 
-//! Weighting at grid interface for DG method
+//! Upwinding type
+struct RichardsDGUpwinding
+{
+  enum Type {
+    fullUpwind, //!< Upwind flux and conductiviy
+    semiUpwind, //!< Only upwind conductivity factor
+    none //!< No upwinding
+  };
+};
+
+//! Gradient weighting
 struct RichardsDGWeights
 {
-  enum Type { weightsOn, //!< Use saturated conductivities at interface for weighting (default)
-              weightsOff //!< No weighting
-            };
+  enum Type {
+    weightsOn, //!< harmonic weighting of flux contributions
+    weightsOff //!< no weighting of flux contributions
+  };
 };
 
 //! Operator internal BC handling
@@ -40,6 +53,52 @@ struct BCType
             };
 };
 
+/// Read experimental operator settings from a parameter file
+/* \return Tuple of RichardsDGMethod, RichardsDGUpwinding, and
+ *    RichardsDGWeights to be inserted into the local operator constructor.
+ */
+inline auto read_experimental_operator_settings(
+  const Dune::ParameterTree &inifile)
+    -> std::tuple<Operator::RichardsDGMethod::Type,
+                  Operator::RichardsDGUpwinding::Type,
+                  Operator::RichardsDGWeights::Type>
+{
+  RichardsDGMethod::Type method;
+  RichardsDGUpwinding::Type upwinding;
+  RichardsDGWeights::Type weights;
+
+  const auto method_str = inifile.get<std::string>("dg.experimental.method");
+  if (method_str == "SIPG")
+    method = RichardsDGMethod::SIPG;
+  else if (method_str == "NIPG")
+    method = RichardsDGMethod::NIPG;
+  else if (method_str == "OOB")
+    method = RichardsDGMethod::OOB;
+  else if (method_str == "IIP")
+    method = RichardsDGMethod::IIP;
+  else
+    DUNE_THROW(Dune::IOError, "Unknown DG method '" + method_str + "'!");
+
+  const auto upwinding_str
+    = inifile.get<std::string>("dg.experimental.upwinding");
+  if (upwinding_str == "none")
+    upwinding = RichardsDGUpwinding::none;
+  else if (upwinding_str == "semiUpwind")
+    upwinding = RichardsDGUpwinding::semiUpwind;
+  else if (upwinding_str == "fullUpwind")
+    upwinding = RichardsDGUpwinding::fullUpwind;
+  else
+    DUNE_THROW(Dune::IOError, "Unknown upwinding '" + upwinding_str + "'!");
+
+  const auto weights_str = inifile.get<bool>("dg.experimental.weights");
+  if (weights_str)
+    weights = RichardsDGWeights::weightsOn;
+  else
+    weights = RichardsDGWeights::weightsOff;
+
+  return std::make_tuple(method, upwinding, weights);
+}
+
 /// A local operator for solving the Richards equation using upwinding.
 /** Written by Ole Klein, IWR, Heidelberg and modified by the
  *  DORiE developers, IUP, Heidelberg \n
@@ -84,14 +143,15 @@ protected:
   const Boundary&   boundary; //!< class providing boundary conditions
   const SourceTerm& sourceTerm; //!< class providing source term information
   const RichardsDGMethod::Type method; //!< Interior penalty type
-  const RichardsDGWeights::Type weights; //!< DG method weights switch
-  const RF penalty_factor; //!< Interior penalty factor
+  const RichardsDGUpwinding::Type upwinding; //!< DG method weights switch
+  const RichardsDGWeights::Type weights; //!< Gradient weighting
   const Mapper mapper; //!< maps grid elements (cells) to natural numbers
 
   const int intorderadd; //!< integration order addend
   const int quadrature_factor; //!< factor to FEM integration order
-  RF time; //!< operator internal time for BC, parameter, and source term queries
-  RF theta; //!< Penalty factor
+  RF penalty_factor; //!< Interior penalty factor
+  RF theta; //!< Symmetry factor
+  RF time;  //!< operator internal time for BC, parameter, and source term queries
   //! Cache for fixing the BC type over one entire time step
   mutable std::map<typename std::size_t,BoundaryCondition::Type> bc_type_cache;
 
@@ -115,28 +175,35 @@ public:
    *  \param view_ Grid view
    *  \param boundary_ Object prividing boundary conditions
    *  \param sourceTerm_ Object providing source term information
-   *  \param method_ DG Penalty type \see RichardsDGMethod::Type
+   *  \param method_ DG method type \see RichardsDGMethod::Type
+   *  \param upwinding_ Upwinding type \see RichardsDGUpwinding::Type
    *  \param weights_ DG weigths switch \see RichardsDGWeights::Type
    *  \param intorderadd_ Addend to integration order
    *  \param quadrature_factor_ Factor to FEM integration order
    */
   RichardsDGSpatialOperator (const Dune::ParameterTree& config, const Parameter& param_, const GridView view_,
       Boundary& boundary_, const SourceTerm& sourceTerm_,
-      //RichardsDGMethod::Type method_ = RichardsDGMethod::NIPG,
       RichardsDGMethod::Type method_ = RichardsDGMethod::SIPG,
+      RichardsDGUpwinding::Type upwinding_ = RichardsDGUpwinding::none,
       RichardsDGWeights::Type weights_ = RichardsDGWeights::weightsOn,
       int intorderadd_ = 2, int quadrature_factor_ = 2)
     : Dune::PDELab::NumericalJacobianVolume<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
     Dune::PDELab::NumericalJacobianSkeleton<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
     Dune::PDELab::NumericalJacobianBoundary<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
-    param(param_), boundary(boundary_), sourceTerm(sourceTerm_), method(method_), weights(weights_),
-    penalty_factor(config.get<RF>("dg.penaltyFactor")), mapper(view_),
+    param(param_), boundary(boundary_), sourceTerm(sourceTerm_), method(method_), upwinding(upwinding_), weights(weights_), mapper(view_),
     intorderadd(intorderadd_), quadrature_factor(quadrature_factor_),
+    penalty_factor(config.get<RF>("dg.penaltyFactor")),
+    time(0.0),
     cache(20)
   {
-    theta = 1.;
-    if (method == RichardsDGMethod::SIPG) theta = -1.;
-    //theta = 0.;
+    theta = 0.; // IIP
+    if (method == RichardsDGMethod::NIPG || method == RichardsDGMethod::OOB)
+      theta = 1.;
+    else if (method == RichardsDGMethod::SIPG)
+      theta = -1.;
+    
+    if (method == RichardsDGMethod::OOB)
+      penalty_factor = 0.;
   }
 
   /// Volume integral depending on test and ansatz functions
@@ -306,16 +373,36 @@ public:
         // compute jump in solution
         const RF jump = u_s - u_n;
 
-        // conductivity at saturation
+        // saturated conductivity (including Miller scaling)
         const RF satCond_s = param.condRefToMiller(param.K(p_global_s), p_global_s);
         const RF satCond_n = param.condRefToMiller(param.K(p_global_n), p_global_n);
 
+        // conductivity factor
+        const RF head_s = param.headMillerToRef(u_s, p_global_s);
+        const RF saturation_s = param.saturation(head_s, p_global_s);
+        const RF relCond_s = param.condFactor(saturation_s, p_global_s);
+
+        const RF head_n = param.headMillerToRef(u_n, p_global_n);
+        const RF saturation_n = param.saturation(head_n, p_global_n);
+        const RF relCond_n = param.condFactor(saturation_n, p_global_n);
+
         // compute weights
         RF omega_s = 0.5, omega_n = 0.5;
-        if(weights == RichardsDGWeights::weightsOn)
+        if (weights == RichardsDGWeights::weightsOn)
         {
-          omega_s = satCond_n / (satCond_s + satCond_n);
-          omega_n = satCond_s / (satCond_s + satCond_n);
+          if (upwinding == RichardsDGUpwinding::none)
+          {
+            omega_s = relCond_n * satCond_n
+              / ((relCond_s * satCond_s) + (relCond_n * satCond_n));
+            omega_n = relCond_s * satCond_s
+              / ((relCond_s * satCond_s) + (relCond_n * satCond_n));
+          }
+          else if (upwinding == RichardsDGUpwinding::semiUpwind
+            || upwinding == RichardsDGUpwinding::fullUpwind)
+          {
+            omega_s = satCond_n / (satCond_s + satCond_n);
+            omega_n = satCond_s / (satCond_s + satCond_n);
+          }
         }
 
         // integration factor
@@ -324,40 +411,88 @@ public:
         const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
 
         // compute numerical flux
-        const RF numFlux_s = satCond_s * ( - (gradu_s * normal) + penalty * jump);
-        const RF numFlux_n = satCond_n * ( - (gradu_n * normal) + penalty * jump);
-        const RF numFlux   = omega_s * numFlux_s + omega_n * numFlux_n;
-
-        // compute upwind conductivity factor
-        RF relCond_upwind;
-        if (numFlux > 0) {
-          const RF head_upwind = param.headMillerToRef(u_s, p_global_s);
-          const RF saturation_upwind = param.saturation(head_upwind, p_global_s);
-          relCond_upwind = param.condFactor(saturation_upwind, p_global_s);
+        RF numFlux_s = satCond_s * ( - (gradu_s * normal) + penalty * jump);
+        RF numFlux_n = satCond_n * ( - (gradu_n * normal) + penalty * jump);
+        if (upwinding == RichardsDGUpwinding::none)
+        {
+          numFlux_s *= relCond_s;
+          numFlux_n *= relCond_n;
         }
-        else {
-          const RF head_upwind = param.headMillerToRef(u_n, p_global_n);
-          const RF saturation_upwind = param.saturation(head_upwind, p_global_n);
-          relCond_upwind = param.condFactor(saturation_upwind, p_global_n);
+        const RF numFlux = omega_s * numFlux_s + omega_n * numFlux_n;
+
+        // compute upwind quantities
+        RF relCond_upwind, satCond_upwind, numFlux_upwind;
+        if (upwinding == RichardsDGUpwinding::semiUpwind
+          || upwinding == RichardsDGUpwinding::fullUpwind)
+        {
+          if (numFlux > 0)
+          {
+            const RF head_upwind = param.headMillerToRef(u_s, p_global_s);
+            const RF saturation_upwind = param.saturation(head_upwind, p_global_s);
+            relCond_upwind = param.condFactor(saturation_upwind, p_global_s);
+            satCond_upwind = satCond_s;
+            numFlux_upwind = numFlux_s;
+          }
+          else
+          {
+            const RF head_upwind = param.headMillerToRef(u_n, p_global_n);
+            const RF saturation_upwind = param.saturation(head_upwind, p_global_n);
+            relCond_upwind = param.condFactor(saturation_upwind, p_global_n);
+            satCond_upwind = satCond_n;
+            numFlux_upwind = numFlux_n;
+          }
         }
 
         // update residual (flux term)
         // diffusion term
         // consistency term
         // + penalty term
-        for (unsigned int i = 0; i<lfsv_s.size(); i++)
-          r_s.accumulate(lfsv_s,i,   relCond_upwind * numFlux * phi_s[i] * factor);
-        for (unsigned int i = 0; i<lfsv_n.size(); i++)
-          r_n.accumulate(lfsv_n,i, - relCond_upwind * numFlux * phi_n[i] * factor);
+        if (upwinding == RichardsDGUpwinding::none)
+        {
+          for (unsigned int i = 0; i < lfsv_s.size(); i++)
+            r_s.accumulate(lfsv_s, i,   numFlux * phi_s[i] * factor);
+          for (unsigned int i = 0; i < lfsv_n.size(); i++)
+            r_n.accumulate(lfsv_n, i, - numFlux * phi_n[i] * factor);
+        }
+        else if (upwinding == RichardsDGUpwinding::semiUpwind)
+        {
+          for (unsigned int i = 0; i<lfsv_s.size(); i++)
+            r_s.accumulate(lfsv_s,i,   relCond_upwind * numFlux * phi_s[i] * factor);
+          for (unsigned int i = 0; i<lfsv_n.size(); i++)
+            r_n.accumulate(lfsv_n,i, - relCond_upwind * numFlux * phi_n[i] * factor);
+        }
+        else if (upwinding == RichardsDGUpwinding::fullUpwind)
+        {
+          for (unsigned int i = 0; i < lfsv_s.size(); i++)
+            r_s.accumulate(lfsv_s, i,   relCond_upwind * numFlux_upwind * phi_s[i] * factor);
+          for (unsigned int i = 0; i < lfsv_n.size(); i++)
+            r_n.accumulate(lfsv_n, i, - relCond_upwind * numFlux_upwind * phi_n[i] * factor);
+        }
 
         // update residual (symmetry term)
         // (non-)symmetric IP term
         // symmetry term
-        for (unsigned int i = 0; i<lfsv_s.size(); i++)
-          r_s.accumulate(lfsv_s,i, theta * omega_s * (tgradphi_s[i] * normal) * relCond_upwind * satCond_s * jump * factor);
-        for (unsigned int i = 0; i<lfsv_n.size(); i++)
-          r_n.accumulate(lfsv_n,i, theta * omega_n * (tgradphi_n[i] * normal) * relCond_upwind * satCond_n * jump * factor);
-
+        if (upwinding == RichardsDGUpwinding::none)
+        {
+          for (unsigned int i = 0; i < lfsv_s.size(); i++)
+            r_s.accumulate(lfsv_s, i, theta * omega_s * (tgradphi_s[i] * normal) * relCond_s * satCond_s * jump * factor);
+          for (unsigned int i = 0; i < lfsv_n.size(); i++)
+            r_n.accumulate(lfsv_n, i, theta * omega_n * (tgradphi_n[i] * normal) * relCond_n * satCond_n * jump * factor);
+        }
+        else if (upwinding == RichardsDGUpwinding::semiUpwind)
+        {
+          for (unsigned int i = 0; i<lfsv_s.size(); i++)
+            r_s.accumulate(lfsv_s,i, theta * omega_s * (tgradphi_s[i] * normal) * relCond_upwind * satCond_s * jump * factor);
+          for (unsigned int i = 0; i<lfsv_n.size(); i++)
+            r_n.accumulate(lfsv_n,i, theta * omega_n * (tgradphi_n[i] * normal) * relCond_upwind * satCond_n * jump * factor);
+        }
+        else if (upwinding == RichardsDGUpwinding::fullUpwind)
+        {
+          for (unsigned int i = 0; i < lfsv_s.size(); i++)
+            r_s.accumulate(lfsv_s, i, theta * (tgradphi_s[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
+          for (unsigned int i = 0; i < lfsv_n.size(); i++)
+            r_n.accumulate(lfsv_n, i, theta * (tgradphi_n[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
+        }
       }
     }
 
@@ -469,27 +604,40 @@ public:
           // conductivity at saturation
           const RF satCond_s = param.condRefToMiller(param.K(p_global_s), p_global_s);
 
-          const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
+          // conductivity factor
+          const RF head_s = param.headMillerToRef(u_s, p_global_s);
+          const RF saturation_s = param.saturation(head_s, p_global_s);
+          const RF relCond_s = param.condFactor(saturation_s, p_global_s);
 
-          // compute numerical flux estimation
-          const RF numFlux = satCond_s * ( - (gradu_s * normal) + penalty * jump);
+          const RF saturation_n = param.saturation(g, p_global_s);
+          const RF relCond_n = param.condFactor(saturation_n, p_global_s);
 
-          RF head_upwind;
-          if (numFlux > 0)
-            head_upwind = param.headMillerToRef(u_s, p_global_s);
-          else
-            head_upwind = g; // Miller here?!
+          // penalty
+          const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
 
-          // compute saturation from matrix head
-          const RF saturation_upwind = param.saturation(head_upwind, p_global_s);
+          // compute numerical flux
+          const RF numFlux = satCond_s * (-(gradu_s * normal) + penalty * jump);
 
-          // compute relative conductivity
-          const RF relCond_upwind = param.condFactor(saturation_upwind, p_global_s);
+          // compute conductivity factor
+          RF relCond;
+          if (upwinding == RichardsDGUpwinding::semiUpwind
+            || upwinding == RichardsDGUpwinding::fullUpwind)
+          {
+            if (numFlux > 0)
+              relCond = relCond_s;
+            else
+              relCond = relCond_n;
+          }
+          else if (upwinding == RichardsDGUpwinding::none)
+          {
+            // harmonic average of conductivity factors
+            relCond = 2 * relCond_s * relCond_n / (relCond_s + relCond_n);
+          }
 
           // switches between Neumann and Dirichlet bc
           // this choice is kept constant for one time step
           normal_flux = boundary.j(ig.intersection(),it.position(),time);
-          if(relCond_upwind * numFlux < normal_flux) {
+          if(relCond * numFlux < normal_flux) {
             bcType = BCType::dirichlet;
             bc_type_cache.emplace(it_id,BoundaryCondition::Dirichlet);
           }
@@ -547,35 +695,48 @@ public:
           // conductivity at saturation
           const RF satCond_s = param.condRefToMiller(param.K(p_global_s), p_global_s);
 
+          // conductivity factor
+          const RF head_s = param.headMillerToRef(u_s, p_global_s);
+          const RF saturation_s = param.saturation(head_s, p_global_s);
+          const RF relCond_s = param.condFactor(saturation_s, p_global_s);
+
+          const RF saturation_n = param.saturation(g, p_global_s);
+          const RF relCond_n = param.condFactor(saturation_n, p_global_s);
+
+          // penalty
           const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
 
           // compute numerical flux
           const RF numFlux = satCond_s * ( - (gradu_s * normal) + penalty * jump);
 
-          RF head_upwind;
-          if (numFlux > 0)
-            head_upwind = param.headMillerToRef(u_s, p_global_s);
-          else
-            head_upwind = g; // Miller here?!
-
-          // compute saturation from matrix head
-          const RF saturation_upwind = param.saturation(head_upwind, p_global_s);
-
-          // compute relative conductivity
-          const RF relCond_upwind = param.condFactor(saturation_upwind, p_global_s);
+          // compute conductivity factor
+          RF relCond;
+          if (upwinding == RichardsDGUpwinding::semiUpwind
+            || upwinding == RichardsDGUpwinding::fullUpwind)
+          {
+            if (numFlux > 0)
+              relCond = relCond_s;
+            else
+              relCond = relCond_n;
+          }
+          else if (upwinding == RichardsDGUpwinding::none)
+          {
+            // harmonic average of conductivity factors
+            relCond = 2 * relCond_s * relCond_n / (relCond_s + relCond_n);
+          }
 
           // update residual (flux term)
           // diffusion term
           // consistency term
           // + penalty term
           for (unsigned int i = 0; i<lfsv_s.size(); i++)
-            r_s.accumulate(lfsv_s,i, relCond_upwind * numFlux * phi_s[i] * factor);
+            r_s.accumulate(lfsv_s,i, relCond * numFlux * phi_s[i] * factor);
 
           // update residual (symmetry term)
           // (non-)symmetric IP term
           // symmetry term
           for (unsigned int i = 0; i<lfsv_s.size(); i++)
-            r_s.accumulate(lfsv_s,i, theta * (tgradphi_s[i] * normal) * relCond_upwind * satCond_s * jump * factor);
+            r_s.accumulate(lfsv_s,i, theta * (tgradphi_s[i] * normal) * relCond * satCond_s * jump * factor);
 
           continue;
         }

dune/dorie/solver/operator_error_indicator.hh

@@ -12,6 +12,7 @@
 
 namespace Dune {
   namespace Dorie {
+  namespace Operator { 
 
     /// Local operator for residual-based error estimation.
     /*
@@ -48,18 +49,23 @@ namespace Dune {
 
       const Parameter& param;
       const Boundary& boundary;
-      const RF penalty_factor;
       const int intorderadd;
       const int quadrature_factor;
+      RF penalty_factor;
       RF time;
 
       FluxErrorEstimator (const Dune::ParameterTree& config,
         const Parameter& param_, const Boundary& boundary_,
+        RichardsDGMethod::Type 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>("dg.penaltyFactor")),
-        intorderadd(intorderadd_), quadrature_factor(quadrature_factor_)
-      { }
+        time(0.0)
+      {
+        if (method_ == RichardsDGMethod::OOB)
+          penalty_factor = 0.0;
+      }
 
       /// skeleton integral depending on test and ansatz functions
       /** each face is only visited ONCE!
@@ -330,7 +336,8 @@ namespace Dune {
 
     };
 
-  }
-}
+  } // namespace Operator
+  } // namespace Dorie
+} // namespace Dune
 
 #endif // DUNE_DORIE_ERROR_INDICATOR_HH
\ No newline at end of file

dune/dorie/solver/param_van_genuchten.hh

@@ -22,8 +22,8 @@ namespace Dune {
       explicit MualemVanGenuchten(const ParameterTree& fieldConfig_, const Dune::MPIHelper& helper_,