From 6b5ef7b0fcf4474795bde1f85aac62446c896b98 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 15:16:40 +0200
Subject: [PATCH] [mpnc][initialhelper] split function
---
.../mpnc/initialconditionhelper.hh | 77 ++++---------------
.../mpnc/2p2ccomparison/problem.hh | 22 +++---
.../porousmediumflow/mpnc/obstacle/problem.hh | 35 ++++-----
.../mpnc/thermalnonequilibrium/problem.hh | 26 +++----
4 files changed, 57 insertions(+), 103 deletions(-)
diff --git a/dumux/porousmediumflow/mpnc/initialconditionhelper.hh b/dumux/porousmediumflow/mpnc/initialconditionhelper.hh
index 9f12971ac2..69c00dbfa4 100644
--- a/dumux/porousmediumflow/mpnc/initialconditionhelper.hh
+++ b/dumux/porousmediumflow/mpnc/initialconditionhelper.hh
@@ -36,10 +36,7 @@ struct NotAllPhasesPresent { int refPhaseIdx; };
} // namespace MPNCInitialConditions
-template<class PrimaryVariables,
- class FluidState,
- class FluidSystem,
- class ModelTraits>
+template<class PrimaryVariables, class FluidSystem, class ModelTraits>
class MPNCInitialConditionHelper
{
using Scalar = typename PrimaryVariables::value_type;
@@ -47,65 +44,28 @@ class MPNCInitialConditionHelper
using NotAllPhasesPresent = MPNCInitialConditions::NotAllPhasesPresent;
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)
+ template<class FluidState>
+ void solve(FluidState& fs, const AllPhasesPresent& allPhases) const
{
- PrimaryVariables priVars(0.0);
-
- // make the fluid state consistent with local thermodynamic equilibrium
typename FluidSystem::ParameterCache paramCache;
- MiscibleMultiPhaseComposition<Scalar, FluidSystem>::solve(fluidState_,
+ MiscibleMultiPhaseComposition<Scalar, FluidSystem>::solve(fs,
paramCache,
allPhases.refPhaseIdx);
-
- static constexpr auto numComponents = FluidSystem::numComponents;
- static constexpr auto numPhases = FluidSystem::numPhases;
- static constexpr auto fug0Idx = ModelTraits::Indices::fug0Idx;
- static constexpr auto s0Idx = ModelTraits::Indices::s0Idx;
- static constexpr auto p0Idx = ModelTraits::Indices::p0Idx;
-
- // all N component fugacities
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- priVars[fug0Idx + compIdx] = fluidState_.fugacity(0, compIdx);
-
- // first M - 1 saturations
- for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
- priVars[s0Idx + phaseIdx] = fluidState_.saturation(phaseIdx);
-
- static constexpr auto pressureFormulation = ModelTraits::pressureFormulation();
- 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 priVars;
}
- PrimaryVariables solveForPrimaryVariables(const NotAllPhasesPresent& notAllPhases)
+ template<class FluidState>
+ void solve(FluidState& fs, const NotAllPhasesPresent& notAllPhases) const
{
- PrimaryVariables priVars(0.0);
-
- // make the fluid state consistent with local thermodynamic equilibrium
typename FluidSystem::ParameterCache paramCache;
- ComputeFromReferencePhase<Scalar, FluidSystem>::solve(fluidState_,
+ ComputeFromReferencePhase<Scalar, FluidSystem>::solve(fs,
paramCache,
notAllPhases.refPhaseIdx);
+ }
+
+ template<class FluidState>
+ PrimaryVariables getPrimaryVariables(const FluidState& fs) const
+ {
+ PrimaryVariables priVars(0.0);
static constexpr auto numComponents = FluidSystem::numComponents;
static constexpr auto numPhases = FluidSystem::numPhases;
@@ -115,25 +75,22 @@ public:
// all N component fugacities
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- priVars[fug0Idx + compIdx] = fluidState_.fugacity(0, compIdx);
+ priVars[fug0Idx + compIdx] = fs.fugacity(0, compIdx);
// first M - 1 saturations
for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
- priVars[s0Idx + phaseIdx] = fluidState_.saturation(phaseIdx);
+ priVars[s0Idx + phaseIdx] = fs.saturation(phaseIdx);
static constexpr auto pressureFormulation = ModelTraits::pressureFormulation();
if (pressureFormulation == MpNcPressureFormulation::mostWettingFirst)
- priVars[p0Idx] = fluidState_.pressure(/*phaseIdx=*/0);
+ priVars[p0Idx] = fs.pressure(/*phaseIdx=*/0);
else if (pressureFormulation == MpNcPressureFormulation::leastWettingFirst)
- priVars[p0Idx] = fluidState_.pressure(numPhases-1);
+ priVars[p0Idx] = fs.pressure(numPhases-1);
else
DUNE_THROW(Dune::InvalidStateException,"unknown pressure formulation");
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 e7bab36b88..a0026d1529 100644
--- a/test/porousmediumflow/mpnc/2p2ccomparison/problem.hh
+++ b/test/porousmediumflow/mpnc/2p2ccomparison/problem.hh
@@ -181,23 +181,23 @@ private:
// the internal method for the initial condition
PrimaryVariables initial_(const GlobalPosition &globalPos) const
{
- 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);
+
+ FluidState fs;
+ fs.setTemperature(this->spatialParams().temperatureAtPos(globalPos));
+ fs.setSaturation(liquidPhaseIdx, liquidPhaseSaturation);
+ fs.setSaturation(gasPhaseIdx, 1.0 - liquidPhaseSaturation);
+ fs.setPressure(gasPhaseIdx, gasPhasePressure);
const auto& fm = this->spatialParams().fluidMatrixInteractionAtPos(globalPos);
const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
- const auto pc = fm.capillaryPressures(helper.fluidState(), wPhaseIdx);
- helper.setPressure(liquidPhaseIdx, gasPhasePressure + pc[liquidPhaseIdx] - pc[gasPhaseIdx]);
+ const auto pc = fm.capillaryPressures(fs, wPhaseIdx);
+ fs.setPressure(liquidPhaseIdx, gasPhasePressure + pc[liquidPhaseIdx] - pc[gasPhaseIdx]);
- return helper.solveForPrimaryVariables(MPNCInitialConditions::AllPhasesPresent{.refPhaseIdx = 0});
+ MPNCInitialConditionHelper<PrimaryVariables, FluidSystem, ModelTraits> helper;
+ helper.solve(fs, MPNCInitialConditions::AllPhasesPresent{.refPhaseIdx = 0});
+ return helper.getPrimaryVariables(fs);
}
bool onInlet_(const GlobalPosition &globalPos) const
diff --git a/test/porousmediumflow/mpnc/obstacle/problem.hh b/test/porousmediumflow/mpnc/obstacle/problem.hh
index 706a99b3ca..159e8b0536 100644
--- a/test/porousmediumflow/mpnc/obstacle/problem.hh
+++ b/test/porousmediumflow/mpnc/obstacle/problem.hh
@@ -215,14 +215,11 @@ private:
// the internal method for the initial condition
PrimaryVariables initial_(const GlobalPosition &globalPos) const
{
- MPNCInitialConditionHelper<
- PrimaryVariables, FluidState, FluidSystem, ModelTraits
- > helper;
-
int refPhaseIdx, otherPhaseIdx;
Scalar refPhasePressure, refPhaseSaturation;
- helper.setTemperature(this->spatialParams().temperatureAtPos(globalPos));
+ FluidState fs;
+ fs.setTemperature(this->spatialParams().temperatureAtPos(globalPos));
if (onInlet_(globalPos))
{
// only liquid on inlet
@@ -231,10 +228,10 @@ private:
refPhasePressure = 2e5;
refPhaseSaturation = 1.0;
- helper.setSaturation(liquidPhaseIdx, refPhaseSaturation);
- helper.setPressure(liquidPhaseIdx, refPhasePressure);
- helper.setMoleFraction(liquidPhaseIdx, N2Idx, 0.0);
- helper.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
+ fs.setSaturation(liquidPhaseIdx, refPhaseSaturation);
+ fs.setPressure(liquidPhaseIdx, refPhasePressure);
+ fs.setMoleFraction(liquidPhaseIdx, N2Idx, 0.0);
+ fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
}
else {
// elsewhere, only gas
@@ -243,21 +240,21 @@ private:
refPhasePressure = 1e5;
refPhaseSaturation = 1.0;
- helper.setSaturation(gasPhaseIdx, refPhaseSaturation);
- helper.setPressure(gasPhaseIdx, refPhasePressure);
- helper.setMoleFraction(gasPhaseIdx, N2Idx, 0.99);
- helper.setMoleFraction(gasPhaseIdx, H2OIdx, 0.01);
+ fs.setSaturation(gasPhaseIdx, refPhaseSaturation);
+ fs.setPressure(gasPhaseIdx, refPhasePressure);
+ fs.setMoleFraction(gasPhaseIdx, N2Idx, 0.99);
+ fs.setMoleFraction(gasPhaseIdx, H2OIdx, 0.01);
}
- helper.setSaturation(otherPhaseIdx, 1.0 - refPhaseSaturation);
+ fs.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(helper.fluidState(), wPhaseIdx);
- helper.setPressure(otherPhaseIdx, refPhasePressure + (pc[otherPhaseIdx] - pc[refPhaseIdx]));
+ const auto pc = fluidMatrixInteraction.capillaryPressures(fs, wPhaseIdx);
+ fs.setPressure(otherPhaseIdx, refPhasePressure + (pc[otherPhaseIdx] - pc[refPhaseIdx]));
- return helper.solveForPrimaryVariables(
- MPNCInitialConditions::NotAllPhasesPresent{.refPhaseIdx = refPhaseIdx}
- );
+ MPNCInitialConditionHelper<PrimaryVariables, FluidSystem, ModelTraits> helper;
+ helper.solve(fs, MPNCInitialConditions::NotAllPhasesPresent{.refPhaseIdx = refPhaseIdx});
+ return helper.getPrimaryVariables(fs);
}
bool onInlet_(const GlobalPosition &globalPos) const
diff --git a/test/porousmediumflow/mpnc/thermalnonequilibrium/problem.hh b/test/porousmediumflow/mpnc/thermalnonequilibrium/problem.hh
index 877b319077..ce771da488 100644
--- a/test/porousmediumflow/mpnc/thermalnonequilibrium/problem.hh
+++ b/test/porousmediumflow/mpnc/thermalnonequilibrium/problem.hh
@@ -288,35 +288,35 @@ private:
if(onRightBoundary_(globalPos))
thisTemperature = TRight_;
- MPNCInitialConditionHelper<PrimaryVariables, FluidState, FluidSystem, ModelTraits> helper;
+ FluidState fs;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
- helper.setSaturation(phaseIdx, S[phaseIdx]);
- helper.setTemperature(thisTemperature);
+ fs.setSaturation(phaseIdx, S[phaseIdx]);
+ fs.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(helper.fluidState(), wettingPhaseIdx);
+ const auto capPress = fm.capillaryPressures(fs, 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]);
+ fs.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);
+ fs.setMoleFraction(wPhaseIdx, wCompIdx, 1.0);
+ fs.setMoleFraction(wPhaseIdx, nCompIdx, 0.0);
+ fs.setMoleFraction(nPhaseIdx, wCompIdx, 1.0);
+ fs.setMoleFraction(nPhaseIdx, nCompIdx, 0.0);
- const int refPhaseIdx = inPM_(globalPos) ? wPhaseIdx : nPhaseIdx;
+ const int refPhaseIdx = inPM_(globalPos) ? wPhaseIdx : nPhaseIdx;
- auto priVars = helper.solveForPrimaryVariables(
- MPNCInitialConditions::NotAllPhasesPresent{.refPhaseIdx = refPhaseIdx}
- );
+ MPNCInitialConditionHelper<PrimaryVariables, FluidSystem, ModelTraits> helper;
+ helper.solve(fs, MPNCInitialConditions::NotAllPhasesPresent{.refPhaseIdx = refPhaseIdx});
+ auto priVars = helper.getPrimaryVariables(fs);
// additionally set the temperature for thermal non-equilibrium for the phases
priVars[energyEq0Idx] = thisTemperature;
--
GitLab