From d276dab9a758bd254b04da46cdae59c9fa80c4a2 Mon Sep 17 00:00:00 2001
From: eike-radeisen <eike.radeisen@bgr.de>
Date: Mon, 4 Dec 2023 14:26:02 +0000
Subject: [PATCH] Update 2 files
- /Tests/MaterialLib/TestMPLSaturationVanGenuchtenWithVolumetricStrain.cpp
- /MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchtenWithVolumetricStrain.cpp
---
...rationVanGenuchtenWithVolumetricStrain.cpp | 16 ++++----
...rationVanGenuchtenWithVolumetricStrain.cpp | 37 ++++++++++++++++++-
2 files changed, 44 insertions(+), 9 deletions(-)
diff --git a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchtenWithVolumetricStrain.cpp b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchtenWithVolumetricStrain.cpp
index b4c08d67dc5..50f30991355 100644
--- a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchtenWithVolumetricStrain.cpp
+++ b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchtenWithVolumetricStrain.cpp
@@ -71,9 +71,9 @@ PropertyDataType SaturationVanGenuchtenWithVolumetricStrain::value(
double const p_M = p_cap / p_b_M;
double const p_to_n_M = std::pow(p_M, n);
double const S_eff_M = std::pow(p_to_n_M + 1., -m_);
- double const S_eff = (((e_m_ + (a_ * d_e)) * (S_eff_mi)) +
- ((e_0_ - e_m_ - (a_ * d_e)) * (S_eff_M))) /
- e_0_;
+ double const S_eff =
+ (e_m_ * S_eff_mi + ((e_0_ - e_m_ - (a_ * d_e)) * S_eff_M)) /
+ (e_0_ - (a_ * d_e));
double const S = S_eff * S_L_max_ - S_eff * S_L_res_ + S_L_res_;
return std::clamp(S, S_L_res_, S_L_max_);
@@ -108,9 +108,9 @@ PropertyDataType SaturationVanGenuchtenWithVolumetricStrain::dValue(
double const p_M = p_cap / p_b_M;
double const p_to_n_M = std::pow(p_M, n);
double const S_eff_M = std::pow(p_to_n_M + 1., -m_);
- double const S_eff = (((e_m_ + (a_ * d_e)) * (S_eff_mi)) +
- ((e_0_ - e_m_ - (a_ * d_e)) * (S_eff_M))) /
- e_0_;
+ double const S_eff =
+ (e_m_ * S_eff_mi + ((e_0_ - e_m_ - (a_ * d_e)) * (S_eff_M))) /
+ (e_0_- (a_ * d_e));
double const S = S_eff * S_L_max_ - S_eff * S_L_res_ + S_L_res_;
if (S < S_L_res_ || S > S_L_max_)
@@ -119,11 +119,11 @@ PropertyDataType SaturationVanGenuchtenWithVolumetricStrain::dValue(
}
double const dS_eff_dp_cap =
- (-(e_m_ + (a_ * d_e)) * n * m_ * p_to_n *
+ (-(e_m_) * n * m_ * p_to_n *
std::pow(p_to_n + 1., -m_ - 1) -
(e_0_ - e_m_ - (a_ * d_e)) * n * m_ * p_to_n_M *
std::pow(p_to_n_M + 1., -m_ - 1)) /
- (e_0_ * p_cap);
+ ((e_0_- (a_ * d_e)) * p_cap);
return dS_eff_dp_cap * (S_L_max_ - S_L_res_);
}
} // namespace MaterialPropertyLib
diff --git a/Tests/MaterialLib/TestMPLSaturationVanGenuchtenWithVolumetricStrain.cpp b/Tests/MaterialLib/TestMPLSaturationVanGenuchtenWithVolumetricStrain.cpp
index 489212b31b7..0215523343b 100644
--- a/Tests/MaterialLib/TestMPLSaturationVanGenuchtenWithVolumetricStrain.cpp
+++ b/Tests/MaterialLib/TestMPLSaturationVanGenuchtenWithVolumetricStrain.cpp
@@ -70,8 +70,43 @@ TEST(MaterialPropertyLib, SaturationVanGenuchtenWithVolumetricStrain)
double const S = pressure_saturation.template value<double>(
variable_array, pos, t, dt);
- double const S_expected = 0.71943039;
+ double const S_expected = 0.73068816;
ASSERT_NEAR(S, S_expected, 1e-5);
}
+ // Test derivative (dValue) of saturation with respect to capillary pressure
+ {
+ double const vol_s = 0.002;
+ double const p_L = 1000000;
+
+ variable_array.capillary_pressure = p_L;
+ variable_array.volumetric_strain = vol_s;
+
+ // Calculate saturation
+ double const S = pressure_saturation.template value<double>(
+ variable_array, pos, t, dt);
+
+ // Calculate the derivative numerically using a small perturbation
+ double const dP = 1.0e-6;
+ variable_array.capillary_pressure = p_L + dP;
+ double const S_plus_dP = pressure_saturation.template value<double>(
+ variable_array, pos, t, dt);
+
+ double const approximated_dS_dP = (S_plus_dP - S) / dP;
+
+ // Calculate the derivative analytically using dValue
+ double const analytic_dS_dP =
+ pressure_saturation.template dValue<double>(
+ variable_array,
+ MaterialPropertyLib::Variable::capillary_pressure, pos, t, dt);
+
+ // Check if the numerical and analytical derivatives are close
+ ASSERT_LE(std::fabs(approximated_dS_dP - analytic_dS_dP), 1e-6)
+ << "for expected derivative of saturation with respect to "
+ "capillary pressure "
+ << approximated_dS_dP
+ << " and for computed derivative of saturation with respect to "
+ "capillary pressure "
+ << analytic_dS_dP;
+ }
}
--
GitLab