From 82b83f75dceb82891be1b779cc3b5b7c806120dd Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Wed, 30 Sep 2020 17:49:11 +0200
Subject: [PATCH] [test][mpnckinetic] Adapt spatialparams
---
.../mpnc/implicit/kinetic/problem.hh | 15 +-
.../mpnc/implicit/kinetic/spatialparams.hh | 249 ++++++++----------
2 files changed, 116 insertions(+), 148 deletions(-)
diff --git a/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh b/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh
index 92a34e8a35..8423989dd1 100644
--- a/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh
+++ b/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh
@@ -51,14 +51,6 @@
#include "spatialparams.hh"
-// material laws for interfacial area
-#include <dumux/material/fluidmatrixinteractions/2pia/efftoabslawia.hh>
-#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfacepolynomial2ndorder.hh>
-#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfacepolynomialedgezero2ndorder.hh>
-#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfaceexpfct.hh>
-#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfacepcmaxfct.hh>
-#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfaceexpswpcto3.hh>
-
namespace Dumux {
/*!
* \ingroup MPNCTests
@@ -389,15 +381,12 @@ private:
equilibriumFluidState.setTemperature(phaseIdx, TInitial_ );
}
- const auto &materialParams =
- this->spatialParams().materialLawParamsAtPos(globalPos);
std::vector<Scalar> capPress(numPhases);
// obtain pc according to saturation
- using MaterialLaw = typename ParentType::SpatialParams::MaterialLaw;
- using MPAdapter = MPAdapter<MaterialLaw, numPhases>;
+ using MPAdapter = FluidMatrix::MPAdapter<numPhases>;
const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
- MPAdapter::capillaryPressures(capPress, materialParams, equilibriumFluidState, wPhaseIdx);
+ MPAdapter::capillaryPressures(capPress, this->spatialParams().fluidMatrixInteractionAtPos(globalPos), equilibriumFluidState, wPhaseIdx);
Scalar p[numPhases];
if (this->spatialParams().inPM_(globalPos)){
diff --git a/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh b/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh
index 644fb6d817..15c3c3863c 100644
--- a/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh
+++ b/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh
@@ -28,18 +28,17 @@
#include <dumux/porousmediumflow/properties.hh>
#include <dumux/material/spatialparams/fvnonequilibrium.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/vangenuchtenoftemperature.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+
+#include <dumux/material/fluidmatrixinteractions/fluidmatrixinteraction.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/brookscorey.hh>
+
#include <dumux/common/parameters.hh>
// material laws for interfacial area
-#include <dumux/material/fluidmatrixinteractions/2pia/efftoabslawia.hh>
-#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfacepolynomial2ndorder.hh>
-#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfacepcmaxfct.hh>
-#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfaceexpswpcto3.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/interfacialarea/polynomial2ndorder.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/interfacialarea/pcmax.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/interfacialarea/exponentialcubic.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/interfacialarea/interfacialarea.hh>
namespace Dumux {
@@ -61,37 +60,33 @@ class EvaporationAtmosphereSpatialParams
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
static constexpr auto dimWorld = GridView::dimensionworld;
+
+ using PcKrSwCurve = FluidMatrix::BrooksCoreyDefault<Scalar>;
+
+ using NonwettingSolidInterfacialArea = FluidMatrix::InterfacialArea<Scalar,
+ FluidMatrix::InterfacialAreaExponentialCubic,
+ FluidMatrix::NonwettingSolidInterfacialAreaPcSw>;
+ using WettingSolidInterfacialArea = FluidMatrix::InterfacialArea<Scalar,
+ FluidMatrix::InterfacialAreaPolynomialSecondOrder,
+ FluidMatrix::WettingSolidInterfacialAreaPcSw>;
+ using WettingNonwettingInterfacialArea = FluidMatrix::InterfacialArea<Scalar,
+ FluidMatrix::InterfacialAreaPcMax,
+ FluidMatrix::WettingNonwettingInterfacialAreaPcSw>;
public:
//! Export the type used for the permeability
using PermeabilityType = Scalar;
- //! Export the material law type used
- using MaterialLaw = EffToAbsLaw<RegularizedBrooksCorey<Scalar>>;
- //! Convenience aliases of the law parameters
- using MaterialLawParams = typename MaterialLaw::Params;
-
- //! Export the types used for interfacial area calculations
- using EffectiveIALawAws = AwnSurfacePolynomial2ndOrder<Scalar>;
- using EffectiveIALawAwn = AwnSurfacePcMaxFct<Scalar>;
- using EffectiveIALawAns = AwnSurfaceExpSwPcTo3<Scalar>;
- using AwnSurface = EffToAbsLawIA<EffectiveIALawAwn, MaterialLawParams>;
- using AwsSurface = EffToAbsLawIA<EffectiveIALawAws, MaterialLawParams>;
- using AnsSurface = EffToAbsLawIA<EffectiveIALawAns, MaterialLawParams>;
-
- using AwnSurfaceParams = typename AwnSurface::Params;
- using AwsSurfaceParams = typename AwsSurface::Params;
- using AnsSurfaceParams = typename AnsSurface::Params;
EvaporationAtmosphereSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
- heightPM_ = getParam<std::vector<Scalar>>("Grid.Positions1")[1];
- heightDomain_ = getParam<std::vector<Scalar>>("Grid.Positions1")[2];
+ heightPM_ = getParam<std::vector<Scalar>>("Grid.Positions1")[1];
+ heightDomain_ = getParam<std::vector<Scalar>>("Grid.Positions1")[2];
- porosityPM_ = getParam<Scalar>("SpatialParams.PorousMedium.porosity");
- intrinsicPermeabilityPM_ = getParam<Scalar>("SpatialParams.PorousMedium.permeability");
+ porosityPM_ = getParam<Scalar>("SpatialParams.PorousMedium.porosity");
+ intrinsicPermeabilityPM_ = getParam<Scalar>("SpatialParams.PorousMedium.permeability");
- porosityFF_ = getParam<Scalar>("SpatialParams.FreeFlow.porosity");
- intrinsicPermeabilityFF_ = getParam<Scalar>("SpatialParams.FreeFlow.permeability");
+ porosityFF_ = getParam<Scalar>("SpatialParams.FreeFlow.porosity");
+ intrinsicPermeabilityFF_ = getParam<Scalar>("SpatialParams.FreeFlow.permeability");
aWettingNonwettingA1_ = getParam<Scalar>("SpatialParams.soil.aWettingNonwettingA1");
aWettingNonwettingA2_ = getParam<Scalar>("SpatialParams.soil.aWettingNonwettingA2");
@@ -101,31 +96,26 @@ public:
aNonwettingSolidA2_ = getParam<Scalar>("SpatialParams.soil.aNonwettingSolidA2");
aNonwettingSolidA3_ = getParam<Scalar>("SpatialParams.soil.aNonwettingSolidA3");
- BCPd_ = getParam<Scalar>("SpatialParams.soil.BCPd");
- BClambda_ = getParam<Scalar>("SpatialParams.soil.BClambda");
- Swr_ = getParam<Scalar>("SpatialParams.soil.Swr");
- Snr_ = getParam<Scalar>("SpatialParams.soil.Snr");
+ BCPd_ = getParam<Scalar>("SpatialParams.soil.BCPd");
+ BClambda_ = getParam<Scalar>("SpatialParams.soil.BClambda");
+ Swr_ = getParam<Scalar>("SpatialParams.soil.Swr");
+ Snr_ = getParam<Scalar>("SpatialParams.soil.Snr");
- characteristicLengthFF_ = getParam<Scalar>("SpatialParams.FreeFlow.meanPoreSize");
- characteristicLengthPM_ = getParam<Scalar>("SpatialParams.PorousMedium.meanPoreSize");
+ characteristicLengthFF_ = getParam<Scalar>("SpatialParams.FreeFlow.meanPoreSize");
+ characteristicLengthPM_ = getParam<Scalar>("SpatialParams.PorousMedium.meanPoreSize");
factorEnergyTransfer_ = getParam<Scalar>("SpatialParams.PorousMedium.factorEnergyTransfer");
factorMassTransfer_ = getParam<Scalar>("SpatialParams.PorousMedium.factorMassTransfer");
- // residual saturations
- materialParamsFF_.setSwr(0.0);
- materialParamsFF_.setSnr(0.00);
-
- materialParamsPM_.setSwr(Swr_);
- materialParamsPM_.setSnr(Snr_);
-
- // pc / kr parameters
- materialParamsPM_.setLambda(BClambda_);
- materialParamsPM_.setPe(BCPd_);
+ // PM parameters for Brooks-Corey
+ typename PcKrSwCurve::BasicParams paramsPM(BCPd_, BClambda_);
+ typename PcKrSwCurve::EffToAbsParams effToAbsParamsPM(Swr_, Snr_);
+ pcKrSwCurvePM_ = std::make_unique<PcKrSwCurve>(paramsPM, effToAbsParamsPM);
- // for making pc == 0 in the FF
- materialParamsFF_.setLambda(42.);
- materialParamsFF_.setPe(0.);
+ // FF parameters for Brooks-Corey
+ typename PcKrSwCurve::BasicParams paramsFF(0/*dummy pe*/, 42/*dummy lambda*/);
+ typename PcKrSwCurve::EffToAbsParams effToAbsParamsFF(0.0/*swr*/, 0.0/*snr*/);
+ pcKrSwCurveFF_ = std::make_unique<PcKrSwCurve>(paramsFF, effToAbsParamsFF);
// determine maximum capillary pressure for wetting-nonwetting surface
/* Of course physically there is no such thing as a maximum capillary pressure.
@@ -136,30 +126,33 @@ public:
* 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>>;
- const auto pcMax = TwoPLaw::pc(materialParamsPM_, /*sw = */0.0);
- aWettingNonwettingSurfaceParams_.setPcMax(pcMax);
+ const auto pcMax = pcKrSwCurvePM_->pc(/*sw = */0.0);
+
+ // non-wetting-solid
+ using NonwettingSolidInterfacialAreaParams = typename NonwettingSolidInterfacialArea::BasicParams;
+ NonwettingSolidInterfacialAreaParams ansParams;
+ ansParams.setA1(aNonwettingSolidA1_);
+ ansParams.setA2(aNonwettingSolidA2_);
+ ansParams.setA3(aNonwettingSolidA3_);
+ aNs_ = std::make_unique<NonwettingSolidInterfacialArea>(ansParams, effToAbsParamsPM);
// wetting-non wetting: surface which goes to zero on the edges, but is a polynomial
- aWettingNonwettingSurfaceParams_.setA1(aWettingNonwettingA1_);
- aWettingNonwettingSurfaceParams_.setA2(aWettingNonwettingA2_);
- aWettingNonwettingSurfaceParams_.setA3(aWettingNonwettingA3_);
-
- // nonwetting-solid
- aNonwettingSolidSurfaceParams_.setA1(aNonwettingSolidA1_);
- aNonwettingSolidSurfaceParams_.setA2(aNonwettingSolidA2_);
- aNonwettingSolidSurfaceParams_.setA3(aNonwettingSolidA3_);
-
- // dummys for free flow: no interface where there is only one phase
- aWettingNonwettingSurfaceParamsFreeFlow_.setA1(0.);
- aWettingNonwettingSurfaceParamsFreeFlow_.setA2(0.);
- aWettingNonwettingSurfaceParamsFreeFlow_.setA3(0.);
- aWettingNonwettingSurfaceParamsFreeFlow_.setPcMax(42.); // not needed because it is anyways zero;
-
- // dummys for free flow: no interface where there is only one phase
- aNonwettingSolidSurfaceParamsFreeFlow_.setA1(0.);
- aNonwettingSolidSurfaceParamsFreeFlow_.setA2(0.);
- aNonwettingSolidSurfaceParamsFreeFlow_.setA3(0.);
+ using WettingNonwettingInterfacialAreaParams = typename WettingNonwettingInterfacialArea::BasicParams;
+ WettingNonwettingInterfacialAreaParams anwParams;
+ anwParams.setPcMax(pcMax);
+ anwParams.setA1(aWettingNonwettingA1_);
+ anwParams.setA2(aWettingNonwettingA2_);
+ anwParams.setA3(aWettingNonwettingA3_);
+ aNw_ = std::make_unique<WettingNonwettingInterfacialArea>(anwParams, effToAbsParamsPM);
+
+ // zero-initialized dummys for free flow: no interface where there is only one phase
+ aNwFreeFlow_ = std::make_unique<WettingNonwettingInterfacialArea>(WettingNonwettingInterfacialAreaParams(), effToAbsParamsFF);
+ aNsFreeFlow_ = std::make_unique<NonwettingSolidInterfacialArea>(NonwettingSolidInterfacialAreaParams(), effToAbsParamsFF);
+
+ // dummy
+ using WettingSolidInterfacialAreaParams = typename WettingSolidInterfacialArea::BasicParams;
+ WettingSolidInterfacialAreaParams awsParams;
+ aWs_ = std::make_unique<WettingSolidInterfacialArea>(awsParams, effToAbsParamsPM); // use ctor without effToAbs
}
template<class ElementSolution>
@@ -200,81 +193,63 @@ public:
DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
}
- template<class ElementSolution>
- const MaterialLawParams& materialLawParams(const Element& element,
- const SubControlVolume& scv,
- const ElementSolution& elemSol) const
- { return materialLawParamsAtPos(scv.dofPosition()); }
-
- const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
+ /*!
+ * \brief Returns the parameters for the material law at a given location
+ * \param globalPos A global coordinate vector
+ */
+ const NonwettingSolidInterfacialArea& nonwettingSolidInterfacialAreaAtPos(const GlobalPosition &globalPos) const
{
- if (inFF_(globalPos))
- return materialParamsFF_;
+ if (inFF_(globalPos) )
+ return *aNsFreeFlow_ ;
else if (inPM_(globalPos))
- return materialParamsPM_;
+ return *aNs_ ;
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 wetting and nonwetting phase.
- *
- * 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
+ /*!
+ * \brief Returns the parameters for the material law at a given location
+ * \param globalPos A global coordinate vector
*/
- template<class ElementSolution>
- const AwnSurfaceParams& aWettingNonwettingSurfaceParams(const Element &element,
- const SubControlVolume &scv,
- const ElementSolution &elemSol) const
+ const WettingSolidInterfacialArea& wettingSolidInterfacialAreaAtPos(const GlobalPosition &globalPos) const
{
- const auto& globalPos = scv.dofPosition();
- if (inFF_(globalPos) )
- return aWettingNonwettingSurfaceParamsFreeFlow_ ;
- else if (inPM_(globalPos))
- return aWettingNonwettingSurfaceParams_ ;
- else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
+ DUNE_THROW(Dune::InvalidStateException, "Should not be called in this test.");
}
- /*!\brief Returns a reference to the container object for the
- * parametrization of the surface between nonwetting and solid phase.
- *
- * 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
+ /*!
+ * \brief Returns the parameters for the material law at a given location
+ * \param globalPos A global coordinate vector
*/
- template<class ElementSolution>
- const AnsSurfaceParams& aNonwettingSolidSurfaceParams(const Element &element,
- const SubControlVolume &scv,
- const ElementSolution &elemSol) const
+ const WettingNonwettingInterfacialArea& wettingNonwettingInterfacialAreaAtPos(const GlobalPosition &globalPos) const
{
- const auto& globalPos = scv.dofPosition();
if (inFF_(globalPos) )
- return aNonwettingSolidSurfaceParamsFreeFlow_ ;
+ return *aNwFreeFlow_ ;
else if (inPM_(globalPos))
- return aNonwettingSolidSurfaceParams_ ;
+ return *aNw_ ;
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 wetting and solid phase.
- *
- * 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
+ /*!
+ * \brief Returns the parameters for the material law at a given location
*/
template<class ElementSolution>
- const AwsSurfaceParams& aWettingSolidSurfaceParams(const Element &element,
- const SubControlVolume &scv,
- const ElementSolution &elemSol) const
+ auto fluidMatrixInteraction(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol) const
+ {
+ return fluidMatrixInteractionAtPos(scv.dofPosition());
+ }
+
+ /*!
+ * \brief Returns the parameters for the material law at a given location
+ */
+ auto fluidMatrixInteractionAtPos(const GlobalPosition &globalPos) const
{
- DUNE_THROW(Dune::NotImplemented, "wetting-solid-interface surface params");
+ if (inFF_(globalPos))
+ return makeFluidMatrixInteraction(*pcKrSwCurveFF_, *aNsFreeFlow_, *aNwFreeFlow_, *aWs_);
+ else if (inPM_(globalPos))
+ return makeFluidMatrixInteraction(*pcKrSwCurvePM_, *aNs_, *aNw_, *aWs_);
+ else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
}
/*!
@@ -368,11 +343,6 @@ private:
static constexpr Scalar eps_ = 1e-6;
Scalar heightDomain_ ;
- AwnSurfaceParams aWettingNonwettingSurfaceParams_;
- AnsSurfaceParams aNonwettingSolidSurfaceParams_ ;
- AwnSurfaceParams aWettingNonwettingSurfaceParamsFreeFlow_;
- AnsSurfaceParams aNonwettingSolidSurfaceParamsFreeFlow_ ;
-
// Porous Medium Domain
Scalar intrinsicPermeabilityPM_ ;
Scalar porosityPM_ ;
@@ -380,13 +350,13 @@ private:
Scalar factorEnergyTransfer_ ;
Scalar factorMassTransfer_ ;
Scalar characteristicLengthPM_ ;
- MaterialLawParams materialParamsPM_ ;
+ std::unique_ptr<PcKrSwCurve> pcKrSwCurvePM_;
// Free Flow Domain
Scalar porosityFF_ ;
Scalar intrinsicPermeabilityFF_ ;
Scalar characteristicLengthFF_ ;
- MaterialLawParams materialParamsFF_ ;
+ std::unique_ptr<PcKrSwCurve> pcKrSwCurveFF_;
// interfacial area parameters
Scalar aWettingNonwettingA1_ ;
@@ -403,6 +373,15 @@ private:
Scalar Swr_ ;
Scalar Snr_ ;
std::vector<Scalar> gridVector_;
+
+ std::unique_ptr<NonwettingSolidInterfacialArea> aNs_;
+ std::unique_ptr<WettingNonwettingInterfacialArea> aNw_;
+
+ std::unique_ptr<WettingNonwettingInterfacialArea> aNwFreeFlow_;
+ std::unique_ptr<NonwettingSolidInterfacialArea> aNsFreeFlow_;
+
+ std::unique_ptr<WettingSolidInterfacialArea> aWs_; // dummy, never used
+
};
} // end namespace Dumux
--
GitLab