Skip to content

GitLab

Commit d92cd40d authored by Martin Schneider's avatar Martin Schneider
Browse files

[test][impes] Update impes test problem

parent 4b774a52
Hide whitespace changes
Inline Side-by-side
test/porousmediumflow/2p/sequential/impes/pressureproblem.hh
View file @ d92cd40d
......@@ -45,10 +45,10 @@ namespace Dumux {
template<class TypeTag> class TwoPTestProblem;
namespace Properties {
NEW_TYPE_TAG(TwoPImpes, INHERITS_FROM(CCMpfaModel,Pressure));
NEW_TYPE_TAG(TwoPImpes, INHERITS_FROM(CCTpfaModel,Pressure));
// Set the grid type
SET_TYPE_PROP(TwoPImpes, Grid, Dune::UGGrid<2>);
SET_TYPE_PROP(TwoPImpes, Grid, Dune::YaspGrid<2>);
// Set the problem type
SET_TYPE_PROP(TwoPImpes, Problem, TwoPTestProblem<TypeTag>);
......@@ -58,7 +58,11 @@ SET_PROP(TwoPImpes, FluidSystem)
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
#if PROBLEM == 2
using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Trichloroethene<Scalar> >;
#else
using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
#endif
using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>;
};
......@@ -74,8 +78,8 @@ public:
// Enable caching
SET_BOOL_PROP(TwoPImpes, EnableGridVolumeVariablesCache, false);
SET_BOOL_PROP(TwoPImpes, EnableGridFluxVariablesCache, true);
SET_BOOL_PROP(TwoPImpes, EnableFVGridGeometryCache, true);
SET_BOOL_PROP(TwoPImpes, EnableGridFluxVariablesCache, false);
SET_BOOL_PROP(TwoPImpes, EnableFVGridGeometryCache, false);
} // end namespace Properties
/*!
......@@ -101,9 +105,22 @@ class TwoPTestProblem : public PorousMediumFlowProblem<TypeTag>
pressureIdx = Indices::pressureIdx
};
static constexpr int dimWorld = GridView::dimensionworld;
public:
TwoPTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry) {}
: ParentType(fvGridGeometry)
{
Scalar inletWidth = getParam<Scalar>("Problem.InletWidth", 1.0);
GlobalPosition inletCenter = this->fvGridGeometry().bBoxMax();
inletCenter[0] *= 0.5;
inletLeftCoord_ = inletCenter;
inletLeftCoord_[0] -=0.5*inletWidth;
inletRightCoord_ = inletCenter;
inletRightCoord_[0] +=0.5*inletWidth;
inFlux_ = getParam<Scalar>("Problem.InjectionFlux", 1e-4);
}
/*!
* \brief Specifies which kind of boundary condition should be
......@@ -114,12 +131,25 @@ public:
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
#if PROBLEM == 2
BoundaryTypes values;
if (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos))
{
values.setAllDirichlet();
}
else
{
values.setAllNeumann();
}
return values;
#else
BoundaryTypes values;
if (onLeftBoundary_(globalPos))
values.setAllDirichlet();
else
values.setAllNeumann();
return values;
#endif
}
/*!
......@@ -133,8 +163,18 @@ public:
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(1.0e5);
if (onLeftBoundary_(globalPos))
values[pressureIdx] = 1.0e5;
// if (onLeftBoundary_(globalPos))
// values[pressureIdx] = 1.0e5;
Scalar pRef = 1.0e5;
using WettingPhase = typename GET_PROP(TypeTag, FluidSystem)::WettingPhase;
if (onUpperBoundary_(globalPos))
{
pRef = 2.0e5;
}
values[pressureIdx] = (pRef - (this->fvGridGeometry().bBoxMax()- globalPos) * this->gravity() * WettingPhase::density(temperature(), pRef));
return values;
}
......@@ -153,8 +193,8 @@ public:
NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
{
NumEqVector values(0.0);
if(onRightBoundary_(globalPos))
values[pressureEqIdx] = 3e-3;
//if(isInlet(globalPos))
// values[pressureEqIdx] = -inFlux_;
return values;
}
......@@ -207,7 +247,25 @@ private:
return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_;
}
bool isInlet(const GlobalPosition& globalPos) const
{
if (!onUpperBoundary_(globalPos))
return false;
for (int i = 0; i < dimWorld; i++)
{
if (globalPos[i] < inletLeftCoord_[i] - eps_)
return false;
if (globalPos[i] > inletRightCoord_[i] + eps_)
return false;
}
return true;
}
static constexpr Scalar eps_ = 1e-6;
Scalar inFlux_;
GlobalPosition inletLeftCoord_;
GlobalPosition inletRightCoord_;
};
} // end namespace Dumux
......
test/porousmediumflow/2p/sequential/impes/pressurespatialparams.hh
View file @ d92cd40d
......@@ -60,24 +60,30 @@ public:
TwoPTestSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft");
lensUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensUpperRight");
materialLawParamsBackground_.setSwr(0.);
materialLawParamsBackground_.setSnr(0.);
// residual saturations
lensMaterialParams_.setSwr(0.0);
lensMaterialParams_.setSnr(0.0);
outerMaterialParams_.setSwr(0.0);
outerMaterialParams_.setSnr(0.0);
materialLawParamsLenses_.setSwr(0.);
materialLawParamsLenses_.setSnr(0.);
// parameters for the Van Genuchten law
// alpha and n
lensMaterialParams_.setLambda(2.0);
lensMaterialParams_.setPe(0);
outerMaterialParams_.setLambda(2.0);
outerMaterialParams_.setPe(0);
//parameters for Brooks-Corey law
lensK_ = 1.0e-7;
outerK_ = 1.0e-7;
// entry pressures
materialLawParamsBackground_.setPe(getParam<Scalar>("SpatialParams.BackgroundEntryPressure", 0.0));
materialLawParamsLenses_.setPe(getParam<Scalar>("SpatialParams.LenseEntryPressure", 0.0));
materialLawParamsBackground_.setLambda(getParam<Scalar>("SpatialParams.BackgroundLambda", 3.0));
materialLawParamsLenses_.setLambda(getParam<Scalar>("SpatialParams.LenseLambda", 2.0));
permBackground_ = getParam<Scalar>("SpatialParams.BackgroundPermeability", 1e-10);
permLenses_ = getParam<Scalar>("SpatialParams.LensPermeability", 1e-10);
lensOneLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensOneLowerLeft", GlobalPosition(0.0));
lensOneUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensOneUpperRight", GlobalPosition(0.0));
lensTwoLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensTwoLowerLeft", GlobalPosition(0.0));
lensTwoUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensTwoUpperRight", GlobalPosition(0.0));
lensThreeLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensThreeLowerLeft", GlobalPosition(0.0));
lensThreeUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensThreeUpperRight", GlobalPosition(0.0));
wettingSaturation_.resize(fvGridGeometry->numDofs(),0.0);
}
......@@ -98,9 +104,15 @@ public:
{
// do not use a less permeable lens in the test with inverted wettability
if (isInLens_(element.geometry().center()))
return lensK_;
return outerK_;
auto globalPos = element.geometry().center();
if (isLensOne(globalPos))
return permLenses_;
else if (isLensTwo(globalPos))
return permLenses_;
else if (isLensThree(globalPos))
return permLenses_;
else
return permBackground_;
}
/*!
......@@ -109,7 +121,7 @@ public:
* \param globalPos The global position
*/
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{ return 0.2; }
{ return 0.4; }
/*!
* \brief Returns the parameter object for the Brooks-Corey material law.
......@@ -125,10 +137,15 @@ public:
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
// do not use different parameters in the test with inverted wettability
if (isInLens_(element.geometry().center()))
return lensMaterialParams_;
return outerMaterialParams_;
auto globalPos = element.geometry().center();
if (isLensOne(globalPos))
return materialLawParamsLenses_;
else if (isLensTwo(globalPos))
return materialLawParamsLenses_;
else if (isLensThree(globalPos))
return materialLawParamsLenses_;
else
return materialLawParamsBackground_;
}
/*!
......@@ -161,22 +178,50 @@ public:
{ wettingSaturation_ = sat;}
private:
bool isInLens_(const GlobalPosition &globalPos) const
bool isLensOne(const GlobalPosition& globalPos) const
{
for (int i = 0; i < dimWorld; ++i) {
if (globalPos[i] < lensLowerLeft_[i] + eps_ || globalPos[i] > lensUpperRight_[i] - eps_)
for (int i = 0; i < dimWorld; i++)
{
if (globalPos[i] < lensOneLowerLeft_[i] - eps_ || globalPos[i] > lensOneUpperRight_[i] + eps_)
{
return false;
}
}
return true;
}
bool isLensTwo(const GlobalPosition& globalPos) const
{
for (int i = 0; i < dimWorld; i++)
{
if (globalPos[i] < lensTwoLowerLeft_[i] - eps_ || globalPos[i] > lensTwoUpperRight_[i] + eps_)
{
return false;
}
}
return true;
}
bool isLensThree(const GlobalPosition& globalPos) const
{
for (int i = 0; i < dimWorld; i++)
{
if (globalPos[i] < lensThreeLowerLeft_[i] - eps_ || globalPos[i] > lensThreeUpperRight_[i] + eps_)
{
return false;
}
}
return true;
}
GlobalPosition lensLowerLeft_;
GlobalPosition lensUpperRight_;
Scalar lensK_;
Scalar outerK_;
MaterialLawParams lensMaterialParams_;
MaterialLawParams outerMaterialParams_;
MaterialLawParams materialLawParamsBackground_;
MaterialLawParams materialLawParamsLenses_;
PermeabilityType permBackground_;
PermeabilityType permLenses_;
GlobalPosition lensOneLowerLeft_;
GlobalPosition lensOneUpperRight_;
GlobalPosition lensTwoLowerLeft_;
GlobalPosition lensTwoUpperRight_;
GlobalPosition lensThreeLowerLeft_;
GlobalPosition lensThreeUpperRight_;
std::vector<Scalar> wettingSaturation_;
......
test/porousmediumflow/2p/sequential/impes/test_impes.input
View file @ d92cd40d
[TimeLoop]
DtInitial = 1e2# [s]
MaxTimeStepSize = 1e2
TEnd = 7e5 # [s]
MaxTimeStepSize = 1
TEnd = 5e4 # [s]
[Grid]
File = ../grids/test_mpfa2p.dgf #grid for buckley-leverett or mcwhorter problem
Refinement = 2
#LowerLeft = 0 0
#UpperRight = 30 3
#Cells = 100 3
#File = ../grids/test_mpfa2p.dgf #grid for buckley-leverett or mcwhorter problem
#Refinement = 0
LowerLeft = 0 0
UpperRight = 20 10
Cells = 100 50
[SpatialParams]
LensLowerLeft = 20.0 0.0 # [m] coordinates of the lower left lens corner
LensUpperRight = 40.0 1.0 # [m] coordinates of the upper right lens corner
BackgroundPermeability = 1e-10
LensPermeability= 1e-14
LensOneLowerLeft = 7 6
LensOneUpperRight = 13 7
LensTwoLowerLeft = 2 4
LensTwoUpperRight = 8 5
LensThreeLowerLeft = 10 2
LensThreeUpperRight = 18 3
BackgroundEntryPressure = 0
LenseEntryPressure = 0
[Problem]
Name = test_buckleyleverett
EnableGravity = false
Name = test_lenses
EnableGravity = true
InletWidth = 2
InjectionFlux = 0.1
[Component]
LiquidDensity = 1000
......
test/porousmediumflow/2p/sequential/impes/test_impes_fv.cc
View file @ d92cd40d
......@@ -21,6 +21,8 @@
*
* \brief test for the two-phase porousmedium flow model
*/
#define PROBLEM 2
#include <config.h>
#include <ctime>
......@@ -195,11 +197,11 @@ int main(int argc, char** argv) try
twoPTransportAssembler->setLinearSystem(As, rs);
// the linear solver
using LinearSolver = AMGBackend<TwoPImpesTT>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, twoPImpesFvGridGeometry->dofMapper());
using LinearSolver = UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>();
//! set some check points for the time loop
timeLoop->setPeriodicCheckPoint(1.0e4);
timeLoop->setPeriodicCheckPoint(4);
//! start the time loop
timeLoop->start();
......@@ -223,8 +225,12 @@ int main(int argc, char** argv) try
{
if (scvf.boundary())
{
const auto idx = scvf.outsideScvIdx();
satVals[idx] = elemVolVars[idx].saturation(0);
const auto bcTypes = twoPImpesProblem->boundaryTypes(element, scvf);
if (bcTypes.hasOnlyDirichlet())
{
const auto idx = scvf.outsideScvIdx();
satVals[idx] = 1.0;
}
}
}
}
......@@ -281,6 +287,10 @@ int main(int argc, char** argv) try
fluxVars.init(*twoPImpesProblem, element, fvGeometry, elemVolVars, scvf, elemFluxVars);
volumeFlux[idx] = fluxVars.advectiveFlux(0, upwindTerm);
}
else
{
//volumeFlux[idx] = -0.1;
}
}
}
}
......
test/porousmediumflow/2p/sequential/impes/transportproblem.hh
View file @ d92cd40d
......@@ -47,7 +47,7 @@ template<class TypeTag> class TwoPTransport;
namespace Properties {
NEW_TYPE_TAG(TwoPTransport, INHERITS_FROM(CCTpfaModel, Transport));
// Set the grid type
SET_TYPE_PROP(TwoPTransport, Grid, Dune::UGGrid<2>);
SET_TYPE_PROP(TwoPTransport, Grid, Dune::YaspGrid<2>);
// Set the problem type
SET_TYPE_PROP(TwoPTransport, Problem, TwoPTransport<TypeTag>);
......@@ -57,7 +57,11 @@ SET_PROP(TwoPTransport, FluidSystem)
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
#if PROBLEM == 2
using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Trichloroethene<Scalar> >;
#else
using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
#endif
using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>;
};
......@@ -73,8 +77,8 @@ public:
// Enable caching
SET_BOOL_PROP(TwoPTransport, EnableGridVolumeVariablesCache, false);
SET_BOOL_PROP(TwoPTransport, EnableGridFluxVariablesCache, true);
SET_BOOL_PROP(TwoPTransport, EnableFVGridGeometryCache, true);
SET_BOOL_PROP(TwoPTransport, EnableGridFluxVariablesCache, false);
SET_BOOL_PROP(TwoPTransport, EnableFVGridGeometryCache, false);
} // end namespace Properties
/*!
......@@ -99,10 +103,23 @@ class TwoPTransport : public PorousMediumFlowProblem<TypeTag>
transportEqIdx = Indices::transportEqIdx,
saturationIdx = Indices::saturationIdx
};
static constexpr int dimWorld = GridView::dimensionworld;
public:
TwoPTransport(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry) {}
: ParentType(fvGridGeometry)
{
Scalar inletWidth = getParam<Scalar>("Problem.InletWidth", 1.0);
GlobalPosition inletCenter = this->fvGridGeometry().bBoxMax();
inletCenter[0] *= 0.5;
inletLeftCoord_ = inletCenter;
inletLeftCoord_[0] -=0.5*inletWidth;
inletRightCoord_ = inletCenter;
inletRightCoord_[0] +=0.5*inletWidth;
inFlux_ = getParam<Scalar>("Problem.InjectionFlux", 1e-4);
}
/*!
* \brief Specifies which kind of boundary condition should be
......@@ -113,12 +130,25 @@ public:
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
#if PROBLEM == 2
BoundaryTypes values;
if (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos))
{
values.setAllDirichlet();
}
else
{
values.setAllNeumann();
}
return values;
#else
BoundaryTypes values;
if (onLeftBoundary_(globalPos))
values.setAllDirichlet();
else
values.setAllNeumann();
return values;
#endif
}
/*!
......@@ -132,8 +162,16 @@ public:
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values;
if (onLeftBoundary_(globalPos))
values[saturationIdx] = 0.8;
// if (onLeftBoundary_(globalPos))
// values[saturationIdx] = 0.8;
using WettingPhase = typename GET_PROP(TypeTag, FluidSystem)::WettingPhase;
values[saturationIdx] = 1.0;
if (isInlet(globalPos))
{
values[saturationIdx] = 0.0;
}
return values;
}
......@@ -166,7 +204,7 @@ public:
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values;
values[saturationIdx] = 0;
values[saturationIdx] = 1.0;
return values;
}
......@@ -204,7 +242,25 @@ private:
return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_;
}
bool isInlet(const GlobalPosition& globalPos) const
{
if (!onUpperBoundary_(globalPos))
return false;
for (int i = 0; i < dimWorld; i++)
{
if (globalPos[i] < inletLeftCoord_[i] - eps_)
return false;
if (globalPos[i] > inletRightCoord_[i] + eps_)
return false;
}
return true;
}
static constexpr Scalar eps_ = 1e-6;
Scalar inFlux_;
GlobalPosition inletLeftCoord_;
GlobalPosition inletRightCoord_;
};
} // end namespace Dumux
......
test/porousmediumflow/2p/sequential/impes/transportspatialparams.hh
View file @ d92cd40d
......@@ -60,24 +60,31 @@ public:
TwoPTransportSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft");
lensUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensUpperRight");
// residual saturations
lensMaterialParams_.setSwr(0.0);
lensMaterialParams_.setSnr(0.0);
outerMaterialParams_.setSwr(0.0);
outerMaterialParams_.setSnr(0.0);
// parameters for the Van Genuchten law
// alpha and n
lensMaterialParams_.setLambda(2.0);
lensMaterialParams_.setPe(0);
outerMaterialParams_.setLambda(2.0);
outerMaterialParams_.setPe(0);
lensK_ = 1.0e-7;
outerK_ = 1.0e-7;
materialLawParamsBackground_.setSwr(0.);
materialLawParamsBackground_.setSnr(0.);
materialLawParamsLenses_.setSwr(0.);
materialLawParamsLenses_.setSnr(0.);
//parameters for Brooks-Corey law
// entry pressures
materialLawParamsBackground_.setPe(getParam<Scalar>("SpatialParams.BackgroundEntryPressure", 0.0));
materialLawParamsLenses_.setPe(getParam<Scalar>("SpatialParams.LenseEntryPressure", 0.0));
materialLawParamsBackground_.setLambda(getParam<Scalar>("SpatialParams.BackgroundLambda", 3.0));
materialLawParamsLenses_.setLambda(getParam<Scalar>("SpatialParams.LenseLambda", 2.0));
permBackground_ = getParam<Scalar>("SpatialParams.BackgroundPermeability", 1e-10);
permLenses_ = getParam<Scalar>("SpatialParams.LensPermeability", 1e-10);
lensOneLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensOneLowerLeft", GlobalPosition(0.0));
lensOneUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensOneUpperRight", GlobalPosition(0.0));
lensTwoLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensTwoLowerLeft", GlobalPosition(0.0));
lensTwoUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensTwoUpperRight", GlobalPosition(0.0));
lensThreeLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensThreeLowerLeft", GlobalPosition(0.0));
lensThreeUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensThreeUpperRight", GlobalPosition(0.0));
}
/*!
......@@ -94,11 +101,16 @@ public:
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
auto globalPos = element.geometry().center();
// do not use a less permeable lens in the test with inverted wettability
if (isInLens_(element.geometry().center()))
return lensK_;
return outerK_;
if (isLensOne(globalPos))
return permLenses_;
else if (isLensTwo(globalPos))
return permLenses_;
else if (isLensThree(globalPos))
return permLenses_;