[next] Implement eval flux method return a single flux residual

dumux/implicit/cellcentered/localassembler.hh

@@ -147,6 +147,7 @@ private:
         const auto globalI = fvGridGeometry.elementMapper().index(element);
 
         // check for boundaries on the element
+        // TODO Do we need them for cell-centered models?
         ElementBoundaryTypes elemBcTypes;
         elemBcTypes.update(problem, element, fvGeometry);
 
@@ -194,15 +195,12 @@ private:
             neighborElements.emplace_back(fvGridGeometry.element(dataJ.globalJ));
             for (const auto scvfIdx : dataJ.scvfsJ)
             {
-                ElementSolutionVector flux(1);
-                localResidual.evalFlux(flux, problem,
-                                       neighborElements.back(),
-                                       fvGeometry,
-                                       curElemVolVars,
-                                       elemBcTypes,
-                                       elemFluxVarsCache,
-                                       fvGeometry.scvf(scvfIdx));
-                origFlux[j] += flux[0];
+                origFlux[j] += localResidual.evalFlux(problem,
+                                                      neighborElements.back(),
+                                                      fvGeometry,
+                                                      curElemVolVars,
+                                                      elemFluxVarsCache,
+                                                      fvGeometry.scvf(scvfIdx));
             }
             // increment neighbor counter
             ++j;
@@ -262,15 +260,12 @@ private:
                 for (std::size_t k = 0; k < numNeighbors; ++k)
                     for (auto scvfIdx : connectivityMap[globalI][k].scvfsJ)
                     {
-                        ElementSolutionVector flux(1);
-                        localResidual.evalFlux(flux, problem,
-                                               neighborElements[k],
-                                               fvGeometry,
-                                               curElemVolVars,
-                                               elemBcTypes,
-                                               elemFluxVarsCache,
-                                               fvGeometry.scvf(scvfIdx));
-                        neighborDeriv[k] += flux[0];
+                        neighborDeriv[k] += localResidual.evalFlux(problem,
+                                                                   neighborElements[k],
+                                                                   fvGeometry,
+                                                                   curElemVolVars,
+                                                                   elemFluxVarsCache,
+                                                                   fvGeometry.scvf(scvfIdx));
                     }
             }
             else
@@ -310,15 +305,12 @@ private:
                 for (std::size_t k = 0; k < numNeighbors; ++k)
                     for (auto scvfIdx : connectivityMap[globalI][k].scvfsJ)
                     {
-                        ElementSolutionVector flux(1);
-                        localResidual.evalFlux(flux, problem,
-                                               neighborElements[k],
-                                               fvGeometry,
-                                               curElemVolVars,
-                                               elemBcTypes,
-                                               elemFluxVarsCache,
-                                               fvGeometry.scvf(scvfIdx));
-                        neighborDeriv[k] += flux[0];
+                        neighborDeriv[k] += localResidual.evalFlux(problem,
+                                                                   neighborElements[k],
+                                                                   fvGeometry,
+                                                                   curElemVolVars,
+                                                                   elemFluxVarsCache,
+                                                                   fvGeometry.scvf(scvfIdx));
                     }
             }
             else

dumux/implicit/cellcentered/localresidual.hh

@@ -47,6 +47,7 @@ class CCLocalResidual : public ImplicitLocalResidual<TypeTag>
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
     using ElementResidualVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
@@ -66,11 +67,22 @@ public:
     {
         const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
         const auto localScvIdx = scv.indexInElement();
+        residual[localScvIdx] += evalFlux(problem, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
+    }
+
+    ResidualVector evalFlux(const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const ElementFluxVariablesCache& elemFluxVarsCache,
+                            const SubControlVolumeFace& scvf) const
+    {
+        ResidualVector flux(0.0);
 
         // inner faces
         if (!scvf.boundary())
         {
-            residual[localScvIdx] += this->asImp_().computeFlux(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+            flux += this->asImp_().computeFlux(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
         }
 
         // boundary faces
@@ -80,7 +92,7 @@ public:
 
             // Dirichlet boundaries
             if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann())
-                residual[localScvIdx] += this->asImp_().computeFlux(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+                flux += this->asImp_().computeFlux(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
 
             // Neumann and Robin ("solution dependent Neumann") boundary conditions
             else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet())
@@ -88,15 +100,17 @@ public:
                 auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, scvf);
 
                 // multiply neumann fluxes with the area and the extrusion factor
+                const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
                 neumannFluxes *= scvf.area()*elemVolVars[scv].extrusionFactor();
 
-                residual[localScvIdx] += neumannFluxes;
+                flux += neumannFluxes;
             }
 
             else
                 DUNE_THROW(Dune::NotImplemented, "Mixed boundary conditions. Use pure boundary conditions by converting Dirichlet BCs to Robin BCs");
         }
 
+        return flux;
     }
 };
 

dumux/implicit/localresidual.hh

@@ -321,6 +321,13 @@ public:
                   const ElementFluxVariablesCache& elemFluxVarsCache,
                   const SubControlVolumeFace& scvf) const {}
 
+    ResidualVector evalFlux(const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const ElementFluxVariablesCache& elemFluxVarsCache,
+                            const SubControlVolumeFace& scvf) const {}
+
     void evalBoundary(ElementResidualVector& residual,
                       const Problem& problem,
                       const Element& element,