From 33cfd00aa53016258cf18524377a62aec116e9ff Mon Sep 17 00:00:00 2001 From: Timo Koch <timo.koch@iws.uni-stuttgart.de> Date: Tue, 27 Oct 2020 18:14:01 +0100 Subject: [PATCH] [box][2p] Use new fluid matrix interface for box interface solver. Update tests. --- CHANGELOG.md | 1 + ...faceparams.hh => boxmaterialinterfaces.hh} | 65 +++++++++---------- .../2p/saturationreconstruction.hh | 26 ++++---- .../2p/implicit/boxdfm/main.cc | 2 +- .../2p/implicit/boxdfm/params.input | 10 +++ .../2p/implicit/boxdfm/spatialparams.hh | 62 +++++++----------- .../2p/implicit/incompressible/main.cc | 2 +- .../2p/implicit/incompressible/params.input | 10 +++ .../implicit/incompressible/spatialparams.hh | 56 ++++++---------- 9 files changed, 114 insertions(+), 120 deletions(-) rename dumux/porousmediumflow/2p/{boxmaterialinterfaceparams.hh => boxmaterialinterfaces.hh} (66%) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1dfadefcaa..14bdc03df9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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. diff --git a/dumux/porousmediumflow/2p/boxmaterialinterfaceparams.hh b/dumux/porousmediumflow/2p/boxmaterialinterfaces.hh similarity index 66% rename from dumux/porousmediumflow/2p/boxmaterialinterfaceparams.hh rename to dumux/porousmediumflow/2p/boxmaterialinterfaces.hh index 5dcb9055d4..3296cd5444 100644 --- a/dumux/porousmediumflow/2p/boxmaterialinterfaceparams.hh +++ b/dumux/porousmediumflow/2p/boxmaterialinterfaces.hh @@ -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()] = ¶ms; + 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()] = ¶ms; + 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 diff --git a/dumux/porousmediumflow/2p/saturationreconstruction.hh b/dumux/porousmediumflow/2p/saturationreconstruction.hh index ef2b49216b..10a49b2644 100644 --- a/dumux/porousmediumflow/2p/saturationreconstruction.hh +++ b/dumux/porousmediumflow/2p/saturationreconstruction.hh @@ -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); } }; diff --git a/test/porousmediumflow/2p/implicit/boxdfm/main.cc b/test/porousmediumflow/2p/implicit/boxdfm/main.cc index 6fe43aa027..e59caef33e 100644 --- a/test/porousmediumflow/2p/implicit/boxdfm/main.cc +++ b/test/porousmediumflow/2p/implicit/boxdfm/main.cc @@ -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>; diff --git a/test/porousmediumflow/2p/implicit/boxdfm/params.input b/test/porousmediumflow/2p/implicit/boxdfm/params.input index 18cb950644..6705416a25 100644 --- a/test/porousmediumflow/2p/implicit/boxdfm/params.input +++ b/test/porousmediumflow/2p/implicit/boxdfm/params.input @@ -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 diff --git a/test/porousmediumflow/2p/implicit/boxdfm/spatialparams.hh b/test/porousmediumflow/2p/implicit/boxdfm/spatialparams.hh index 5d0455254f..8183314ecb 100644 --- a/test/porousmediumflow/2p/implicit/boxdfm/spatialparams.hh +++ b/test/porousmediumflow/2p/implicit/boxdfm/spatialparams.hh @@ -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; }; diff --git a/test/porousmediumflow/2p/implicit/incompressible/main.cc b/test/porousmediumflow/2p/implicit/incompressible/main.cc index 0d7c8192bb..3824474246 100644 --- a/test/porousmediumflow/2p/implicit/incompressible/main.cc +++ b/test/porousmediumflow/2p/implicit/incompressible/main.cc @@ -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>; diff --git a/test/porousmediumflow/2p/implicit/incompressible/params.input b/test/porousmediumflow/2p/implicit/incompressible/params.input index 336df6cf83..17b87d2465 100644 --- a/test/porousmediumflow/2p/implicit/incompressible/params.input +++ b/test/porousmediumflow/2p/implicit/incompressible/params.input @@ -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 diff --git a/test/porousmediumflow/2p/implicit/incompressible/spatialparams.hh b/test/porousmediumflow/2p/implicit/incompressible/spatialparams.hh index 1232ca5d2b..fdb750dcb2 100644 --- a/test/porousmediumflow/2p/implicit/incompressible/spatialparams.hh +++ b/test/porousmediumflow/2p/implicit/incompressible/spatialparams.hh @@ -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; }; -- GitLab