Finalize DG scheme options in local operator

- fixed full upwdining: only use satCond for weighting
- re-implemented option to disable weights in skeleton term
- moved experimental settings function to Operator namespace

dune/dorie/interface/simulation.cc

@@ -32,9 +32,15 @@ Simulation<Traits>::Simulation (Dune::MPIHelper& _helper, std::shared_ptr<Grid>
 
 	// --- Local Operators ---
 #ifdef EXPERIMENTAL_DG_FEATURES
-	const auto settings = read_experimental_operator_settings(inifile);
+	// 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,
-																settings.first, settings.second);
+																method, upwinding, weights);
 #else
 	slop = std::make_unique<SLOP>(inifile,*param,gv,*fboundary,*fsource);
 #endif // EXPERIMENTAL_DG_FEATURES
@@ -135,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
@@ -147,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/interface/simulation.hh

@@ -127,46 +127,6 @@ protected:
 
 	const int verbose;
 
-private:
-
-	/// Read experimental operator settings from a parameter file
-	/* \return Pair of RichardsDGMethod and RichardsDGUpwinding
-	 * 		to be inserted into the local operator constructor.
-	 */
-	auto read_experimental_operator_settings (const Dune::ParameterTree& inifile)
-	-> std::pair<Operator::RichardsDGMethod::Type,
-							Operator::RichardsDGUpwinding::Type>
-	{
-		using namespace Operator;
-		RichardsDGMethod::Type method;
-		RichardsDGUpwinding::Type upwinding;
-
-		const auto method_str = inifile.get<std::string>("dg.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.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 upwdinding '" + upwinding_str + "'!");
-		
-		return std::make_pair(method, upwinding);
-	}	
-	
-
 public:
 
 	/// Execute the simulation until tEnd is reached.

dune/dorie/solver/operator_DG.hh

@@ -35,6 +35,15 @@ struct RichardsDGUpwinding
   };
 };
 
+//! Gradient weighting
+struct RichardsDGWeights
+{
+  enum Type {
+    weightsOn, //!< harmonic weighting of flux contributions
+    weightsOff //!< no weighting of flux contributions
+  };
+};
+
 //! Operator internal BC handling
 struct BCType
 {
@@ -44,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
@@ -87,8 +142,9 @@ protected:
   const Parameter&  param; //!< class providing equation parameters
   const Boundary&   boundary; //!< class providing boundary conditions
   const SourceTerm& sourceTerm; //!< class providing source term information
-  RichardsDGMethod::Type method; //!< Interior penalty type
-  RichardsDGUpwinding::Type upwinding; //!< DG method weights switch
+  const RichardsDGMethod::Type method; //!< Interior penalty type
+  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
@@ -120,6 +176,7 @@ public:
    *  \param boundary_ Object prividing boundary conditions
    *  \param sourceTerm_ Object providing source term information
    *  \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
@@ -128,11 +185,12 @@ public:
       Boundary& boundary_, const SourceTerm& sourceTerm_,
       RichardsDGMethod::Type method_ = RichardsDGMethod::SIPG,
       RichardsDGUpwinding::Type upwinding_ = RichardsDGUpwinding::semiUpwind,
+      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_), upwinding(upwinding_), 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),
@@ -330,18 +388,21 @@ public:
 
         // compute weights
         RF omega_s = 0.5, omega_n = 0.5;
-        if (upwinding == RichardsDGUpwinding::none
-          || upwinding == RichardsDGUpwinding::fullUpwind)
-        {
-          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)
+        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
@@ -352,16 +413,15 @@ public:
         // compute numerical flux
         RF numFlux_s = satCond_s * ( - (gradu_s * normal) + penalty * jump);
         RF numFlux_n = satCond_n * ( - (gradu_n * normal) + penalty * jump);
-        if (upwinding == RichardsDGUpwinding::none
-          || upwinding == RichardsDGUpwinding::fullUpwind)
+        if (upwinding == RichardsDGUpwinding::none)
         {
           numFlux_s *= relCond_s;
           numFlux_n *= relCond_n;
         }
         const RF numFlux = omega_s * numFlux_s + omega_n * numFlux_n;
 
-        // compute upwind conductivity factor
-        RF relCond_upwind = 1.0;
+        // compute upwind quantities
+        RF relCond_upwind, satCond_upwind, numFlux_upwind;
         if (upwinding == RichardsDGUpwinding::semiUpwind
           || upwinding == RichardsDGUpwinding::fullUpwind)
         {
@@ -370,26 +430,16 @@ public:
             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);
-          }
-        }
-
-        // compute upwind flux
-        RF numFlux_upwind;
-        RF satCond_upwind;
-        if (upwinding == RichardsDGUpwinding::fullUpwind) {
-          if (numFlux > 0) {
-            numFlux_upwind = numFlux_s;
-            satCond_upwind = satCond_s;
-          }
-          else {
-            numFlux_upwind = numFlux_n;
             satCond_upwind = satCond_n;
+            numFlux_upwind = numFlux_n;
           }
         }
 
@@ -414,9 +464,9 @@ public:
         else if (upwinding == RichardsDGUpwinding::fullUpwind)
         {
           for (unsigned int i = 0; i < lfsv_s.size(); i++)
-            r_s.accumulate(lfsv_s, i,   numFlux_upwind * phi_s[i] * factor);
+            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, - numFlux_upwind * phi_n[i] * factor);
+            r_n.accumulate(lfsv_n, i, - relCond_upwind * numFlux_upwind * phi_n[i] * factor);
         }
 
         // update residual (symmetry term)