From 224e6b4ec1fe3841b4cfc2c3b2b0962002afcf6f Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 1 Dec 2016 20:22:08 +0100 Subject: [PATCH] [upwind] Define the upwing quantities instead of the rule in local residual Now we just tell the flux variables which quantities should be receive upwind treatment. The flux variables will decide how to upwind. This way upwind gets more flexible. This is necessary to treat e.g. bifurcations differently where we don't have only two participating volvars but still need to know the physical quantities to upwind from the local residual. --- .../3p3c/implicit/localresidual.hh | 23 ++------ .../compositional/localresidual.hh | 44 ++++----------- .../immiscible/localresidual.hh | 28 ++-------- .../implicit/fluxvariables.hh | 56 +++++++++++++++---- .../nonisothermal/implicit/localresidual.hh | 15 ++--- 5 files changed, 69 insertions(+), 97 deletions(-) diff --git a/dumux/porousmediumflow/3p3c/implicit/localresidual.hh b/dumux/porousmediumflow/3p3c/implicit/localresidual.hh index 1a766f7457..bd1712ffea 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 6a674c5d57..65d24df2d4 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 5eebf3a4a6..ae195f67d5 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 d37d1c24cf..c0fb4691a9 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 1a1bb7d77c..30f8e53a2b 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 -- GitLab