From bb04056f75c7a5e041cbd92289c5114166c44f65 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Thu, 29 Oct 2020 20:13:20 +0100
Subject: [PATCH] [test][3p/infiltration] Adapt spatialParams
---
.../3p/implicit/infiltration/params.input | 11 ++-
.../3p/implicit/infiltration/problem.hh | 23 ++---
.../3p/implicit/infiltration/spatialparams.hh | 85 ++++++++++---------
3 files changed, 64 insertions(+), 55 deletions(-)
diff --git a/test/porousmediumflow/3p/implicit/infiltration/params.input b/test/porousmediumflow/3p/implicit/infiltration/params.input
index 2fc97e822c..b313b40818 100644
--- a/test/porousmediumflow/3p/implicit/infiltration/params.input
+++ b/test/porousmediumflow/3p/implicit/infiltration/params.input
@@ -12,8 +12,15 @@ Name = infiltration3pcc
[SpatialParams]
permeability = 1.e-11 # m^2
porosity = 0.40
-vanGenuchtenAlpha = 0.0005
-vanGenuchtenN = 4.0
+
+ParkerVanGenuchtenAlpha = 0.0005
+ParkerVanGenuchtenN = 4
+ParkerVanGenuchtenBetaNw = 1
+ParkerVanGenuchtenBetaGw = 1
+ParkerVanGenuchtenBetaGn = 1
+Swr = 0.12
+Snr = 0.07
+Sgr = 0.03
[Output]
PlotFluidMatrixInteractions = false
diff --git a/test/porousmediumflow/3p/implicit/infiltration/problem.hh b/test/porousmediumflow/3p/implicit/infiltration/problem.hh
index 79cb2e8974..8dfc7cd19c 100644
--- a/test/porousmediumflow/3p/implicit/infiltration/problem.hh
+++ b/test/porousmediumflow/3p/implicit/infiltration/problem.hh
@@ -291,35 +291,36 @@ private:
void initial_(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
- Scalar y = globalPos[1];
- Scalar x = globalPos[0];
+ const Scalar y = globalPos[1];
+ const Scalar x = globalPos[0];
Scalar sw, swr=0.12, sgr=0.03;
- if(y > (-1.E-3*x+5) - eps_)
+ if (y > (-1e-3*x + 5) - eps_)
{
- Scalar pc = 9.81 * 1000.0 * (y - (-5E-4*x+5));
+ Scalar pc = 9.81 * 1000.0 * (y - (-5e-4*x + 5));
if (pc < 0.0) pc = 0.0;
sw = invertPcgw_(pc,
- this->spatialParams().materialLawParamsAtPos(globalPos));
+ this->spatialParams().fluidMatrixInteractionAtPos(globalPos));
if (sw < swr) sw = swr;
if (sw > 1.-sgr) sw = 1.-sgr;
values[pressureIdx] = 1e5 ;
values[swIdx] = sw;
values[snIdx] = 0.;
- }else {
- values[pressureIdx] = 1e5 + 9.81 * 1000.0 * ((-5E-4*x+5) - y);
+ }
+ else
+ {
+ values[pressureIdx] = 1e5 + 9.81 * 1000.0 * ((-5e-4*x + 5) - y);
values[swIdx] = 1.-sgr;
values[snIdx] = 0.;
}
}
// small solver inverting the pc curve
- template<class MaterialLawParams>
- static Scalar invertPcgw_(Scalar pcIn, const MaterialLawParams &pcParams)
+ template<class FMInteraction>
+ static Scalar invertPcgw_(Scalar pcIn, const FMInteraction& fluidMatrixInteraction)
{
- using MaterialLaw = typename ParentType::SpatialParams::MaterialLaw;
Scalar lower,upper;
int k;
int maxIt = 50;
@@ -329,7 +330,7 @@ private:
for (k=1; k<=25; k++)
{
sw = 0.5*(upper+lower);
- pcgw = MaterialLaw::pcgw(pcParams, sw);
+ pcgw = fluidMatrixInteraction.pcgw(sw, 0.0);
Scalar delta = pcgw-pcIn;
if (delta<0.) delta*=-1.;
if (delta<bisLimit)
diff --git a/test/porousmediumflow/3p/implicit/infiltration/spatialparams.hh b/test/porousmediumflow/3p/implicit/infiltration/spatialparams.hh
index 6fe052f318..42c0ccf0da 100644
--- a/test/porousmediumflow/3p/implicit/infiltration/spatialparams.hh
+++ b/test/porousmediumflow/3p/implicit/infiltration/spatialparams.hh
@@ -28,11 +28,9 @@
#include <dumux/porousmediumflow/properties.hh>
#include <dumux/material/spatialparams/fv.hh>
-#include <dumux/material/fluidmatrixinteractions/3p/regularizedparkervangen3p.hh>
-#include <dumux/material/fluidmatrixinteractions/3p/regularizedparkervangen3pparams.hh>
-#include <dumux/material/fluidmatrixinteractions/3p/efftoabslaw.hh>
+#include <dumux/material/fluidmatrixinteractions/3p/parkervangenuchten.hh>
#include <dumux/io/gnuplotinterface.hh>
-#include <dumux/io/plotmateriallaw3p.hh>
+#include <dumux/common/enumerate.hh>
namespace Dumux {
/*!
@@ -53,18 +51,15 @@ class InfiltrationThreePSpatialParams
using GlobalPosition = typename SubControlVolume::GlobalPosition;
- using EffectiveLaw = RegularizedParkerVanGen3P<Scalar>;
+ using ThreePhasePcKrSw = FluidMatrix::ParkerVanGenuchten3PDefault<Scalar>;
public:
//! Export permeability type
using PermeabilityType = Scalar;
- //! Get the material law from the property system
- using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
- using MaterialLawParams = typename MaterialLaw::Params;
-
InfiltrationThreePSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
+ , pcKrSwCurve_("SpatialParams")
{
// intrinsic permeabilities
fineK_ = getParam<Scalar>("SpatialParams.permeability");
@@ -72,22 +67,6 @@ public:
// porosities
porosity_ = getParam<Scalar>("SpatialParams.porosity");
- vGAlpha_ = getParam<Scalar>("SpatialParams.vanGenuchtenAlpha");
- vGN_ = getParam<Scalar>("SpatialParams.vanGenuchtenN");
- // residual saturations
- materialParams_.setSwr(0.12);
- materialParams_.setSnr(0.07);
- materialParams_.setSgr(0.03);
-
- // parameters for the 3phase van Genuchten law
- materialParams_.setVgAlpha(vGAlpha_);
- materialParams_.setVgn(vGN_);
- materialParams_.setKrRegardsSnr(false);
-
- // parameters for adsorption
- materialParams_.setKdNAPL(0.);
- materialParams_.setRhoBulk(1500.);
-
plotFluidMatrixInteractions_ = getParam<bool>("Output.PlotFluidMatrixInteractions");
}
@@ -99,15 +78,40 @@ public:
{
GnuplotInterface<Scalar> gnuplot(plotFluidMatrixInteractions_);
gnuplot.setOpenPlotWindow(plotFluidMatrixInteractions_);
- PlotMaterialLaw<Scalar, MaterialLaw> plotMaterialLaw(plotFluidMatrixInteractions_);
+
+ const Scalar sg = 0.2; // assume a fixed gas saturation
+ auto swRange = linspace(0.2, 1.0, 1000);
+
+ // assume fixed sn = 0.2 for pcgw curve
+ auto pcgw = swRange;
+ std::transform(swRange.begin(), swRange.end(), pcgw.begin(), [&](auto x){ return this->pcKrSwCurve_.pcgw(x, 0.2); });
gnuplot.resetAll();
- plotMaterialLaw.addpc(gnuplot, materialParams_);
- gnuplot.plot("pc");
+ gnuplot.setXlabel("Sw");
+ gnuplot.setYlabel("pcgw");
+ gnuplot.addDataSetToPlot(swRange, pcgw, "pcgw", "w l");
+ gnuplot.plot("pcgw-sw");
+
+ // plot kr
+ swRange = linspace(0.2, 0.8, 1000);
+ auto krw = swRange;
+ auto krn = swRange;
+ auto krg = swRange;
+ for (const auto& [i, sw] : enumerate(swRange))
+ {
+ const Scalar sn = 1.0 - sg - sw;
+ krw[i] = pcKrSwCurve_.krw(sw, sn);
+ krn[i] = pcKrSwCurve_.krn(sw, sn);
+ krg[i] = pcKrSwCurve_.krg(sw, sn);
+ }
gnuplot.resetAll();
- plotMaterialLaw.addkr(gnuplot, materialParams_);
- gnuplot.plot("kr");
+ gnuplot.setXlabel("Sw");
+ gnuplot.setYlabel("kr");
+ gnuplot.addDataSetToPlot(swRange, krw, "krw", "w l");
+ gnuplot.addDataSetToPlot(swRange, krn, "krn", "w l");
+ gnuplot.addDataSetToPlot(swRange, krg, "krg", "w l");
+ gnuplot.plot("kr-sw");
}
/*!
@@ -139,29 +143,26 @@ public:
}
/*!
- * \brief Returns the parameter object for the Brooks-Corey material law
+ * \brief Returns the parameters for the material law at a given location
*
- * \param globalPos The global position
+ * \param globalPos The global coordinates for the given location
*/
- const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
+ auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const
{
- return materialParams_;
+ return makeFluidMatrixInteraction(pcKrSwCurve_);
}
+
private:
bool isFineMaterial_(const GlobalPosition &globalPos) const
- { return
- 70. - eps_ <= globalPos[0] && globalPos[0] <= 85. + eps_ &&
- 7.0 - eps_ <= globalPos[1] && globalPos[1] <= 7.50 + eps_;
+ {
+ return 70. - eps_ <= globalPos[0] && globalPos[0] <= 85. + eps_ &&
+ 7.0 - eps_ <= globalPos[1] && globalPos[1] <= 7.50 + eps_;
}
Scalar fineK_;
Scalar coarseK_;
-
Scalar porosity_;
- Scalar vGN_;
- Scalar vGAlpha_;
-
- MaterialLawParams materialParams_;
+ const ThreePhasePcKrSw pcKrSwCurve_;
bool plotFluidMatrixInteractions_;
--
GitLab