From cc9ee64379c63d80c7afdcc9a93fd17f086d368e Mon Sep 17 00:00:00 2001
From: heck <katharina.heck@iws.uni-stuttgart.de>
Date: Wed, 29 Jun 2022 15:06:51 +0200
Subject: [PATCH] [mpnc] add an initial helper that makes the choice of the
constraintsolver depending on your saturation
---
.../mpnc/initialconditionhelper.hh | 72 +++++++++++++++++++
.../mpnc/2p2ccomparison/problem.hh | 8 +--
test/porousmediumflow/mpnc/kinetic/problem.hh | 12 ++--
.../porousmediumflow/mpnc/obstacle/problem.hh | 11 ++-
4 files changed, 87 insertions(+), 16 deletions(-)
create mode 100644 dumux/porousmediumflow/mpnc/initialconditionhelper.hh
diff --git a/dumux/porousmediumflow/mpnc/initialconditionhelper.hh b/dumux/porousmediumflow/mpnc/initialconditionhelper.hh
new file mode 100644
index 0000000000..6328de46e2
--- /dev/null
+++ b/dumux/porousmediumflow/mpnc/initialconditionhelper.hh
@@ -0,0 +1,72 @@
+// -*- 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 NavierStokesModel
+ * \brief Navier Stokes scalar boundary flux helper
+ */
+#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 {
+
+/*!
+ * \ingroup MPNCModel
+ * \brief A helper function to get the correct initial conditions by updating the fluidstate for the MPNC model
+ */
+template <class Scalar, class FluidSystem>
+struct MPNCInitialConditionHelper
+{
+
+ /*!
+ * \brief Return the area-specific outflow fluxes for all scalar balance equations.
+ * This should only be used of flow reversal does never occur.
+ * A (deactivable) warning is emitted otherwise.
+ */
+ template<class FluidState, class ParameterCache>
+ static void solveFluidStateForMPNCInitialCondition(FluidState &fluidState,
+ ParameterCache ¶mCache,
+ int refPhaseIdx)
+ {
+ if (fluidState.saturation(0) < 1.0 && fluidState.saturation(0) > 0)
+ {
+ // make the fluid state consistent with local thermodynamic
+ // equilibrium
+ using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
+
+ ParameterCache paramCache;
+ MiscibleMultiPhaseComposition::solve(fluidState, paramCache);
+ }
+ else
+ {
+ using ComputeFromReferencePhase = ComputeFromReferencePhase<Scalar, FluidSystem>;
+
+ ParameterCache paramCache;
+ ComputeFromReferencePhase::solve(fluidState,
+ paramCache,
+ refPhaseIdx);
+ }
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/mpnc/2p2ccomparison/problem.hh b/test/porousmediumflow/mpnc/2p2ccomparison/problem.hh
index 2980df54f9..0d134b8daa 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 {
@@ -200,12 +200,10 @@ private:
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>;
+ using InitialHelper = MPNCInitialConditionHelper<Scalar, FluidSystem>;
ParameterCache paramCache;
- MiscibleMultiPhaseComposition::solve(fs, paramCache);
+ InitialHelper::solveFluidStateForMPNCInitialCondition(fs, paramCache, 0);
///////////
// assign the primary variables
diff --git a/test/porousmediumflow/mpnc/kinetic/problem.hh b/test/porousmediumflow/mpnc/kinetic/problem.hh
index d055479993..ed0300193a 100644
--- a/test/porousmediumflow/mpnc/kinetic/problem.hh
+++ b/test/porousmediumflow/mpnc/kinetic/problem.hh
@@ -42,9 +42,11 @@
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/porousmediumflow/mpnc/pressureformulation.hh>
+#include <dumux/porousmediumflow/mpnc/initialconditionhelper.hh>
#include <dumux/material/binarycoefficients/h2o_n2.hh>
#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
+
namespace Dumux {
/*!
@@ -95,7 +97,7 @@ class EvaporationAtmosphereProblem: public PorousMediumFlowProblem<TypeTag>
enum { numEnergyEqSolid = ModelTraits::numEnergyEqSolid() };
static constexpr bool enableChemicalNonEquilibrium = getPropValue<TypeTag, Properties::EnableChemicalNonEquilibrium>();
- using ConstraintSolver = MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
+ using InitialHelper = MPNCInitialConditionHelper<Scalar, FluidSystem>;
// formulations
static constexpr auto pressureFormulation = ModelTraits::pressureFormulation();
@@ -205,9 +207,11 @@ public:
fluidState.setTemperature(liquidPhaseIdx, TInitial_ ); // this value is a good one, TInject does not work
// This solves the system of equations defining x=x(p,T)
+ using ConstraintSolver = Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
ConstraintSolver::solve(fluidState,
dummyCache) ;
+
// Now let's make the air phase less than fully saturated with water
fluidState.setMoleFraction(gasPhaseIdx, wCompIdx, fluidState.moleFraction(gasPhaseIdx, wCompIdx)*percentOfEquil_ ) ;
fluidState.setMoleFraction(gasPhaseIdx, nCompIdx, 1.-fluidState.moleFraction(gasPhaseIdx, wCompIdx) ) ;
@@ -352,10 +356,8 @@ private:
for (int phaseIdx=0; phaseIdx<numPhases; phaseIdx++)
equilibriumFluidState.setPressure(phaseIdx, p[phaseIdx]);
- // This solves the system of equations defining x=x(p,T)
- ParameterCache dummyCache;
- ConstraintSolver::solve(equilibriumFluidState,
- dummyCache) ;
+ ParameterCache paramCache;
+ InitialHelper::solveFluidStateForMPNCInitialCondition(equilibriumFluidState, paramCache, 0);
FluidState dryFluidState(equilibriumFluidState);
// Now let's make the air phase less than fully saturated with vapor
diff --git a/test/porousmediumflow/mpnc/obstacle/problem.hh b/test/porousmediumflow/mpnc/obstacle/problem.hh
index b32cf5f06f..5059fa02f7 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 {
@@ -269,13 +269,12 @@ private:
+ (pc[otherPhaseIdx] - pc[refPhaseIdx]));
// make the fluid state consistent with local thermodynamic
- // equilibrium
- using ComputeFromReferencePhase = ComputeFromReferencePhase<Scalar, FluidSystem>;
+ // equilibrium using the initialhelper that selects the correct constraintsolver
+ using InitialHelper = MPNCInitialConditionHelper<Scalar, FluidSystem>;
ParameterCache paramCache;
- ComputeFromReferencePhase::solve(fs,
- paramCache,
- refPhaseIdx);
+ InitialHelper::solveFluidStateForMPNCInitialCondition(fs, paramCache, refPhaseIdx);
+
///////////
// assign the primary variables
--
GitLab