From cb4e03e1fd2e42222d566d150bffe0ac6060a7c0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= <dennis.glaeser@iws.uni-stuttgart.de>
Date: Wed, 27 Jul 2022 08:46:51 +0200
Subject: [PATCH] [mpnc] let initialhelper have a state
---
.../mpnc/initialconditionhelper.hh | 92 +++++++++++--------
.../mpnc/2p2ccomparison/problem.hh | 40 ++++----
.../porousmediumflow/mpnc/obstacle/problem.hh | 66 +++++--------
.../mpnc/thermalnonequilibrium/problem.hh | 57 ++++--------
4 files changed, 113 insertions(+), 142 deletions(-)
diff --git a/dumux/porousmediumflow/mpnc/initialconditionhelper.hh b/dumux/porousmediumflow/mpnc/initialconditionhelper.hh
index 9911984132..9f12971ac2 100644
--- a/dumux/porousmediumflow/mpnc/initialconditionhelper.hh
+++ b/dumux/porousmediumflow/mpnc/initialconditionhelper.hh
@@ -29,27 +29,48 @@
namespace Dumux {
-namespace MPNCInitialConditions{
+namespace MPNCInitialConditions {
struct AllPhasesPresent { int refPhaseIdx; };
struct NotAllPhasesPresent { int refPhaseIdx; };
} // namespace MPNCInitialConditions
-template <class Scalar, class PrimaryVariables, class FluidSystem, class ModelTraits>
-struct MPNCInitialConditionHelper
+template<class PrimaryVariables,
+ class FluidState,
+ class FluidSystem,
+ class ModelTraits>
+class MPNCInitialConditionHelper
{
+ using Scalar = typename PrimaryVariables::value_type;
+ using AllPhasesPresent = MPNCInitialConditions::AllPhasesPresent;
+ using NotAllPhasesPresent = MPNCInitialConditions::NotAllPhasesPresent;
- template<class FluidState, class ParamCache>
- static auto solveForPrimaryVariables(FluidState& fluidState,
- ParamCache& paramCache,
- const MPNCInitialConditions::AllPhasesPresent& allPhases)
+public:
+ const FluidState& fluidState() const
+ { return fluidState_; }
+
+ void setPressure(int phaseIdx, Scalar value)
+ { fluidState_.setPressure(phaseIdx, value); }
+
+ void setSaturation(int phaseIdx, Scalar value)
+ { fluidState_.setSaturation(phaseIdx, value); }
+
+ void setTemperature(Scalar value)
+ { fluidState_.setTemperature(value); }
+
+ void setMoleFraction(int phaseIdx, int compIdx, Scalar value)
+ { fluidState_.setMoleFraction(phaseIdx, compIdx, value); }
+
+ PrimaryVariables solveForPrimaryVariables(const AllPhasesPresent& allPhases)
{
+ PrimaryVariables priVars(0.0);
// make the fluid state consistent with local thermodynamic equilibrium
- using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
-
- MiscibleMultiPhaseComposition::solve(fluidState, paramCache, allPhases.refPhaseIdx);
+ typename FluidSystem::ParameterCache paramCache;
+ MiscibleMultiPhaseComposition<Scalar, FluidSystem>::solve(fluidState_,
+ paramCache,
+ allPhases.refPhaseIdx);
static constexpr auto numComponents = FluidSystem::numComponents;
static constexpr auto numPhases = FluidSystem::numPhases;
@@ -57,40 +78,34 @@ struct MPNCInitialConditionHelper
static constexpr auto s0Idx = ModelTraits::Indices::s0Idx;
static constexpr auto p0Idx = ModelTraits::Indices::p0Idx;
- PrimaryVariables values(0.0);
// all N component fugacities
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- values[fug0Idx + compIdx] = fluidState.fugacity(0, compIdx);
+ priVars[fug0Idx + compIdx] = fluidState_.fugacity(0, compIdx);
// first M - 1 saturations
for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
- values[s0Idx + phaseIdx] = fluidState.saturation(phaseIdx);
+ priVars[s0Idx + phaseIdx] = fluidState_.saturation(phaseIdx);
static constexpr auto pressureFormulation = ModelTraits::pressureFormulation();
-
- // first pressure
- if(pressureFormulation == MpNcPressureFormulation::mostWettingFirst)
- values[p0Idx] = fluidState.pressure(/*phaseIdx=*/0);
- else if(pressureFormulation == MpNcPressureFormulation::leastWettingFirst)
- values[p0Idx] = fluidState.pressure(numPhases-1);
+ if (pressureFormulation == MpNcPressureFormulation::mostWettingFirst)
+ priVars[p0Idx] = fluidState_.pressure(/*phaseIdx=*/0);
+ else if (pressureFormulation == MpNcPressureFormulation::leastWettingFirst)
+ priVars[p0Idx] = fluidState_.pressure(numPhases-1);
else
DUNE_THROW(Dune::InvalidStateException,"unknown pressure formulation");
- return values;
+ return priVars;
}
- template<class FluidState, class ParamCache>
- static auto solveForPrimaryVariables(FluidState& fluidState,
- ParamCache& paramCache,
- const MPNCInitialConditions::NotAllPhasesPresent& notAllPhases)
+ PrimaryVariables solveForPrimaryVariables(const NotAllPhasesPresent& notAllPhases)
{
+ PrimaryVariables priVars(0.0);
// make the fluid state consistent with local thermodynamic equilibrium
- using ComputeFromReferencePhase = ComputeFromReferencePhase<Scalar, FluidSystem>;
-
- ComputeFromReferencePhase::solve(fluidState,
- paramCache,
- notAllPhases.refPhaseIdx);
+ typename FluidSystem::ParameterCache paramCache;
+ ComputeFromReferencePhase<Scalar, FluidSystem>::solve(fluidState_,
+ paramCache,
+ notAllPhases.refPhaseIdx);
static constexpr auto numComponents = FluidSystem::numComponents;
static constexpr auto numPhases = FluidSystem::numPhases;
@@ -98,26 +113,27 @@ struct MPNCInitialConditionHelper
static constexpr auto s0Idx = ModelTraits::Indices::s0Idx;
static constexpr auto p0Idx = ModelTraits::Indices::p0Idx;
- PrimaryVariables values(0.0);
// all N component fugacities
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- values[fug0Idx + compIdx] = fluidState.fugacity(0, compIdx);
+ priVars[fug0Idx + compIdx] = fluidState_.fugacity(0, compIdx);
// first M - 1 saturations
for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
- values[s0Idx + phaseIdx] = fluidState.saturation(phaseIdx);
+ priVars[s0Idx + phaseIdx] = fluidState_.saturation(phaseIdx);
static constexpr auto pressureFormulation = ModelTraits::pressureFormulation();
- // first pressure
- if(pressureFormulation == MpNcPressureFormulation::mostWettingFirst)
- values[p0Idx] = fluidState.pressure(/*phaseIdx=*/0);
- else if(pressureFormulation == MpNcPressureFormulation::leastWettingFirst)
- values[p0Idx] = fluidState.pressure(numPhases-1);
+ if (pressureFormulation == MpNcPressureFormulation::mostWettingFirst)
+ priVars[p0Idx] = fluidState_.pressure(/*phaseIdx=*/0);
+ else if (pressureFormulation == MpNcPressureFormulation::leastWettingFirst)
+ priVars[p0Idx] = fluidState_.pressure(numPhases-1);
else
DUNE_THROW(Dune::InvalidStateException,"unknown pressure formulation");
- return values;
+ return priVars;
}
+
+private:
+ FluidState fluidState_;
};
} // end namespace Dumux
diff --git a/test/porousmediumflow/mpnc/2p2ccomparison/problem.hh b/test/porousmediumflow/mpnc/2p2ccomparison/problem.hh
index aa458b3426..e7bab36b88 100644
--- a/test/porousmediumflow/mpnc/2p2ccomparison/problem.hh
+++ b/test/porousmediumflow/mpnc/2p2ccomparison/problem.hh
@@ -61,7 +61,7 @@ class MPNCComparisonProblem
using Element = typename GridView::template Codim<0>::Entity;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FluidState = GetPropType<TypeTag, Properties::FluidState>;
- using ParameterCache = typename FluidSystem::ParameterCache;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
static constexpr auto numPhases = GetPropType<TypeTag, Properties::ModelTraits>::numFluidPhases();
static constexpr auto numComponents = GetPropType<TypeTag, Properties::ModelTraits>::numFluidComponents();
@@ -181,31 +181,23 @@ private:
// the internal method for the initial condition
PrimaryVariables initial_(const GlobalPosition &globalPos) const
{
- PrimaryVariables values(0.0);
- FluidState fs;
-
- // set the fluid temperatures
- fs.setTemperature(this->spatialParams().temperatureAtPos(globalPos));
-
- // set water saturation
- fs.setSaturation(liquidPhaseIdx, 0.8);
- fs.setSaturation(gasPhaseIdx, 1.0 - fs.saturation(liquidPhaseIdx));
- // set pressure of the gas phase
- fs.setPressure(gasPhaseIdx, 1e5);
- // calulate the capillary pressure
- const auto& fm =
- this->spatialParams().fluidMatrixInteractionAtPos(globalPos);
+ MPNCInitialConditionHelper<
+ PrimaryVariables, FluidState, FluidSystem, ModelTraits
+ > helper;
+
+ const Scalar liquidPhaseSaturation = 0.8;
+ const Scalar gasPhasePressure = 1e5;
+ helper.setTemperature(this->spatialParams().temperatureAtPos(globalPos));
+ helper.setSaturation(liquidPhaseIdx, liquidPhaseSaturation);
+ helper.setSaturation(gasPhaseIdx, 1.0 - liquidPhaseSaturation);
+ helper.setPressure(gasPhaseIdx, gasPhasePressure);
+
+ const auto& fm = this->spatialParams().fluidMatrixInteractionAtPos(globalPos);
const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
- const auto pc = fm.capillaryPressures(fs, wPhaseIdx);
- fs.setPressure(liquidPhaseIdx,
- fs.pressure(gasPhaseIdx) + pc[liquidPhaseIdx] - pc[gasPhaseIdx]);
+ const auto pc = fm.capillaryPressures(helper.fluidState(), wPhaseIdx);
+ helper.setPressure(liquidPhaseIdx, gasPhasePressure + pc[liquidPhaseIdx] - pc[gasPhaseIdx]);
- ParameterCache paramCache;
-
- using InitialHelper = MPNCInitialConditionHelper<Scalar, PrimaryVariables, FluidSystem, GetPropType<TypeTag, Properties::ModelTraits>>;
- values = InitialHelper::solveForPrimaryVariables(fs, paramCache, MPNCInitialConditions::AllPhasesPresent{.refPhaseIdx = 0});
-
- return values;
+ return helper.solveForPrimaryVariables(MPNCInitialConditions::AllPhasesPresent{.refPhaseIdx = 0});
}
bool onInlet_(const GlobalPosition &globalPos) const
diff --git a/test/porousmediumflow/mpnc/obstacle/problem.hh b/test/porousmediumflow/mpnc/obstacle/problem.hh
index efc685f936..706a99b3ca 100644
--- a/test/porousmediumflow/mpnc/obstacle/problem.hh
+++ b/test/porousmediumflow/mpnc/obstacle/problem.hh
@@ -80,7 +80,6 @@ class ObstacleProblem
using Element = typename GridView::template Codim<0>::Entity;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FluidState = GetPropType<TypeTag, Properties::FluidState>;
- using ParameterCache = typename FluidSystem::ParameterCache;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using Indices = typename ModelTraits::Indices;
@@ -216,66 +215,49 @@ private:
// the internal method for the initial condition
PrimaryVariables initial_(const GlobalPosition &globalPos) const
{
- PrimaryVariables values(0.0);
- FluidState fs;
+ MPNCInitialConditionHelper<
+ PrimaryVariables, FluidState, FluidSystem, ModelTraits
+ > helper;
- int refPhaseIdx;
- int otherPhaseIdx;
-
- // set the fluid temperatures
- fs.setTemperature(this->spatialParams().temperatureAtPos(globalPos));
+ int refPhaseIdx, otherPhaseIdx;
+ Scalar refPhasePressure, refPhaseSaturation;
+ helper.setTemperature(this->spatialParams().temperatureAtPos(globalPos));
if (onInlet_(globalPos))
{
// only liquid on inlet
refPhaseIdx = liquidPhaseIdx;
otherPhaseIdx = gasPhaseIdx;
+ refPhasePressure = 2e5;
+ refPhaseSaturation = 1.0;
- // set liquid saturation
- fs.setSaturation(liquidPhaseIdx, 1.0);
-
- // set pressure of the liquid phase
- fs.setPressure(liquidPhaseIdx, 2e5);
-
- // set the liquid composition to pure water
- fs.setMoleFraction(liquidPhaseIdx, N2Idx, 0.0);
- fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
+ helper.setSaturation(liquidPhaseIdx, refPhaseSaturation);
+ helper.setPressure(liquidPhaseIdx, refPhasePressure);
+ helper.setMoleFraction(liquidPhaseIdx, N2Idx, 0.0);
+ helper.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
}
else {
// elsewhere, only gas
refPhaseIdx = gasPhaseIdx;
otherPhaseIdx = liquidPhaseIdx;
+ refPhasePressure = 1e5;
+ refPhaseSaturation = 1.0;
- // set gas saturation
- fs.setSaturation(gasPhaseIdx, 1.0);
-
- // set pressure of the gas phase
- fs.setPressure(gasPhaseIdx, 1e5);
-
- // set the gas composition to 99% nitrogen and 1% steam
- fs.setMoleFraction(gasPhaseIdx, N2Idx, 0.99);
- fs.setMoleFraction(gasPhaseIdx, H2OIdx, 0.01);
+ helper.setSaturation(gasPhaseIdx, refPhaseSaturation);
+ helper.setPressure(gasPhaseIdx, refPhasePressure);
+ helper.setMoleFraction(gasPhaseIdx, N2Idx, 0.99);
+ helper.setMoleFraction(gasPhaseIdx, H2OIdx, 0.01);
}
- // set the other saturation
- fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));
-
- // calculate the capillary pressure
+ helper.setSaturation(otherPhaseIdx, 1.0 - refPhaseSaturation);
const auto fluidMatrixInteraction = this->spatialParams().fluidMatrixInteractionAtPos(globalPos);
const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
- const auto pc = fluidMatrixInteraction.capillaryPressures(fs, wPhaseIdx);
- fs.setPressure(otherPhaseIdx,
- fs.pressure(refPhaseIdx)
- + (pc[otherPhaseIdx] - pc[refPhaseIdx]));
-
- // make the fluid state consistent with local thermodynamic
- // equilibrium using the initialhelper that selects the correct constraintsolver
-
- ParameterCache paramCache;
+ const auto pc = fluidMatrixInteraction.capillaryPressures(helper.fluidState(), wPhaseIdx);
+ helper.setPressure(otherPhaseIdx, refPhasePressure + (pc[otherPhaseIdx] - pc[refPhaseIdx]));
- using InitialHelper = MPNCInitialConditionHelper<Scalar, PrimaryVariables, FluidSystem, ModelTraits>;
- values = InitialHelper::solveForPrimaryVariables(fs, paramCache, MPNCInitialConditions::NotAllPhasesPresent{.refPhaseIdx = refPhaseIdx});
- return values;
+ return helper.solveForPrimaryVariables(
+ MPNCInitialConditions::NotAllPhasesPresent{.refPhaseIdx = refPhaseIdx}
+ );
}
bool onInlet_(const GlobalPosition &globalPos) const
diff --git a/test/porousmediumflow/mpnc/thermalnonequilibrium/problem.hh b/test/porousmediumflow/mpnc/thermalnonequilibrium/problem.hh
index 28365b4f4d..877b319077 100644
--- a/test/porousmediumflow/mpnc/thermalnonequilibrium/problem.hh
+++ b/test/porousmediumflow/mpnc/thermalnonequilibrium/problem.hh
@@ -274,72 +274,53 @@ private:
// the internal method for the initial condition
PrimaryVariables initial_(const GlobalPosition &globalPos) const
{
- PrimaryVariables priVars(0.0);
const Scalar curPos = globalPos[0];
- const Scalar slope = (SwBoundary_-SwOneComponentSys_) / (this->spatialParams().lengthPM());
+ const Scalar slope = (SwBoundary_ - SwOneComponentSys_)/this->spatialParams().lengthPM();
Scalar S[numPhases];
const Scalar thisSaturation = SwOneComponentSys_ + curPos * slope;
S[wPhaseIdx] = SwBoundary_;
- if (inPM_(globalPos) ) {
+ if (inPM_(globalPos) )
S[wPhaseIdx] = thisSaturation;
- }
-
S[nPhaseIdx] = 1. - S[wPhaseIdx];
- FluidState fluidState;
-
Scalar thisTemperature = TInitial_;
if(onRightBoundary_(globalPos))
- thisTemperature = TRight_;
-
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
- fluidState.setSaturation(phaseIdx, S[phaseIdx]);
-
- fluidState.setTemperature(thisTemperature );
+ thisTemperature = TRight_;
+ MPNCInitialConditionHelper<PrimaryVariables, FluidState, FluidSystem, ModelTraits> helper;
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ helper.setSaturation(phaseIdx, S[phaseIdx]);
+ helper.setTemperature(thisTemperature);
}
//obtain pc according to saturation
const int wettingPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
const auto& fm = this->spatialParams().fluidMatrixInteractionAtPos(globalPos);
- const auto capPress = fm.capillaryPressures(fluidState, wettingPhaseIdx);
-
+ const auto capPress = fm.capillaryPressures(helper.fluidState(), wettingPhaseIdx);
Scalar p[numPhases];
using std::abs;
p[wPhaseIdx] = pnInitial_ - abs(capPress[wPhaseIdx]);
p[nPhaseIdx] = p[wPhaseIdx] + abs(capPress[wPhaseIdx]);
+ for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
+ helper.setPressure(phaseIdx, p[phaseIdx]);
- for (int phaseIdx=0; phaseIdx<numPhases; phaseIdx++)
- fluidState.setPressure(phaseIdx, p[phaseIdx]);
+ helper.setMoleFraction(wPhaseIdx, wCompIdx, 1.0);
+ helper.setMoleFraction(wPhaseIdx, nCompIdx, 0.0);
+ helper.setMoleFraction(nPhaseIdx, wCompIdx, 1.0);
+ helper.setMoleFraction(nPhaseIdx, nCompIdx, 0.0);
- fluidState.setMoleFraction(wPhaseIdx, wCompIdx, 1.0);
- fluidState.setMoleFraction(wPhaseIdx, nCompIdx, 0.0);
-
- fluidState.setMoleFraction(nPhaseIdx, wCompIdx, 1.0);
- fluidState.setMoleFraction(nPhaseIdx, nCompIdx, 0.0);
-
- int refPhaseIdx;
+ const int refPhaseIdx = inPM_(globalPos) ? wPhaseIdx : nPhaseIdx;
- // on right boundary: reference is gas
- refPhaseIdx = nPhaseIdx;
+ auto priVars = helper.solveForPrimaryVariables(
+ MPNCInitialConditions::NotAllPhasesPresent{.refPhaseIdx = refPhaseIdx}
+ );
- if(inPM_(globalPos)) {
- refPhaseIdx = wPhaseIdx;
- }
-
- ParameterCache paramCache;
- // obtain fugacities and the primary variables for chemical equilibrium
- using InitialHelper = MPNCInitialConditionHelper<Scalar, PrimaryVariables, FluidSystem, ModelTraits>;
- priVars = InitialHelper::solveForPrimaryVariables(fluidState, paramCache, MPNCInitialConditions::NotAllPhasesPresent{.refPhaseIdx = refPhaseIdx});
-
- //////////////////////////////////////
// additionally set the temperature for thermal non-equilibrium for the phases
- //////////////////////////////////////
priVars[energyEq0Idx] = thisTemperature;
priVars[energyEqSolidIdx] = thisTemperature;
-
return priVars;
}
--
GitLab