Commit 2cf6ffb8 authored by Timo Koch's avatar Timo Koch Committed by Kilian Weishaupt
Browse files

[test][mpnc] Use new fluidmatrixinteraction interface

parent 0dc56527
...@@ -6,9 +6,13 @@ TEnd = 1e4 # [s] ...@@ -6,9 +6,13 @@ TEnd = 1e4 # [s]
UpperRight = 60 40 UpperRight = 60 40
Cells = 24 16 Cells = 24 16
[SpatialParams]
Swr = 0.2
BrooksCoreyPcEntry = 1e4
BrooksCoreyLambda = 2.0
[LinearSolver] [LinearSolver]
ResidualReduction = 1e-12 ResidualReduction = 1e-12
[Problem] [Problem]
Name = obstacle Name = obstacle
...@@ -261,20 +261,16 @@ private: ...@@ -261,20 +261,16 @@ private:
// set pressure of the gas phase // set pressure of the gas phase
fs.setPressure(gasPhaseIdx, 1e5); fs.setPressure(gasPhaseIdx, 1e5);
// calulate the capillary pressure // calulate the capillary pressure
const auto& matParams = const auto& fm =
this->spatialParams().materialLawParamsAtPos(globalPos); this->spatialParams().fluidMatrixInteractionAtPos(globalPos);
PhaseVector pc;
using MaterialLaw = typename ParentType::SpatialParams::MaterialLaw;
using MPAdapter = MPAdapter<MaterialLaw, numPhases>;
const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos); const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
MPAdapter::capillaryPressures(pc, matParams, fs, wPhaseIdx); const auto pc = fm.capillaryPressures(fs, wPhaseIdx);
fs.setPressure(liquidPhaseIdx, fs.setPressure(liquidPhaseIdx,
fs.pressure(gasPhaseIdx) + pc[liquidPhaseIdx] - pc[gasPhaseIdx]); fs.pressure(gasPhaseIdx) + pc[liquidPhaseIdx] - pc[gasPhaseIdx]);
// make the fluid state consistent with local thermodynamic // make the fluid state consistent with local thermodynamic
// equilibrium // equilibrium
using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem>; using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
ParameterCache paramCache; ParameterCache paramCache;
MiscibleMultiPhaseComposition::solve(fs, paramCache); MiscibleMultiPhaseComposition::solve(fs, paramCache);
......
...@@ -27,10 +27,8 @@ ...@@ -27,10 +27,8 @@
#include <dumux/porousmediumflow/properties.hh> #include <dumux/porousmediumflow/properties.hh>
#include <dumux/material/spatialparams/fv.hh> #include <dumux/material/spatialparams/fv.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh> #include <dumux/material/fluidmatrixinteractions/2p/brookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh> #include <dumux/material/fluidmatrixinteractions/mp/mpadapter.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
namespace Dumux { namespace Dumux {
...@@ -53,14 +51,16 @@ class MPNCComparisonSpatialParams ...@@ -53,14 +51,16 @@ class MPNCComparisonSpatialParams
MPNCComparisonSpatialParams<GridGeometry, Scalar>>; MPNCComparisonSpatialParams<GridGeometry, Scalar>>;
using GlobalPosition = typename SubControlVolume::GlobalPosition; using GlobalPosition = typename SubControlVolume::GlobalPosition;
using EffectiveLaw = RegularizedBrooksCorey<Scalar>;
using PcKrSwCurve = FluidMatrix::BrooksCoreyDefault<Scalar>;
using MPAdapter = Dumux::FluidMatrix::MPAdapter<PcKrSwCurve, 2>;
public: public:
using PermeabilityType = Scalar; using PermeabilityType = Scalar;
using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
using MaterialLawParams = typename MaterialLaw::Params;
MPNCComparisonSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry) : ParentType(gridGeometry) MPNCComparisonSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
, pcKrSw_("SpatialParams")
{ {
// intrinsic permeabilities // intrinsic permeabilities
coarseK_ = 1e-12; coarseK_ = 1e-12;
...@@ -68,18 +68,6 @@ public: ...@@ -68,18 +68,6 @@ public:
// the porosity // the porosity
porosity_ = 0.3; porosity_ = 0.3;
// residual saturations
fineMaterialParams_.setSwr(0.2);
fineMaterialParams_.setSnr(0.0);
coarseMaterialParams_.setSwr(0.2);
coarseMaterialParams_.setSnr(0.0);
// parameters for the Brooks-Corey law
fineMaterialParams_.setPe(1e4);
coarseMaterialParams_.setPe(1e4);
fineMaterialParams_.setLambda(2.0);
coarseMaterialParams_.setLambda(2.0);
} }
template<class ElementSolution> template<class ElementSolution>
...@@ -95,7 +83,6 @@ public: ...@@ -95,7 +83,6 @@ public:
/*! /*!
* \brief Defines the porosity \f$[-]\f$ of the soil * \brief Defines the porosity \f$[-]\f$ of the soil
*
* \param globalPos The global Position * \param globalPos The global Position
*/ */
Scalar porosityAtPos(const GlobalPosition& globalPos) const Scalar porosityAtPos(const GlobalPosition& globalPos) const
...@@ -105,23 +92,16 @@ public: ...@@ -105,23 +92,16 @@ public:
/*! /*!
* \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.). * \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.).
*
* \param globalPos The global position of the sub-control volume. * \param globalPos The global position of the sub-control volume.
* \return The material parameters object
*/ */
const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const
{ {
if (isFineMaterial_(globalPos)) return makeFluidMatrixInteraction(MPAdapter(pcKrSw_));
return fineMaterialParams_;
else
return coarseMaterialParams_;
} }
/*! /*!
* \brief Function for defining which phase is to be considered as the wetting phase. * \brief Function for defining which phase is to be considered as the wetting phase.
*
* \param globalPos The global position * \param globalPos The global position
* \return The wetting phase index
*/ */
template<class FluidSystem> template<class FluidSystem>
int wettingPhaseAtPos(const GlobalPosition& globalPos) const int wettingPhaseAtPos(const GlobalPosition& globalPos) const
...@@ -144,8 +124,7 @@ private: ...@@ -144,8 +124,7 @@ private:
Scalar coarseK_; Scalar coarseK_;
Scalar fineK_; Scalar fineK_;
Scalar porosity_; Scalar porosity_;
MaterialLawParams fineMaterialParams_; PcKrSwCurve pcKrSw_;
MaterialLawParams coarseMaterialParams_;
static constexpr Scalar eps_ = 1e-6; static constexpr Scalar eps_ = 1e-6;
}; };
......
...@@ -381,12 +381,10 @@ private: ...@@ -381,12 +381,10 @@ private:
equilibriumFluidState.setTemperature(phaseIdx, TInitial_ ); equilibriumFluidState.setTemperature(phaseIdx, TInitial_ );
} }
std::vector<Scalar> capPress(numPhases);
// obtain pc according to saturation // obtain pc according to saturation
using MPAdapter = FluidMatrix::MPAdapter<numPhases>;
const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos); const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
MPAdapter::capillaryPressures(capPress, this->spatialParams().fluidMatrixInteractionAtPos(globalPos), equilibriumFluidState, wPhaseIdx); const auto& fm = this->spatialParams().fluidMatrixInteractionAtPos(globalPos);
const auto capPress = fm.capillaryPressures(equilibriumFluidState, wPhaseIdx);
Scalar p[numPhases]; Scalar p[numPhases];
if (this->spatialParams().inPM_(globalPos)){ if (this->spatialParams().inPM_(globalPos)){
......
...@@ -31,6 +31,7 @@ ...@@ -31,6 +31,7 @@
#include <dumux/material/fluidmatrixinteractions/fluidmatrixinteraction.hh> #include <dumux/material/fluidmatrixinteractions/fluidmatrixinteraction.hh>
#include <dumux/material/fluidmatrixinteractions/2p/brookscorey.hh> #include <dumux/material/fluidmatrixinteractions/2p/brookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/mp/mpadapter.hh>
#include <dumux/common/parameters.hh> #include <dumux/common/parameters.hh>
...@@ -62,6 +63,7 @@ class EvaporationAtmosphereSpatialParams ...@@ -62,6 +63,7 @@ class EvaporationAtmosphereSpatialParams
static constexpr auto dimWorld = GridView::dimensionworld; static constexpr auto dimWorld = GridView::dimensionworld;
using PcKrSwCurve = FluidMatrix::BrooksCoreyDefault<Scalar>; using PcKrSwCurve = FluidMatrix::BrooksCoreyDefault<Scalar>;
using MPAdapter = Dumux::FluidMatrix::MPAdapter<PcKrSwCurve, 2>;
using NonwettingSolidInterfacialArea = FluidMatrix::InterfacialArea<Scalar, using NonwettingSolidInterfacialArea = FluidMatrix::InterfacialArea<Scalar,
FluidMatrix::InterfacialAreaExponentialCubic, FluidMatrix::InterfacialAreaExponentialCubic,
...@@ -246,9 +248,9 @@ public: ...@@ -246,9 +248,9 @@ public:
auto fluidMatrixInteractionAtPos(const GlobalPosition &globalPos) const auto fluidMatrixInteractionAtPos(const GlobalPosition &globalPos) const
{ {
if (inFF_(globalPos)) if (inFF_(globalPos))
return makeFluidMatrixInteraction(*pcKrSwCurveFF_, *aNsFreeFlow_, *aNwFreeFlow_, *aWs_); return makeFluidMatrixInteraction(MPAdapter(*pcKrSwCurveFF_), *aNsFreeFlow_, *aNwFreeFlow_, *aWs_);
else if (inPM_(globalPos)) else if (inPM_(globalPos))
return makeFluidMatrixInteraction(*pcKrSwCurvePM_, *aNs_, *aNw_, *aWs_); return makeFluidMatrixInteraction(MPAdapter(*pcKrSwCurvePM_), *aNs_, *aNw_, *aWs_);
else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]); else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
} }
......
...@@ -331,11 +331,8 @@ private: ...@@ -331,11 +331,8 @@ private:
// calculate the capillary pressure // calculate the capillary pressure
const auto fluidMatrixInteraction = this->spatialParams().fluidMatrixInteractionAtPos(globalPos); const auto fluidMatrixInteraction = this->spatialParams().fluidMatrixInteractionAtPos(globalPos);
PhaseVector pc;
using MPAdapter = FluidMatrix::MPAdapter<numPhases>;
const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos); const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
MPAdapter::capillaryPressures(pc, fluidMatrixInteraction, fs, wPhaseIdx); const auto pc = fluidMatrixInteraction.capillaryPressures(fs, wPhaseIdx);
fs.setPressure(otherPhaseIdx, fs.setPressure(otherPhaseIdx,
fs.pressure(refPhaseIdx) fs.pressure(refPhaseIdx)
+ (pc[otherPhaseIdx] - pc[refPhaseIdx])); + (pc[otherPhaseIdx] - pc[refPhaseIdx]));
......
...@@ -29,6 +29,7 @@ ...@@ -29,6 +29,7 @@
#include <dumux/material/spatialparams/fv.hh> #include <dumux/material/spatialparams/fv.hh>
#include <dumux/material/fluidmatrixinteractions/fluidmatrixinteraction.hh> #include <dumux/material/fluidmatrixinteractions/fluidmatrixinteraction.hh>
#include <dumux/material/fluidmatrixinteractions/2p/smoothedlinearlaw.hh> #include <dumux/material/fluidmatrixinteractions/2p/smoothedlinearlaw.hh>
#include <dumux/material/fluidmatrixinteractions/mp/mpadapter.hh>
namespace Dumux { namespace Dumux {
...@@ -53,6 +54,7 @@ class ObstacleSpatialParams ...@@ -53,6 +54,7 @@ class ObstacleSpatialParams
using GlobalPosition = typename SubControlVolume::GlobalPosition; using GlobalPosition = typename SubControlVolume::GlobalPosition;
using PcKrSwCurve = FluidMatrix::SmoothedLinearLaw<Scalar>; using PcKrSwCurve = FluidMatrix::SmoothedLinearLaw<Scalar>;
using MPAdapter = Dumux::FluidMatrix::MPAdapter<PcKrSwCurve, 2>;
public: public:
//! Export the type used for the permeability //! Export the type used for the permeability
...@@ -91,7 +93,7 @@ public: ...@@ -91,7 +93,7 @@ public:
*/ */
auto fluidMatrixInteractionAtPos(const GlobalPosition &globalPos) const auto fluidMatrixInteractionAtPos(const GlobalPosition &globalPos) const
{ {
return makeFluidMatrixInteraction(pcKrSwCurve_); return makeFluidMatrixInteraction(MPAdapter(pcKrSwCurve_));
} }
/*! /*!
......
...@@ -449,12 +449,11 @@ private: ...@@ -449,12 +449,11 @@ private:
////////////////////////////////////// //////////////////////////////////////
priVars[energyEq0Idx] = thisTemperature; priVars[energyEq0Idx] = thisTemperature;
priVars[energyEqSolidIdx] = thisTemperature; priVars[energyEqSolidIdx] = thisTemperature;
std::array<Scalar, numPhases> capPress;
//obtain pc according to saturation //obtain pc according to saturation
const int wettingPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos); const int wettingPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
using MPAdapter = FluidMatrix::MPAdapter<numPhases>; const auto& fm = this->spatialParams().fluidMatrixInteractionAtPos(globalPos);
MPAdapter::capillaryPressures(capPress, this->spatialParams().fluidMatrixInteractionAtPos(globalPos), fluidState, wettingPhaseIdx); const auto capPress = fm.capillaryPressures(fluidState, wettingPhaseIdx);
Scalar p[numPhases]; Scalar p[numPhases];
......
...@@ -31,6 +31,7 @@ ...@@ -31,6 +31,7 @@
#include <dumux/material/spatialparams/fvnonequilibrium.hh> #include <dumux/material/spatialparams/fvnonequilibrium.hh>
#include <dumux/material/fluidmatrixinteractions/2p/heatpipelaw.hh> #include <dumux/material/fluidmatrixinteractions/2p/heatpipelaw.hh>
#include <dumux/material/fluidmatrixinteractions/mp/mpadapter.hh>
#include <dumux/material/fluidmatrixinteractions/1pia/fluidsolidinterfacialareashiwang.hh> #include <dumux/material/fluidmatrixinteractions/1pia/fluidsolidinterfacialareashiwang.hh>
#include <dumux/porousmediumflow/properties.hh> #include <dumux/porousmediumflow/properties.hh>
#include <dumux/material/spatialparams/fv.hh> #include <dumux/material/spatialparams/fv.hh>
...@@ -56,6 +57,7 @@ class CombustionSpatialParams ...@@ -56,6 +57,7 @@ class CombustionSpatialParams
using GlobalPosition = typename SubControlVolume::GlobalPosition; using GlobalPosition = typename SubControlVolume::GlobalPosition;
using PcKrSwCurve = FluidMatrix::HeatPipeLaw<Scalar>; using PcKrSwCurve = FluidMatrix::HeatPipeLaw<Scalar>;
using MPAdapter = Dumux::FluidMatrix::MPAdapter<PcKrSwCurve, 2>;
public: public:
//! Export the type used for the permeability //! Export the type used for the permeability
...@@ -192,7 +194,7 @@ public: ...@@ -192,7 +194,7 @@ public:
*/ */
auto fluidMatrixInteractionAtPos(const GlobalPosition &globalPos) const auto fluidMatrixInteractionAtPos(const GlobalPosition &globalPos) const
{ {
return makeFluidMatrixInteraction(*pcKrSwCurve_); return makeFluidMatrixInteraction(MPAdapter(*pcKrSwCurve_));
} }
private: private:
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment