complete saturation derivative by outer derivative

MaterialLib/MPL/Properties/SaturationBrooksCorey.cpp

@@ -61,6 +61,8 @@ PropertyDataType SaturationBrooksCorey::dValue(
            "SaturationBrooksCorey::dValue is implemented for "
            " derivatives with respect to capillary pressure only.");
 
+    const double s_L_res = _residual_liquid_saturation;
+    const double s_L_max = 1.0 - _residual_gas_saturation;
     const double p_b = _entry_pressure;
     const double p_cap = std::get<double>(
         variable_array[static_cast<int>(Variable::capillary_pressure)]);
@@ -72,8 +74,9 @@ PropertyDataType SaturationBrooksCorey::dValue(
                          .template value<double>(variable_array, pos, t);
 
     const double lambda = _exponent;
+    const double ds_L_d_s_eff = 1. / (s_L_max - s_L_res);
 
-    return -lambda / p_cap * s_L;
+    return -lambda / p_cap * s_L * ds_L_d_s_eff;
 }
 
 PropertyDataType SaturationBrooksCorey::d2Value(

Tests/MaterialLib/TestMPLSaturationBrooksCorey.cpp

@@ -86,10 +86,12 @@ TEST(MaterialPropertyLib, SaturationBrooksCorey)
         const double s_ref =
             s_eff * (max_saturation - ref_residual_liquid_saturation) +
             ref_residual_liquid_saturation;
-        const double ds_dpc =
+        const double ds_eff_dpc =
             (p_cap <= ref_entry_pressure) ? 0. : -ref_lambda / p_cap * s_ref;
+            const double ds_L_ds_eff = 1. / (max_saturation - ref_residual_liquid_saturation);
+            const double ds_L_dpc = ds_L_ds_eff * ds_eff_dpc;
 
         ASSERT_NEAR(s_L, s_ref, 1.e-10);
-        ASSERT_NEAR(ds_L_dp_cap, ds_dpc, 1.e-10);
+        ASSERT_NEAR(ds_L_dp_cap, ds_L_dpc, 1.e-10);
     }
 }
\ No newline at end of file