Added new unit tests to check the effective stresses

tests/robertson_cpt_interpretation/test_robertson_cpt_interpretation.py

@@ -1137,3 +1137,138 @@ def test_pwp_level_calc_user_already_inputted_pwp_error(self):
         # call test functions
         with pytest.raises(ValueError):
             interpreter.pwp_level_calc()
+
+    def test_stress_calc_total_stress(self, setup_interpreter):
+        interpreter = RobertsonCptInterpretation()
+        interpreter.data = type('test', (object,), {})()
+        # case that we are bellow the surface
+
+        interpreter.data.depth = np.array([0, 1, 2, 3, 4, 5])
+        interpreter.data.depth_to_reference = np.array([0, 1, 2, 3, 4, 5])
+        interpreter.data.pwp = 2.0
+        interpreter.data.g = 9.81
+        interpreter.data.gamma = np.array([18, 18, 18, 18, 18, 18])
+        interpreter = setup_interpreter
+        interpreter.stress_calc()
+        expected_total_stress = np.array([0, 18, 36, 54, 72, 90])
+        np.testing.assert_array_almost_equal(interpreter.data.total_stress, expected_total_stress)
+
+
+    def test_stress_calc_no_pwp_effect(self):
+        """
+        Test when the computed pore water pressure (pwp) is zero,
+        so that effective stress equals total stress.
+        """
+        instance = RobertsonCptInterpretation()
+        instance.data = type('test', (object,), {})()
+        # Set a simple depth array (in meters)
+        instance.data.depth = np.array([0, 1, 2, 3, 4], dtype=float)
+        # set the pwp to a very low value so that the pwp is negative which will be set to zero
+        instance.data.pwp = -10  # Arbitrary; z_aux will still be 0.
+        instance.data.depth_to_reference = instance.data.depth - 1
+        instance.data.g = 9.81
+        # Use a constant unit weight (gamma) of 18 kN/m3
+        instance.gamma = np.full(instance.data.depth.shape, 18.0)
+
+        instance.stress_calc()
+
+        # Compute expected total stress:
+        # z = diff([0,1,2,3,4]) = [1,1,1,1] then appended becomes [1,1,1,1,1].
+        # Cumulative sum of (18 * 1): [18, 36, 54, 72, 90].
+        expected_total = np.array([18, 36, 54, 72, 90], dtype=float)
+        # Since pwp is zero, effective stress equals total stress.
+        expected_effective = expected_total.copy()
+
+        np.testing.assert_allclose(instance.data.total_stress, expected_total, rtol=1e-5)
+        np.testing.assert_allclose(instance.data.effective_stress, expected_effective, rtol=1e-5)
+
+    def test_stress_calc_positive_pwp(self):
+        """
+        Test a scenario where the computed pwp is positive and reduces effective stress.
+        """
+        instance = RobertsonCptInterpretation()
+        instance.data = type('test', (object,), {})()
+        # Depth starting at 5 m
+        instance.data.depth = np.array([5, 6, 7, 8, 9], dtype=float)
+        # Let depth_to_reference be constant at -2 m (typical for below ground surface)
+        instance.data.depth_to_reference = np.full(instance.data.depth.shape, -2.0)
+        instance.data.pwp = 10
+        instance.data.g = 9.81
+        instance.gamma = np.full(instance.data.depth.shape, 18.0)
+
+        instance.stress_calc()
+
+        # Compute expected total stress:
+        # z = diff([5,6,7,8,9]) = [1,1,1,1] appended to [1,1,1,1,1]
+        # Cumulative sum: [18, 36, 54, 72, 90] then add 5 * mean(gamma) = 5*18 = 90,
+        # so total stress = [108, 126, 144, 162, 180]
+        expected_total = np.array([108, 126, 144, 162, 180], dtype=float)
+        # Compute expected pwp:
+        # z_aux = min(pwp, abs(depth_to_reference[0]) + depth[0]) = min(10, 2 + 5) = 7.
+        # Then, pwp = (7 - (-2)) * 9.81 = 9 * 9.81 ≈ 88.29 (for all depths).
+        pwp_val = 9 * 9.81
+        expected_effective = expected_total - pwp_val
+
+        np.testing.assert_allclose(instance.data.total_stress, expected_total, rtol=1e-5)
+        np.testing.assert_allclose(instance.data.effective_stress, expected_effective, rtol=1e-5)
+
+    def test_stress_calc_negative_pwp_clipped(self):
+        """
+        Test a case where the computed pwp would be negative (indicating suction)
+        but is then clipped to zero.
+        """
+        instance = RobertsonCptInterpretation()
+        instance.data = type('test', (object,), {})()
+        instance.data.depth = np.array([0, 1, 2, 3, 4], dtype=float)
+        # Set depth_to_reference such that its values are higher than z_aux will be.
+        instance.data.depth_to_reference = np.array([2, 3, 3, 3, 3], dtype=float)
+        instance.data.pwp = 10
+        instance.data.g = 9.81
+        instance.gamma = np.full(instance.data.depth.shape, 18.0)
+
+        instance.stress_calc()
+
+        # Expected total stress is as in the first test: [18, 36, 54, 72, 90]
+        expected_total = np.array([18, 36, 54, 72, 90], dtype=float)
+        # Compute expected pwp:
+        # z_aux = min(10, abs(2)+0) = min(10,2) = 2.
+        # pwp = (2 - depth_to_reference)*9.81 = [0, -9.81, -9.81, -9.81, -9.81] then clipped to zeros.
+        # Therefore, effective stress = total stress.
+        expected_effective = expected_total.copy()
+
+        np.testing.assert_allclose(instance.data.total_stress, expected_total, rtol=1e-5)
+        np.testing.assert_allclose(instance.data.effective_stress, expected_effective, rtol=1e-5)
+
+
+    def test_stress_calc_negative_effective_clipped(self):
+        """
+        Test a case where the pore water pressure is so high that the computed effective
+        stress would be negative, and then it is clipped to zero.
+        """
+        instance = RobertsonCptInterpretation()
+        instance.data = type('test', (object,), {})()
+        # Use a shallow depth so that total stress is relatively low.
+        instance.data.depth = np.array([0, 0.5, 1, 1.5, 2], dtype=float)
+        # Let depth_to_reference be -1 m.
+        instance.data.depth_to_reference = np.full(instance.data.depth.shape, -1.0, dtype=float)
+        # Set pwp high relative to the depth: here, pwp=2 leads to z_aux = min(2, 1+0)=1.
+        instance.data.pwp = 2
+        instance.data.g = 9.81
+        instance.gamma = np.full(instance.data.depth.shape, 18.0)
+
+        instance.stress_calc()
+
+        # Compute expected total stress:
+        # z = diff([0, 0.5, 1, 1.5, 2]) = [0.5, 0.5, 0.5, 0.5] then appended: [0.5, 0.5, 0.5, 0.5, 0.5]
+        # Cumulative sum: 18*0.5 = 9, so: [9, 18, 27, 36, 45]
+        expected_total = np.array([9, 18, 27, 36, 45], dtype=float)
+        # Compute pwp:
+        # z_aux = min(2, abs(-1)+0)=min(2,1)=1.
+        # pwp = (1 - (-1))*9.81 = 2*9.81 = 19.62 for every depth.
+        effective = expected_total - 19.62
+        # Clip negative effective stresses to zero.
+        effective[effective < 0] = 0
+        expected_effective = effective
+
+        np.testing.assert_allclose(instance.data.total_stress, expected_total, rtol=1e-5)
+        np.testing.assert_allclose(instance.data.effective_stress, expected_effective, rtol=1e-5)