diff --git a/dumux/porousmediumflow/nonequilibrium/volumevariables.hh b/dumux/porousmediumflow/nonequilibrium/volumevariables.hh
index f4d454e5f8bbfd4af15e695b223b97d0829e1d7b..3cbee9316c2982aad6e44f20a94edb8ca9c4a4ff 100644
--- a/dumux/porousmediumflow/nonequilibrium/volumevariables.hh
+++ b/dumux/porousmediumflow/nonequilibrium/volumevariables.hh
@@ -180,41 +180,40 @@ public:
//obtain parameters for interfacial area constitutive relations
const auto& aWettingNonWettingSurfaceParams = problem.spatialParams().aWettingNonWettingSurfaceParams(element, scv, elemSol);
- const auto& aNonWettingSolidSurfaceParams = problem.spatialParams().aNonWettingSolidSurfaceParams(element, scv, elemSol);
const Scalar pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);
const Scalar Sw = fluidState.saturation(phase0Idx);
- Scalar awn;
-
using AwnSurface = typename Problem::SpatialParams::AwnSurface;
- awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc );
+ const auto awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc);
interfacialArea_[phase0Idx][phase1Idx] = awn;
interfacialArea_[phase1Idx][phase0Idx] = interfacialArea_[phase0Idx][phase1Idx];
interfacialArea_[phase0Idx][phase0Idx] = 0.;
using AnsSurface = typename Problem::SpatialParams::AnsSurface;
- Scalar ans = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams,Sw, pc);
+ const auto& aNonWettingSolidSurfaceParams = problem.spatialParams().aNonWettingSolidSurfaceParams(element, scv, elemSol);
+ const auto ans = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, Sw, pc);
- // Switch for using a a_{wn} relations that has some "maximum capillary pressure" as parameter // That value is obtained by regularization of the pc(Sw) function.
-#if USE_PCMAX
- const Scalar pcMax = problem.spatialParams().pcMax(element, scv, elemSol);
- //I know the solid surface from the pore network. But it is more consistent to use the fitvalue.
- using AnsSurface = typename Problem::SpatialParams::AnsSurface;
- solidSurface_ = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, /*Sw=*/0., pcMax);
+ // Switch for using a a_{wn} relations that has some "maximum capillary pressure" as parameter
+ // That value is obtained by regularization of the pc(Sw) function.
+ static const bool computeAwsFromAnsAndPcMax = getParam<bool>("SpatialParams.ComputeAwsFromAnsAndPcMax", true);
+ if (computeAwsFromAnsAndPcMax)
+ {
+ // I know the solid surface from the pore network. But it is more consistent to use the fit value.
+ const Scalar pcMax = aWettingNonWettingSurfaceParams.pcMax();
+ const auto solidSurface = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, /*Sw=*/0., pcMax);
+ interfacialArea_[phase0Idx][sPhaseIdx] = solidSurface - ans;
+ }
+ else
+ {
+ using AwsSurface = typename Problem::SpatialParams::AwsSurface;
+ const auto& aWettingSolidSurfaceParams = problem.spatialParams().aWettingSolidSurfaceParams(element, scv, elemSol);
+ interfacialArea_[phase0Idx][sPhaseIdx] = AwsSurface::interfacialArea(aWettingSolidSurfaceParams, materialParams, Sw, pc);
+ }
- const Scalar aws = solidSurface_ - ans;
- interfacialArea_[phase0Idx][sPhaseIdx] = aws;
- interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
- interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.;
-#else
- using AwsSurface = typename Problem::SpatialParams::AwsSurface;
- const auto& aWettingSolidSurfaceParams = problem.spatialParams().aWettingSolidSurfaceParams(element, scv, elemSol);
- const auto aws = AwsSurface::interfacialArea(aWettingSolidSurfaceParams,materialParams, Sw, pc);
- interfacialArea_[phase0Idx][sPhaseIdx] = aws ;
interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.;
-#endif
+
interfacialArea_[phase1Idx][sPhaseIdx] = ans;
interfacialArea_[sPhaseIdx][phase1Idx] = interfacialArea_[phase1Idx][sPhaseIdx];
interfacialArea_[phase1Idx][phase1Idx] = 0.;
@@ -259,7 +258,6 @@ private:
Scalar characteristicLength_;
Scalar factorEnergyTransfer_;
- Scalar solidSurface_ ;
InterfacialAreasArray interfacialArea_;
};
@@ -720,46 +718,43 @@ public:
const Scv& scv)
{
// obtain (standard) material parameters (needed for the residual saturations)
- const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
+ const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
//obtain parameters for interfacial area constitutive relations
const auto& aWettingNonWettingSurfaceParams = problem.spatialParams().aWettingNonWettingSurfaceParams(element, scv, elemSol);
- const auto& aNonWettingSolidSurfaceParams = problem.spatialParams().aNonWettingSolidSurfaceParams(element, scv, elemSol);
const Scalar pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);
const Scalar Sw = fluidState.saturation(phase0Idx);
- Scalar awn;
-
using AwnSurface = typename Problem::SpatialParams::AwnSurface;
- awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc );
+ const auto awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc);
interfacialArea_[phase0Idx][phase1Idx] = awn;
interfacialArea_[phase1Idx][phase0Idx] = interfacialArea_[phase0Idx][phase1Idx];
interfacialArea_[phase0Idx][phase0Idx] = 0.;
using AnsSurface = typename Problem::SpatialParams::AnsSurface;
- Scalar ans = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams,Sw, pc);
+ const auto& aNonWettingSolidSurfaceParams = problem.spatialParams().aNonWettingSolidSurfaceParams(element, scv, elemSol);
+ const auto ans = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, Sw, pc);
// Switch for using a a_{wn} relations that has some "maximum capillary pressure" as parameter.
// That value is obtained by regularization of the pc(Sw) function.
-#if USE_PCMAX
- const Scalar pcMax = problem.spatialParams().pcMax(element, scv, elemSol);
- // I know the solid surface from the pore network. But it is more consistent to use the fit value.
- using AnsSurface = typename Problem::SpatialParams::AnsSurface;
- solidSurface_ = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, /*Sw=*/0., pcMax);
+ static const bool computeAwsFromAnsAndPcMax = getParam<bool>("SpatialParams.ComputeAwsFromAnsAndPcMax", true);
+ if (computeAwsFromAnsAndPcMax)
+ {
+ // I know the solid surface from the pore network. But it is more consistent to use the fit value.
+ const Scalar pcMax = aWettingNonWettingSurfaceParams.pcMax();
+ const auto solidSurface = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, /*Sw=*/0., pcMax);
+ interfacialArea_[phase0Idx][sPhaseIdx] = solidSurface - ans;
+ }
+ else
+ {
+ using AwsSurface = typename Problem::SpatialParams::AwsSurface;
+ const auto& aWettingSolidSurfaceParams = problem.spatialParams().aWettingSolidSurfaceParams(element, scv, elemSol);
+ interfacialArea_[phase0Idx][sPhaseIdx] = AwsSurface::interfacialArea(aWettingSolidSurfaceParams, materialParams, Sw, pc);
+ }
- const Scalar aws = solidSurface_ - ans;
- interfacialArea_[phase0Idx][sPhaseIdx] = aws;
interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.;
-#else
- using AwsSurface = typename Problem::SpatialParams::AwsSurface;
- const auto& aWettingSolidSurfaceParams = problem.spatialParams().aWettingSolidSurfaceParams();
- const auto aws = AwsSurface::interfacialArea(aWettingSolidSurfaceParams,materialParams, Sw, pc );
- interfacialArea_[phase0Idx][sPhaseIdx] = aws ;
- interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
- interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.;
-#endif
interfacialArea_[phase1Idx][sPhaseIdx] = ans;
interfacialArea_[sPhaseIdx][phase1Idx] = interfacialArea_[phase1Idx][sPhaseIdx];
@@ -819,7 +814,6 @@ private:
Scalar characteristicLength_;
Scalar factorEnergyTransfer_;
Scalar factorMassTransfer_;
- Scalar solidSurface_ ;
InterfacialAreasArray interfacialArea_;
};
diff --git a/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh b/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh
index e539043254e9d4ed5786c8056dc35ccf595e74cd..89c79251d680b15e1a748d69f8488fc2b41853c5 100644
--- a/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh
+++ b/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh
@@ -35,12 +35,6 @@
#ifndef DUMUX_EVAPORATION_ATMOSPHERE_PROBLEM_HH
#define DUMUX_EVAPORATION_ATMOSPHERE_PROBLEM_HH
-// this sets that the relation using pc_max is used.
-// i.e. - only parameters for awn, ans are given,
-// - the fit for ans involves the maximum value for pc, where Sw, awn are zero.
-// setting it here, because it impacts volume variables and spatialparameters
-#define USE_PCMAX 1
-
#include <dune/grid/yaspgrid.hh>
#include <dumux/common/properties.hh>
diff --git a/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh b/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh
index b546524a802ea95d4884f827109f6613f661d93a..092879eea3079fbb2dfe8bfd9e37a9ea678d5652 100644
--- a/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh
+++ b/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh
@@ -128,9 +128,17 @@ public:
materialParamsFF_.setPe(0.);
// determine maximum capillary pressure for wetting-nonwetting surface
+ /* Of course physically there is no such thing as a maximum capillary pressure.
+ * The parametrization (VG/BC) goes to infinity and physically there is only one pressure
+ * for single phase conditions.
+ * Here, this is used for fitting the interfacial area surface: the capillary pressure,
+ * where the interfacial area is zero.
+ * Technically this value is obtained as the capillary pressure of saturation zero.
+ * This value of course only exists for the case of a regularized pc-Sw relation.
+ */
using TwoPLaw = EffToAbsLaw<RegularizedBrooksCorey<Scalar>>;
- pcMax_ = TwoPLaw::pc(materialParamsPM_, /*sw = */0.0);
- aWettingNonWettingSurfaceParams_.setPcMax(pcMax_);
+ const auto pcMax = TwoPLaw::pc(materialParamsPM_, /*sw = */0.0);
+ aWettingNonWettingSurfaceParams_.setPcMax(pcMax);
// wetting-non wetting: surface which goes to zero on the edges, but is a polynomial
aWettingNonWettingSurfaceParams_.setA1(aWettingNonWettingA1_);
@@ -228,7 +236,7 @@ public:
else if (inPM_(globalPos))
return aWettingNonWettingSurfaceParams_ ;
else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
- }
+ }
/*!\brief Returns a reference to the container object for the
* parametrization of the surface between non-Wetting and solid phase.
@@ -250,26 +258,24 @@ public:
else if (inPM_(globalPos))
return aNonWettingSolidSurfaceParams_ ;
else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
- }
+ }
- /*!\brief Returns the maximum capillary pressure for the given pc-Sw curve
+ /*!\brief Returns a reference to the container object for the
+ * parametrization of the surface between wetting and solid phase.
*
- * Of course physically there is no such thing as a maximum capillary pressure.
- * The parametrization (VG/BC) goes to infinity and physically there is only one pressure
- * for single phase conditions.
- * Here, this is used for fitting the interfacial area surface: the capillary pressure,
- * where the interfacial area is zero.
- * Technically this value is obtained as the capillary pressure of saturation zero.
- * This value of course only exists for the case of a regularized pc-Sw relation.
+ * The position is determined based on the coordinate of
+ * the vertex belonging to the considered sub-control volume.
* \param element The finite element
* \param scv The sub-control volume
* \param elemSol The element solution
*/
template<class ElementSolution>
- const Scalar pcMax(const Element &element,
- const SubControlVolume &scv,
- const ElementSolution &elemSol) const
- { return aWettingNonWettingSurfaceParams_.pcMax() ; }
+ const AwsSurfaceParams& aWettingSolidSurfaceParams(const Element &element,
+ const SubControlVolume &scv,
+ const ElementSolution &elemSol) const
+ {
+ DUNE_THROW(Dune::NotImplemented, "wetting-solid-interface surface params");
+ }
/*!
* \brief Returns the characteristic length for the mass transfer.
@@ -367,8 +373,6 @@ private:
AwnSurfaceParams aWettingNonWettingSurfaceParamsFreeFlow_;
AnsSurfaceParams aNonWettingSolidSurfaceParamsFreeFlow_ ;
- Scalar pcMax_ ;
-
// Porous Medium Domain
Scalar intrinsicPermeabilityPM_ ;
Scalar porosityPM_ ;
@@ -401,6 +405,6 @@ private:
std::vector<Scalar> gridVector_;
};
-}
+} // end namespace Dumux
#endif // GUARDIAN