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

(cherry picked from commit dac2a4e1)

d760a799 Add new options for symmetry/penalty in local operator
6cdbdb13 Added option to change upwinding in local operator
99db5435 Merge branch 'bugix/conductivity-in-dirichlet-bc' into feature/local-operator-schemes
57441a3a Consider upwinding schemes in Dirichlet BC
6432d739 softened vanGenuchten parameter limits
7d47d613 Fix an error in the DG skeleton flux term
9128001f Add CMake option to control DG method with the config file
fe7a21a4 Finalize DG scheme options in local operator
5eb13c57 Added docs for experimental DG keys
047b6016 Moved experimental DG parmeters to the end of the parameter list
e09f0b28 Merge branch 'feature/41-enable-vtk-vertex-plotting' into 'master'
5576913a Set default upwinding method to 'none'
ad0ef464 Apply CPP definition for experimental features only to main dorie exe
678a9a75 Add section on experimental features to Readme

CMakeLists.txt

@@ -11,6 +11,12 @@ if(CMAKE_BUILD_TYPE_UPPER MATCHES DEBUG)
 	set(CMAKE_VERBOSE_MAKEFILE ON)
 endif()
 
+# 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

@@ -3,4 +3,9 @@ add_subdirectory(solver)
 add_subdirectory(impl)
 
 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