Consider upwinding schemes in Dirichlet BC

For no upwinding, the harmonic average of theconductivity factors is computed.

dune/dorie/solver/operator_DG.hh

@@ -538,27 +538,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);
           }
@@ -616,35 +629,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;
         }