Commit 33cfd00a authored by Timo Koch's avatar Timo Koch Committed by Kilian Weishaupt
Browse files

[box][2p] Use new fluid matrix interface for box interface solver. Update tests.

parent c199687e
......@@ -71,6 +71,7 @@ An additional new option is `Vtk.CoordPrecision` which changes the precision of
- __Van Genuchten__: Corrected VanGenuchten-Mualem exponent in the nonwetting saturation formula (`1/3` instead of `1/2` (or `l`, see above))
- __Van Genuchten__: Corrected VanGenuchten-Mualem implementation of `dkrn/dSw`
- __Brooks-Corey__: Corrected Brooks-Corey implementation of `dkrn/dSw` and added the derivatives for the regularized version
- __Box interface solver__: The box interface solver immediately requires the new material law interface without deprecation period. Use the new class `BoxMaterialInterfaces` and update your spatial params to use the new fluidmatrixinteraction interface to be able to use the box interface solver in version 3.3.
- __AMGBackend__: The internal structure of the AMGBackend and the ParallelISTLHelper has been overhauled, as only used by the AMG, we did not make the changes backwards-compatible
- The global default parameters for linear solvers have been removed and moved to the class `LinearSolver`.
This only affects users that directly obtain this parameter via `getParam` somewhere in the code.
......
......@@ -22,8 +22,8 @@
* \copydoc Dumux::BoxMaterialInterfaceParams
*/
#ifndef DUMUX_2P_BOX_MATERIAL_INTERFACE_PARAMS_HH
#define DUMUX_2P_BOX_MATERIAL_INTERFACE_PARAMS_HH
#ifndef DUMUX_2P_BOX_MATERIAL_INTERFACES_HH
#define DUMUX_2P_BOX_MATERIAL_INTERFACES_HH
#include <dune/common/exceptions.hh>
......@@ -42,11 +42,19 @@ namespace Dumux {
* of freedom. On the basis of these parameters, the saturations in the
* remaining sub-control volumes connected to the vertex can be reconstructed.
*/
template<class SpatialParams>
class BoxMaterialInterfaceParams
template<class GridGeometry, class PcKrSw>
class BoxMaterialInterfaces
{
using SubControlVolume = typename GridGeometry::SubControlVolume;
public:
using MaterialLawParams = typename SpatialParams::MaterialLaw::Params;
template<class SpatialParams, class SolutionVector>
BoxMaterialInterfaces(const GridGeometry& gridGeometry,
const SpatialParams& spatialParams,
const SolutionVector& x)
{
update(gridGeometry, spatialParams, x);
}
/*!
* \brief Updates the scv -> dofparameter map
......@@ -55,30 +63,18 @@ public:
* \param spatialParams Class encapsulating the spatial parameters
* \param x The current state of the solution vector
*/
template<class GridGeometry, class SolutionVector>
template<class SpatialParams, class SolutionVector>
void update(const GridGeometry& gridGeometry,
const SpatialParams& spatialParams,
const SolutionVector& x)
{
using MaterialLaw = typename SpatialParams::MaterialLaw;
// Make sure the spatial params return a const ref and no copy!
using Elem = typename GridGeometry::GridView::template Codim<0>::Entity;
using ElemSol = decltype( elementSolution(Elem(), x, gridGeometry) );
using Scv = typename GridGeometry::SubControlVolume;
using ReturnType = decltype(spatialParams.materialLawParams(Elem(), Scv(), ElemSol()));
static_assert(std::is_lvalue_reference<ReturnType>::value,
"In order to use the box-interface solver please provide access "
"to the material law parameters via returning (const) references");
// make sure this is only called for geometries of the box method!
if (GridGeometry::discMethod != DiscretizationMethod::box)
DUNE_THROW(Dune::InvalidStateException, "Determination of the interface material parameters with "
"this class only makes sense when using the box method!");
isUpdated_ = true;
isOnMaterialInterface_.resize(gridGeometry.numDofs(), false);
dofParams_.resize(gridGeometry.numDofs(), nullptr);
pcSwAtDof_.resize(gridGeometry.numDofs(), nullptr);
for (const auto& element : elements(gridGeometry.gridView()))
{
const auto elemSol = elementSolution(element, x, gridGeometry);
......@@ -87,40 +83,43 @@ public:
fvGeometry.bind(element);
for (const auto& scv : scvs(fvGeometry))
{
const auto& params = spatialParams.materialLawParams(element, scv, elemSol);
const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol);
const auto& pcKrSw = fluidMatrixInteraction.pcSwCurve();
// assert current preconditions on requiring the spatial params to store the pckrSw curve
static_assert(std::is_lvalue_reference<typename std::decay_t<decltype(fluidMatrixInteraction)>::PcKrSwType>::value,
"In order to use the box-interface solver please provide access "
"to the material law parameters via returning (const) references");
// if no parameters had been set, set them now
if (dofParams_[scv.dofIndex()] == nullptr)
dofParams_[scv.dofIndex()] = &params;
if (!pcSwAtDof_[scv.dofIndex()])
pcSwAtDof_[scv.dofIndex()] = &pcKrSw;
// otherwise only use the current ones if endPointPc (e.g. Brooks-Corey entry pressure) is lower
else if (MaterialLaw::endPointPc( params ) < MaterialLaw::endPointPc( *(dofParams_[scv.dofIndex()]) ))
else if (pcKrSw.endPointPc() < pcSwAtDof_[scv.dofIndex()]->endPointPc())
{
dofParams_[scv.dofIndex()] = &params;
pcSwAtDof_[scv.dofIndex()] = &pcKrSw;
isOnMaterialInterface_[scv.dofIndex()] = true;
}
// keep track of material interfaces in any case
else if ( !(params == *(dofParams_[scv.dofIndex()])) )
else if ( !(pcKrSw == *(pcSwAtDof_[scv.dofIndex()])) )
isOnMaterialInterface_[scv.dofIndex()] = true;
}
}
}
//! Returns if this scv is connected to a material interface
template<class Scv>
bool isOnMaterialInterface(const Scv& scv) const
{ assert(isUpdated_); return isOnMaterialInterface_[scv.dofIndex()]; }
bool isOnMaterialInterface(const SubControlVolume& scv) const
{ return isOnMaterialInterface_[scv.dofIndex()]; }
//! Returns the material parameters associated with the dof
template<class Scv>
const MaterialLawParams& getDofParams(const Scv& scv) const
{ assert(isUpdated_); return *(dofParams_[scv.dofIndex()]); }
const PcKrSw& pcSwAtDof(const SubControlVolume& scv) const
{ return *(pcSwAtDof_[scv.dofIndex()]); }
private:
bool isUpdated_{false};
std::vector<bool> isOnMaterialInterface_;
std::vector<const MaterialLawParams*> dofParams_;
std::vector<const PcKrSw*> pcSwAtDof_;
};
} // end namespace Dumux
......
......@@ -61,8 +61,8 @@ public:
const Element& element,
const Scv& scv,
const ElemSol& elemSol,
typename ElemSol::PrimaryVariables::value_type Sn)
{ return Sn; }
typename ElemSol::PrimaryVariables::value_type sn)
{ return sn; }
};
//! Specialization for the box scheme with the interface solver enabled
......@@ -85,23 +85,25 @@ public:
const Element& element,
const Scv& scv,
const ElemSol& elemSol,
typename ElemSol::PrimaryVariables::value_type Sn)
typename ElemSol::PrimaryVariables::value_type sn)
{
// if this dof doesn't lie on a material interface, simply return Sn
if (!spatialParams.materialInterfaceParams().isOnMaterialInterface(scv))
return Sn;
const auto& materialInterfaces = spatialParams.materialInterfaces();
if (!materialInterfaces.isOnMaterialInterface(scv))
return sn;
using MaterialLaw = typename SpatialParams::MaterialLaw;
// compute capillary pressure using material parameters associated with the dof
const auto& ifMaterialParams = spatialParams.materialInterfaceParams().getDofParams(scv);
const auto pc = MaterialLaw::pc(ifMaterialParams, /*Sw=*/1.0 - Sn);
const auto& interfacePcSw = materialInterfaces.pcSwAtDof(scv);
const auto pc = interfacePcSw.pc(/*ww=*/1.0 - sn);
// reconstruct by inverting the pc-sw curve
const auto& materialLawParams = spatialParams.materialLawParams(element, scv, elemSol);
const auto pcMin = MaterialLaw::endPointPc(materialLawParams);
const auto& pcSw = spatialParams.fluidMatrixInteraction(element, scv, elemSol).pcSwCurve();
const auto pcMin = pcSw.endPointPc();
if (pc < pcMin && pcMin > 0.0) return 0.0;
else return 1.0 - MaterialLaw::sw(materialLawParams, pc);
if (pc < pcMin && pcMin > 0.0)
return 0.0;
else
return 1.0 - pcSw.sw(pc);
}
};
......
......@@ -131,7 +131,7 @@ int main(int argc, char** argv)
auto xOld = x;
// update the interface parameters
problem->spatialParams().updateMaterialInterfaceParams(x);
problem->spatialParams().updateMaterialInterfaces(x);
// the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
......
......@@ -9,6 +9,16 @@ File = grids/durlofsky.msh
[SpatialParams]
FractureAperture = 1e-3
[SpatialParams.Matrix]
Swr = 0.18
BrooksCoreyPcEntry = 1e4
BrooksCoreyLambda = 2
[SpatialParams.Fracture]
Swr = 0.05
BrooksCoreyPcEntry = 1e3
BrooksCoreyLambda = 2
[Problem]
Name = 2pboxdfm
EnableGravity = false
......@@ -28,11 +28,9 @@
#include <dumux/discretization/method.hh>
#include <dumux/material/spatialparams/fv.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#include <dumux/material/fluidmatrixinteractions/2p/brookscorey.hh>
#include <dumux/porousmediumflow/2p/boxmaterialinterfaceparams.hh>
#include <dumux/porousmediumflow/2p/boxmaterialinterfaces.hh>
namespace Dumux {
......@@ -41,7 +39,8 @@ namespace Dumux {
* \brief The spatial params for the incompressible 2p test.
*/
template<class GridGeometry, class Scalar>
class TwoPTestSpatialParams : public FVSpatialParams< GridGeometry, Scalar, TwoPTestSpatialParams<GridGeometry, Scalar> >
class TwoPTestSpatialParams
: public FVSpatialParams< GridGeometry, Scalar, TwoPTestSpatialParams<GridGeometry, Scalar> >
{
using ThisType = TwoPTestSpatialParams<GridGeometry, Scalar>;
using ParentType = FVSpatialParams<GridGeometry, Scalar, ThisType>;
......@@ -55,25 +54,16 @@ class TwoPTestSpatialParams : public FVSpatialParams< GridGeometry, Scalar, TwoP
static constexpr int dimWorld = GridView::dimensionworld;
using PcKrSw = FluidMatrix::BrooksCoreyDefault<Scalar>;
using MaterialInterfaces = BoxMaterialInterfaces<GridGeometry, PcKrSw>;
public:
using MaterialLaw = EffToAbsLaw< RegularizedBrooksCorey<Scalar> >;
using MaterialLawParams = typename MaterialLaw::Params;
using PermeabilityType = Scalar;
TwoPTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry) : ParentType(gridGeometry)
{
// residual saturations
matrixMaterialParams_.setSwr(0.18);
matrixMaterialParams_.setSnr(0.0);
fractureMaterialParams_.setSwr(0.05);
fractureMaterialParams_.setSnr(0.0);
// parameters for the Brooks-Corey law
matrixMaterialParams_.setPe(1e4);
matrixMaterialParams_.setLambda(2);
fractureMaterialParams_.setPe(1e3);
fractureMaterialParams_.setLambda(2);
}
TwoPTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
, pcKrSwMatrix_("SpatialParams.Matrix")
, pcKrSwFracture_("SpatialParams.Fracture")
{}
/*!
* \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
......@@ -115,24 +105,23 @@ public:
}
/*!
* \brief Returns the parameter object for the Brooks-Corey material law.
* \brief Returns the fluid-matrix interaction law.
*
* In this test, we use element-wise distributed material parameters.
*
* \param element The current element
* \param scv The sub-control volume inside the element.
* \param elemSol The solution at the dofs connected to the element.
* \return The material parameters object
*/
template<class ElementSolution>
const MaterialLawParams& materialLawParams(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
auto fluidMatrixInteraction(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
if (scv.isOnFracture())
return fractureMaterialParams_;
return makeFluidMatrixInteraction(pcKrSwFracture_);
else
return matrixMaterialParams_;
return makeFluidMatrixInteraction(pcKrSwMatrix_);
}
/*!
......@@ -147,21 +136,18 @@ public:
//! Updates the map of which material parameters are associated with a nodal dof.
template<class SolutionVector>
void updateMaterialInterfaceParams(const SolutionVector& x)
{
materialInterfaceParams_.update(this->gridGeometry(), *this, x);
}
void updateMaterialInterfaces(const SolutionVector& x)
{ materialInterfaces_ = std::make_unique<MaterialInterfaces>(this->gridGeometry(), *this, x); }
//! Returns the material parameters associated with a nodal dof
const BoxMaterialInterfaceParams<ThisType>& materialInterfaceParams() const
{ return materialInterfaceParams_; }
const MaterialInterfaces& materialInterfaces() const
{ return *materialInterfaces_; }
private:
MaterialLawParams matrixMaterialParams_;
MaterialLawParams fractureMaterialParams_;
PcKrSw pcKrSwMatrix_;
PcKrSw pcKrSwFracture_;
// Determines the parameters associated with the dofs at material interfaces
BoxMaterialInterfaceParams<ThisType> materialInterfaceParams_;
std::unique_ptr<MaterialInterfaces> materialInterfaces_;
static constexpr Scalar eps_ = 1.5e-7;
};
......
......@@ -151,7 +151,7 @@ int main(int argc, char** argv)
// maybe update the interface parameters
if (ENABLEINTERFACESOLVER)
problem->spatialParams().updateMaterialInterfaceParams(x);
problem->spatialParams().updateMaterialInterfaces(x);
// the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
......
......@@ -11,6 +11,16 @@ Cells = 48 32
LensLowerLeft = 1.0 2.0 # [m] coordinates of the lower left lens corner
LensUpperRight = 4.0 3.0 # [m] coordinates of the upper right lens corner
[SpatialParams.Lens]
Swr = 0.18
VanGenuchtenAlpha = 0.00045
VanGenuchtenN = 7.3
[SpatialParams.Outer]
Swr = 0.05
VanGenuchtenAlpha = 0.0037
VanGenuchtenN = 4.7
[Problem]
Name = 2p
EnableGravity = true
......
......@@ -26,10 +26,8 @@
#define DUMUX_INCOMPRESSIBLE_TWOP_TEST_SPATIAL_PARAMS_HH
#include <dumux/material/spatialparams/fv.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#include <dumux/porousmediumflow/2p/boxmaterialinterfaceparams.hh>
#include <dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh>
#include <dumux/porousmediumflow/2p/boxmaterialinterfaces.hh>
namespace Dumux {
......@@ -51,36 +49,24 @@ class TwoPTestSpatialParams
static constexpr int dimWorld = GridView::dimensionworld;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using EffectiveLaw = RegularizedVanGenuchten<Scalar>;
using PcKrSw = FluidMatrix::VanGenuchtenDefault<Scalar>;
using MaterialInterfaces = BoxMaterialInterfaces<GridGeometry, PcKrSw>;
public:
using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
using MaterialLawParams = typename MaterialLaw::Params;
using PermeabilityType = Scalar;
TwoPTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
, lensPcKrSw_("SpatialParams.Lens")
, outerPcKrSw_("SpatialParams.Outer")
{
lensIsOilWet_ = getParam<bool>("SpatialParams.LensIsOilWet", false);
lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft");
lensUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensUpperRight");
// residual saturations
lensMaterialParams_.setSwr(0.18);
lensMaterialParams_.setSnr(0.0);
outerMaterialParams_.setSwr(0.05);
outerMaterialParams_.setSnr(0.0);
// parameters for the Van Genuchten law
// alpha and n
lensMaterialParams_.setVgAlpha(0.00045);
lensMaterialParams_.setVgn(7.3);
outerMaterialParams_.setVgAlpha(0.0037);
outerMaterialParams_.setVgn(4.7);
lensK_ = getParam<Scalar>("SpatialParams.lensK", 9.05e-12);
outerK_ = getParam<Scalar>("SpatialParams.outerK", 4.6e-10);
lensK_ = getParam<Scalar>("SpatialParams.Lens.Permeability", 9.05e-12);
outerK_ = getParam<Scalar>("SpatialParams.Outer.Permeability", 4.6e-10);
}
/*!
......@@ -123,14 +109,14 @@ public:
* \return The material parameters object
*/
template<class ElementSolution>
const MaterialLawParams& materialLawParams(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
auto fluidMatrixInteraction(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
// do not use different parameters in the test with inverted wettability
if (isInLens_(element.geometry().center()) && !lensIsOilWet_)
return lensMaterialParams_;
return outerMaterialParams_;
return makeFluidMatrixInteraction(lensPcKrSw_);
return makeFluidMatrixInteraction(outerPcKrSw_);
}
/*!
......@@ -149,15 +135,15 @@ public:
//! Updates the map of which material parameters are associated with a nodal dof.
template<class SolutionVector>
void updateMaterialInterfaceParams(const SolutionVector& x)
void updateMaterialInterfaces(const SolutionVector& x)
{
if (GridGeometry::discMethod == DiscretizationMethod::box)
materialInterfaceParams_.update(this->gridGeometry(), *this, x);
materialInterfaces_ = std::make_unique<MaterialInterfaces>(this->gridGeometry(), *this, x);
}
//! Returns the material parameters associated with a nodal dof
const BoxMaterialInterfaceParams<ThisType>& materialInterfaceParams() const
{ return materialInterfaceParams_; }
const MaterialInterfaces& materialInterfaces() const
{ return *materialInterfaces_; }
//! Returns whether or not the lens is oil wet
bool lensIsOilWet() const { return lensIsOilWet_; }
......@@ -178,11 +164,11 @@ private:
Scalar lensK_;
Scalar outerK_;
MaterialLawParams lensMaterialParams_;
MaterialLawParams outerMaterialParams_;
// Determines the parameters associated with the dofs at material interfaces
BoxMaterialInterfaceParams<ThisType> materialInterfaceParams_;
const PcKrSw lensPcKrSw_;
const PcKrSw outerPcKrSw_;
std::unique_ptr<MaterialInterfaces> materialInterfaces_;
static constexpr Scalar eps_ = 1.5e-7;
};
......
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