From 71e3e763c2c7c87caeb0ec5014e3ebc2b8eaedfb Mon Sep 17 00:00:00 2001
From: DennisGlaeser <dennis.glaeser@iws.uni-stuttgart.de>
Date: Tue, 5 Jul 2016 11:55:32 +0200
Subject: [PATCH] [cc][localres] change boundary treatment

the boundary treatment is now performed after all the interior
fluxes have been calculated. This corresponds more to the previous
implementation and is assumend to be more correct for mixed BCs.
---
 dumux/implicit/cellcentered/localresidual.hh | 47 ++++++++++++++------
 1 file changed, 33 insertions(+), 14 deletions(-)

diff --git a/dumux/implicit/cellcentered/localresidual.hh b/dumux/implicit/cellcentered/localresidual.hh
index c9e0b970b9..0d04f04f67 100644
--- a/dumux/implicit/cellcentered/localresidual.hh
+++ b/dumux/implicit/cellcentered/localresidual.hh
@@ -53,7 +53,7 @@ class CCLocalResidual : public ImplicitLocalResidual<TypeTag>
 
     enum { constantBC = GET_PROP_VALUE(TypeTag, ConstantBoundaryConditions) };
 
-    enum { enableVolVarsCache = GET_PROP_VALUE(TypeTag, EnableVolumeVariablesCache) };
+    enum { enableVolVarsCache = GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache) };
 
     typedef typename GridView::template Codim<0>::Entity Element;
     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
@@ -70,23 +70,44 @@ public:
 
 protected:
 
+    void evalFluxes_()
+    {
+        // calculate the mass flux over the scv faces
+        const auto& fvGeometry = this->fvGeometry_();
+        for (const auto& scvf : fvGeometry.scvfs())
+        {
+            this->residual_[0] += this->asImp_().computeFlux_(scvf);
+        }
+    }
+
     PrimaryVariables computeFlux_(const SubControlVolumeFace &scvf)
     {
         if (!scvf.boundary() /*TODO: || GET_PROP_VALUE(TypeTag, BoundaryReconstruction)*/)
             return this->asImp_().computeFlux(scvf);
         else
-            return this->asImp_().evalBoundary_(scvf);
+            return PrimaryVariables(0.0);
+
     }
 
-    PrimaryVariables evalBoundary_(const SubControlVolumeFace &scvf)
+    void evalBoundary_()
     {
-        PrimaryVariables flux = evalBoundaryFluxes_(scvf);
+        const auto& fvGeometry = this->fvGeometry_();
+
+        for (const auto& scvf : fvGeometry.scvfs())
+        {
+            if (scvf.boundary())
+                this->residual_[0] += evalBoundaryFluxes_(scvf);
+        }
 
         // additionally treat mixed D/N conditions in a strong sense
         if (this->bcTypes_().hasDirichlet())
-            this->asImp_().evalDirichlet_(flux, scvf);
-
-        return flux;
+        {
+            for (const auto& scvf : fvGeometry.scvfs())
+            {
+                if (scvf.boundary())
+                    this->asImp_().evalDirichlet_(scvf);
+            }
+        }
     }
 
     /*!
@@ -99,11 +120,11 @@ protected:
 
         // evaluate the Neumann conditions at the boundary face
         PrimaryVariables flux(0);
-        if (bcTypes.hasNeumann())
+        if (bcTypes.hasNeumann() /*TODO: && !GET_PROP_VALUE(TypeTag, BoundaryReconstruction)*/)
             flux += this->asImp_().evalNeumannSegment_(scvf, bcTypes);
 
         // TODO: evaluate the outflow conditions at the boundary face
-        //if (bcTypes.hasOutflow())
+        //if (bcTypes.hasOutflow() /*TODO: && !GET_PROP_VALUE(TypeTag, BoundaryReconstruction)*/)
         //    flux += this->asImp_().evalOutflowSegment_(&intersection, bcTypes);
 
         // evaluate the pure Dirichlet conditions at the boundary face
@@ -117,12 +138,12 @@ protected:
      * \brief Evaluate Dirichlet conditions on faces that have mixed
      *        Dirichlet/Neumann boundary conditions.
      */
-    void evalDirichlet_(PrimaryVariables &flux, const SubControlVolumeFace &scvf)
+    void evalDirichlet_(const SubControlVolumeFace &scvf)
     {
         BoundaryTypes bcTypes = this->problem_().boundaryTypes(this->element_(), scvf);
 
         if (bcTypes.hasDirichlet() && bcTypes.hasNeumann())
-            this->asImp_().evalDirichletSegmentMixed_(flux, scvf, bcTypes);
+            this->asImp_().evalDirichletSegmentMixed_(scvf, bcTypes);
     }
 
     /*!
@@ -240,8 +261,7 @@ protected:
      * \brief Treat Dirichlet boundary conditions in a strong sense for a
      *        single intersection that has mixed D/N boundary conditions
      */
-    void evalDirichletSegmentMixed_(PrimaryVariables &flux,
-                                    const SubControlVolumeFace &scvf,
+    void evalDirichletSegmentMixed_(const SubControlVolumeFace &scvf,
                                     const BoundaryTypes &bcTypes)
     {
         // temporary vector to store the Dirichlet values
@@ -257,7 +277,6 @@ protected:
             {
                 auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
                 this->residual_[0][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
-                flux[eqIdx] = 0;
             }
         }
     }
-- 
GitLab