diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1dfadefcaade3aaa0500d26ed3512b7443da3b8e..311f823640d6235ab5144044dc254f993f85f27e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -15,14 +15,18 @@ Differences Between DuMux 3.3 and DuMux 3.2
In order to use this, include ``. `format`, `format_to`, `format_to_n`, `formatted_size` are available in the `Dumux::Fmt` namespace.
The string formatting is documented [here](https://en.cppreference.com/w/cpp/utility/format/formatter#Standard_format_specification) and follows the Python string formatting rules.
The function are documented on [cppreference](https://en.cppreference.com/w/cpp/utility/format).
+ - The usage of laws for pc-Sw and kr-Sw has been completely revised. A caller does not have to pass a `parameters` object to the laws anymore. The `spatialParams` now provide a `fluidMatrixInteraction` function which bundles an arbitrary number of
+ different interaction laws such as a pc-Sw and kr-Sw curve and interfacial areas.
+ New pre-cached spline laws were added which can help to increase efficiency. The usage of the old interface is deprecated and warnings will be raised. The old interface will be removed after the release of 3.3.
### Immediate interface changes not allowing/requiring a deprecation period:
-
+- __Flash/Constraintsolver__: The flashes depending on material laws are immediately required to use new-style material laws (fluidMatrixInteraction interface in spatialparams)
+- __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.
- For the "sequential" models, the property `BoundaryTypes` has been simply renamed to `SequentialBoundaryTypes`
- __Quadmath__: Dumux::Quad has been removed without deprecation. Use Dune::Float128 instead.
- Within the RANS group, two additional runtime parameters have been included 'IsFlatWallBounded' and 'WriteFlatWallBoundedFields'.
For both the K-Epsilon and Zero-eq RANS models the 'IsFlatWallBounded' runtime parameter should be set as True,
-as wall topology is not supported for these models with our geometric contraints. If not set as true, the geometry
+as wall topology is not supported for these models with our geometric constraints. If not set as true, the geometry
will be checked before the model is run. If either the runtime parameter or the geometry check indicate non-flat walls,
the model will terminate. To add FlatWallBounded specific output to the vtk output, WriteFlatWallBoundedFields can be set as True.
- __1d3d coupling__: The kernel coupling manager has been replaced with the one from Koch et al (2020) JCP https://doi.org/10.1016/j.jcp.2020.109370
diff --git a/doc/handbook/dumux-handbook.bib b/doc/handbook/dumux-handbook.bib
index 9b5588c9ac46388b98565dc4233d2b0c85e0192e..a3b6a3ecb47963c6146aa25bee52984a00febc0c 100644
--- a/doc/handbook/dumux-handbook.bib
+++ b/doc/handbook/dumux-handbook.bib
@@ -426,6 +426,17 @@
url = {http://elib.uni-stuttgart.de/opus/volltexte/2013/8141}
}
+@PhdThesis{nuske2014,
+ author = {Philipp Nuske},
+ title = {{Beyond local equilibrium : relaxing local equilibrium assumptions in multiphase flow in porous media}},
+ school = {Universit\"at Stuttgart},
+ year = {2014},
+ address = {Holzgartenstr. 16, 70174 Stuttgart},
+ isbn = {978-3-942036-41-2},
+ language = {eng},
+ url = {https://elib.uni-stuttgart.de/handle/11682/614}
+}
+
@ARTICLE{A3:Braun:2002,
author = {Braun, C. and Helmig, R. and Manthey, S.},
title = {{Determination of constitutive relationships for two-phase flow processes
diff --git a/dumux/common/deprecated.hh b/dumux/common/deprecated.hh
index 3e191a6891295a010468ba4c9be6dc07b0b3f6d2..9164b15ce1c0fcf98e0a82bb7ba063a85e8f84c8 100644
--- a/dumux/common/deprecated.hh
+++ b/dumux/common/deprecated.hh
@@ -26,6 +26,9 @@
#define DUMUX_COMMON_DEPRECATED_HH
#include
+#include
+#include
+#include
namespace Dumux {
@@ -43,9 +46,521 @@ namespace Deprecated {
////////////////////////////////////////////
///////////////////////////////////////////////////////////////
-// Deprecation warnings for effective diffusion coefficients //
+// Deprecation warnings for the new material law //
///////////////////////////////////////////////////////////////
+// support old interface of the effective thermal conductivity laws
+template
+struct HasNewFIAIF
+{
+ template
+ auto operator()(S&& sp)
+ -> decltype(sp.fluidMatrixInteraction(std::declval(),
+ std::declval(),
+ std::declval())) {}
+};
+
+template
+struct HasNewFIAIFAtPos
+{
+ template
+ auto operator()(S&& sp)
+ -> decltype(sp.fluidMatrixInteractionAtPos(std::declval())) {}
+};
+
+
+template
+class PcKrSwHelper : public FluidMatrix::Adapter, FluidMatrix::PcKrSw>
+{
+public:
+ using Scalar = ScalarT;
+
+ // pass scalar so template arguments can all be deduced
+ PcKrSwHelper(const Scalar& scalar,
+ const SpatialParams& sp,
+ const Element& element,
+ const Scv& scv,
+ const ElemSol& elemSol)
+ : spatialParams_(sp), element_(element), scv_(scv), elemSol_(elemSol)
+ {}
+
+ Scalar krw(const Scalar sw) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::krw(params, sw);
+ }
+
+ Scalar krn(const Scalar sw) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::krn(params, sw);
+ }
+
+ Scalar pc(const Scalar sw) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::pc(params, sw);
+ }
+
+ Scalar dpc_dsw(const Scalar sw) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::dpc_dsw(params, sw);
+ }
+
+ Scalar endPointPc() const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::endPointPc(params);
+ }
+
+ Scalar sw(const Scalar pc) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::sw(params, pc);
+ }
+
+ Scalar dsw_dpc(const Scalar pc) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::dsw_dpc(params, pc);
+ }
+
+ Scalar dkrw_dsw(const Scalar sw) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::dkrw_dsw(params, sw);
+ }
+
+ Scalar dkrn_dsw(const Scalar sw) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::dkrn_dsw(params, sw);
+ }
+
+ const auto& basicParams() const
+ { return spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_); }
+
+ const auto& effToAbsParams() const
+ { return spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_); }
+
+private:
+ const SpatialParams& spatialParams_;
+ const Element& element_;
+ const Scv& scv_;
+ const ElemSol& elemSol_;
+};
+
+template
+class PcKrSwThreePHelper : public FluidMatrix::Adapter, FluidMatrix::ThreePhasePcKrSw>
+{
+public:
+ using Scalar = ScalarT;
+
+ // pass scalar so template arguments can all be deduced
+ PcKrSwThreePHelper(const Scalar& scalar,
+ const SpatialParams& sp,
+ const Element& element,
+ const Scv& scv,
+ const ElemSol& elemSol)
+ : spatialParams_(sp), element_(element), scv_(scv), elemSol_(elemSol)
+ {}
+
+ Scalar pcgw(const Scalar sw, const Scalar /*dummySn*/) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::pcgw(params, sw);
+ }
+
+ Scalar pcnw(const Scalar sw, const Scalar /*dummySn*/) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::pcnw(params, sw);
+ }
+
+ Scalar pcgn(const Scalar sw, const Scalar sn) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::pcgn(params, sw + sn);
+ }
+
+ Scalar pcAlpha(const Scalar /*dummySw*/, const Scalar sn) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::pcAlpha(params, sn);
+ }
+
+ Scalar krw(const Scalar sw, const Scalar sn) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::krw(params, sw, sn);
+ }
+
+ Scalar krn(const Scalar sw, const Scalar sn) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::krw(params, sw, sn);
+ }
+
+ Scalar krg(const Scalar sw, const Scalar sn) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::krg(params, sw, sn);
+ }
+
+ Scalar kr(const int phaseIdx, const Scalar sw, const Scalar sn) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ return SpatialParams::MaterialLaw::kr(params, phaseIdx, sw, sn, 1 - sw - sn);
+ }
+
+ const auto& basicParams() const
+ { return spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_); }
+
+ const auto& effToAbsParams() const
+ { return spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_); }
+
+private:
+ const SpatialParams& spatialParams_;
+ const Element& element_;
+ const Scv& scv_;
+ const ElemSol& elemSol_;
+};
+
+// for implicit models
+template
+auto makePcKrSw(const Scalar& scalar,
+ const SpatialParams& sp,
+ const Element& element,
+ const Scv& scv,
+ const ElemSol& elemSol)
+{
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ constexpr bool hasNew = decltype(isValid(HasNewFIAIF()).template check())::value;
+ constexpr bool hasNewAtPos = decltype(isValid(HasNewFIAIFAtPos()).template check())::value;
+ if constexpr (hasNew)
+ return sp.fluidMatrixInteraction(element, scv, elemSol);
+ else if constexpr (hasNewAtPos)
+ return sp.fluidMatrixInteractionAtPos(scv.center());
+ else
+ {
+ if constexpr (numPhases == 2)
+ return makeFluidMatrixInteraction(PcKrSwHelper(scalar, sp, element, scv, elemSol));
+ else
+ return makeFluidMatrixInteraction(PcKrSwThreePHelper(scalar, sp, element, scv, elemSol));
+ }
+}
+
+// for sequential models
+template
+auto makePcKrSw(const Scalar& scalar,
+ const SpatialParams& sp,
+ const Element& element)
+{
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ constexpr bool hasNewAtPos = decltype(isValid(HasNewFIAIFAtPos()).template check())::value;
+ if constexpr (hasNewAtPos)
+ return sp.fluidMatrixInteractionAtPos(element.geometry().center());
+ else
+ {
+ using DummyScv = int;
+ using DummyElemSol = int;
+ return makeFluidMatrixInteraction(PcKrSwHelper(scalar, sp, element, DummyScv(), DummyElemSol()));
+ }
+}
+
+
+///////////////////////////////////////////////////////////////
+// Deprecation warnings for the mp material law stuff //
+///////////////////////////////////////////////////////////////
+
+template
+class PcKrSwMPHelper : public FluidMatrix::Adapter, FluidMatrix::MultiPhasePcKrSw>
+{
+ using MaterialLaw = typename SpatialParams::MaterialLaw;
+ using MPAdapter = Dumux::MPAdapter;
+public:
+ using Scalar = ScalarT;
+
+ // pass scalar so template arguments can all be deduced
+ PcKrSwMPHelper(const Scalar& scalar,
+ const SpatialParams& sp,
+ const Element& element,
+ const Scv& scv,
+ const ElemSol& elemSol,
+ const NumPhases& np)
+ : spatialParams_(sp), element_(element), scv_(scv), elemSol_(elemSol)
+ {}
+
+ template
+ auto capillaryPressures(const FluidState& fs, int wPhaseIdx) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ Dune::FieldVector pc;
+ MPAdapter::capillaryPressures(pc, params, fs, wPhaseIdx);
+ return pc;
+ }
+
+ template
+ auto relativePermeabilities(const FluidState& fs, int wPhaseIdx) const
+ {
+ const auto& params = spatialParams_.materialLawParamsDeprecated(element_, scv_, elemSol_);
+ Dune::FieldVector kr;
+ MPAdapter::capillaryPressures(kr, params, fs, wPhaseIdx);
+ return kr;
+ }
+
+private:
+ const SpatialParams& spatialParams_;
+ const Element& element_;
+ const Scv& scv_;
+ const ElemSol& elemSol_;
+};
+
+template
+auto makeMPPcKrSw(const Scalar& scalar,
+ const SpatialParams& sp,
+ const Element& element,
+ const Scv& scv,
+ const ElemSol& elemSol,
+ const NumPhases& np)
+{
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ constexpr bool hasNew = decltype(isValid(HasNewFIAIF()).template check())::value;
+ constexpr bool hasNewAtPos = decltype(isValid(HasNewFIAIFAtPos()).template check())::value;
+ if constexpr (hasNew)
+ return sp.fluidMatrixInteraction(element, scv, elemSol);
+ else if constexpr (hasNewAtPos)
+ return sp.fluidMatrixInteractionAtPos(scv.center());
+ else
+ return makeFluidMatrixInteraction(PcKrSwMPHelper(scalar, sp, element, scv, elemSol, np));
+}
+
+///////////////////////////////////////////////////////////////
+// Deprecation warnings for the kinetic surface areas //
+///////////////////////////////////////////////////////////////
+
+// support old interface of surface area
+template
+struct HasNewAns
+{
+ template
+ auto operator()(S&& sp)
+ -> decltype(sp.nonwettingSolidInterfacialArea(std::declval(),
+ std::declval(),
+ std::declval())) {}
+};
+
+template
+struct HasNewAnsAtPos
+{
+ template
+ auto operator()(S&& sp)
+ -> decltype(sp.nonwettingSolidInterfacialAreaAtPos(std::declval())) {}
+};
+
+// support old interface of surface area
+template
+struct HasNewAnw
+{
+ template
+ auto operator()(S&& sp)
+ -> decltype(sp.wettingNonwettingInterfacialArea(std::declval(),
+ std::declval(),
+ std::declval())) {}
+};
+
+template
+struct HasNewAnwAtPos
+{
+ template
+ auto operator()(S&& sp)
+ -> decltype(sp.wettingNonwettingInterfacialAreaAtPos(std::declval())) {}
+};
+
+// support old interface of surface area
+template
+struct HasNewAws
+{
+ template
+ auto operator()(S&& sp)
+ -> decltype(sp.wettingSolidInterfacialArea(std::declval(),
+ std::declval(),
+ std::declval())) {}
+};
+
+template
+struct HasNewAwsAtPos
+{
+ template
+ auto operator()(S&& sp)
+ -> decltype(sp.wettingSolidInterfacialAreaAtPos(std::declval())) {}
+};
+
+template
+class WettingNonwettingInterfacialArea : public FluidMatrix::Adapter, FluidMatrix::WettingNonwettingInterfacialAreaPcSw>
+{
+public:
+ using Scalar = ScalarT;
+
+ WettingNonwettingInterfacialArea(const Scalar& scalar,
+ const SpatialParams& sp,
+ const Element& element,
+ const Scv& scv,
+ const ElemSol& elemSol)
+ : spatialParams_(sp), element_(element), scv_(scv), elemSol_(elemSol)
+ {}
+
+ const auto& basicParams() const
+ {
+ return spatialParams_.aWettingNonWettingSurfaceParams(element_, scv_, elemSol_);
+ }
+
+ Scalar area(const Scalar sw, const Scalar pc) const
+ {
+ const auto& surfaceParams = spatialParams_.aWettingNonWettingSurfaceParams(element_, scv_, elemSol_);
+ const auto& materialParams = spatialParams_.materialLawParams(element_, scv_, elemSol_);
+ using AwnSurface = typename SpatialParams::AwnSurface;
+ return AwnSurface::interfacialArea(surfaceParams, materialParams, sw, pc);
+ }
+
+ Scalar darea_dpc(const Scalar sw, const Scalar pc)
+ {
+ const auto& surfaceParams = spatialParams_.aWettingNonWettingSurfaceParams(element_, scv_, elemSol_);
+ using AwnSurface = typename SpatialParams::AwnSurface;
+ return AwnSurface::dawn_dpc(surfaceParams, sw, pc);
+ }
+
+ Scalar darea_dsw(const Scalar sw, const Scalar pc)
+ {
+ const auto& surfaceParams = spatialParams_.aWettingNonWettingSurfaceParams(element_, scv_, elemSol_);
+ using AwnSurface = typename SpatialParams::AwnSurface;
+ return AwnSurface::dawn_dsw(surfaceParams, sw, pc);
+ }
+
+private:
+ const SpatialParams& spatialParams_;
+ const Element& element_;
+ const Scv& scv_;
+ const ElemSol& elemSol_;
+};
+
+template
+class NonwettingSolidInterfacialArea : public FluidMatrix::Adapter, FluidMatrix::NonwettingSolidInterfacialAreaPcSw>
+{
+public:
+ using Scalar = ScalarT;
+
+ NonwettingSolidInterfacialArea(const Scalar& scalar,
+ const SpatialParams& sp,
+ const Element& element,
+ const Scv& scv,
+ const ElemSol& elemSol)
+ : spatialParams_(sp), element_(element), scv_(scv), elemSol_(elemSol)
+ {}
+
+ const auto& basicParams() const
+ {
+ return spatialParams_.aNonWettingSolidSurfaceParams(element_, scv_, elemSol_);
+ }
+
+ Scalar area(const Scalar sw, const Scalar pc) const
+ {
+ const auto& surfaceParams = spatialParams_.aNonWettingSolidSurfaceParams(element_, scv_, elemSol_);
+ const auto& materialParams = spatialParams_.materialLawParams(element_, scv_, elemSol_);
+ using AnsSurface = typename SpatialParams::AnsSurface;
+ return AnsSurface::interfacialArea(surfaceParams, materialParams, sw, pc);
+ }
+
+ Scalar darea_dpc(const Scalar sw, const Scalar pc)
+ {
+ const auto& surfaceParams = spatialParams_.aNonWettingSolidSurfaceParams(element_, scv_, elemSol_);
+ using AnsSurface = typename SpatialParams::AnsSurface;
+ return AnsSurface::dawn_dpc(surfaceParams, sw, pc);
+ }
+
+ Scalar darea_dsw(const Scalar sw, const Scalar pc)
+ {
+ const auto& surfaceParams = spatialParams_.aNonWettingSolidSurfaceParams(element_, scv_, elemSol_);
+ using AnsSurface = typename SpatialParams::AnsSurface;
+ return AnsSurface::dawn_dsw(surfaceParams, sw, pc);
+ }
+
+private:
+ const SpatialParams& spatialParams_;
+ const Element& element_;
+ const Scv& scv_;
+ const ElemSol& elemSol_;
+};
+
+template
+class WettingSolidInterfacialArea : public FluidMatrix::Adapter, FluidMatrix::WettingSolidInterfacialAreaPcSw>
+{
+public:
+ using Scalar = ScalarT;
+
+ WettingSolidInterfacialArea(const Scalar& scalar,
+ const SpatialParams& sp,
+ const Element& element,
+ const Scv& scv,
+ const ElemSol& elemSol)
+ : spatialParams_(sp), element_(element), scv_(scv), elemSol_(elemSol)
+ {}
+
+ const auto& basicParams() const
+ {
+ return spatialParams_.aWettingSolidSurfaceParams(element_, scv_, elemSol_);
+ }
+
+ Scalar area(const Scalar sw, const Scalar pc) const
+ {
+ const auto& surfaceParams = spatialParams_.aWettingSolidSurfaceParams(element_, scv_, elemSol_);
+ const auto& materialParams = spatialParams_.materialLawParams(element_, scv_, elemSol_);
+ using AwsSurface = typename SpatialParams::AwsSurface;
+ return AwsSurface::interfacialArea(surfaceParams, materialParams, sw, pc);
+ }
+
+ Scalar darea_dpc(const Scalar sw, const Scalar pc)
+ {
+ const auto& surfaceParams = spatialParams_.aWettingSolidSurfaceParams(element_, scv_, elemSol_);
+ using AwsSurface = typename SpatialParams::AwsSurface;
+ return AwsSurface::dawn_dpc(surfaceParams, sw, pc);
+ }
+
+ Scalar darea_dsw(const Scalar sw, const Scalar pc)
+ {
+ const auto& surfaceParams = spatialParams_.aWettingSolidSurfaceParams(element_, scv_, elemSol_);
+ using AwsSurface = typename SpatialParams::AwsSurface;
+ return AwsSurface::dawn_dsw(surfaceParams, sw, pc);
+ }
+
+private:
+ const SpatialParams& spatialParams_;
+ const Element& element_;
+ const Scv& scv_;
+ const ElemSol& elemSol_;
+};
+
+template
+auto makeInterfacialArea(const Scalar& scalar,
+ const SpatialParams& sp,
+ const Element& element,
+ const Scv& scv,
+ const ElemSol& elemSol)
+{
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ constexpr bool hasNew = decltype(isValid(HasNewFIAIF()).template check())::value;
+ constexpr bool hasNewAtPos = decltype(isValid(HasNewFIAIFAtPos()).template check())::value;
+ if constexpr (hasNew)
+ return sp.fluidMatrixInteraction(element, scv, elemSol);
+ else if constexpr (hasNewAtPos)
+ return sp.fluidMatrixInteractionAtPos(scv.center());
+ else
+ return makeFluidMatrixInteraction(WettingNonwettingInterfacialArea(scalar, sp, element, scv, elemSol),
+ NonwettingSolidInterfacialArea(scalar, sp, element, scv, elemSol),
+ WettingSolidInterfacialArea(scalar, sp, element, scv, elemSol));
+}
+
} // end namespace Deprecated
#endif
diff --git a/dumux/common/monotonecubicspline.hh b/dumux/common/monotonecubicspline.hh
index 9fa462e63486c9059a5b1107842ede1dd9694ea9..75abcfd61f099df959a3b3512d763a5824f664ff 100644
--- a/dumux/common/monotonecubicspline.hh
+++ b/dumux/common/monotonecubicspline.hh
@@ -77,7 +77,7 @@ public:
// check some requirements
assert (x.size() == y.size());
assert (x.size() >=2);
- assert (std::is_sorted(x.begin(), x.end()));
+ assert (std::is_sorted(x.begin(), x.end()) || std::is_sorted(x.rbegin(), x.rend()));
// save a copy of the control points
x_ = x;
@@ -87,7 +87,8 @@ public:
numPoints_ = x.size();
// whether we are increasing
- increasing_ = y_.back() > y_.front();
+ increasingX_ = x_.back() > x_.front();
+ increasingY_ = y_.back() > y_.front();
// the slope at every control point
m_.resize(numPoints_);
@@ -114,21 +115,12 @@ public:
*/
Scalar eval(const Scalar x) const
{
- if (x <= x_.front())
+ if ((x <= x_.front() && increasingX_) || (x >= x_.front() && !increasingX_))
return y_.front() + m_.front()*(x - x_.front());
- else if (x > x_.back())
+ else if ((x > x_.back() && increasingX_) || (x < x_.back() && !increasingX_))
return y_.back() + m_.back()*(x - x_.back());
- else
- {
- const auto lookUpIndex = std::distance(x_.begin(), std::lower_bound(x_.begin(), x_.end(), x));
- assert(lookUpIndex != 0);
- // interpolate parametrization parameter t in [0,1]
- const auto h = (x_[lookUpIndex] - x_[lookUpIndex-1]);
- const auto t = (x - x_[lookUpIndex-1])/h;
- return y_[lookUpIndex-1]*Basis::h00(t) + h*m_[lookUpIndex-1]*Basis::h10(t)
- + y_[lookUpIndex]*Basis::h01(t) + h*m_[lookUpIndex]*Basis::h11(t);
- }
+ return eval_(x);
}
/*!
@@ -138,22 +130,12 @@ public:
*/
Scalar evalDerivative(const Scalar x) const
{
- if (x <= x_.front())
+ if ((x <= x_.front() && increasingX_) || (x >= x_.front() && !increasingX_))
return m_.front();
- else if (x > x_.back())
+ else if ((x > x_.back() && increasingX_) || (x < x_.back() && !increasingX_))
return m_.back();
- else
- {
- const auto lookUpIndex = std::distance(x_.begin(), std::lower_bound(x_.begin(), x_.end(), x));
- assert(lookUpIndex != 0);
- // interpolate parametrization parameter t in [0,1]
- const auto h = (x_[lookUpIndex] - x_[lookUpIndex-1]);
- const auto t = (x - x_[lookUpIndex-1])/h;
- const auto dtdx = 1.0/h;
- return y_[lookUpIndex-1]*Basis::dh00(t)*dtdx + m_[lookUpIndex-1]*Basis::dh10(t)
- + y_[lookUpIndex]*Basis::dh01(t)*dtdx + m_[lookUpIndex]*Basis::dh11(t);
- }
+ return evalDerivative_(x);
}
/*!
@@ -164,40 +146,39 @@ public:
*/
Scalar evalInverse(const Scalar y) const
{
- if (increasing_)
- {
- if (y <= y_.front())
- return x_.front() + (y - y_.front())/m_.front();
- else if (y > y_.back())
- return x_.back() + (y - y_.back())/m_.back();
- else
- {
- const auto lookUpIndex = std::distance(y_.begin(), std::lower_bound(y_.begin(), y_.end(), y));
- assert(lookUpIndex != 0);
-
- return evalInverse_(y, lookUpIndex);
- }
- }
+ if ((y <= y_.front() && increasingY_) || (y >= y_.front() && !increasingY_))
+ return x_.front() + (y - y_.front())/m_.front();
+ else if ((y > y_.back() && increasingY_) || (y < y_.back() && !increasingY_))
+ return x_.back() + (y - y_.back())/m_.back();
- else
- {
- if (y >= y_.front())
- return x_.front() + (y - y_.front())/m_.front();
- else if (y < y_.back())
- return x_.back() + (y - y_.back())/m_.back();
- else
- {
- const auto lookUpIndex = y_.size() - std::distance(y_.rbegin(), std::lower_bound(y_.rbegin(), y_.rend(), y));
- assert(lookUpIndex != 0);
-
- return evalInverse_(y, lookUpIndex);
- }
- }
+ return evalInverse_(y);
}
private:
- Scalar evalInverse_(const Scalar y, const std::size_t lookUpIndex) const
+ Scalar eval_(const Scalar x) const
{
+ // interpolate parametrization parameter t in [0,1]
+ const auto lookUpIndex = lookUpIndex_(x_, x, increasingX_);
+ const auto h = (x_[lookUpIndex] - x_[lookUpIndex-1]);
+ const auto t = (x - x_[lookUpIndex-1])/h;
+ return y_[lookUpIndex-1]*Basis::h00(t) + h*m_[lookUpIndex-1]*Basis::h10(t)
+ + y_[lookUpIndex]*Basis::h01(t) + h*m_[lookUpIndex]*Basis::h11(t);
+ }
+
+ Scalar evalDerivative_(const Scalar x) const
+ {
+ // interpolate parametrization parameter t in [0,1]
+ const auto lookUpIndex = lookUpIndex_(x_, x, increasingX_);
+ const auto h = (x_[lookUpIndex] - x_[lookUpIndex-1]);
+ const auto t = (x - x_[lookUpIndex-1])/h;
+ const auto dtdx = 1.0/h;
+ return y_[lookUpIndex-1]*Basis::dh00(t)*dtdx + m_[lookUpIndex-1]*Basis::dh10(t)
+ + y_[lookUpIndex]*Basis::dh01(t)*dtdx + m_[lookUpIndex]*Basis::dh11(t);
+ }
+
+ Scalar evalInverse_(const Scalar y) const
+ {
+ const auto lookUpIndex = lookUpIndex_(y_, y, increasingY_);
auto localPolynomial = [&](const auto x) {
// interpolate parametrization parameter t in [0,1]
const auto h = (x_[lookUpIndex] - x_[lookUpIndex-1]);
@@ -211,11 +192,31 @@ private:
return findScalarRootBrent(x_[lookUpIndex-1]-eps, x_[lookUpIndex]+eps, localPolynomial);
}
+ auto lookUpIndex_(const std::vector& vec, const Scalar v, bool increasing) const
+ {
+ return increasing ? lookUpIndexIncreasing_(vec, v) : lookUpIndexDecreasing_(vec, v);
+ }
+
+ auto lookUpIndexIncreasing_(const std::vector& vec, const Scalar v) const
+ {
+ const auto lookUpIndex = std::distance(vec.begin(), std::lower_bound(vec.begin(), vec.end(), v));
+ assert(lookUpIndex != 0);
+ return lookUpIndex;
+ }
+
+ auto lookUpIndexDecreasing_(const std::vector& vec, const Scalar v) const
+ {
+ const auto lookUpIndex = vec.size() - std::distance(vec.rbegin(), std::lower_bound(vec.rbegin(), vec.rend(), v));
+ assert(lookUpIndex != 0);
+ return lookUpIndex;
+ }
+
std::vector x_; //!< the x-coordinates
std::vector y_; //!< the y-coordinates
std::vector m_; //!< the slope for each control point
std::size_t numPoints_; //!< the number of control points
- bool increasing_; //!< if we are increasing monotone or not
+ bool increasingX_; //!< if we are increasing monotone or not
+ bool increasingY_; //!< if we are increasing monotone or not
};
} // end namespace Dumux
diff --git a/dumux/io/CMakeLists.txt b/dumux/io/CMakeLists.txt
index 6032f37994f30212e869b1f36b982781e8679254..2a491c8b784acd8c59e957097192194b075a8343 100644
--- a/dumux/io/CMakeLists.txt
+++ b/dumux/io/CMakeLists.txt
@@ -14,6 +14,7 @@ name.hh
ploteffectivediffusivitymodel.hh
plotmateriallaw.hh
plotmateriallaw3p.hh
+plotpckrsw.hh
plotthermalconductivitymodel.hh
pointcloudvtkwriter.hh
rasterimagereader.hh
diff --git a/dumux/io/plotmateriallaw.hh b/dumux/io/plotmateriallaw.hh
index 92488748d5c28cca2f29f880ab747e8baf84a772..cf2ff58006ef61d1f759481b0ee7649edb34fb8d 100644
--- a/dumux/io/plotmateriallaw.hh
+++ b/dumux/io/plotmateriallaw.hh
@@ -24,6 +24,8 @@
#ifndef DUMUX_PLOT_FLUID_MATRIX_LAW_HH
#define DUMUX_PLOT_FLUID_MATRIX_LAW_HH
+#warning "This header is deprecated and will be removed after 3.3. Use new 2p material laws and plot tools from io/plotpckrsw.hh"
+
#include
#include
#include
diff --git a/dumux/io/plotpckrsw.hh b/dumux/io/plotpckrsw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..75bd4f2119229ef216b4787b8073fc17c9328c04
--- /dev/null
+++ b/dumux/io/plotpckrsw.hh
@@ -0,0 +1,223 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup InputOutput
+ * \brief Interface for plotting the two-phase fluid-matrix-interaction laws
+ */
+#ifndef DUMUX_IO_PLOT_PC_KR_SW_HH
+#define DUMUX_IO_PLOT_PC_KR_SW_HH
+
+#include
+#include
+#include
+#include
+
+namespace Dumux {
+
+namespace Detail {
+template
+Range evalFunctionForRange(const Function& f, const Range& range)
+{
+ Range result = range;
+ std::transform(range.begin(), range.end(), result.begin(), [&](auto x){ return f(x); });
+ return result;
+}
+} // end namespace Detail
+
+/*!
+ * \ingroup InputOutput
+ * \brief sample the pc-sw curve
+ */
+template
+auto samplePcSw(const PcKrSw& curve, const V& sw)
+{
+ return Detail::evalFunctionForRange([&](const auto s){ return curve.pc(s); }, sw);
+}
+
+/*!
+ * \ingroup InputOutput
+ * \brief sample the pc-sw curve derivative wrt sw
+ */
+template
+auto samplePcSwDerivative(const PcKrSw& curve, const V& sw)
+{
+ return Detail::evalFunctionForRange([&](const auto s){ return curve.dpc_dsw(s); }, sw);
+}
+
+/*!
+ * \ingroup InputOutput
+ * \brief sample the sw-pc curve derivative wrt pc
+ */
+template
+auto samplePcSwInverseDerivative(const PcKrSw& curve, const V& pc)
+{
+ return Detail::evalFunctionForRange([&](const auto p){ return curve.dsw_dpc(p); }, pc);
+}
+
+/*!
+ * \ingroup InputOutput
+ * \brief sample sw-pc curve but return the log10 of the capillary pressure
+ */
+template
+auto sampleLog10PcSw(const PcKrSw& curve, const V& sw)
+{
+ return Detail::evalFunctionForRange([&](const auto s){ using std::log10; return log10(curve.pc(s)); }, sw);
+}
+
+/*!
+ * \ingroup InputOutput
+ * \brief sample krw-sw and krn-sw curves
+ */
+template
+auto sampleRelPerms(const PcKrSw& curve, const V& sw)
+{
+ return std::make_tuple(
+ Detail::evalFunctionForRange([&](const auto s){ return curve.krw(s); }, sw),
+ Detail::evalFunctionForRange([&](const auto s){ return curve.krn(s); }, sw)
+ );
+}
+
+/*!
+ * \ingroup InputOutput
+ * \brief sample the derivatives of the krw-sw and krn-sw curves
+ */
+template
+auto sampleRelPermDerivatives(const PcKrSw& curve, const V& sw)
+{
+ return std::make_tuple(
+ Detail::evalFunctionForRange([&](const auto s){ return curve.dkrw_dsw(s); }, sw),
+ Detail::evalFunctionForRange([&](const auto s){ return curve.dkrn_dsw(s); }, sw)
+ );
+}
+
+// forward declaration
+template class GnuplotInterface;
+
+namespace Detail {
+template
+void addDataSetToGnuplot(GnuplotInterface& gnuplot,
+ const V& x, const V& y,
+ const std::string& curveName,
+ const std::string& curveOptions,
+ const std::string& xLabel,
+ const std::string& yLabel)
+{
+ gnuplot.setXlabel(xLabel);
+ gnuplot.setYlabel(yLabel);
+ gnuplot.addDataSetToPlot(x, y, curveName, curveOptions);
+}
+} // end namespace Detail
+
+/*!
+ * \ingroup InputOutput
+ * \brief Helper functions related to gnuplot
+ */
+namespace Gnuplot {
+
+/*!
+ * \ingroup InputOutput
+ * \brief Convenience function for adding material law quantities to gnuplot
+ */
+template
+void addPcSw(GnuplotInterface& gnuplot, const V& sw, const V& pc,
+ const std::string& curveName = "pc-sw",
+ const std::string& curveOptions = "w l",
+ const std::string& xLabel = "wetting phase saturation [-]",
+ const std::string& yLabel = "capillary pressure [Pa]")
+{
+ Detail::addDataSetToGnuplot(gnuplot, sw, pc, curveName, curveOptions, xLabel, yLabel);
+}
+
+/*!
+ * \ingroup InputOutput
+ * \brief Convenience function for adding material law quantities to gnuplot
+ */
+template
+void addPcSwDerivative(GnuplotInterface& gnuplot, const V& sw, const V& dpc_dsw,
+ const std::string& curveName = "dpc-dsw",
+ const std::string& curveOptions = "w l",
+ const std::string& xLabel = "wetting phase saturation [-]",
+ const std::string& yLabel = "derivative of capillary pressure [Pa]")
+{
+ Detail::addDataSetToGnuplot(gnuplot, sw, dpc_dsw, curveName, curveOptions, xLabel, yLabel);
+}
+
+/*!
+ * \ingroup InputOutput
+ * \brief Convenience function for adding material law quantities to gnuplot
+ */
+template
+void addPcSwInverseDerivative(GnuplotInterface& gnuplot, const V& sw, const V& dpc_dsw,
+ const std::string& curveName = "dsw-dpc",
+ const std::string& curveOptions = "w l",
+ const std::string& xLabel = "capillary pressure [Pa]",
+ const std::string& yLabel = "derivative of saturation [1/Pa]")
+{
+ Detail::addDataSetToGnuplot(gnuplot, sw, dpc_dsw, curveName, curveOptions, xLabel, yLabel);
+}
+
+/*!
+ * \ingroup InputOutput
+ * \brief Convenience function for adding material law quantities to gnuplot
+ */
+template
+void addLog10PcSw(GnuplotInterface& gnuplot, const V& sw, const V& log10pc,
+ const std::string& curveName = "log10-pc-sw",
+ const std::string& curveOptions = "w l",
+ const std::string& xLabel = "wetting phase saturation [-]",
+ const std::string& yLabel = "log10 of capillary pressure [Pa]")
+{
+ Detail::addDataSetToGnuplot(gnuplot, sw, log10pc, curveName, curveOptions, xLabel, yLabel);
+}
+
+/*!
+ * \ingroup InputOutput
+ * \brief Convenience function for adding material law quantities to gnuplot
+ */
+template
+void addRelPerms(GnuplotInterface& gnuplot, const V& sw, const V& krw, const V& krn,
+ const std::string& curveName = "relperm",
+ const std::string& curveOptions = "w l",
+ const std::string& xLabel = "wetting phase saturation [-]",
+ const std::string& yLabel = "relative permeability [-]")
+{
+ Detail::addDataSetToGnuplot(gnuplot, sw, krw, curveName + "_krw", curveOptions, xLabel, yLabel);
+ Detail::addDataSetToGnuplot(gnuplot, sw, krn, curveName + "_krn", curveOptions, xLabel, yLabel);
+}
+
+/*!
+ * \ingroup InputOutput
+ * \brief Convenience function for adding material law quantities to gnuplot
+ */
+template
+void addRelPermDerivatives(GnuplotInterface& gnuplot, const V& sw, const V& krw, const V& krn,
+ const std::string& curveName = "relperm_dsw",
+ const std::string& curveOptions = "w l",
+ const std::string& xLabel = "wetting phase saturation [-]",
+ const std::string& yLabel = "derivative of the relative permeability [-]")
+{
+ Detail::addDataSetToGnuplot(gnuplot, sw, krw, curveName + "_krw", curveOptions, xLabel, yLabel);
+ Detail::addDataSetToGnuplot(gnuplot, sw, krn, curveName + "_krn", curveOptions, xLabel, yLabel);
+}
+
+} // end namespace Gnuplot
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/material/constraintsolvers/immiscibleflash.hh b/dumux/material/constraintsolvers/immiscibleflash.hh
index b04e524534c0a833df16ecdc4a29513e8675ecd1..1c3a934be8061bf0177315b143c18395354ed4e1 100644
--- a/dumux/material/constraintsolvers/immiscibleflash.hh
+++ b/dumux/material/constraintsolvers/immiscibleflash.hh
@@ -29,6 +29,7 @@
#include
#include
+#include
namespace Dumux {
@@ -117,15 +118,15 @@ public:
* \param fluidState Thermodynamic state of the fluids
* \param paramCache Container for cache parameters
* \param globalMolarities
- * \param matParams The material law object
+ * \param material The material law object
*
* The phase's fugacities must already be set.
*/
template
- void solve(FluidState &fluidState,
- ParameterCache ¶mCache,
- const typename MaterialLaw::Params &matParams,
- const ComponentVector &globalMolarities)
+ void solve(FluidState& fluidState,
+ ParameterCache& paramCache,
+ const MaterialLaw& material,
+ const ComponentVector& globalMolarities)
{
/////////////////////////
// Check if all fluid phases are incompressible
@@ -164,12 +165,12 @@ public:
// right hand side
Vector b;
- completeFluidState_(fluidState, paramCache, matParams);
+ completeFluidState_(fluidState, paramCache, material);
const int nMax = 50; // <- maximum number of newton iterations
for (int nIdx = 0; nIdx < nMax; ++nIdx) {
// calculate Jacobian matrix and right hand side
- linearize_(J, b, fluidState, paramCache, matParams, globalMolarities);
+ linearize_(J, b, fluidState, paramCache, material, globalMolarities);
// Solve J*x = b
deltaX = 0;
@@ -207,7 +208,7 @@ public:
*/
// update the fluid quantities.
- Scalar relError = update_(fluidState, paramCache, matParams, deltaX);
+ Scalar relError = update_(fluidState, paramCache, material, deltaX);
if (relError < 1e-9)
return;
@@ -263,7 +264,7 @@ protected:
Vector &b,
FluidState &fluidState,
ParameterCache ¶mCache,
- const typename MaterialLaw::Params &matParams,
+ const MaterialLaw& material,
const ComponentVector &globalMolarities)
{
FluidState origFluidState(fluidState);
@@ -287,7 +288,7 @@ protected:
// deviate the mole fraction of the i-th component
Scalar x_i = getQuantity_(fluidState, pvIdx);
const Scalar eps = 1e-10/quantityWeight_(fluidState, pvIdx);
- setQuantity_(fluidState, paramCache, matParams, pvIdx, x_i + eps);
+ setQuantity_(fluidState, paramCache, material, pvIdx, x_i + eps);
// compute derivative of the defect
calculateDefect_(tmp, origFluidState, fluidState, globalMolarities);
@@ -325,7 +326,7 @@ protected:
template
Scalar update_(FluidState &fluidState,
ParameterCache ¶mCache,
- const typename MaterialLaw::Params &matParams,
+ const MaterialLaw& material,
const Vector &deltaX)
{
Scalar relError = 0;
@@ -351,7 +352,7 @@ protected:
setQuantityRaw_(fluidState, pvIdx, tmp - delta);
}
- completeFluidState_(fluidState, paramCache, matParams);
+ completeFluidState_(fluidState, paramCache, material);
return relError;
}
@@ -359,7 +360,7 @@ protected:
template
void completeFluidState_(FluidState &fluidState,
ParameterCache ¶mCache,
- const typename MaterialLaw::Params &matParams)
+ const MaterialLaw& material)
{
// calculate the saturation of the last phase as a function of
// the other saturations
@@ -384,8 +385,7 @@ protected:
// update the pressures using the material law (saturations
// and first pressure are already set because it is implicitly
// solved for.)
- ComponentVector pc;
- MaterialLaw::capillaryPressures(pc, matParams, fluidState, wettingPhaseIdx_);
+ const auto pc = material.capillaryPressures(fluidState, wettingPhaseIdx_);
for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
fluidState.setPressure(phaseIdx,
fluidState.pressure(0)
@@ -432,7 +432,7 @@ protected:
template
void setQuantity_(FluidState &fs,
ParameterCache ¶mCache,
- const typename MaterialLaw::Params &matParams,
+ const MaterialLaw& material,
int pvIdx,
Scalar value)
{
@@ -470,8 +470,7 @@ protected:
// update all fluid pressures using the capillary pressure
// law
- ComponentVector pc;
- MaterialLaw::capillaryPressures(pc, matParams, fs, wettingPhaseIdx_);
+ const auto pc = material.capillaryPressures(fs, wettingPhaseIdx_);
for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
fs.setPressure(phaseIdx,
fs.pressure(0)
diff --git a/dumux/material/constraintsolvers/ncpflash.hh b/dumux/material/constraintsolvers/ncpflash.hh
index 652ba0a9a96689f29c0388908ae882ef3001ea62..f5567a629ecfd488b847d6f40f4912bd6a54cde7 100644
--- a/dumux/material/constraintsolvers/ncpflash.hh
+++ b/dumux/material/constraintsolvers/ncpflash.hh
@@ -29,6 +29,7 @@
#include
#include
+#include
namespace Dumux {
@@ -143,14 +144,14 @@ public:
* \param fluidState Thermodynamic state of the fluids
* \param paramCache Container for cache parameters
* \param globalMolarities
- * \param matParams The material law object
+ * \param material The material law object
*
* The phase's fugacities must already be set.
*/
template
void solve(FluidState &fluidState,
ParameterCache ¶mCache,
- const typename MaterialLaw::Params &matParams,
+ const MaterialLaw& material,
const ComponentVector &globalMolarities)
{
/////////////////////////
@@ -165,9 +166,7 @@ public:
Vector b;
// make the fluid state consistent with the fluid system.
- completeFluidState_(fluidState,
- paramCache,
- matParams);
+ completeFluidState_(fluidState, paramCache, material);
/*
std::cout << "--------------------\n";
@@ -180,7 +179,7 @@ public:
const int nMax = 50; // <- maximum number of newton iterations
for (int nIdx = 0; nIdx < nMax; ++nIdx) {
// calculate Jacobian matrix and right hand side
- linearize_(J, b, fluidState, paramCache, matParams, globalMolarities);
+ linearize_(J, b, fluidState, paramCache, material, globalMolarities);
// Solve J*x = b
deltaX = 0;
@@ -236,7 +235,7 @@ public:
*/
// update the fluid quantities.
- Scalar relError = update_(fluidState, paramCache, matParams, deltaX);
+ Scalar relError = update_(fluidState, paramCache, material, deltaX);
if (relError < 1e-9)
return;
@@ -313,7 +312,7 @@ protected:
Vector &b,
FluidState &fluidState,
ParameterCache ¶mCache,
- const typename MaterialLaw::Params &matParams,
+ const MaterialLaw& material,
const ComponentVector &globalMolarities)
{
FluidState origFluidState(fluidState);
@@ -337,7 +336,7 @@ protected:
// deviate the mole fraction of the i-th component
Scalar x_i = getQuantity_(fluidState, pvIdx);
const Scalar eps = 1e-8/quantityWeight_(fluidState, pvIdx);
- setQuantity_(fluidState, paramCache, matParams, pvIdx, x_i + eps);
+ setQuantity_(fluidState, paramCache, material, pvIdx, x_i + eps);
// compute derivative of the defect
calculateDefect_(tmp, origFluidState, fluidState, globalMolarities);
@@ -418,7 +417,7 @@ protected:
template
Scalar update_(FluidState &fluidState,
ParameterCache ¶mCache,
- const typename MaterialLaw::Params &matParams,
+ const MaterialLaw& material,
const Vector &deltaX)
{
Scalar relError = 0;
@@ -486,7 +485,7 @@ protected:
}
}
- completeFluidState_(fluidState, paramCache, matParams);
+ completeFluidState_(fluidState, paramCache, material);
return relError;
}
@@ -494,7 +493,7 @@ protected:
template
void completeFluidState_(FluidState &fluidState,
ParameterCache ¶mCache,
- const typename MaterialLaw::Params &matParams)
+ const MaterialLaw& material)
{
// calculate the saturation of the last phase as a function of
// the other saturations
@@ -506,8 +505,7 @@ protected:
// update the pressures using the material law (saturations
// and first pressure are already set because it is implicitly
// solved for.)
- ComponentVector pc;
- MaterialLaw::capillaryPressures(pc, matParams, fluidState, wettingPhaseIdx_);
+ const auto pc = material.capillaryPressures(fluidState, wettingPhaseIdx_);
for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
fluidState.setPressure(phaseIdx,
fluidState.pressure(0)
@@ -569,7 +567,7 @@ protected:
template
void setQuantity_(FluidState &fs,
ParameterCache ¶mCache,
- const typename MaterialLaw::Params &matParams,
+ const MaterialLaw& material,
int pvIdx,
Scalar value)
{
@@ -608,8 +606,7 @@ protected:
// update all fluid pressures using the capillary pressure
// law
- ComponentVector pc;
- MaterialLaw::capillaryPressures(pc, matParams, fs, wettingPhaseIdx_);
+ const auto pc = material.capillaryPressures(fs, wettingPhaseIdx_);
for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
fs.setPressure(phaseIdx,
fs.pressure(0)
diff --git a/dumux/material/fluidmatrixinteractions/2p/CMakeLists.txt b/dumux/material/fluidmatrixinteractions/2p/CMakeLists.txt
index d2a4aacc9616dfdca638ab2e6b651cffb0814625..9c39155a73f5d0a8d07c8a007656bb68826e1b08 100644
--- a/dumux/material/fluidmatrixinteractions/2p/CMakeLists.txt
+++ b/dumux/material/fluidmatrixinteractions/2p/CMakeLists.txt
@@ -1,20 +1,27 @@
+add_subdirectory(interfacialarea)
add_subdirectory(thermalconductivity)
install(FILES
brookscorey.hh
brookscoreyparams.hh
+datasplinemateriallaw.hh
+efftoabsdefaultpolicy.hh
efftoabslaw.hh
efftoabslawparams.hh
heatpipelaw.hh
heatpipelawparams.hh
linearmaterial.hh
linearmaterialparams.hh
+materiallaw.hh
+noregularization.hh
regularizedbrookscorey.hh
regularizedbrookscoreyparams.hh
regularizedlinearmaterial.hh
regularizedlinearmaterialparams.hh
regularizedvangenuchten.hh
regularizedvangenuchtenparams.hh
+smoothedlinearlaw.hh
+splinemateriallaw.hh
thermalconductivityjohansen.hh
thermalconductivitysimplefluidlumping.hh
thermalconductivitysomerton.hh
diff --git a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh b/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
index 19cfc930bc2ca985bfe607a918765857c2a91145..a69963c302d342c3b284904148381461a2261abe 100644
--- a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
@@ -22,9 +22,10 @@
* \brief Implementation of the capillary pressure and
* relative permeability <-> saturation relations according to Brooks and Corey.
*/
-#ifndef DUMUX_BROOKS_COREY_HH
-#define DUMUX_BROOKS_COREY_HH
+#ifndef DUMUX_MATERIAL_FLUIDMATRIX_BROOKS_COREY_HH
+#define DUMUX_MATERIAL_FLUIDMATRIX_BROOKS_COREY_HH
+// remove from here after release 3.3 /////////////
#include "brookscoreyparams.hh"
#include
@@ -45,7 +46,7 @@ namespace Dumux {
*\see BrooksCoreyParams
*/
template >
-class BrooksCorey
+class [[deprecated("Use new material laws and FluidMatrix::BrooksCorey instead!")]] BrooksCorey
{
public:
using Params = ParamsT;
@@ -272,5 +273,507 @@ public:
};
} // end namespace Dumux
+// remove until here after release 3.3 /////////////
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+namespace Dumux::FluidMatrix {
+
+/*!
+ * \ingroup Fluidmatrixinteractions
+ *
+ * \brief Implementation of the Brooks-Corey capillary pressure <->
+ * saturation relation. This class bundles the "raw" curves
+ * as static members and doesn't concern itself converting
+ * absolute to effective saturations and vice versa.
+ *
+ * For general info: EffToAbsLaw
+ *
+ *\see BrooksCoreyParams
+ */
+class BrooksCorey
+{
+public:
+ /*!
+ * \brief The parameter type
+ * \tparam Scalar The scalar type
+ * \note The Brooks Corey laws are parameterized with two parameters: \f$\mathrm{p_{ce}, \lambda}\f$,
+ * the capillary entry pressure in \f$\mathrm{[Pa]}\f$] and a dimensionless shape parameter, respectively.
+ */
+ template
+ struct Params
+ {
+ Params(Scalar pcEntry, Scalar lambda)
+ : pcEntry_(pcEntry), lambda_(lambda)
+ {}
+
+ Scalar pcEntry() const{ return pcEntry_; }
+ void setPcEntry(Scalar pcEntry){ pcEntry_ = pcEntry; }
+
+ Scalar lambda() const { return lambda_; }
+ void setLambda(Scalar lambda) { lambda_ = lambda; }
+
+ bool operator== (const Params& p) const
+ {
+ return Dune::FloatCmp::eq(pcEntry(), p.pcEntry(), 1e-6)
+ && Dune::FloatCmp::eq(lambda(), p.lambda(), 1e-6);
+ }
+
+ private:
+ Scalar pcEntry_, lambda_;
+ };
+
+ /*!
+ * \brief Construct from a subgroup from the global parameter tree
+ * \note This will give you nice error messages if a mandatory parameter is missing
+ */
+ template
+ static Params makeParams(const std::string& paramGroup)
+ {
+ const auto pcEntry = getParamFromGroup(paramGroup, "BrooksCoreyPcEntry");
+ const auto lambda = getParamFromGroup(paramGroup, "BrooksCoreyLambda");
+ return {pcEntry, lambda};
+ }
+
+ /*!
+ * \brief The capillary pressure-saturation curve according to Brooks & Corey.
+ *
+ * The Brooks-Corey empirical capillary pressure <-> saturation
+ * function is given by
+ *
+ * \f$\mathrm{ p_c = p_{ce}\overline{S}_w^{-1/\lambda}
+ * }\f$
+ *
+ * \param swe Effective saturation of the wetting phase \f$\mathrm{[\overline{S}_w]}\f$
+ * \param params A container object that is populated with the appropriate coefficients for the respective law.
+ * Therefore, in the (problem specific) spatialParameters first, the material law is chosen,
+ and then the params container is constructed accordingly. Afterwards the values are set there, too.
+ * \return Capillary pressure calculated by Brooks & Corey constitutive relation.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
+ */
+ template
+ static Scalar pc(Scalar swe, const Params& params)
+ {
+ using std::pow;
+ using std::clamp;
+
+ swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+
+ return params.pcEntry()*pow(swe, -1.0/params.lambda());
+ }
+
+ /*!
+ * \brief The saturation-capillary pressure curve according to Brooks & Corey.
+ *
+ * This is the inverse of the capillary pressure-saturation curve:
+ * \f$\mathrm{ \overline{S}_w = (\frac{p_c}{p_{ce}})^{-\lambda}}\f$
+ *
+ * \param pc Capillary pressure \f$\mathrm{[p_c]}\f$ in \f$\mathrm{[Pa]}\f$.
+ * \param params A container object that is populated with the appropriate coefficients for the respective law.
+ * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
+ * is constructed accordingly. Afterwards the values are set there, too.
+ * \return Effective wetting phase saturation calculated as inverse of BrooksCorey constitutive relation.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
+ */
+ template
+ static Scalar swe(Scalar pc, const Params& params)
+ {
+ using std::pow;
+ using std::max;
+
+ pc = max(pc, 0.0); // the equation below is undefined for negative pcs
+
+ return pow(pc/params.pcEntry(), -params.lambda());
+ }
+
+ /*!
+ * \brief The capillary pressure at Swe = 1.0 also called end point capillary pressure
+ *
+ * \param params A container object that is populated with the appropriate coefficients for the respective law.
+ * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
+ * is constructed accordingly. Afterwards the values are set there, too.
+ */
+ template
+ static Scalar endPointPc(const Params& params)
+ { return params.pcEntry(); }
+
+ /*!
+ * \brief The partial derivative of the capillary
+ * pressure w.r.t. the effective saturation according to Brooks & Corey.
+ *
+ * This is equivalent to
+ * \f$\mathrm{\frac{\partial p_c}{\partial \overline{S}_w} =
+ * -\frac{p_{ce}}{\lambda} \overline{S}_w^{-1/\lambda - 1}
+ * }\f$
+ *
+ * \param swe Effective saturation of the wetting phase \f$\mathrm{[\overline{S}_w]}\f$
+ * \param params A container object that is populated with the appropriate coefficients for the respective law.
+ * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
+ * is constructed accordingly. Afterwards the values are set there, too.
+ * \return Partial derivative of \f$\mathrm{[p_c]}\f$ w.r.t. effective saturation according to Brooks & Corey.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
+ */
+ template
+ static Scalar dpc_dswe(Scalar swe, const Params& params)
+ {
+ using std::pow;
+ using std::clamp;
+
+ swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+
+ return -params.pcEntry()/params.lambda() * pow(swe, -1.0/params.lambda() - 1.0);
+ }
+
+ /*!
+ * \brief The partial derivative of the effective
+ * saturation w.r.t. the capillary pressure according to Brooks & Corey.
+ *
+ * \param pc Capillary pressure \f$\mathrm{[p_c]}\f$ in \f$\mathrm{[Pa]}\f$.
+ * \param params A container object that is populated with the appropriate coefficients for the respective law.
+ * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
+ * is constructed accordingly. Afterwards the values are set there, too.
+ * \return Partial derivative of effective saturation w.r.t. \f$\mathrm{[p_c]}\f$ according to Brooks & Corey.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
+ */
+ template
+ static Scalar dswe_dpc(Scalar pc, const Params& params)
+ {
+ using std::pow;
+ using std::max;
+
+ pc = max(pc, 0.0); // the equation below is undefined for negative pcs
+
+ return -params.lambda()/params.pcEntry() * pow(pc/params.pcEntry(), - params.lambda() - 1.0);
+ }
+
+ /*!
+ * \brief The relative permeability for the wetting phase of
+ * the medium implied by the Brooks-Corey
+ * parameterization.
+ *
+ * \param swe The mobile saturation of the wetting phase.
+ * \param params A container object that is populated with the appropriate coefficients for the respective law.
+ * Therefore, in the (problem specific) spatialParameters first, the material law is chosen,
+ * and then the params container is constructed accordingly. Afterwards the values are set there, too.
+ * \return Relative permeability of the wetting phase calculated as implied by Brooks & Corey.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
+ */
+ template
+ static Scalar krw(Scalar swe, const Params& params)
+ {
+ using std::pow;
+ using std::clamp;
+
+ swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+
+ return pow(swe, 2.0/params.lambda() + 3.0);
+ }
+
+ /*!
+ * \brief The derivative of the relative permeability for the
+ * wetting phase with regard to the wetting saturation of the
+ * medium implied by the Brooks-Corey parameterization.
+ *
+ * \param swe The mobile saturation of the wetting phase.
+ * \param params A container object that is populated with the appropriate coefficients for the respective law.
+ * Therefore, in the (problem specific) spatialParameters first, the material law is chosen,
+ * and then the params container is constructed accordingly. Afterwards the values are set there, too.
+ * \return Derivative of the relative permeability of the wetting phase w.r.t. effective wetting phase
+ * saturation calculated as implied by Brooks & Corey.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
+ */
+ template
+ static Scalar dkrw_dswe(Scalar swe, const Params& params)
+ {
+ using std::pow;
+ using std::clamp;
+
+ swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+
+ return (2.0/params.lambda() + 3.0)*pow(swe, 2.0/params.lambda() + 2.0);
+ }
+
+ /*!
+ * \brief The relative permeability for the non-wetting phase of
+ * the medium as implied by the Brooks-Corey
+ * parameterization.
+ *
+ * \param swe The mobile saturation of the wetting phase.
+ * \param params A container object that is populated with the appropriate coefficients for the respective law.
+ * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
+ * is constructed accordingly. Afterwards the values are set there, too.
+ * \return Relative permeability of the non-wetting phase calculated as implied by Brooks & Corey.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
+ */
+ template
+ static Scalar krn(Scalar swe, const Params& params)
+ {
+ using std::pow;
+ using std::clamp;
+
+ swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+
+ const Scalar exponent = 2.0/params.lambda() + 1.0;
+ const Scalar sne = 1.0 - swe;
+ return sne*sne*(1.0 - pow(swe, exponent));
+ }
+
+ /*!
+ * \brief The derivative of the relative permeability for the
+ * non-wetting phase in regard to the wetting saturation of
+ * the medium as implied by the Brooks-Corey
+ * parameterization.
+ *
+ * \param swe The mobile saturation of the wetting phase.
+ * \param params A container object that is populated with the appropriate coefficients for the respective law.
+ * Therefore, in the (problem specific) spatialParameters first, the material law is chosen,
+ * and then the params container is constructed accordingly. Afterwards the values are set there, too.
+ * \return Derivative of the relative permeability of the non-wetting phase w.r.t. effective wetting phase
+ * saturation calculated as implied by Brooks & Corey.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
+ */
+ template
+ static Scalar dkrn_dswe(Scalar swe, const Params& params)
+ {
+ using std::pow;
+ using std::clamp;
+
+ swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+
+ const auto lambdaInv = 1.0/params.lambda();
+ const auto swePow = pow(swe, 2*lambdaInv);
+ return 2.0*(swe - 1.0)*(1.0 + (0.5 + lambdaInv)*swePow - (1.5 + lambdaInv)*swePow*swe);
+ }
+};
+
+/*!
+ * \ingroup Fluidmatrixinteractions
+ * \brief A regularization for the BrooksCorey material law
+ * \note Regularization can be turned of by setting the threshold parameters
+ * out of range (runtime) or by replacing
+ * this class by NoRegularization (compile time).
+ */
+template
+class BrooksCoreyRegularization
+{
+public:
+ //! Regularization parameters
+ template
+ struct Params
+ {
+ Params(S pcLowSwe = 0.01) : pcLowSwe_(pcLowSwe) {}
+
+ S pcLowSwe() const { return pcLowSwe_; }
+ void setPcLowSwe(S pcLowSwe) { pcLowSwe_ = pcLowSwe; }
+
+ private:
+ S pcLowSwe_ = 0.01;
+ };
+
+ //! Initialize the spline
+ template
+ void init(const MaterialLaw* m, const std::string& paramGroup)
+ {
+ pcLowSwe_ = getParamFromGroup(paramGroup, "BrooksCoreyPcLowSweThreshold", 0.01);
+ entryPressure_ = getParamFromGroup(paramGroup, "BrooksCoreyPcEntry");
+
+ initPcParameters_(m, pcLowSwe_);
+ }
+
+ template
+ void init(const MaterialLaw* m, const BaseParams& bp, const EffToAbsParams& etap, const Params& p)
+ {
+ pcLowSwe_ = p.pcLowSwe();
+ entryPressure_ = bp.pcEntry();
+
+ initPcParameters_(m, pcLowSwe_);
+ }
+
+ /*!
+ * \brief Equality comparison with another instance
+ */
+ bool operator== (const BrooksCoreyRegularization& o) const
+ {
+ return Dune::FloatCmp::eq(pcLowSwe_, o.pcLowSwe_, 1e-6)
+ && Dune::FloatCmp::eq(entryPressure_, o.entryPressure_, 1e-6);
+ }
+
+ /*!
+ * \brief The regularized capillary pressure-saturation curve
+ * regularized part:
+ * - low saturation: extend the \f$\mathrm{p_c(S_w)}\f$ curve with the slope at the regularization point (i.e. no kink).
+ * - high saturation: connect the high regularization point with \f$\mathrm{\overline{S}_w =1}\f$
+ * with a spline and continue linearly for \f$\mathrm{\overline{S}_w > 1}\f$
+ */
+ OptionalScalar pc(const Scalar swe) const
+ {
+ // make sure that the capillary pressure observes a derivative
+ // != 0 for 'illegal' saturations. This is favourable for the
+ // newton solver (if the derivative is calculated numerically)
+ // in order to get the saturation moving to the right
+ // direction if it temporarily is in an 'illegal' range.
+ if (swe <= pcLowSwe_)
+ return pcLowSwePcValue_ + pcDerivativeLowSw_*(swe - pcLowSwe_);
+
+ else if (swe > 1.0)
+ return pcDerivativeHighSwEnd_*(swe - 1.0) + entryPressure_;
+
+ else
+ return {}; // no regularization
+ }
+
+ /*!
+ * \brief The regularized partial derivative of the capillary pressure w.r.t. the saturation
+ */
+ OptionalScalar dpc_dswe(const Scalar swe) const
+ {
+ if (swe < pcLowSwe_)
+ return pcDerivativeLowSw_;
+
+ else if (swe > 1.0)
+ return pcDerivativeHighSwEnd_;
+
+ else
+ return {}; // no regularization
+ }
+
+ /*!
+ * \brief The regularized saturation-capillary pressure curve
+ */
+ OptionalScalar swe(const Scalar pc) const
+ {
+ if (pc <= entryPressure_)
+ return 1.0 + (pc - entryPressure_)/pcDerivativeHighSwEnd_;
+
+ else if (pc >= pcLowSwePcValue_)
+ return (pc - pcLowSwePcValue_)/pcDerivativeLowSw_ + pcLowSwe_;
+
+ else
+ return {}; // no regularization
+ }
+
+ /*!
+ * \brief The regularized partial derivative of the saturation to the capillary pressure
+ */
+ OptionalScalar dswe_dpc(const Scalar pc) const
+ {
+ if (pc <= entryPressure_)
+ return 1.0/pcDerivativeHighSwEnd_;
+
+ else if (pc >= pcLowSwePcValue_)
+ return 1.0/pcDerivativeLowSw_;
+
+ else
+ return {}; // no regularization
+ }
+
+ /*!
+ * \brief The regularized relative permeability for the wetting phase
+ */
+ OptionalScalar krw(const Scalar swe) const
+ {
+ if (swe <= 0.0)
+ return 0.0;
+ else if (swe >= 1.0)
+ return 1.0;
+ else
+ return {}; // no regularization
+ }
+
+ /*!
+ * \brief The regularized derivative of the relative permeability for the wetting phase w.r.t. saturation
+ */
+ OptionalScalar dkrw_dswe(const Scalar swe) const
+ {
+ if (swe < 0.0)
+ return 0.0;
+ else if (swe >= 1.0)
+ return 0.0;
+ else
+ return {}; // no regularization
+ }
+
+ /*!
+ * \brief The regularized relative permeability for the non-wetting phase
+ */
+ OptionalScalar krn(const Scalar swe) const
+ {
+ if (swe <= 0.0)
+ return 1.0;
+ else if (swe >= 1.0)
+ return 0.0;
+ else
+ return {}; // no regularization
+ }
+
+ /*!
+ * \brief The regularized derivative of the relative permeability for the non-wetting phase w.r.t. saturation
+ */
+ OptionalScalar dkrn_dswe(const Scalar swe) const
+ {
+ if (swe <= 0.0)
+ return 0.0;
+ else if (swe >= 1.0)
+ return 0.0;
+ else
+ return {}; // no regularization
+ }
+
+private:
+ template
+ void initPcParameters_(const MaterialLaw* m, const Scalar lowSwe)
+ {
+ const auto lowSw = MaterialLaw::EffToAbs::sweToSw(lowSwe, m->effToAbsParams());
+ const auto highSw = MaterialLaw::EffToAbs::sweToSw(1.0, m->effToAbsParams());
+ const auto dsw_dswe = MaterialLaw::EffToAbs::dsw_dswe(m->effToAbsParams());
+ pcDerivativeLowSw_ = m->template dpc_dsw(lowSw)*dsw_dswe;
+ pcDerivativeHighSwEnd_ = m->template dpc_dsw(highSw)*dsw_dswe;
+ pcLowSwePcValue_ = m->template pc(lowSw);
+ }
+
+ Scalar pcLowSwe_;
+ Scalar pcLowSwePcValue_;
+ Scalar entryPressure_;
+ Scalar pcDerivativeLowSw_;
+ Scalar pcDerivativeHighSwEnd_;
+};
+
+/*!
+ * \ingroup Fluidmatrixinteractions
+ * \brief A default configuration for using the VanGenuchten material law
+ */
+template
+using BrooksCoreyDefault = TwoPMaterialLaw, TwoPEffToAbsDefaultPolicy>;
+
+/*!
+ * \ingroup Fluidmatrixinteractions
+ * \brief A default configuration without regularization for using the VanGenuchten material law
+ */
+template
+using BrooksCoreyNoReg = TwoPMaterialLaw;
+
+} // end namespace Dumux::FluidMatrix
#endif
diff --git a/dumux/material/fluidmatrixinteractions/2p/brookscoreyparams.hh b/dumux/material/fluidmatrixinteractions/2p/brookscoreyparams.hh
index ad0f28a6e77020cdb47d0b02ac2393c8e8ce01fc..a8100a145857da7142ed213210117d3a627fcf43 100644
--- a/dumux/material/fluidmatrixinteractions/2p/brookscoreyparams.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/brookscoreyparams.hh
@@ -36,6 +36,7 @@ namespace Dumux {
* for the Brooks Corey constitutive relations.
* \see BrooksCorey
*/
+// "Use new material laws! Removal after 3.3")
template
class BrooksCoreyParams
{
diff --git a/dumux/material/fluidmatrixinteractions/2p/datasplinemateriallaw.hh b/dumux/material/fluidmatrixinteractions/2p/datasplinemateriallaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..990e43fd09648b7e20d08365fc8126d05a079637
--- /dev/null
+++ b/dumux/material/fluidmatrixinteractions/2p/datasplinemateriallaw.hh
@@ -0,0 +1,201 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Fluidmatrixinteractions
+ * \brief Pc- and Kr-sw curves based on monotone splines through given data points
+ */
+#ifndef DUMUX_MATERIAL_FLUIDMATRIX_DATA_SPLINE_MATERIAL_HH
+#define DUMUX_MATERIAL_FLUIDMATRIX_DATA_SPLINE_MATERIAL_HH
+
+#include
+#include
+#include
+#include
+#include
+
+namespace Dumux::FluidMatrix {
+
+/*!
+ * \ingroup Fluidmatrixinteractions
+ * \brief Pc- and Kr-sw curves based on monotone splines through given data points
+ * \tparam S the type for scalar numbers
+ * \tparam approximatePcSwInverse if this is set true, the
+ * spline approximates sw(pc) and evaluating pc(sw) needs spline inversion.
+ * if this is false, the spline approximates pc(sw) and evaluating
+ * sw(pc) needs spline inversion. Spline inversion is rather expensive
+ * since it has to be done numerically.
+ */
+template
+class DataSplineTwoPMaterialLaw
+: public Adapter, PcKrSw>
+{
+public:
+
+ using Scalar = S;
+
+ /*!
+ * \brief Return the number of fluid phases
+ */
+ static constexpr int numFluidPhases()
+ { return 2; }
+
+ static constexpr bool isRegularized()
+ { return false; }
+
+ /*!
+ * \brief Deleted default constructor (so we are never in an undefined state)
+ * \note store owning pointers to laws instead if you need default-constructible objects
+ */
+ DataSplineTwoPMaterialLaw() = delete;
+
+ /*!
+ * \brief Construct from a subgroup from the global parameter tree
+ */
+ explicit DataSplineTwoPMaterialLaw(const std::string& paramGroup)
+ {
+ using V = std::vector;
+ const std::string swPcGroup = paramGroup.empty() ? "Pc" : paramGroup + ".Pc";
+ const auto swPc = getParamFromGroup(swPcGroup, "SwData");
+ const std::string krwPcGroup = paramGroup.empty() ? "Krw" : paramGroup + ".Krw";
+ const auto swKrw = getParamFromGroup(krwPcGroup, "SwData");
+ const std::string krnPcGroup = paramGroup.empty() ? "Krn" : paramGroup + ".Krn";
+ const auto swKrn = getParamFromGroup(krnPcGroup, "SwData");
+
+ const auto pc = getParamFromGroup(paramGroup, "PcData");
+ const auto krw = getParamFromGroup(paramGroup, "KrwData");
+ const auto krn = getParamFromGroup(paramGroup, "KrnData");
+
+ updateData_(swPc, pc, swKrw, krw, swKrn, krn);
+ }
+
+ /*!
+ * \brief Construct from user data
+ */
+ DataSplineTwoPMaterialLaw(const std::vector& swPc,
+ const std::vector& pc,
+ const std::vector& swKrw,
+ const std::vector& krw,
+ const std::vector& swKrn,
+ const std::vector& krn)
+ {
+ updateData_(swPc, pc, swKrw, krw, swKrn, krn);
+ }
+
+ /*!
+ * \brief The capillary pressure-saturation curve
+ */
+ Scalar pc(const Scalar sw) const
+ {
+ if constexpr (approximatePcSwInverse)
+ return pcSpline_.evalInverse(sw);
+ else
+ return pcSpline_.eval(sw);
+ }
+
+ /*!
+ * \brief The partial derivative of the capillary pressure w.r.t. the saturation
+ */
+ Scalar dpc_dsw(const Scalar sw) const
+ {
+ if constexpr (approximatePcSwInverse)
+ return 1.0/pcSpline_.evalDerivative(pcSpline_.evalInverse(sw));
+ else
+ return pcSpline_.evalDerivative(sw);
+ }
+
+ /*!
+ * \brief The saturation-capillary pressure curve
+ */
+ Scalar sw(const Scalar pc) const
+ {
+ if constexpr (approximatePcSwInverse)
+ return pcSpline_.eval(pc);
+ else
+ return pcSpline_.evalInverse(pc);
+ }
+
+ /*!
+ * \brief The partial derivative of the saturation to the capillary pressure
+ */
+ Scalar dsw_dpc(const Scalar pc) const
+ {
+ if constexpr (approximatePcSwInverse)
+ return pcSpline_.evalDerivative(pc);
+ else
+ return 1.0/pcSpline_.evalDerivative(pcSpline_.evalInverse(pc));
+ }
+
+ /*!
+ * \brief The relative permeability for the wetting phase
+ */
+ Scalar krw(const Scalar sw) const
+ {
+ return krwSpline_.eval(sw);
+ }
+
+ /*!
+ * \brief The derivative of the relative permeability for the wetting phase w.r.t. saturation
+ */
+ Scalar dkrw_dsw(const Scalar sw) const
+ {
+ return krwSpline_.evalDerivative(sw);
+ }
+
+ /*!
+ * \brief The relative permeability for the non-wetting phase
+ */
+ Scalar krn(const Scalar sw) const
+ {
+ return krnSpline_.eval(sw);
+ }
+
+ /*!
+ * \brief The derivative of the relative permeability for the non-wetting phase w.r.t. saturation
+ */
+ Scalar dkrn_dsw(const Scalar sw) const
+ {
+ return krnSpline_.evalDerivative(sw);
+ }
+
+private:
+ void updateData_(const std::vector& swPc,
+ const std::vector& pc,
+ const std::vector& swKrw,
+ const std::vector& krw,
+ const std::vector& swKrn,
+ const std::vector& krn)
+ {
+ if constexpr (approximatePcSwInverse)
+ pcSpline_.updatePoints(pc, swPc);
+ else
+ pcSpline_.updatePoints(swPc, pc);
+
+ krwSpline_.updatePoints(swKrw, krw);
+ krnSpline_.updatePoints(swKrn, krn);
+ }
+
+ MonotoneCubicSpline pcSpline_;
+ MonotoneCubicSpline krwSpline_;
+ MonotoneCubicSpline krnSpline_;
+};
+
+} // end namespace Dumux::FluidMatrix
+
+#endif
diff --git a/dumux/material/fluidmatrixinteractions/2p/efftoabsdefaultpolicy.hh b/dumux/material/fluidmatrixinteractions/2p/efftoabsdefaultpolicy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c0af48e267bc868a5014d97509e28c8553f8006a
--- /dev/null
+++ b/dumux/material/fluidmatrixinteractions/2p/efftoabsdefaultpolicy.hh
@@ -0,0 +1,166 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Fluidmatrixinteractions
+ * \brief This is a policy for 2p material laws how to convert absolute to relative
+ * saturations and vice versa.
+ *
+ */
+#ifndef DUMUX_MATERIAL_FLUIDMATRIX_TWOP_EFF_TO_ABS_DEFAULT_POLICY_HH
+#define DUMUX_MATERIAL_FLUIDMATRIX_TWOP_EFF_TO_ABS_DEFAULT_POLICY_HH
+
+namespace Dumux::FluidMatrix {
+
+/*!
+ * \ingroup Fluidmatrixinteractions
+ *
+ * \brief This is a policy for 2p material laws how to convert absolute to relative
+ * saturations and vice versa.
+ *
+ * Material laws (like VanGenuchten or BrooksCorey) are defined for effective saturations.
+ * The numeric calculations however are performed with absolute saturations. The policy class converts
+ * the saturations. This allows for changing the calculation of the effective
+ * saturations easily, as this is subject of discussion / may be problem specific.
+ */
+class TwoPEffToAbsDefaultPolicy
+{
+public:
+ /*!
+ * \brief The parameter type
+ * \tparam Scalar The scalar type
+ * \note The efftoabs policy need two parameters: \f$\mathrm{S_{w,r}}, \mathrm{S_{n,r}}\f$.
+ * For the respective formulas check out the description of the free function.
+ */
+ template
+ struct Params
+ {
+ Params(const Scalar swr = 0.0, const Scalar snr = 0.0)
+ : swr_(swr), snr_(snr)
+ {}
+
+ /*!
+ * \brief Return the residual wetting saturation.
+ */
+ Scalar swr() const
+ { return swr_; }
+
+ /*!
+ * \brief Set the residual wetting saturation.
+ */
+ void setSwr(Scalar v)
+ { swr_ = v; }
+
+ /*!
+ * \brief Return the residual non-wetting saturation.
+ */
+ Scalar snr() const
+ { return snr_; }
+
+ /*!
+ * \brief Set the residual non-wetting saturation.
+ */
+ void setSnr(Scalar v)
+ { snr_ = v; }
+
+ bool operator== (const Params& p) const
+ {
+ return Dune::FloatCmp::eq(swr(), p.swr(), 1e-6)
+ && Dune::FloatCmp::eq(snr(), p.snr(), 1e-6);
+ }
+ private:
+ Scalar swr_;
+ Scalar snr_;
+ };
+
+ /*!
+ * \brief Construct from a subgroup from the global parameter tree
+ * \note This will give you nice error messages if a mandatory parameter is missing
+ */
+ template
+ static Params makeParams(const std::string& paramGroup)
+ {
+ Params params;
+ params.setSwr(getParamFromGroup(paramGroup, "Swr", 0.0));
+ params.setSnr(getParamFromGroup(paramGroup, "Snr", 0.0));
+ return params;
+ }
+
+ /*!
+ * \brief Convert an absolute wetting saturation to an effective one.
+ *
+ * \param sw Absolute saturation of the wetting phase \f$\mathrm{[{S}_w]}\f$.
+ * \param params A container object that is populated with the appropriate coefficients for the respective law.
+ * Therefore, in the (problem specific) spatialParameters first, the material law is chosen,
+ * and then the params container is constructed accordingly. Afterwards the values are set there, too.
+ * \return Effective saturation of the wetting phase.
+ */
+ template
+ static Scalar swToSwe(const Scalar sw, const Params& params)
+ {
+ return (sw - params.swr())/(1.0 - params.swr() - params.snr());
+ }
+
+ /*!
+ * \brief Convert an effective wetting saturation to an absolute one.
+ *
+ * \param swe Effective saturation of the non-wetting phase \f$\mathrm{[\overline{S}_n]}\f$.
+ * \param params A container object that is populated with the appropriate coefficients for the respective law.
+ * Therefore, in the (problem specific) spatialParameters first, the material law is chosen,
+ * and then the params container is constructed accordingly. Afterwards the values are set there, too.
+ * \return Absolute saturation of the non-wetting phase.
+ */
+ template
+ static Scalar sweToSw(const Scalar swe, const Params& params)
+ {
+ return swe*(1.0 - params.swr() - params.snr()) + params.swr();
+ }
+
+ /*!
+ * \brief Derivative of the effective saturation w.r.t. the absolute saturation.
+ *
+ * \param params A container object that is populated with the appropriate coefficients for the respective law.
+ * Therefore, in the (problem specific) spatialParameters first, the material law is chosen,
+ * and then the params container is constructed accordingly. Afterwards the values are set there, too.
+ * \return Derivative of the effective saturation w.r.t. the absolute saturation.
+ */
+ template
+ static Scalar dswe_dsw(const Params& params)
+ {
+ return 1.0/(1.0 - params.swr() - params.snr());
+ }
+
+ /*!
+ * \brief Derivative of the absolute saturation w.r.t. the effective saturation.
+ *
+ * \param params A container object that is populated with the appropriate coefficients for the respective law.
+ * Therefore, in the (problem specific) spatialParameters first, the material law is chosen,
+ * and then the params container is constructed accordingly. Afterwards the values are set there, too.
+ * \return Derivative of the absolute saturation w.r.t. the effective saturation.
+ */
+ template
+ static Scalar dsw_dswe(const Params& params)
+ {
+ return 1.0 - params.swr() - params.snr();
+ }
+};
+
+} // end namespace Dumux::FluidMatrix
+
+#endif
diff --git a/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh b/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh
index 6afd11a0fa45033379384250574a1790fc1d956d..41a49ea7a172f76ea9dfd3d8fa85127e7022c796 100644
--- a/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh
@@ -26,6 +26,8 @@
#ifndef DUMUX_EFF_TO_ABS_LAW_HH
#define DUMUX_EFF_TO_ABS_LAW_HH
+#warning "This header is deprecated. Removal after 3.3. Use new material laws."
+
#include "efftoabslawparams.hh"
namespace Dumux {
diff --git a/dumux/material/fluidmatrixinteractions/2p/efftoabslawparams.hh b/dumux/material/fluidmatrixinteractions/2p/efftoabslawparams.hh
index bd83efb5462a6673701deccb68e465526600f0fc..a9f03997e288a0420dca20ab5752af2f5e726134 100644
--- a/dumux/material/fluidmatrixinteractions/2p/efftoabslawparams.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/efftoabslawparams.hh
@@ -26,6 +26,8 @@
#ifndef DUMUX_EFF_TO_ABS_LAW_PARAMS_HH
#define DUMUX_EFF_TO_ABS_LAW_PARAMS_HH
+#warning "This header is deprecated. Removal after 3.3. Use new material laws."
+
#include
namespace Dumux {
@@ -37,7 +39,7 @@ namespace Dumux {
* saturations.
*/
template
-class EffToAbsLawParams : public EffLawParamsT
+class [[deprecated("Use new material laws! Removal after 3.3")]] EffToAbsLawParams : public EffLawParamsT
{
using EffLawParams = EffLawParamsT;
public:
diff --git a/dumux/material/fluidmatrixinteractions/2p/heatpipelaw.hh b/dumux/material/fluidmatrixinteractions/2p/heatpipelaw.hh
index b21b2ebe5f4ef1c2f57e2749ef4d1a3e3d3292b2..43a96a6170754b17cc21d909f5c0075b177fa2fb 100644
--- a/dumux/material/fluidmatrixinteractions/2p/heatpipelaw.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/heatpipelaw.hh
@@ -22,9 +22,10 @@
* \brief Implementation of the capillary pressure <-> saturation relation
* for the heatpipe problem.
*/
-#ifndef HEATPIPELAW_HH
-#define HEATPIPELAW_HH
+#ifndef DUMUX_MATERIAL_FLUIDMATRIX_TWOP_HEATPIPELAW_HH
+#define DUMUX_MATERIAL_FLUIDMATRIX_TWOP_HEATPIPELAW_HH
+// remove from here after release 3.3 /////////////
#include "heatpipelawparams.hh"
#include
@@ -45,7 +46,7 @@ namespace Dumux {
* converting absolute to effective saturations and vince versa.
*/
template >
-class HeatPipeLaw
+class [[deprecated("Use new material laws and FluidMatrix::HeatPipeLaw instead!")]] HeatPipeLaw
{
public:
using Params = ParamsT;
@@ -166,5 +167,215 @@ private:
};
} // end namespace Dumux
+// remove until here after release 3.3 /////////////
+
+#include
+#include
+
+#include
+#include
+#include
+
+namespace Dumux::FluidMatrix {
+
+/*!
+ * \ingroup Fluidmatrixinteractions
+ * \brief Implementation of the capillary pressure <-> saturation
+ * relation for the heatpipe problem.
+ */
+template
+class HeatPipeLaw : Adapter, PcKrSw>
+{
+
+public:
+ using Scalar = ScalarType;
+ using EffToAbsParams = typename EffToAbsPolicy::template Params;
+ using EffToAbs = EffToAbsPolicy;
+
+ /*!
+ * \brief Return the number of fluid phases
+ */
+ static constexpr int numFluidPhases()
+ { return 2; }
+
+ /*!
+ * \brief Return whether this law is regularized
+ */
+ static constexpr bool isRegularized()
+ { return false; }
+
+ /*!
+ * \brief The parameter type
+ * \tparam Scalar The scalar type
+ */
+ struct Params
+ {
+ Params(Scalar gamma, Scalar p0)
+ : gamma_(gamma), p0_(p0)
+ {}
+
+ Scalar gamma() const { return gamma_; }
+ void setGamma(Scalar gamma) { gamma_ = gamma; }
+
+ Scalar p0() const { return p0_; }
+ void setP0(Scalar p0) { p0_ = p0; }
+
+ bool operator== (const Params& p) const
+ {
+ return Dune::FloatCmp::eq(gamma(), p.gamma(), 1e-6)
+ && Dune::FloatCmp::eq(p0(), p.p0(), 1e-6);
+ }
+
+ private:
+ Scalar gamma_, p0_;
+ };
+
+ /*!
+ * \brief Deleted default constructor (so we are never in an undefined state)
+ * \note store owning pointers to laws instead if you need default-constructible objects
+ */
+ HeatPipeLaw() = delete;
+
+ /*!
+ * \brief Construct from a subgroup from the global parameter tree
+ * \note This will give you nice error messages if a mandatory parameter is missing
+ */
+ explicit HeatPipeLaw(const std::string& paramGroup)
+ {
+ params_ = makeParams(paramGroup);
+ effToAbsParams_ = EffToAbs::template makeParams(paramGroup);
+ }
+
+ /*!
+ * \brief Construct from parameter structs
+ * \note More efficient constructor but you need to ensure all parameters are initialized
+ */
+ HeatPipeLaw(const Params& params,
+ const EffToAbsParams& effToAbsParams = {})
+ : params_(params)
+ , effToAbsParams_(effToAbsParams)
+ , spline_(eps_, 1.0, // x1, x2
+ eps_*eps_*eps_, 1.0, // y1, y2
+ 3.0*eps_*eps_, 0.0) // m1, m2
+ {}
+
+ /*!
+ * \brief Construct from a subgroup from the global parameter tree
+ * \note This will give you nice error messages if a mandatory parameter is missing
+ */
+ static Params makeParams(const std::string& paramGroup)
+ {
+ const auto gamma = getParamFromGroup(paramGroup, "HeatpipeLawGamma");
+ const auto p0 = getParamFromGroup(paramGroup, "HeatpipeLawP0");
+ return {gamma, p0};
+ }
+
+ /*!
+ * \brief The capillary pressure-saturation curve.
+ *
+ * \param swe Effective saturation of the wetting phase \f$\mathrm{\overline{S}_w}\f$
+ */
+ template
+ Scalar pc(Scalar swe) const
+ {
+ const Scalar sne = 1 - swe;
+ const Scalar p0Gamma = params_.p0()*params_.gamma();
+
+ // regularization
+ if (sne >= 1.0)
+ {
+ const Scalar y = p0Gamma*( (1.263*1.0 - 2.120)*1.0 + 1.417)*1.0;
+ const Scalar m = p0Gamma*((3*1.263*1.0 - 2*2.120)*1.0 + 1.417);
+ return (sne - 1)*m + y;
+ }
+ else if (sne <= 0.0)
+ return sne * p0Gamma*1.417;
+ else
+ return p0Gamma*((1.263*sne - 2.120)*sne + 1.417) * sne;
+ }
+
+ /*!
+ * \brief The capillary pressure at Swe = 1.0 also called end point capillary pressure
+ */
+ template
+ Scalar endPointPc() const
+ { return 0.0; }
+
+ /*!
+ * \brief The partial derivative of the capillary
+ * pressure w.r.t. the effective saturation.
+ *
+ * \param swe Effective saturation of the wetting phase \f$\mathrm{\overline{S}_w}\f$
+ */
+ template
+ Scalar dpc_dswe(Scalar swe) const
+ {
+ const Scalar sne = 1 - swe;
+ const Scalar p0Gamma = params_.p0()*params_.gamma();
+ if (sne > 1.0)
+ sne = 1.0;
+ else if (sne <= 0.0)
+ return -p0Gamma*1.417;
+ else
+ return - p0Gamma*((3*1.263*sne - 2*2.120)*sne + 1.417);
+ }
+
+ /*!
+ * \brief The relative permeability for the wetting phase of
+ * the medium.
+ *
+ * \param swe The mobile saturation of the wetting phase.
+ */
+ template
+ Scalar krw(Scalar swe) const
+ {
+ return kr_(swe);
+ }
+
+ /*!
+ * \brief The relative permeability for the non-wetting phase
+ * of the medium.
+ *
+ * \param swe The mobile saturation of the wetting phase.
+ */
+ template
+ Scalar krn(Scalar swe) const
+ {
+ const Scalar sne = 1 - swe; // TODO does this make sense?
+ return kr_(sne);
+ }
+
+ /*!
+ * \brief Equality comparison with another instance
+ */
+ bool operator== (const HeatPipeLaw& o) const
+ {
+ return params_ == o.params_
+ && effToAbsParams_ == o.effToAbsParams_;
+ }
+
+ const EffToAbsParams& effToAbsParams() const
+ { return effToAbsParams_; }
+
+private:
+
+ Scalar kr_(Scalar s) const
+ {
+ if (s >= 1.0)
+ return 1;
+ else if (s <= 0.0)
+ return 0;
+ else if (s > eps_)
+ return spline_.eval(s); // regularize
+ else
+ return s*s*s;
+ }
+
+ Params params_;
+ EffToAbsParams effToAbsParams_;
+ static constexpr Scalar eps_ = 0.95;
+ Dumux::Spline spline_;
+};
+} // end namespace Dumux::FluidMatrix
#endif
diff --git a/dumux/material/fluidmatrixinteractions/2p/heatpipelawparams.hh b/dumux/material/fluidmatrixinteractions/2p/heatpipelawparams.hh
index 4f127d89c4f6a2a64556451ce51822d2027419f3..b51e786235b0ab8f64793c86bf7941ad1ac1bfe4 100644
--- a/dumux/material/fluidmatrixinteractions/2p/heatpipelawparams.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/heatpipelawparams.hh
@@ -32,6 +32,7 @@ namespace Dumux {
* \brief Reference implementation of a params for the heat pipe's
* material law
*/
+// "Use new material laws! Removal after 3.3")
template
class HeatPipeLawParams
{
diff --git a/dumux/material/fluidmatrixinteractions/2p/interfacialarea/CMakeLists.txt b/dumux/material/fluidmatrixinteractions/2p/interfacialarea/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f908e650e54f0f5f52649911c2473d76d69d6dbd
--- /dev/null
+++ b/dumux/material/fluidmatrixinteractions/2p/interfacialarea/CMakeLists.txt
@@ -0,0 +1,8 @@
+install(FILES
+exponential.hh
+exponentialcubic.hh
+interfacialarea.hh
+pcmax.hh
+polynomial2ndorder.hh
+polynomialedgezero2ndorder.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/material/fluidmatrixinteractions/2p/interfacialarea)
diff --git a/dumux/material/fluidmatrixinteractions/2p/interfacialarea/exponential.hh b/dumux/material/fluidmatrixinteractions/2p/interfacialarea/exponential.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4d3b6b2b81439865bbda73ee94bbe833314ada5e
--- /dev/null
+++ b/dumux/material/fluidmatrixinteractions/2p/interfacialarea/exponential.hh
@@ -0,0 +1,147 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see