vangenuchten.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_MATERIAL_FLUIDMATRIX_VAN_GENUCHTEN_HH
#define DUMUX_MATERIAL_FLUIDMATRIX_VAN_GENUCHTEN_HH
 
#include <cmath>
#include <algorithm>
 
#include <dumux/common/parameters.hh>
#include <dumux/common/spline.hh>
#include <dumux/common/optionalscalar.hh>
#include <dumux/material/fluidmatrixinteractions/2p/materiallaw.hh>
 
namespace Dumux::FluidMatrix {
 
class VanGenuchten
{
 
public:
    template<class Scalar>
    struct Params
    {
        Params(Scalar alpha, Scalar n, Scalar l = 0.5)
        : alpha_(alpha), n_(n), m_(1.0 - 1.0/n), l_(l)
        {}
 
        Scalar alpha() const { return alpha_; }
        void setAlpha(Scalar alpha) { alpha_ = alpha; }
 
        Scalar m() const { return m_; }
        void setM(Scalar m) { m_ = m; n_ = 1.0/(1.0 - m); }
 
        Scalar n() const{ return n_; }
        void setN(Scalar n){ n_ = n; m_ = 1.0 - 1.0/n; }
 
        Scalar l() const { return l_; }
        void setL(Scalar l) { l_ = l; }
 
        bool operator== (const Params& p) const
        {
            return Dune::FloatCmp::eq(alpha_, p.alpha_, 1e-6)
                   && Dune::FloatCmp::eq(n_, p.n_, 1e-6)
                   && Dune::FloatCmp::eq(m_, p.m_, 1e-6)
                   && Dune::FloatCmp::eq(l_, p.l_, 1e-6);
        }
 
    private:
        Scalar alpha_, n_, m_, l_;
    };
 
    template<class Scalar = double>
    static Params<Scalar> makeParams(const std::string& paramGroup)
    {
        const auto n = getParamFromGroup<Scalar>(paramGroup, "VanGenuchtenN");
        const auto alpha = getParamFromGroup<Scalar>(paramGroup, "VanGenuchtenAlpha");
        // l is usually chosen to be 0.5 (according to Mualem (1976), WRR)
        const auto l = getParamFromGroup<Scalar>(paramGroup, "VanGenuchtenL", 0.5);
        return Params<Scalar>(alpha, n, l);
    }
 
    template<class Scalar>
    static Scalar pc(Scalar swe, const Params<Scalar>& 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 pc = pow(pow(swe, -1.0/params.m()) - 1, 1.0/params.n())/params.alpha();
        return pc;
    }
 
    template<class Scalar>
    static Scalar swe(Scalar pc, const Params<Scalar>& params)
    {
        using std::pow;
        using std::max;
 
        pc = max(pc, 0.0); // the equation below is undefined for negative pcs
 
        const Scalar sw = pow(pow(params.alpha()*pc, params.n()) + 1, -params.m());
        return sw;
    }
 
    template<class Scalar>
    static Scalar endPointPc(const Params<Scalar>& params)
    { return 0.0; }
 
    template<class Scalar>
    static Scalar dpc_dswe(Scalar swe, const Params<Scalar>& 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 powSwe = pow(swe, -1/params.m());
        return - 1.0/params.alpha() * pow(powSwe - 1, 1.0/params.n() - 1)/params.n()
                                  * powSwe/swe/params.m();
    }
 
    template<class Scalar>
    static Scalar dswe_dpc(Scalar pc, const Params<Scalar>& params)
    {
        using std::pow;
        using std::max;
 
        pc = max(pc, 0.0); // the equation below is undefined for negative pcs
 
        const Scalar powAlphaPc = pow(params.alpha()*pc, params.n());
        return -pow(powAlphaPc + 1, -params.m()-1)*params.m()*powAlphaPc/pc*params.n();
    }
 
    template<class Scalar>
    static Scalar krw(Scalar swe, const Params<Scalar>& 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 r = 1.0 - pow(1.0 - pow(swe, 1.0/params.m()), params.m());
        return pow(swe, params.l())*r*r;
    }
 
    template<class Scalar>
    static Scalar dkrw_dswe(Scalar swe, const Params<Scalar>& 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 x = 1.0 - pow(swe, 1.0/params.m());
        const Scalar xToM = pow(x, params.m());
        return (1.0 - xToM)*pow(swe, params.l()-1) * ( (1.0 - xToM)*params.l() + 2*xToM*(1.0-x)/x );
    }
 
    template<class Scalar>
    static Scalar krn(Scalar swe, const Params<Scalar>& 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(1 - swe, params.l()) * pow(1 - pow(swe, 1.0/params.m()), 2*params.m());
    }
 
    template<class Scalar>
    static Scalar dkrn_dswe(Scalar swe, const Params<Scalar>& 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 sne = 1.0 - swe;
        const auto x = 1.0 - pow(swe, 1.0/params.m());
        return -pow(sne, params.l()-1.0) * pow(x, 2*params.m() - 1.0) * ( params.l()*x + 2.0*sne/swe*(1.0 - x) );
    }
};
 
template <class Scalar>
class VanGenuchtenRegularization
{
public:
    template<class S>
    struct Params
    {
        void setPcLowSwe(Scalar pcLowSwe)
        { pcLowSwe_ = pcLowSwe; }
 
