diff --git a/dumux/material/constraintsolvers/compositionfromfugacities.hh b/dumux/material/constraintsolvers/compositionfromfugacities.hh
index 6b2644e472b36b5704a29d9c44f5448c2d309152..5aade25c60cc970f4029bfacb34bda4b4f6cfeb9 100644
--- a/dumux/material/constraintsolvers/compositionfromfugacities.hh
+++ b/dumux/material/constraintsolvers/compositionfromfugacities.hh
@@ -51,8 +51,7 @@ namespace Dumux {
template <class Scalar, class FluidSystem>
class CompositionFromFugacities
{
- enum { numComponents = FluidSystem::numComponents };
-
+ static constexpr int numComponents = FluidSystem::numComponents;
using ParameterCache = typename FluidSystem::ParameterCache;
public:
diff --git a/dumux/porousmediumflow/mpnc/initialconditionhelper.hh b/dumux/porousmediumflow/mpnc/initialconditionhelper.hh
new file mode 100644
index 0000000000000000000000000000000000000000..69c00dbfa4e9ed08f96647ca820b8b1da95d07c2
--- /dev/null
+++ b/dumux/porousmediumflow/mpnc/initialconditionhelper.hh
@@ -0,0 +1,98 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup MPNCModel
+ * \brief A helper function to get the correct initial conditions by updating the fluidstate and defining the primary variables needed for equilibrium mpnc models for the MPNC model
+ */
+#ifndef DUMUX_MPNC_INITIALCONDITION_HELPER_HH
+#define DUMUX_MPNC_INITIALCONDITION_HELPER_HH
+
+#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
+#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
+
+namespace Dumux {
+
+namespace MPNCInitialConditions {
+
+struct AllPhasesPresent { int refPhaseIdx; };
+struct NotAllPhasesPresent { int refPhaseIdx; };
+
+} // namespace MPNCInitialConditions
+
+template<class PrimaryVariables, class FluidSystem, class ModelTraits>
+class MPNCInitialConditionHelper
+{
+ using Scalar = typename PrimaryVariables::value_type;
+ using AllPhasesPresent = MPNCInitialConditions::AllPhasesPresent;
+ using NotAllPhasesPresent = MPNCInitialConditions::NotAllPhasesPresent;
+
+public:
+ template<class FluidState>
+ void solve(FluidState& fs, const AllPhasesPresent& allPhases) const
+ {
+ typename FluidSystem::ParameterCache paramCache;
+ MiscibleMultiPhaseComposition<Scalar, FluidSystem>::solve(fs,
+ paramCache,
+ allPhases.refPhaseIdx);
+ }
+
+ template<class FluidState>
+ void solve(FluidState& fs, const NotAllPhasesPresent& notAllPhases) const
+ {
+ typename FluidSystem::ParameterCache paramCache;
+ 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;
+ 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] = fs.fugacity(0, compIdx);
+
+ // first M - 1 saturations
+ for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
+ priVars[s0Idx + phaseIdx] = fs.saturation(phaseIdx);
+
+ static constexpr auto pressureFormulation = ModelTraits::pressureFormulation();
+ if (pressureFormulation == MpNcPressureFormulation::mostWettingFirst)
+ priVars[p0Idx] = fs.pressure(/*phaseIdx=*/0);
+ else if (pressureFormulation == MpNcPressureFormulation::leastWettingFirst)
+ priVars[p0Idx] = fs.pressure(numPhases-1);
+ else
+ DUNE_THROW(Dune::InvalidStateException,"unknown pressure formulation");
+
+ return priVars;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/mpnc/2p2ccomparison/problem.hh b/test/porousmediumflow/mpnc/2p2ccomparison/problem.hh
index 2980df54f95139e6ef0ad2ce770f8a7dfe03d033..a0026d152960baffaf31de66001ca85da4ec8985 100644
--- a/test/porousmediumflow/mpnc/2p2ccomparison/problem.hh
+++ b/test/porousmediumflow/mpnc/2p2ccomparison/problem.hh
@@ -32,7 +32,7 @@
#include <dumux/common/numeqvector.hh>
#include <dumux/porousmediumflow/problem.hh>
-#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
+#include <dumux/porousmediumflow/mpnc/initialconditionhelper.hh>
namespace Dumux {
@@ -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,47 +181,23 @@ private:
// the internal method for the initial condition
PrimaryVariables initial_(const GlobalPosition &globalPos) const
{
- PrimaryVariables values(0.0);
- FluidState fs;
+ const Scalar liquidPhaseSaturation = 0.8;
+ const Scalar gasPhasePressure = 1e5;
- // set the fluid temperatures
+ FluidState fs;
fs.setTemperature(this->spatialParams().temperatureAtPos(globalPos));
+ fs.setSaturation(liquidPhaseIdx, liquidPhaseSaturation);
+ fs.setSaturation(gasPhaseIdx, 1.0 - liquidPhaseSaturation);
+ fs.setPressure(gasPhaseIdx, gasPhasePressure);
- // 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);
+ 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]);
-
- // make the fluid state consistent with local thermodynamic
- // equilibrium
- using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
-
- ParameterCache paramCache;
- MiscibleMultiPhaseComposition::solve(fs, paramCache);
+ fs.setPressure(liquidPhaseIdx, gasPhasePressure + pc[liquidPhaseIdx] - pc[gasPhaseIdx]);
- ///////////
- // assign the primary variables
- ///////////
-
- // all N component fugacities
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- values[fug0Idx + compIdx] = fs.fugacity(gasPhaseIdx, compIdx);
-
- // first M - 1 saturations
- for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
- values[s0Idx + phaseIdx] = fs.saturation(phaseIdx);
-
- // first pressure
- values[p0Idx] = fs.pressure(/*phaseIdx=*/0);
- return values;
+ 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 b32cf5f06f9cba00495cb56431b9d2c5b3941488..159e8b0536524a7b852972ae541e68732e93845e 100644
--- a/test/porousmediumflow/mpnc/obstacle/problem.hh
+++ b/test/porousmediumflow/mpnc/obstacle/problem.hh
@@ -34,7 +34,7 @@
#include <dumux/common/numeqvector.hh>
#include <dumux/porousmediumflow/problem.hh>
-#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
+#include <dumux/porousmediumflow/mpnc/initialconditionhelper.hh>
namespace Dumux {
@@ -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,28 +215,21 @@ private:
// the internal method for the initial condition
PrimaryVariables initial_(const GlobalPosition &globalPos) const
{
- PrimaryVariables values(0.0);
- FluidState fs;
-
- int refPhaseIdx;
- int otherPhaseIdx;
+ int refPhaseIdx, otherPhaseIdx;
+ Scalar refPhasePressure, refPhaseSaturation;
- // set the fluid temperatures
+ FluidState fs;
fs.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.setSaturation(liquidPhaseIdx, refPhaseSaturation);
+ fs.setPressure(liquidPhaseIdx, refPhasePressure);
fs.setMoleFraction(liquidPhaseIdx, N2Idx, 0.0);
fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
}
@@ -245,53 +237,24 @@ private:
// 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.setSaturation(gasPhaseIdx, refPhaseSaturation);
+ fs.setPressure(gasPhaseIdx, refPhasePressure);
fs.setMoleFraction(gasPhaseIdx, N2Idx, 0.99);
fs.setMoleFraction(gasPhaseIdx, H2OIdx, 0.01);
}
- // set the other saturation
- fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));
-
- // calculate the capillary pressure
+ 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(fs, wPhaseIdx);
- fs.setPressure(otherPhaseIdx,
- fs.pressure(refPhaseIdx)
- + (pc[otherPhaseIdx] - pc[refPhaseIdx]));
-
- // make the fluid state consistent with local thermodynamic
- // equilibrium
- using ComputeFromReferencePhase = ComputeFromReferencePhase<Scalar, FluidSystem>;
-
- ParameterCache paramCache;
- ComputeFromReferencePhase::solve(fs,
- paramCache,
- refPhaseIdx);
-
- ///////////
- // assign the primary variables
- ///////////
-
- // all N component fugacities
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- values[fug0Idx + compIdx] = fs.fugacity(refPhaseIdx, compIdx);
-
- // first M - 1 saturations
- for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
- values[s0Idx + phaseIdx] = fs.saturation(phaseIdx);
-
- // first pressure
- values[p0Idx] = fs.pressure(/*phaseIdx=*/0);
- return values;
+ fs.setPressure(otherPhaseIdx, refPhasePressure + (pc[otherPhaseIdx] - pc[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 3e6bfa6167e478eff631e62074993ad33e8629b8..ce771da48821ecb6bffeadf412cd3487465c7499 100644
--- a/test/porousmediumflow/mpnc/thermalnonequilibrium/problem.hh
+++ b/test/porousmediumflow/mpnc/thermalnonequilibrium/problem.hh
@@ -36,8 +36,7 @@
#include <dumux/common/numeqvector.hh>
#include <dumux/porousmediumflow/problem.hh>
-#include <dumux/porousmediumflow/mpnc/pressureformulation.hh>
-#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
+#include <dumux/porousmediumflow/mpnc/initialconditionhelper.hh>
namespace Dumux {
@@ -73,27 +72,22 @@ class CombustionProblemOneComponent: public PorousMediumFlowProblem<TypeTag>
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using Indices = typename ModelTraits::Indices;
- enum {dimWorld = GridView::dimensionworld};
- enum {numComponents = ModelTraits::numFluidComponents()};
- enum {s0Idx = Indices::s0Idx};
- enum {p0Idx = Indices::p0Idx};
- enum {conti00EqIdx = Indices::conti0EqIdx};
- enum {energyEq0Idx = Indices::energyEqIdx};
- enum {numEnergyEqFluid = ModelTraits::numEnergyEqFluid()};
- enum {numEnergyEqSolid = ModelTraits::numEnergyEqSolid()};
- enum {energyEqSolidIdx = energyEq0Idx + numEnergyEqFluid + numEnergyEqSolid - 1};
- enum {wPhaseIdx = FluidSystem::wPhaseIdx};
- enum {nPhaseIdx = FluidSystem::nPhaseIdx};
- enum {wCompIdx = FluidSystem::H2OIdx};
- enum {nCompIdx = FluidSystem::N2Idx};
+ static constexpr int dimWorld = GridView::dimensionworld;
+ static constexpr int numComponents = ModelTraits::numFluidComponents();
+ static constexpr int s0Idx = Indices::s0Idx;
+ static constexpr int p0Idx = Indices::p0Idx;
+ static constexpr int conti00EqIdx = Indices::conti0EqIdx;
+ static constexpr int energyEq0Idx = Indices::energyEqIdx;
+ static constexpr int numEnergyEqFluid = ModelTraits::numEnergyEqFluid();
+ static constexpr int numEnergyEqSolid = ModelTraits::numEnergyEqSolid();
+ static constexpr int energyEqSolidIdx = energyEq0Idx + numEnergyEqFluid + numEnergyEqSolid - 1;
+ static constexpr int wPhaseIdx = FluidSystem::wPhaseIdx;
+ static constexpr int nPhaseIdx = FluidSystem::nPhaseIdx;
+ static constexpr int wCompIdx = FluidSystem::H2OIdx;
+ static constexpr int nCompIdx = FluidSystem::N2Idx;
static constexpr auto numPhases = ModelTraits::numFluidPhases();
- // formulations
- static constexpr auto pressureFormulation = ModelTraits::pressureFormulation();
- static constexpr auto mostWettingFirst = MpNcPressureFormulation::mostWettingFirst;
- static constexpr auto leastWettingFirst = MpNcPressureFormulation::leastWettingFirst;
-
public:
CombustionProblemOneComponent(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
@@ -280,102 +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];
- //////////////////////////////////////
- // Set saturation
- //////////////////////////////////////
- for (int i = 0; i < numPhases - 1; ++i) {
- priVars[s0Idx + i] = S[i];
- }
-
- 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_;
+ FluidState fs;
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ fs.setSaturation(phaseIdx, S[phaseIdx]);
+ fs.setTemperature(thisTemperature);
}
- //////////////////////////////////////
- // Set temperature
- //////////////////////////////////////
- priVars[energyEq0Idx] = thisTemperature;
- priVars[energyEqSolidIdx] = 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(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++)
+ fs.setPressure(phaseIdx, p[phaseIdx]);
- for (int phaseIdx=0; phaseIdx<numPhases; phaseIdx++)
- fluidState.setPressure(phaseIdx, p[phaseIdx]);
+ fs.setMoleFraction(wPhaseIdx, wCompIdx, 1.0);
+ fs.setMoleFraction(wPhaseIdx, nCompIdx, 0.0);
+ fs.setMoleFraction(nPhaseIdx, wCompIdx, 1.0);
+ fs.setMoleFraction(nPhaseIdx, nCompIdx, 0.0);
- //////////////////////////////////////
- // Set pressure
- //////////////////////////////////////
- if(pressureFormulation == mostWettingFirst) {
- // This means that the pressures are sorted from the most wetting to the least wetting-1 in the primary variables vector.
- // For two phases this means that there is one pressure as primary variable: pw
- priVars[p0Idx] = p[wPhaseIdx];
- }
- else if(pressureFormulation == leastWettingFirst) {
- // This means that the pressures are sorted from the least wetting to the most wetting-1 in the primary variables vector.
- // For two phases this means that there is one pressure as primary variable: pn
- priVars[p0Idx] = p[nPhaseIdx];
- }
- else
- DUNE_THROW(Dune::InvalidStateException, "CombustionProblemOneComponent does not support the chosen pressure formulation.");
-
- fluidState.setMoleFraction(wPhaseIdx, wCompIdx, 1.0);
- fluidState.setMoleFraction(wPhaseIdx, nCompIdx, 0.0);
-
- fluidState.setMoleFraction(nPhaseIdx, wCompIdx, 1.0);
- fluidState.setMoleFraction(nPhaseIdx, nCompIdx, 0.0);
+ const int refPhaseIdx = inPM_(globalPos) ? wPhaseIdx : nPhaseIdx;
- int refPhaseIdx;
+ MPNCInitialConditionHelper<PrimaryVariables, FluidSystem, ModelTraits> helper;
+ helper.solve(fs, MPNCInitialConditions::NotAllPhasesPresent{.refPhaseIdx = refPhaseIdx});
+ auto priVars = helper.getPrimaryVariables(fs);
- // on right boundary: reference is gas
- refPhaseIdx = nPhaseIdx;
-
- if(inPM_(globalPos)) {
- refPhaseIdx = wPhaseIdx;
- }
-
- // obtain fugacities
- using ComputeFromReferencePhase = ComputeFromReferencePhase<Scalar, FluidSystem>;
- ParameterCache paramCache;
- ComputeFromReferencePhase::solve(fluidState,
- paramCache,
- refPhaseIdx);
-
- //////////////////////////////////////
- // Set fugacities
- //////////////////////////////////////
- for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
- priVars[conti00EqIdx + compIdx] = fluidState.fugacity(refPhaseIdx,compIdx);
- }
+ // additionally set the temperature for thermal non-equilibrium for the phases
+ priVars[energyEq0Idx] = thisTemperature;
+ priVars[energyEqSolidIdx] = thisTemperature;
return priVars;
}