Update 2 files

- /Tests/MaterialLib/TestMPLSaturationVanGenuchtenWithVolumetricStrain.cpp
- /MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchtenWithVolumetricStrain.cpp

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

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;
+    }
 }