        Scalar pcLowSwe() const
        { return pcLowSwe_; }
 
        void setPcHighSwe(Scalar pcHighSwe)
        { pcHighSwe_ = pcHighSwe; }
 
        Scalar pcHighSwe() const
        { return pcHighSwe_; }
 
        void setKrnLowSwe(Scalar krnLowSwe)
        { krnLowSwe_ = krnLowSwe; }
 
        Scalar krnLowSwe() const
        { return krnLowSwe_; }
 
        void setKrwHighSwe(Scalar krwHighSwe)
        { krwHighSwe_ = krwHighSwe; }
 
        Scalar krwHighSwe() const
        { return krwHighSwe_; }
 
private:
        S pcLowSwe_ = 0.01;
        S pcHighSwe_ = 0.99;
        S krnLowSwe_ = 0.1;
        S krwHighSwe_ = 0.9;
    };
 
    template<class MaterialLaw>
    void init(const MaterialLaw* m, const std::string& paramGroup)
    {
        pcLowSwe_ = getParamFromGroup<Scalar>(paramGroup, "VanGenuchtenPcLowSweThreshold", 0.01);
        pcHighSwe_ = getParamFromGroup<Scalar>(paramGroup, "VanGenuchtenPcHighSweThreshold", 0.99);
        krwHighSwe_ = getParamFromGroup<Scalar>(paramGroup, "VanGenuchtenKrwHighSweThreshold", 0.9);
        krnLowSwe_ = getParamFromGroup<Scalar>(paramGroup, "VanGenuchtenKrnLowSweThreshold", 0.1);
 
        initPcParameters_(m, pcLowSwe_, pcHighSwe_);
        initKrParameters_(m, krnLowSwe_, krwHighSwe_);
    }
 
    template<class MaterialLaw, class BaseParams, class EffToAbsParams>
    void init(const MaterialLaw* m, const BaseParams& bp, const EffToAbsParams& etap, const Params<Scalar>& p)
    {
        pcLowSwe_ = p.pcLowSwe();
        pcHighSwe_ = p.pcHighSwe();
        krwHighSwe_ = p.krwHighSwe();
        krnLowSwe_ = p.krnLowSwe();
 
        initPcParameters_(m, pcLowSwe_, pcHighSwe_);
        initKrParameters_(m, krnLowSwe_, krwHighSwe_);
    }
 
    bool operator== (const VanGenuchtenRegularization& o) const
    {
        return Dune::FloatCmp::eq(pcLowSwe_, o.pcLowSwe_, 1e-6)
               && Dune::FloatCmp::eq(pcHighSwe_, o.pcHighSwe_, 1e-6)
               && Dune::FloatCmp::eq(krwHighSwe_, o.krwHighSwe_, 1e-6)
               && Dune::FloatCmp::eq(krnLowSwe_, o.krnLowSwe_, 1e-6);
    }
 
    OptionalScalar<Scalar> 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 pcDerivativeHighSweEnd_*(swe - 1.0);
 
        else if (swe > pcHighSwe_)
            return pcSpline_.eval(swe);
 
        else
            return {}; // no regularization
    }
 
    OptionalScalar<Scalar> dpc_dswe(const Scalar swe) const
    {
        if (swe <= pcLowSwe_)
            return pcDerivativeLowSw_;
 
        else if (swe >= 1.0)
            return pcDerivativeHighSweEnd_;
 
        else if (swe > pcHighSwe_)
            return pcSpline_.evalDerivative(swe);
 
        else
            return {}; // no regularization
    }
 
    OptionalScalar<Scalar> swe(const Scalar pc) const
    {
        if (pc <= 0.0)
        {
            if (pcHighSwe_ >= 1.0)
                return 1.0;
            else
                return pc/pcDerivativeHighSweEnd_ + 1.0;
        }
 
        // invert spline
        else if (pc <= pcHighSwePcValue_)
            return pcSpline_.intersectInterval(pcHighSwe_, 1.0, 0.0, 0.0, 0.0, pc);
 
        else if (pc >= pcLowSwePcValue_)
            return (pc - pcLowSwePcValue_)/pcDerivativeLowSw_ + pcLowSwe_;
 
        else
            return {}; // no regularization
    }
 
    OptionalScalar<Scalar> dswe_dpc(const Scalar pc) const
    {
        if (pc <= 0.0)
        {
            if (pcHighSwe_ >= 1.0)
                return 0.0;
            else
                return 1.0/pcDerivativeHighSweEnd_;
        }
 
        // derivative of the inverse of the function is one over derivative of the function
        else if (pc <= pcHighSwePcValue_)
            return 1.0/pcSpline_.evalDerivative(pcSpline_.intersectInterval(pcHighSwe_, 1.0, 0.0, 0.0, 0.0, pc));
 
        else if (pc >= pcLowSwePcValue_)
            return 1.0/pcDerivativeLowSw_;
 
        else
            return {}; // no regularization
    }
 
    OptionalScalar<Scalar> krw(const Scalar swe) const
    {
        if (swe <= 0.0)
            return 0.0;
        else if (swe >= 1.0)
            return 1.0;
        else if (swe >= krwHighSwe_)
            return krwSpline_.eval(swe);
        else
            return {}; // no regularization
    }
 
    OptionalScalar<Scalar> dkrw_dswe(const Scalar swe) const
    {
        if (swe <= 0.0)
            return 0.0;
        else if (swe >= 1.0)
            return 0.0;
        else if (swe >= krwHighSwe_)
            return krwSpline_.evalDerivative(swe);
        else
            return {}; // no regularization
    }
 
    OptionalScalar<Scalar> krn(const Scalar swe) const
    {
        if (swe <= 0.0)
            return 1.0;
        else if (swe >= 1.0)
            return 0.0;
        else if (swe <= krnLowSwe_)
            return krnSpline_.eval(swe);
        else
            return {}; // no regularization
    }
 
    OptionalScalar<Scalar> dkrn_dswe(const Scalar swe) const
    {
        if (swe <= 0.0)
            return 0.0;
        else if (swe >= 1.0)
            return 0.0;
        else if (swe <= krnLowSwe_)
            return krnSpline_.evalDerivative(swe);
        else
            return {}; // no regularization
    }
 
private:
    template<class MaterialLaw>
    void initPcParameters_(const MaterialLaw* m, const Scalar lowSwe, const Scalar highSwe)
    {
        const auto lowSw = MaterialLaw::EffToAbs::sweToSw(lowSwe, m->effToAbsParams());
        const auto highSw = MaterialLaw::EffToAbs::sweToSw(highSwe, m->effToAbsParams());
        const auto dsw_dswe = MaterialLaw::EffToAbs::dsw_dswe(m->effToAbsParams());
 
        pcDerivativeLowSw_ = m->template dpc_dsw<false>(lowSw)*dsw_dswe;
 
        pcDerivativeHighSweThreshold_ = m->template dpc_dsw<false>(highSw)*dsw_dswe;
        pcDerivativeHighSweEnd_ = 2.0*(0.0 - m->template pc<false>(highSw))/(1.0 - highSwe);
 
        pcLowSwePcValue_ = m->template pc<false>(lowSw);
        pcHighSwePcValue_ = m->template pc<false>(highSw);
 
        // Only initialize regularization spline if given parameters are in
        // the admissible range. When constructing with non-sensible parameters
        // the spline construction might fail (e.g. highSwe == 1.0)
        if (highSwe < 1.0)
            pcSpline_ = Spline<Scalar>(highSwe, 1.0, // x0, x1
                                       pcHighSwePcValue_, 0, // y0, y1
                                       pcDerivativeHighSweThreshold_, pcDerivativeHighSweEnd_); // m0, m1
    }
 
    template<class MaterialLaw>
    void initKrParameters_(const MaterialLaw* m, const Scalar lowSwe, const Scalar highSwe)
    {
        const auto lowSw = MaterialLaw::EffToAbs::sweToSw(lowSwe, m->effToAbsParams());
        const auto highSw = MaterialLaw::EffToAbs::sweToSw(highSwe, m->effToAbsParams());
        const auto dsw_dswe = MaterialLaw::EffToAbs::dsw_dswe(m->effToAbsParams());
 
        const auto krwHighSw = m->template krw<false>(highSw);
        const auto dkrwHighSw = m->template dkrw_dsw<false>(highSw)*dsw_dswe;
 
        const auto krnLowSw = m->template krn<false>(lowSw);
        const auto dkrnLowSw = m->template dkrn_dsw<false>(lowSw)*dsw_dswe;
 
        if (highSwe < 1.0)
            krwSpline_ = Spline<Scalar>(highSwe, 1.0, // x0, x1
                                        krwHighSw, 1.0, // y0, y1
                                        dkrwHighSw, 0.0); // m0, m1
 
        if (lowSwe > 0.0)
            krnSpline_ = Spline<Scalar>(0.0, lowSwe, // x0, x1
                                        1.0, krnLowSw, // y0, y1
                                        0.0, dkrnLowSw); // m0, m1
    }
 
    Scalar pcLowSwe_, pcHighSwe_;
    Scalar pcLowSwePcValue_, pcHighSwePcValue_;
    Scalar krwHighSwe_, krnLowSwe_;
    Scalar pcDerivativeLowSw_;
    Scalar pcDerivativeHighSweThreshold_, pcDerivativeHighSweEnd_;
 
    Spline<Scalar> pcSpline_;
    Spline<Scalar> krwSpline_;
    Spline<Scalar> krnSpline_;
};
 
template<typename Scalar = double>
using VanGenuchtenDefault = TwoPMaterialLaw<Scalar, VanGenuchten, VanGenuchtenRegularization<Scalar>, TwoPEffToAbsDefaultPolicy>;
 
template<typename Scalar = double>
using VanGenuchtenNoReg = TwoPMaterialLaw<Scalar, VanGenuchten, NoRegularization, TwoPEffToAbsDefaultPolicy>;
 
} // end namespace Dumux::FluidMatrix
 
#endif