van_genuchten.h

Go to the documentation of this file.

/*
 * richards.h
 *
 *  Created on: 21.12.2019
 *      Author: anaegel
 */
 
 
#ifndef RICHARDS__VAN_GENUCHTEN_H_
#define RICHARDS__VAN_GENUCHTEN_H_
 
// Standard lib.
#include <cmath>
 
// Use ADOL-C for automatic differentiation.
// #include <adolc/adtl.h>
 
// Use autodiff for automatic differentiation.
#include <autodiff/forward/dual.hpp>
 
 
 
// UG4 lib.
#include "common/util/smart_pointer.h"
#include "common/assert.h"
#include "registry/class.h"
 
// My libs.
#include <nlohmann/json.hpp>
#include "json_basics.hh"
 
 
 
namespace ug{
namespace Richards{
 
typedef autodiff::dual dual;
 
// const double one = 1.0;
 
inline void UGCheckValues(const dual &number)
{
  if (!std::isfinite(number.val) || std::isnan(number.val))
  {
    std::cerr << "WARNING: Invalid value:" << number.val << std::endl;
    UG_ASSERT(std::isfinite(number.val), "Function is not bounded");
  }
 
  if (!std::isfinite(number.grad) || std::isnan(number.grad))
  {
    std::cerr << "WARNING: Invalid value:" << number.grad << std::endl;
    UG_ASSERT(std::isfinite(number.grad), "Derivative is not bounded");
  }
 
}
inline void UGCheckValues(const dual &number) {…}
 
 
 
 
// Parameters for a van Genuchten model.
struct VanGenuchtenParameters
{
  double alpha=1.0;
 
  double n=2.0;
  double m=0.5;     // default: 1.0 - (1.0/n);}
 
  double thetaR=0.0;  // default: 0.0
  double thetaS=1.0;  // default: 1.0
 
  double Ksat=1.0;    // saturated conductivity K_s (optional)
 
};
struct VanGenuchtenParameters {…};
 
#ifdef UG_JSON
  NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE_WITH_DEFAULT(VanGenuchtenParameters, alpha, n, m, thetaS, thetaR, Ksat);
#endif
 
 
#ifdef UG_JSON
 void CreateJSONMap(const JSONType &array, std::map<std::string, JSONType> &map);
#endif
 
/*
class BrooksCorey
{
public:
  BrooksCorrey()
  : m_pd(1.0), m_lambda(1.0)
  {}
 
  double capillary_pressure(double Sw_eff)
  {
    return pow(m_pd*Sw_eff, 1.0/m_lambda);
  }
 
 
 
protected:
  double m_pd;      // entry pressure
  double m_lambda;    // ~pore size distribution
 
 
};
*/
 
 
 
 
 
 
template <typename TDerived>
class IRichardsModel
{
 
protected:
  // CRTP functions.
  const TDerived* me() const
  { return static_cast<const TDerived*>(this); }
  const TDerived* me() const {…}
 
  TDerived* me()
  { return static_cast<TDerived*>(this); }
  TDerived* me() {…}
 
  template <typename TFunc>
  void get_value_and_deriv(TFunc F, double H, double &f, double &df) const
  {
    dual h= H;
    /*f = F(h).val; df = derivative(F, wrt(h), at(h));*/
    auto [F0, Df0] = derivatives(F, wrt(h), at(h));
    f = F0; df = Df0;
  };
  void get_value_and_deriv(TFunc F, double H, double &f, double &df) const {…}
 
public:
  void get_saturation(double H, double &S, double &dSdH) const
  {
    auto sat = [this](dual h){ return me()->Saturation_(h);};
    get_value_and_deriv(sat, H, S, dSdH);
  }
  void get_saturation(double H, double &S, double &dSdH) const {…}
 
  void get_saturations(const double *H, double *S, double *dSdH, size_t n) const
  {
    for (size_t i=0; i<n; ++i)
    { this->get_saturation(H[i], S[i], dSdH[i]); }
  }
  void get_saturations(const double *H, double *S, double *dSdH, size_t n) const {…}
 
