Loading [MathJax]/extensions/tex2jax.js
OGS
RelPermNonWettingPhaseVanGenuchtenMualem.cpp
Go to the documentation of this file.
1
13
14#include <algorithm>
15#include <cmath>
16#include <limits>
17
18#include "BaseLib/Error.h"
22
23namespace MaterialPropertyLib
24{
25double computeVanGenuchtenMualemValue(const double S_L, const double S_L_r,
26 const double S_L_max, const double m)
27{
28 const double Se = (S_L - S_L_r) / (S_L_max - S_L_r);
29 return std::sqrt(1.0 - Se) * std::pow(1.0 - std::pow(Se, 1.0 / m), 2.0 * m);
30}
31
34 const double S_L_r,
35 const double S_n_r,
36 const double m,
37 const double krel_min,
38 const double a)
39 : S_L_r_(S_L_r),
40 S_L_max_(1. - S_n_r),
41 m_(m),
42 krel_min_(krel_min),
43 S_L_for_krel_min_(computeSaturationForMinimumRelativePermeability()),
44 a_(a)
45{
46 name_ = std::move(name);
48}
49
51 VariableArray const& variable_array,
52 ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
53 double const /*dt*/) const
54{
55 const double S_L =
56 std::clamp(variable_array.liquid_saturation, S_L_r_, S_L_max_);
57
58 const double krel =
60 return std::max(krel_min_, a_ * krel);
61}
62
64 VariableArray const& variable_array, Variable const variable,
65 ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
66 double const /*dt*/) const
67{
68 if (variable != Variable::liquid_saturation)
69 {
71 "RelPermNonWettingPhaseVanGenuchtenMualem::dValue is implemented "
72 "for the derivative with respect to liquid saturation only.");
73 }
74
75 const double S_L = variable_array.liquid_saturation;
76 if (S_L < S_L_r_ || S_L > S_L_for_krel_min_)
77 {
78 return 0.0;
79 }
80
81 if (std::fabs(S_L - S_L_max_) < std::numeric_limits<double>::epsilon())
82 {
83 return 0.0;
84 }
85
86 const double Se = (S_L - S_L_r_) / (S_L_max_ - S_L_r_);
87
88 const double val1 = std::sqrt(1.0 - Se);
89 const double val2 = 1.0 - std::pow(Se, 1.0 / m_);
90
91 return a_ *
92 (-0.5 * std::pow(val2, 2.0 * m_) / val1 -
93 2.0 * std::pow(Se, -1.0 + 1.0 / m_) * val1 *
94 std::pow(val2, 2.0 * m_ - 1.0)) /
95 (S_L_max_ - S_L_r_);
96}
97
100{
101 auto f = [this](double S_L) -> double
102 {
103 return computeVanGenuchtenMualemValue(S_L, this->S_L_r_, this->S_L_max_,
104 this->m_) -
105 this->krel_min_;
106 };
107
108 auto root_finder =
110 f, S_L_r_, S_L_max_);
111
112 root_finder.step(1000);
113
114 return root_finder.getResult();
115}
116} // namespace MaterialPropertyLib
#define OGS_FATAL(...)
Definition Error.h:26
virtual PropertyDataType value() const
Definition Property.cpp:76
RelPermNonWettingPhaseVanGenuchtenMualem(std::string name, const double S_L_r, const double S_n_r, const double m, const double krel_min, const double a)
PropertyDataType dValue(VariableArray const &variable_array, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const override
double computeVanGenuchtenMualemValue(const double S_L, const double S_L_r, const double S_L_max, const double m)
void checkVanGenuchtenExponentRange(const double m)
std::variant< double, Eigen::Matrix< double, 2, 1 >, Eigen::Matrix< double, 3, 1 >, Eigen::Matrix< double, 2, 2 >, Eigen::Matrix< double, 3, 3 >, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 >, Eigen::MatrixXd > PropertyDataType
Definition Property.h:31
RegulaFalsi< SubType, Function > makeRegulaFalsi(Function &&f, double const a, double const b)
Definition Root1D.h:137