diff --git a/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/c_SaturationVanGenuchten.md b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/c_SaturationVanGenuchten.md
index 2a8f10769050f001e0532878f196e583db65c309..8acbdb5f1d2891cacb5dc94a11ba6f060e993c5a 100644
--- a/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/c_SaturationVanGenuchten.md
+++ b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/c_SaturationVanGenuchten.md
@@ -1 +1,6 @@
\copydoc MaterialPropertyLib::SaturationVanGenuchten
+
+Either the `exponent` (which is the pressure exponent \f$m\f$ in the equations)
+must be set and then the saturation exponent is \f$n = 1 / (1 - m)\f$.
+Or both values the `pressure_exponent` and the `saturation_exponent` must be set
+independently.
diff --git a/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_exponent.md b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_exponent.md
index 260b1dd015e7c9cb8b04a65e4df7587dc81b9224..7456dd0d43f8fe8a61723d09ef9bbfd47811c550 100644
--- a/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_exponent.md
+++ b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_exponent.md
@@ -1 +1,5 @@
-The exponent of the van Genuchten saturation function.
+The "pressure" exponent \f$m\f$ of the van Genuchten saturation function.
+The "saturation" exponent is then set to \f$n = 1 / (1 - m)\f$.
+
+See `pressure_exponent` and `saturation_exponent` to set the values
+independently.
diff --git a/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_pressure_exponent.md b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_pressure_exponent.md
new file mode 100644
index 0000000000000000000000000000000000000000..abbc870763d012f4e45b91226844c4c0489f136d
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_pressure_exponent.md
@@ -0,0 +1,2 @@
+The "pressure" exponent \f$m\f$ of the van Genuchten saturation function.
+The "saturation" exponent \f$n\f$ must be set as well.
diff --git a/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_saturation_exponent.md b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_saturation_exponent.md
new file mode 100644
index 0000000000000000000000000000000000000000..cfbb9085e467f847aca9c1924e29338e3d9228ab
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_saturation_exponent.md
@@ -0,0 +1,2 @@
+The "saturation" exponent \f$n\f$ of the van Genuchten saturation function.
+The "pressure" exponent \f$m\f$ must be set as well.
diff --git a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationVanGenuchten.cpp b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationVanGenuchten.cpp
index cd874d71d229af667f5ab4f9e2227f3e00277716..e929134f0e291b81660c35189e50d4af7dec09be 100644
--- a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationVanGenuchten.cpp
+++ b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationVanGenuchten.cpp
@@ -31,14 +31,31 @@ std::unique_ptr<SaturationVanGenuchten> createSaturationVanGenuchten(
auto const residual_gas_saturation =
//! \ogs_file_param{properties__property__SaturationVanGenuchten__residual_gas_saturation}
config.getConfigParameter<double>("residual_gas_saturation");
- auto const exponent =
- //! \ogs_file_param{properties__property__SaturationVanGenuchten__exponent}
- config.getConfigParameter<double>("exponent");
+
+ double pressure_exponent;
+ double saturation_exponent;
+ if (auto const optional_exponent =
+ //! \ogs_file_param{properties__property__SaturationVanGenuchten__exponent}
+ config.getConfigParameterOptional<double>("exponent"))
+ {
+ pressure_exponent = *optional_exponent;
+ saturation_exponent = 1.0 / (1.0 - pressure_exponent);
+ }
+ else
+ {
+ pressure_exponent =
+ //! \ogs_file_param{properties__property__SaturationVanGenuchten__pressure_exponent}
+ config.getConfigParameter<double>("pressure_exponent");
+ saturation_exponent =
+ //! \ogs_file_param{properties__property__SaturationVanGenuchten__saturation_exponent}
+ config.getConfigParameter<double>("saturation_exponent");
+ }
+
//! \ogs_file_param{properties__property__SaturationVanGenuchten__p_b}
auto const p_b = config.getConfigParameter<double>("p_b");
return std::make_unique<SaturationVanGenuchten>(
std::move(property_name), residual_liquid_saturation,
- residual_gas_saturation, exponent, p_b);
+ residual_gas_saturation, pressure_exponent, saturation_exponent, p_b);
}
} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchten.cpp b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchten.cpp
index 0d6bfd75fd1cf05d400a1b1f4e05eb8206456291..8fd236cf90cc5c2a0af24b90109b8c1c238a9298 100644
--- a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchten.cpp
+++ b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchten.cpp
@@ -21,11 +21,13 @@ SaturationVanGenuchten::SaturationVanGenuchten(
std::string name,
double const residual_liquid_saturation,
double const residual_gas_saturation,
- double const exponent,
+ double const pressure_exponent,
+ double const saturation_exponent,
double const p_b)
: S_L_res_(residual_liquid_saturation),
S_L_max_(1. - residual_gas_saturation),
- m_(exponent),
+ m_(pressure_exponent),
+ n_(saturation_exponent),
p_b_(p_b)
{
name_ = std::move(name);
@@ -33,10 +35,18 @@ SaturationVanGenuchten::SaturationVanGenuchten(
if (!(m_ > 0 && m_ < 1))
{
OGS_FATAL(
- "The exponent value m = {:g} of van Genuchten saturation model, is "
- "out of its range of (0, 1)",
+ "The pressure exponent value m = {:g} of van Genuchten saturation "
+ "model, is out of its range of (0, 1)",
m_);
}
+
+ if (n_ < 1)
+ {
+ OGS_FATAL(
+ "The saturation exponent value n = {:g} of van Genuchten "
+ "saturation model, is out of its range of [1, +inf)",
+ n_);
+ }
}
PropertyDataType SaturationVanGenuchten::value(
@@ -52,8 +62,7 @@ PropertyDataType SaturationVanGenuchten::value(
}
double const p = p_cap / p_b_;
- double const n = 1. / (1. - m_);
- double const p_to_n = std::pow(p, n);
+ double const p_to_n = std::pow(p, n_);
double const S_eff = std::pow(p_to_n + 1., -m_);
double const S = S_eff * S_L_max_ - S_eff * S_L_res_ + S_L_res_;
@@ -80,8 +89,7 @@ PropertyDataType SaturationVanGenuchten::dValue(
}
double const p = p_cap / p_b_;
- double const n = 1. / (1. - m_);
- double const p_to_n = std::pow(p, n);
+ double const p_to_n = std::pow(p, n_);
double const S_eff = std::pow(p_to_n + 1., -m_);
double const S = S_eff * S_L_max_ - S_eff * S_L_res_ + S_L_res_;
@@ -91,9 +99,8 @@ PropertyDataType SaturationVanGenuchten::dValue(
return 0.;
}
- double const dS_eff_dp_cap = -m_ * std::pow(p, n - 1) *
- std::pow(1 + p_to_n, -1. - m_) /
- (p_b_ * (1. - m_));
+ double const dS_eff_dp_cap =
+ -m_ * n_ * p_to_n * S_eff / (p_cap * (p_to_n + 1.));
return dS_eff_dp_cap * (S_L_max_ - S_L_res_);
}
@@ -117,8 +124,7 @@ PropertyDataType SaturationVanGenuchten::d2Value(
}
double const p = p_cap / p_b_;
- double const n = 1. / (1. - m_);
- double const p_to_n = std::pow(p, n);
+ double const p_to_n = std::pow(p, n_);
double const S_eff = std::pow(p_to_n + 1., -m_);
double const S = S_eff * S_L_max_ - S_eff * S_L_res_ + S_L_res_;
@@ -129,8 +135,8 @@ PropertyDataType SaturationVanGenuchten::d2Value(
}
double const d2S_eff_dp_cap2 =
- m_ * p_to_n * std::pow(p_to_n + 1., -m_ - 2.) * (p_to_n - m_) /
- (p_cap * p_cap * (m_ - 1.) * (m_ - 1.));
+ m_ * n_ * n_ * p_to_n * S_eff * (p_to_n - m_) /
+ ((p_cap * p_to_n + p_cap) * (p_cap * p_to_n + p_cap));
return d2S_eff_dp_cap2 * (S_L_max_ - S_L_res_);
}
} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchten.h b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchten.h
index dbe3a3a7d18089a5815e628d3f15e627a61e21e7..d084061255122bbaac232bc21a735659e15a5918 100644
--- a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchten.h
+++ b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchten.h
@@ -24,13 +24,12 @@ class Component;
* \f[S_\text{eff}=\frac{S-S_r}{S_{\text{max}}-S_r}.\f]
* Above, \f$S_r\f$ and \f$S_{\text{max}}\f$ are the residual and the maximum
* saturations.
- * The exponent \f$m \in (0,1)\f$ and the pressure scaling parameter \f$p_b\f$
- * (it is equal to \f$\rho g/\alpha\f$ in original publication) are given by the
- * user.
- * The scaling parameter \f$p_b\f$ is given in same units as pressure.
+ * The (pressure) exponent \f$m \in (0,1)\f$ and the pressure scaling parameter
+ * \f$p_b\f$ (it is equal to \f$\rho g/\alpha\f$ in original publication) are
+ * given by the user. The scaling parameter \f$p_b\f$ is given in same units as
+ * pressure.
*
- * In the original work another exponent \f$n\f$ is used, but usually set to
- * \f$n = 1 / (1 - m)\f$, and also in this implementation.
+ * Another (saturaton) exponent \f$n\f$ is usually set to \f$n = 1 / (1 - m)\f$.
*
* The saturation is computed from the capillary pressure as follows:
* \f[S(p_c)=
@@ -49,7 +48,8 @@ public:
SaturationVanGenuchten(std::string name,
double const residual_liquid_saturation,
double const residual_gas_saturation,
- double const exponent,
+ double const pressure_exponent,
+ double const saturation_exponent,
double const p_b);
void checkScale() const override
@@ -83,6 +83,7 @@ private:
double const S_L_res_;
double const S_L_max_;
double const m_;
+ double const n_;
double const p_b_;
};
} // namespace MaterialPropertyLib
diff --git a/Tests/MaterialLib/TestMPLCapillaryPressureVanGenuchten.cpp b/Tests/MaterialLib/TestMPLCapillaryPressureVanGenuchten.cpp
index c5cfb8f6518b4ee9c2efa7f7711ec8fbe6a1dddb..be7876014529568bc9023da9d9d5fd126dc79008 100644
--- a/Tests/MaterialLib/TestMPLCapillaryPressureVanGenuchten.cpp
+++ b/Tests/MaterialLib/TestMPLCapillaryPressureVanGenuchten.cpp
@@ -21,21 +21,20 @@ TEST(MaterialPropertyLib, CapillaryPressureVanGenuchten)
{
double const residual_liquid_saturation = 0.1;
double const residual_gas_saturation = 0.05;
- double const exponent = 0.79;
+ double const pressure_exponent = 0.79;
+ double const saturation_exponent = 1. / (1. - pressure_exponent);
double const p_b = 10000;
double const maximum_capillary_pressure = 20000;
- MPL::Property const& pressure =
- MPL::CapillaryPressureVanGenuchten{"capillary_pressure",
- residual_liquid_saturation,
- residual_gas_saturation,
- exponent,
- p_b,
- maximum_capillary_pressure};
-
- MPL::Property const& saturation =
- MPL::SaturationVanGenuchten{"saturation", residual_liquid_saturation,
- residual_gas_saturation, exponent, p_b};
+ MPL::Property const& pressure = MPL::CapillaryPressureVanGenuchten{
+ "capillary_pressure", residual_liquid_saturation,
+ residual_gas_saturation, pressure_exponent,
+ // TODO? saturation_exponent,
+ p_b, maximum_capillary_pressure};
+
+ MPL::Property const& saturation = MPL::SaturationVanGenuchten{
+ "saturation", residual_liquid_saturation, residual_gas_saturation,
+ pressure_exponent, saturation_exponent, p_b};
MPL::VariableArray variable_array;
ParameterLib::SpatialPosition const pos;
diff --git a/Tests/MaterialLib/TestMPLSaturationVanGenuchten.cpp b/Tests/MaterialLib/TestMPLSaturationVanGenuchten.cpp
index 32b0e4bdd6a184fbf583423ffbb84965787d927e..bba6aa5c14b380ba32cea9dcbb4bd105785f8f3d 100644
--- a/Tests/MaterialLib/TestMPLSaturationVanGenuchten.cpp
+++ b/Tests/MaterialLib/TestMPLSaturationVanGenuchten.cpp
@@ -20,12 +20,13 @@ TEST(MaterialPropertyLib, SaturationVanGenuchten)
{
double const residual_liquid_saturation = 0.1;
double const residual_gas_saturation = 0.05;
- double const exponent = 0.79;
+ double const pressure_exponent = 0.79;
+ double const saturation_exponent = 1. / (1. - pressure_exponent);
double const p_b = 5000;
- MPL::Property const& pressure_saturation =
- MPL::SaturationVanGenuchten{"saturation", residual_liquid_saturation,
- residual_gas_saturation, exponent, p_b};
+ MPL::Property const& pressure_saturation = MPL::SaturationVanGenuchten{
+ "saturation", residual_liquid_saturation, residual_gas_saturation,
+ pressure_exponent, saturation_exponent, p_b};
MPL::VariableArray variable_array;
ParameterLib::SpatialPosition const pos;
diff --git a/Tests/ProcessLib/RichardsMechanics/MicroporosityComputation.cpp b/Tests/ProcessLib/RichardsMechanics/MicroporosityComputation.cpp
index e7b58c5eec00387e72ddd58e76b75fd3ca4de1a1..17061d63e6304b43634a6ae39c8cd4e64202fc81 100644
--- a/Tests/ProcessLib/RichardsMechanics/MicroporosityComputation.cpp
+++ b/Tests/ProcessLib/RichardsMechanics/MicroporosityComputation.cpp
@@ -78,11 +78,12 @@ TEST(RichardsMechanics, computeMicroPorosity)
// Saturation
constexpr double residual_liquid_saturation = 0;
constexpr double residual_gas_saturation = 0;
- constexpr double exponent = 0.4;
+ constexpr double pressure_exponent = 0.4;
+ constexpr double saturation_exponent = 1. / (1. - pressure_exponent);
constexpr double p_b = 300e6;
MPL::Property const& saturation_micro = MPL::SaturationVanGenuchten{
"saturation_micro", residual_liquid_saturation, residual_gas_saturation,
- exponent, p_b};
+ pressure_exponent, saturation_exponent, p_b};
// Swelling
constexpr std::array swelling_pressures{1e7, 1e7, 1e7};