  void get_relative_permeability(double H, double &Kr, double &dKrdH) const
  {
    auto perm = [this](dual h){ return me()->RelativePermeability_(h);};
    get_value_and_deriv(perm, H, Kr, dKrdH);
  }
  void get_relative_permeability(double H, double &Kr, double &dKrdH) const {…}
 
  void get_relative_permeabilities(const double *H, double *K, double *dKdH, size_t n) const
  {
    for (size_t i=0; i<n; ++i)
    { this->get_relative_permeability(H[i], K[i], dKdH[i]); }
  }
  void get_relative_permeabilities(const double *H, double *K, double *dKdH, size_t n) const {…}
 
  void get_conductivity(double H, double &K, double &dKdH) const
  {
    auto cond = [this](dual h){ return me()->Conductivity_(h);};
    get_value_and_deriv(cond, H, K, dKdH);
 
  }
  void get_conductivity(double H, double &K, double &dKdH) const {…}
 
  void get_conductivities(const double *H, double *K, double *dKdH, size_t n) const
  {
    for (size_t i=0; i<n; ++i)
    { this->get_conductivity(H[i], K[i], dKdH[i]); }
  }
  void get_conductivities(const double *H, double *K, double *dKdH, size_t n) const {…}
 
 
public: // The following functions can be used from LUA
 
  // Saturation
  double Saturation(double H) const
  {
    dual h=H;
    auto func = [this](dual h){ return me()->Saturation_(h);};
    return func(h).val;
  }
  double Saturation(double H) const {…}
  double dSaturation_dH(double H) const
  {
    dual h=H;
    auto func = [this](dual h){ return me()->Saturation_(h);};
    return derivative(func, wrt(h), at(h));
  }
  double dSaturation_dH(double H) const {…}
 
  double Conductivity(double H)
  {
    dual h=H;
    auto func = [this](dual h){ return me()->Conductivity_(h);};
    return func(h).val;
  }
  double Conductivity(double H) {…}
 
  double dConductivity_dH(double H)
  {
    dual h=H;
    auto func = [this](dual h){ return me()->Conductivity_(h);};
    return derivative(func, wrt(h), at(h));
  }
  double dConductivity_dH(double H) {…}
 
};
class IRichardsModel {…};
 
 
 
struct BrooksCoreyFunctions
{
  static dual ComputeEffectiveSaturation(dual pc, double lambda, double pb)
  {
    dual val = pow(pb/pc, lambda);
    UGCheckValues(val);
    return val;
  };
  static dual ComputeEffectiveSaturation(dual pc, double lambda, double pb) {…}
};
struct BrooksCoreyFunctions {…};
 
 
template <typename TParameter>
struct IParameterizedModel
{
  IParameterizedModel() {}
 
  IParameterizedModel(const TParameter &p) : m_param(p) {}
 
  void set_parameters(const TParameter &p)
  { m_param = p; }
  void set_parameters(const TParameter &p) {…}
 
  const TParameter& get_parameters() const
  { return m_param; }
  const TParameter& get_parameters() const {…}
 
  std::string config_string() const
  {
    nlohmann::json jaux = m_param;
    return jaux.dump();
  }
  std::string config_string() const {…}
 
  TParameter m_param;
};
struct IParameterizedModel {…};
 
 
/*********************************************
 * Exponential model.
 *********************************************/
 
struct ExponentialModelParameters
{
  double pentry = 1.0;
  double alpha = 1.0;
  double beta = 1.0;
 
  double thetaR = 0.0;  // default: 0.0
  double thetaS = 1.0;  // default: 1.0
 
