Added option to change upwinding in local operator

**Note: This only affects the skeleton not the boundary term!**
Upwinding can now be turned off (`none`)
or set to full upwinding (`fullUpwind`).
Removed 'weights' option.

dune/dorie/solver/operator_DG.hh

@@ -25,12 +25,14 @@ struct RichardsDGMethod
             };
 };
 
-//! Weighting at grid interface for DG method
-struct RichardsDGWeights
+//! Upwinding type
+struct RichardsDGUpwinding
 {
-  enum Type { weightsOn, //!< Use saturated conductivities at interface for weighting (default)
-              weightsOff //!< No weighting
-            };
+  enum Type {
+    fullUpwind, //!< Upwind flux and conductiviy
+    semiUpwind, //!< Only upwind conductivity factor
+    none //!< No upwinding
+  };
 };
 
 //! Operator internal BC handling
@@ -85,8 +87,8 @@ protected:
   const Parameter&  param; //!< class providing equation parameters
   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
+  RichardsDGMethod::Type method; //!< Interior penalty type
+  RichardsDGUpwinding::Type upwinding; //!< DG method weights switch
   const Mapper mapper; //!< maps grid elements (cells) to natural numbers
 
   const int intorderadd; //!< integration order addend
@@ -125,12 +127,12 @@ public:
   RichardsDGSpatialOperator (const Dune::ParameterTree& config, const Parameter& param_, const GridView view_,
       Boundary& boundary_, const SourceTerm& sourceTerm_,
       RichardsDGMethod::Type method_ = RichardsDGMethod::SIPG,
-      RichardsDGWeights::Type weights_ = RichardsDGWeights::weightsOn,
+      RichardsDGUpwinding::Type upwinding_ = RichardsDGUpwinding::semiUpwind,
       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_), mapper(view_),
+    param(param_), boundary(boundary_), sourceTerm(sourceTerm_), method(method_), upwinding(upwinding_), mapper(view_),
     intorderadd(intorderadd_), quadrature_factor(quadrature_factor_),
     penalty_factor(config.get<RF>("dg.penaltyFactor")),
     time(0.0),
@@ -313,13 +315,30 @@ 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 (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)
         {
           omega_s = satCond_n / (satCond_s + satCond_n);
           omega_n = satCond_s / (satCond_s + satCond_n);
@@ -331,21 +350,47 @@ 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;
+        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)
+        {
+          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;
-        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 relCond_upwind = 1.0;
+        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);
+          }
+          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);
+          }
         }
-        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 = numFlux;
+        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;
+          }
         }
 
         // update residual (flux term)
@@ -353,17 +398,34 @@ public:
         // 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_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 * phi_n[i] * factor);
+          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);
+        }
 
       }
     }