[mpnc] Fix localResidual

* set phaseNcp in method evalFluxSource as this one should be called
  in all possible cases for assembly (implicit/explicit/numeric/analytic)

dumux/porousmediumflow/mpnc/localresidual.hh

@@ -45,16 +45,10 @@ template<class TypeTag>
 class MPNCLocalResidual : public CompositionalLocalResidual<TypeTag>
 {
     using ParentType = CompositionalLocalResidual<TypeTag>;
-    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
-    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
-    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
-    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
-    using ElementResidualVector = Dune::BlockVector<typename GET_PROP_TYPE(TypeTag, NumEqVector)>;
     using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
     using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
-    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
 
@@ -64,102 +58,22 @@ class MPNCLocalResidual : public CompositionalLocalResidual<TypeTag>
 public:
     using ParentType::ParentType;
 
+    using typename ParentType::ElementResidualVector;
 
-    /*!
-     * \brief Compute the local residual, i.e. the deviation of the
-     *        equations from zero for instationary problems.
-     * \param problem The problem to solve
-
-     * \param element The DUNE Codim<0> entity for which the residual
-     *                ought to be calculated
-     * \param fvGeometry The finite-volume geometry of the element
-     * \param prevVolVars The volume averaged variables for all
-     *                    sub-control volumes of the element at the previous
-     *                    time level
-     * \param curVolVars The volume averaged variables for all
-     *                   sub-control volumes of the element at the current
-     *                   time level
-     * \param bcTypes The types of the boundary conditions for all
-     *                vertices of the element
-     */
-    ElementResidualVector eval(const Problem& problem,
-                               const Element& element,
-                               const FVElementGeometry& fvGeometry,
-                               const ElementVolumeVariables& prevElemVolVars,
-                               const ElementVolumeVariables& curElemVolVars,
-                               const ElementBoundaryTypes &bcTypes,
-                               const ElementFluxVariablesCache& elemFluxVarsCache) const
-    {
-        // initialize the residual vector for all scvs in this element
-        ElementResidualVector residual(fvGeometry.numScv());
-        residual = 0.0;
-
-        // evaluate the volume terms (storage + source terms)
-        for (auto&& scv : scvs(fvGeometry))
-        {
-            //! foward to the local residual specialized for the discretization methods
-            this->asImp().evalStorage(residual, problem, element, fvGeometry, prevElemVolVars, curElemVolVars, scv);
-            this->asImp().evalSource(residual, problem, element, fvGeometry, curElemVolVars, scv);
-             //here we need to set the constraints of the mpnc model into the residual
-                const auto localScvIdx = scv.indexInElement();
-                for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-                {
-                    residual[localScvIdx][phase0NcpIdx + phaseIdx] = curElemVolVars[scv].phaseNcp(phaseIdx);
-                }
-        }
-
-        for (auto&& scvf : scvfs(fvGeometry))
-        {
-            //! foward to the local residual specialized for the discretization methods
-            this->asImp().evalFlux(residual, problem, element, fvGeometry, curElemVolVars, bcTypes, elemFluxVarsCache, scvf);
-        }
-
-        return residual;
-    }
-
-    /*!
-     * \brief Compute the local residual, i.e. the deviation of the
-     *        equations from zero for stationary problem.
-     * \param problem The problem to solve
-
-     * \param element The DUNE Codim<0> entity for which the residual
-     *                ought to be calculated
-     * \param fvGeometry The finite-volume geometry of the element
-     * \param curVolVars The volume averaged variables for all
-     *                   sub-control volumes of the element at the current
-     *                   time level
-     * \param bcTypes The types of the boundary conditions for all
-     *                vertices of the element
-     */
-    ElementResidualVector eval(const Problem& problem,
-                               const Element& element,
-                               const FVElementGeometry& fvGeometry,
-                               const ElementVolumeVariables& curElemVolVars,
-                               const ElementBoundaryTypes &bcTypes,
-                               const ElementFluxVariablesCache& elemFluxVarsCache) const
+    ElementResidualVector evalFluxSource(const Element& element,
+                                         const FVElementGeometry& fvGeometry,
+                                         const ElementVolumeVariables& elemVolVars,
+                                         const ElementFluxVariablesCache& elemFluxVarsCache,
+                                         const ElementBoundaryTypes &bcTypes) const
     {
-        // initialize the residual vector for all scvs in this element
-        ElementResidualVector residual(fvGeometry.numScv());
-        residual = 0.0;
+        ElementResidualVector residual = ParentType::evalFluxSource(element, fvGeometry, elemVolVars, elemFluxVarsCache, bcTypes);
 
-        // evaluate the volume terms (storage + source terms)
         for (auto&& scv : scvs(fvGeometry))
         {
-            //! foward to the local residual specialized for the discretization methods
-            this->asImp().evalSource(residual, problem, element, fvGeometry, curElemVolVars, scv);
-
-               //here we need to set the constraints of the mpnc model into the residual
-                const auto localScvIdx = scv.indexInElement();
-                for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-                {
-                    residual[localScvIdx][phase0NcpIdx + phaseIdx] = curElemVolVars[scv].phaseNcp(phaseIdx);
-                }
-        }
-
-        for (auto&& scvf : scvfs(fvGeometry))
-        {
-            //! foward to the local residual specialized for the discretization methods
-            this->asImp().evalFlux(residual, problem, element, fvGeometry, curElemVolVars, bcTypes, elemFluxVarsCache, scvf);
+            //here we need to set the constraints of the mpnc model into the residual
+            const auto localScvIdx = scv.indexInElement();
+            for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                residual[localScvIdx][phase0NcpIdx + phaseIdx] = elemVolVars[scv].phaseNcp(phaseIdx);
         }
 
         return residual;