  double Ksat = 1.0;    // saturated conductivity K_s (optional)
};
struct ExponentialModelParameters {…};
 
#ifdef UG_JSON
  NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE_WITH_DEFAULT(ExponentialModelParameters, pentry, alpha, beta, thetaR, thetaS, Ksat)
#endif
 
 
class ExponentialModel :
  public IRichardsModel<ExponentialModel>,
  public IParameterizedModel<ExponentialModelParameters>
{
public:
  typedef IRichardsModel<ExponentialModel> base_type;
  typedef ExponentialModelParameters parameter_type;
  typedef IParameterizedModel<ExponentialModelParameters> parameterized_model_type;
 
  friend base_type;
  friend parameterized_model_type;
 
  ExponentialModel(const parameter_type &p) : parameterized_model_type(p) {}
 
protected:
 
  // p_c = p_a - p_w         : Kapillardruck >=0
  // psi = p_c / (rho_w * g) : Kapillardruckhoehe >=0
 
  inline dual EffSaturation_(dual pc) const
  {
    const parameter_type &p = this->get_parameters();
    return exp(-p.alpha * (pc/p.pentry));
  }
  inline dual EffSaturation_(dual pc) const {…}
 
  inline dual Saturation_(dual pc) const
  {
    const parameter_type &p = this->get_parameters();
    UGCheckValues(pc);
 
    if (pc.val <= 0) return p.thetaS; // Medium is saturated
    const dual Seff = EffSaturation_(pc); UGCheckValues(Seff);
    //UG_LOG("seff:" << (pc.val/p.pentry) << " -> " << Seff.val << std::endl);
    return p.thetaR + (p.thetaS-p.thetaR)*Seff;
 
  }
  inline dual Saturation_(dual pc) const {…}
 
  inline dual RelativePermeability_(dual pc) const
  {
    const parameter_type &p = this->get_parameters();
    UGCheckValues(pc);
 
    if (pc.val<= 0) return 1.0;
    const dual kval = exp(-p.beta * (pc/p.pentry));
 
    //UG_LOG("k:" << (pc.val/p.pentry) << " -> " << kval.val << std::endl);
    return kval;
 
  }
  inline dual RelativePermeability_(dual pc) const {…}
 
  // Conductivity C:=K_{sat}*k_r
  inline dual Conductivity_(dual pc) const
  {
    return get_parameters().Ksat*RelativePermeability_(pc);
  }
  inline dual Conductivity_(dual pc) const {…}
 
};
class ExponentialModel : {…};
 
/*********************************************
 * Van Genuchten-Mualem
 *********************************************/
 
struct VanGenuchtenFunctions
{
 
  static dual ComputeEffectiveSaturation(dual psi_, double alpha, double n, double m)
  {
        UGCheckValues(psi_);
        dual apn = alpha*psi_;
        UGCheckValues(apn);
 
        apn = pow(apn, n);
        UGCheckValues(apn);
 
        dual val = pow(1.0/(1.0+apn), m);
        UGCheckValues(val);
        return val;
  }
  static dual ComputeEffectiveSaturation(dual psi_, double alpha, double n, double m) {…}
 
  static dual ComputeEffectiveSaturation(dual psi_, double alpha, double n)
  { return ComputeEffectiveSaturation(psi_, alpha, n, 1.0-(1.0/n)); }
  static dual ComputeEffectiveSaturation(dual psi_, double alpha, double n) {…}
};
struct VanGenuchtenFunctions {…};
 
 
 
class VanGenuchtenModel
: public IRichardsModel<VanGenuchtenModel>, public IParameterizedModel<VanGenuchtenParameters>
{
public:
  typedef IRichardsModel<VanGenuchtenModel> base_type;
  typedef IParameterizedModel<VanGenuchtenParameters> parameterized_model_type;
  typedef VanGenuchtenParameters parameter_type;
 
  friend base_type;
  friend parameterized_model_type;
 
  VanGenuchtenModel(const parameter_type &p) :  parameterized_model_type(p) {}
 
 
protected:
 
  // p_c = p_a - p_w         : Kapillardruck
  // psi = p_c / (rho_w * g) : Kapillardruckhoehe
 
  inline dual EffSaturation_(dual psi_) const
  {
    const VanGenuchtenParameters &p = m_param;
    return VanGenuchtenFunctions::ComputeEffectiveSaturation(psi_, p.alpha, p.n, p.m);
  }
  inline dual EffSaturation_(dual psi_) const {…}
 
  inline dual Saturation_(dual psi_) const
  {
    const VanGenuchtenParameters &p = m_param;
    UGCheckValues(psi_);
 
    if (psi_ <= 0) return p.thetaS;
    else
    {
      dual Seff = EffSaturation_(psi_);
      UGCheckValues(Seff);
      return p.thetaR + (p.thetaS-p.thetaR)*Seff;
    }
  }
  inline dual Saturation_(dual psi_) const {…}
 
  inline dual RelativePermeability_(dual psi_) const
  {
    const VanGenuchtenParameters &p = m_param;
    UGCheckValues(psi_);
 
    if (psi_<= 0) return 1.0;
 
    dual Seff = EffSaturation_(psi_);   // Seff -> 1 is dangerous!
    if (fabs(Seff.val-1.0)<1e-15) return 1.0;
    UGCheckValues(Seff);
 
    double m = p.m;
    dual brackS = 1.0-pow(Seff,1.0/m);  // => brackS -> 0.0
    UGCheckValues(brackS);
 
    dual auxS = pow(brackS,m);      // => this derivative may explode!
    UGCheckValues(auxS);
 
    return sqrt(Seff)*(1.0-auxS)*(1.0-auxS); // check!
  }
  inline dual RelativePermeability_(dual psi_) const {…}
 
  // Conductivity C:=K_{sat}*k_r
  inline dual Conductivity_(dual psi_) const
  {
    const VanGenuchtenParameters &p = m_param;
    return p.Ksat*RelativePermeability_(psi_);
  }
  inline dual Conductivity_(dual psi_) const {…}
 
};
class VanGenuchtenModel {…};
 
 
/*********************************************
 * Haverkamp model.
 *********************************************/
 
 
 
// Parameters for Haverkamp model.
struct HaverkampParameters
{
  double alpha;
  double n;
 
