From d22f7c0b9971773c2801046b15ec8876a7e6cd74 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Fri, 12 Jan 2018 21:37:20 +0100
Subject: [PATCH] [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 | 108 ++-----------------
1 file changed, 11 insertions(+), 97 deletions(-)
diff --git a/dumux/porousmediumflow/mpnc/localresidual.hh b/dumux/porousmediumflow/mpnc/localresidual.hh
index 5bf07dc6ea..d982de82f8 100644
--- a/dumux/porousmediumflow/mpnc/localresidual.hh
+++ b/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;
--
GitLab