diff --git a/dumux/porousmediumflow/3p3c/implicit/localresidual.hh b/dumux/porousmediumflow/3p3c/implicit/localresidual.hh index 1a766f745793ed9f5383a5e615439aa5324a0dc8..bd1712ffea438f11a0d327197e068bb4c8b79cd6 100644 --- a/dumux/porousmediumflow/3p3c/implicit/localresidual.hh +++ b/dumux/porousmediumflow/3p3c/implicit/localresidual.hh @@ -76,10 +76,6 @@ class ThreePThreeCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual public: - ThreePThreeCLocalResidual() : ParentType() - { - upwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight); - } /*! * \brief Evaluate the amount of all conservation quantities * (e.g. phase mass) within a sub-control volume. @@ -142,7 +138,6 @@ public: fluxVarsCache); // get upwind weights into local scope - auto w = upwindWeight_; PrimaryVariables flux(0.0); // advective fluxes @@ -150,23 +145,16 @@ public: { for (int compIdx = 0; compIdx < numComponents; ++compIdx) { - auto upwindRule = [w, phaseIdx, compIdx](const VolumeVariables& up, const VolumeVariables& dn) - { - return ( w )*up.molarDensity(phaseIdx) - *up.moleFraction(phaseIdx, compIdx) - *up.mobility(phaseIdx) - +(1-w)*dn.molarDensity(phaseIdx) - *dn.moleFraction(phaseIdx, compIdx) - *dn.mobility(phaseIdx); - }; + auto upwindTerm = [phaseIdx, compIdx](const VolumeVariables& volVars) + { return volVars.molarDensity(phaseIdx)*volVars.moleFraction(phaseIdx, compIdx)*volVars.mobility(phaseIdx); }; // get equation index auto eqIdx = conti0EqIdx + compIdx; - flux[eqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindRule); + flux[eqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm); } //! Add advective phase energy fluxes. For isothermal model the contribution is zero. - EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx, w); + EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx); } //! Add diffusive energy fluxes. For isothermal model the contribution is zero. @@ -203,9 +191,6 @@ protected: const Implementation *asImp_() const { return static_cast (this); } - -private: - Scalar upwindWeight_; }; } // end namespace diff --git a/dumux/porousmediumflow/compositional/localresidual.hh b/dumux/porousmediumflow/compositional/localresidual.hh index 6a674c5d572fd92fab3e2f3a9a3633cf5233b4b9..65d24df2d41ef8ae40f7bdcb6f5b02ef642a8f66 100644 --- a/dumux/porousmediumflow/compositional/localresidual.hh +++ b/dumux/porousmediumflow/compositional/localresidual.hh @@ -70,14 +70,6 @@ class CompositionalLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidua public: - CompositionalLocalResidual() - : ParentType() - { - // retrieve the upwind weight for the mass conservation equations. Use the value - // specified via the property system as default, and overwrite - // it by the run-time parameter from the Dune::ParameterTree - upwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight); - } /*! * \brief Evaluate the amount of all conservation quantities * (e.g. phase mass) within a sub-control volume. @@ -181,7 +173,6 @@ public: fluxVarsCache); // get upwind weights into local scope - auto w = upwindWeight_; PrimaryVariables flux(0.0); // formulation with mole balances @@ -194,17 +185,12 @@ public: { // get equation index auto eqIdx = conti0EqIdx + compIdx; - auto upwindRule = [w, phaseIdx, compIdx](const VolumeVariables& up, const VolumeVariables& dn) - { - return ( w )*up.molarDensity(phaseIdx) - *up.moleFraction(phaseIdx, compIdx) - *up.mobility(phaseIdx) - +(1-w)*dn.molarDensity(phaseIdx) - *dn.moleFraction(phaseIdx, compIdx) - *dn.mobility(phaseIdx); - }; - const auto advFlux = fluxVars.advectiveFlux(phaseIdx, upwindRule); + // the physical quantities for which we perform upwinding + auto upwindTerm = [phaseIdx, compIdx](const VolumeVariables& volVars) + { return volVars.molarDensity(phaseIdx)*volVars.moleFraction(phaseIdx, compIdx)*volVars.mobility(phaseIdx); }; + + const auto advFlux = fluxVars.advectiveFlux(phaseIdx, upwindTerm); if (eqIdx != replaceCompEqIdx) flux[eqIdx] += advFlux; @@ -225,7 +211,7 @@ public: } //! Add advective phase energy fluxes. For isothermal model the contribution is zero. - EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx, w); + EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx); } //! Add diffusive energy fluxes. For isothermal model the contribution is zero. @@ -241,17 +227,12 @@ public: { // get equation index auto eqIdx = conti0EqIdx + compIdx; - auto upwindRule = [w, phaseIdx, compIdx](const VolumeVariables& up, const VolumeVariables& dn) - { - return ( w )*up.density(phaseIdx) - *up.massFraction(phaseIdx, compIdx) - *up.mobility(phaseIdx) - +(1-w)*dn.density(phaseIdx) - *dn.massFraction(phaseIdx, compIdx) - *dn.mobility(phaseIdx); - }; - const auto advFlux = fluxVars.advectiveFlux(phaseIdx, upwindRule); + // the physical quantities for which we perform upwinding + auto upwindTerm = [phaseIdx, compIdx](const VolumeVariables& volVars) + { return volVars.density(phaseIdx)*volVars.massFraction(phaseIdx, compIdx)*volVars.mobility(phaseIdx); }; + + const auto advFlux = fluxVars.advectiveFlux(phaseIdx, upwindTerm); if (eqIdx != replaceCompEqIdx) flux[eqIdx] += advFlux; @@ -288,9 +269,6 @@ protected: const Implementation *asImp_() const { return static_cast (this); } - -private: - Scalar upwindWeight_; }; } // end namespace Dumux diff --git a/dumux/porousmediumflow/immiscible/localresidual.hh b/dumux/porousmediumflow/immiscible/localresidual.hh index 5eebf3a4a653082d4f4c5f9bd3df197dc1aca937..ae195f67d53bd34c70bf8456586fae733f0dc310 100644 --- a/dumux/porousmediumflow/immiscible/localresidual.hh +++ b/dumux/porousmediumflow/immiscible/localresidual.hh @@ -59,18 +59,6 @@ class ImmiscibleLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual) public: - /*! - * \brief Constructor. Gets the upwind weight. - */ - ImmiscibleLocalResidual() - : ParentType() - { - // retrieve the upwind weight for the mass conservation equations. Use the value - // specified via the property system as default, and overwrite - // it by the run-time parameter from the Dune::ParameterTree - upwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight); - } - /*! * \brief Evaluate the rate of change of all conservation * quantites (e.g. phase mass) within a sub-control @@ -122,22 +110,18 @@ public: scvf, fluxVarsCache); - // copy weight to local scope for use in lambda expression - auto w = upwindWeight_; - PrimaryVariables flux; for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - // the upwinding scheme - auto upwindRule = [w, phaseIdx](const VolumeVariables& up, const VolumeVariables& dn) - { return (up.density(phaseIdx)*up.mobility(phaseIdx))*(w) - + (dn.density(phaseIdx)*dn.mobility(phaseIdx))*(1-w); }; + // the physical quantities for which we perform upwinding + auto upwindTerm = [phaseIdx](const VolumeVariables& volVars) + { return volVars.density(phaseIdx)*volVars.mobility(phaseIdx); }; auto eqIdx = conti0EqIdx + phaseIdx; - flux[eqIdx] = fluxVars.advectiveFlux(phaseIdx, upwindRule); + flux[eqIdx] = fluxVars.advectiveFlux(phaseIdx, upwindTerm); //! Add advective phase energy fluxes. For isothermal model the contribution is zero. - EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx, w); + EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx); } //! Add diffusive energy fluxes. For isothermal model the contribution is zero. @@ -152,8 +136,6 @@ private: const Implementation *asImp_() const { return static_cast (this); } - - Scalar upwindWeight_; }; } diff --git a/dumux/porousmediumflow/implicit/fluxvariables.hh b/dumux/porousmediumflow/implicit/fluxvariables.hh index d37d1c24cfd0f51d44ce11fed245718163216058..c0fb4691a9d06707c7aca93a3e70f1b46301b9b3 100644 --- a/dumux/porousmediumflow/implicit/fluxvariables.hh +++ b/dumux/porousmediumflow/implicit/fluxvariables.hh @@ -37,6 +37,7 @@ NEW_PROP_TAG(NumPhases); NEW_PROP_TAG(EnableAdvection); NEW_PROP_TAG(EnableMolecularDiffusion); NEW_PROP_TAG(EnableEnergyBalance); +NEW_PROP_TAG(ImplicitMassUpwindWeight); } // forward declaration @@ -84,10 +85,14 @@ public: const FluxVariablesCache& fluxVarsCache) { ParentType::init(problem, element, fvGeometry, elemVolVars, scvFace, fluxVarsCache); + // retrieve the upwind weight for the mass conservation equations. Use the value + // specified via the property system as default, and overwrite + // it by the run-time parameter from the Dune::ParameterTree + upwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight); } template - Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindFunction) + Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindTerm) { Scalar flux = AdvectionType::flux(this->problem(), this->element(), @@ -122,10 +127,13 @@ public: } } + // upwind scheme if (std::signbit(flux)) - return flux*upwindFunction(outsideVolVars, insideVolVars); + return flux*(upwindWeight_*upwindTerm(outsideVolVars) + + (1.0 - upwindWeight_)*upwindTerm(insideVolVars)); else - return flux*upwindFunction(insideVolVars, outsideVolVars); + return flux*(upwindWeight_*upwindTerm(insideVolVars) + + (1.0 - upwindWeight_)*upwindTerm(outsideVolVars)); } Stencil computeStencil(const Problem& problem, @@ -133,6 +141,9 @@ public: const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvFace) { return AdvectionType::stencil(problem, element, fvGeometry, scvFace); } + +private: + Scalar upwindWeight_; }; @@ -169,10 +180,14 @@ public: { advFluxCached_.reset(); ParentType::init(problem, element, fvGeometry, elemVolVars, scvFace, fluxVarsCache); + // retrieve the upwind weight for the mass conservation equations. Use the value + // specified via the property system as default, and overwrite + // it by the run-time parameter from the Dune::ParameterTree + upwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight); } template - Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindFunction) + Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindTerm) { if (!advFluxCached_[phaseIdx]) { @@ -194,9 +209,11 @@ public: const auto& outsideVolVars = this->elemVolVars()[this->scvFace().outsideScvIdx()]; if (std::signbit(advPreFlux_[phaseIdx])) - return advPreFlux_[phaseIdx]*upwindFunction(outsideVolVars, insideVolVars); + return advPreFlux_[phaseIdx]*(upwindWeight_*upwindTerm(outsideVolVars) + + (1.0 - upwindWeight_)*upwindTerm(insideVolVars)); else - return advPreFlux_[phaseIdx]*upwindFunction(insideVolVars, outsideVolVars); + return advPreFlux_[phaseIdx]*(upwindWeight_*upwindTerm(insideVolVars) + + (1.0 - upwindWeight_)*upwindTerm(outsideVolVars)); } Scalar molecularDiffusionFlux(const int phaseIdx, const int compIdx) @@ -237,6 +254,7 @@ private: //! simple caching if advection flux is used twice with different upwind function std::bitset advFluxCached_; std::array advPreFlux_; + Scalar upwindWeight_; }; // specialization for non-isothermal advective flow (e.g. non-isothermal one-phase darcy equation) @@ -272,10 +290,14 @@ public: { advFluxCached_.reset(); ParentType::init(problem, element, fvGeometry, elemVolVars, scvFace, fluxVarsCache); + // retrieve the upwind weight for the mass conservation equations. Use the value + // specified via the property system as default, and overwrite + // it by the run-time parameter from the Dune::ParameterTree + upwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight); } template - Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindFunction) + Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindTerm) { if (!advFluxCached_[phaseIdx]) { @@ -297,9 +319,11 @@ public: const auto& outsideVolVars = this->elemVolVars()[this->scvFace().outsideScvIdx()]; if (std::signbit(advPreFlux_[phaseIdx])) - return advPreFlux_[phaseIdx]*upwindFunction(outsideVolVars, insideVolVars); + return advPreFlux_[phaseIdx]*(upwindWeight_*upwindTerm(outsideVolVars) + + (1.0 - upwindWeight_)*upwindTerm(insideVolVars)); else - return advPreFlux_[phaseIdx]*upwindFunction(insideVolVars, outsideVolVars); + return advPreFlux_[phaseIdx]*(upwindWeight_*upwindTerm(insideVolVars) + + (1.0 - upwindWeight_)*upwindTerm(outsideVolVars)); } Scalar heatConductionFlux() @@ -340,6 +364,7 @@ private: //! simple caching if advection flux is used twice with different upwind function std::bitset advFluxCached_; std::array advPreFlux_; + Scalar upwindWeight_; }; // specialization for non-isothermal advection and difussion equations (e.g. non-isothermal three-phase three-component flow) @@ -376,10 +401,14 @@ public: { advFluxCached_.reset(); ParentType::init(problem, element, fvGeometry, elemVolVars, scvFace, fluxVarsCache); + // retrieve the upwind weight for the mass conservation equations. Use the value + // specified via the property system as default, and overwrite + // it by the run-time parameter from the Dune::ParameterTree + upwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight); } template - Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindFunction) + Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindTerm) { if (!advFluxCached_[phaseIdx]) { @@ -401,9 +430,11 @@ public: const auto& outsideVolVars = this->elemVolVars()[this->scvFace().outsideScvIdx()]; if (std::signbit(advPreFlux_[phaseIdx])) - return advPreFlux_[phaseIdx]*upwindFunction(outsideVolVars, insideVolVars); + return advPreFlux_[phaseIdx]*(upwindWeight_*upwindTerm(outsideVolVars) + + (1.0 - upwindWeight_)*upwindTerm(insideVolVars)); else - return advPreFlux_[phaseIdx]*upwindFunction(insideVolVars, outsideVolVars); + return advPreFlux_[phaseIdx]*(upwindWeight_*upwindTerm(insideVolVars) + + (1.0 - upwindWeight_)*upwindTerm(outsideVolVars)); } Scalar molecularDiffusionFlux(const int phaseIdx, const int compIdx) @@ -458,6 +489,7 @@ private: //! simple caching if advection flux is used twice with different upwind function std::bitset advFluxCached_; std::array advPreFlux_; + Scalar upwindWeight_; }; } // end namespace diff --git a/dumux/porousmediumflow/nonisothermal/implicit/localresidual.hh b/dumux/porousmediumflow/nonisothermal/implicit/localresidual.hh index 1a1bb7d77c73db4235ee5641a20a9beedb8091c8..30f8e53a2b6a1152590bc3cbbbe0c7f1de759d99 100644 --- a/dumux/porousmediumflow/nonisothermal/implicit/localresidual.hh +++ b/dumux/porousmediumflow/nonisothermal/implicit/localresidual.hh @@ -74,8 +74,7 @@ public: //! The advective phase energy fluxes static void heatConvectionFlux(PrimaryVariables& flux, FluxVariables& fluxVars, - int phaseIdx, - Scalar upwindWeight) + int phaseIdx) {} //! The diffusive energy fluxes @@ -124,16 +123,12 @@ public: //! The advective phase energy fluxes static void heatConvectionFlux(PrimaryVariables& flux, FluxVariables& fluxVars, - int phaseIdx, - Scalar upwindWeight) + int phaseIdx) { - auto upwindRule = [upwindWeight, phaseIdx](const VolumeVariables& up, const VolumeVariables& dn) - { - return (up.density(phaseIdx)*up.mobility(phaseIdx)*up.enthalpy(phaseIdx))*(upwindWeight) - + (dn.density(phaseIdx)*dn.mobility(phaseIdx)*dn.enthalpy(phaseIdx))*(1-upwindWeight); - }; + auto upwindTerm = [phaseIdx](const VolumeVariables& volVars) + { return volVars.density(phaseIdx)*volVars.mobility(phaseIdx)*volVars.enthalpy(phaseIdx); }; - flux[energyEqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindRule); + flux[energyEqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm); } //! The diffusive energy fluxes