  double beta;
  double m;
 
  double thetaR=0.0;  // default: 0.0
  double thetaS=1.0;  // default: 1.0
 
  double Ksat=1.0;    // saturated conductivity K_s (optional)
};
struct HaverkampParameters {…};
 
 
 
#ifdef UG_JSON
 
 void to_json(JSONType& j, const HaverkampParameters& p);
 
  void from_json(const JSONType& j, HaverkampParameters& p);
#endif
 
 
 
 
class HaverkampModel  : public IRichardsModel<HaverkampModel>, public IParameterizedModel<HaverkampParameters>
{
public:
  typedef IRichardsModel<HaverkampModel> base_type;
  typedef IParameterizedModel<HaverkampParameters> parameterized_model_type;
  typedef HaverkampParameters parameter_type;
 
  friend base_type;
  friend parameterized_model_type;
 
  HaverkampModel(const parameter_type &p) : parameterized_model_type(p)
  {}
  HaverkampModel(const parameter_type &p) : parameterized_model_type(p) {…}
 
#ifdef UG_JSON
  // TODO: Move to Factory
  HaverkampModel(const char* json)
  {
    nlohmann::json j = nlohmann::json::parse(json);
    parameter_type p = j.get<parameter_type>();
    set_parameters(p);
  }
#endif
 
 
 
protected:
 
 
  // p_c = p_a - p_w         : Kapillardruck
  // psi = p_c / (rho_w * g) : Kapillardruckhoehe
 
 
  dual Saturation_(dual psi_) const
  {
    const parameter_type &p = get_parameters();
    if (psi_ <= 0.0) return p.thetaS;
 
    dual Seff = VanGenuchtenFunctions::ComputeEffectiveSaturation(psi_, p.alpha, p.n, 1.0);
    return p.thetaR + (p.thetaS-p.thetaR)*Seff;
 
  }
  dual Saturation_(dual psi_) const {…}
 
  dual RelativePermeability_(dual psi_) const
  {
    const parameter_type &p = get_parameters();
    if (psi_<= 0.0) return 1.0;
 
    return VanGenuchtenFunctions::ComputeEffectiveSaturation(psi_, p.beta, p.m, 1.0);
  }
  dual RelativePermeability_(dual psi_) const {…}
 
  // Conductivity K:=K_{sat}*k_r
  dual Conductivity_(dual psi_) const
  {
    const parameter_type &p = get_parameters();
    return p.Ksat*RelativePermeability_(psi_);
  }
  dual Conductivity_(dual psi_) const {…}
 
};
class HaverkampModel  : public IRichardsModel<HaverkampModel>, public IParameterizedModel<HaverkampParameters> {…};
 
 
class GardnerModel :
    public IParameterizedModel<HaverkampParameters>,
    public IRichardsModel<GardnerModel>
{
public:
  typedef IRichardsModel<GardnerModel> base_type;
  typedef IParameterizedModel<HaverkampParameters> parameterized_model_type;
  typedef VanGenuchtenParameters parameter_type;
 
  GardnerModel(const VanGenuchtenParameters &p) : m_param(p) {}
 
  void get_saturation(double H, double &S, double &dSdH) const
  {
    auto sat = [this](dual h){ return me()->Saturation_(h, m_param);};
    base_type::get_value_and_deriv(sat, H, S, dSdH);
  }
  void get_saturation(double H, double &S, double &dSdH) const {…}
 
  void get_conductivity(double H, double &K, double &dKdH) const
  {
    auto cond = [this](dual h){ return Conductivity_(h, m_param);};
    base_type::get_value_and_deriv(cond, H, K, dKdH);
  }
  void get_conductivity(double H, double &K, double &dKdH) const {…}
 
  std::string config_string() const;
 
protected:
  static double one;
  VanGenuchtenParameters m_param;
 
 
  // p_c = p_a - p_w         : Kapillardruck
  // hpc = p_c / (rho_w * g) : Kapillardruckhoehe
 
 
  static dual EffSaturation_(dual psi_, const VanGenuchtenParameters &p)
  { return exp(-p.alpha*psi_); }
  static dual EffSaturation_(dual psi_, const VanGenuchtenParameters &p) {…}
 
  static dual Saturation_(dual psi_, const VanGenuchtenParameters &p)
  {
    if (psi_ <= 0) return p.thetaS;
    else
    {
      dual Seff = EffSaturation_(psi_, p);
      return p.thetaR + (p.thetaS-p.thetaR)*Seff;
    }
  }
  static dual Saturation_(dual psi_, const VanGenuchtenParameters &p) {…}
 
  // Conductivity K_s*k_r
  static dual Conductivity_(dual psi_ , const VanGenuchtenParameters &p)
  { return p.Ksat* EffSaturation_(psi_, p); }
  static dual Conductivity_(dual psi_ , const VanGenuchtenParameters &p) {…}
 
};
class GardnerModel : {…};
 
 
 
 
#ifdef UG_JSON
 
 
template <typename TModel, typename TParameter = typename TModel::parameter_type>
SmartPtr<TModel> CreateModel(const char *jstring) {
 
  SmartPtr<TModel> inst = SPNULL;
  try
  {
      nlohmann::json j = nlohmann::json::parse(jstring);
      TParameter p = j.get<TParameter>();
      inst = make_sp(new TModel(p));
  }
  catch (...)
  {
      UG_THROW ("Construction failed: " << jstring);
  }
  return inst;
};
 
 
 
SmartPtr<VanGenuchtenModel> CreateVanGenuchtenModel(const char *json);
 
 
struct RichardsModelFactory {
  SmartPtr<ExponentialModel> create_exponential(const char *json);
  SmartPtr<VanGenuchtenModel> create_van_genuchten(const char *json);
  SmartPtr<HaverkampModel> create_haverkamp(const char *json);
};
 
#endif
 
 
 
 
}
 
 
#ifdef UG_JSON
  template <> struct is_json_constructible<typename Richards::ExponentialModelParameters> { const static bool value = true;};
#endif
}
 
 
 
#ifdef UG_JSON
// Finally, we need to define how to construct from JSON.
namespace nlohmann {
 
  template <>
  struct adl_serializer<ug::Richards::VanGenuchtenModel>
  {
       static ug::Richards::VanGenuchtenModel from_json(const json& j)
       { return j.get<typename ug::Richards::VanGenuchtenParameters>(); }; // initialize from model }
 
       static void to_json(json& j, typename ug::Richards::VanGenuchtenModel m)
       { j = m.get_parameters();}
   };
};
#endif
 
#endif /* RICHARDS__VAN_GENUCHTEN_H_ */