From e3207cba331f304975e92ee186386ec0e0bc9eb4 Mon Sep 17 00:00:00 2001
From: Richard Faasse <56549273+rfaasse@users.noreply.github.com>
Date: Fri, 10 Jan 2025 15:59:50 +0100
Subject: [PATCH] [GeoMechanicsApplication] Create exit criterion for c phi
 process for too small reduction increments (#12980)

---
 .../apply_c_phi_reduction_process.cpp         |  2 ++
 .../test_apply_c_phi_reduction_process.cpp    | 29 +++++++++++++++++++
 2 files changed, 31 insertions(+)

diff --git a/applications/GeoMechanicsApplication/custom_processes/apply_c_phi_reduction_process.cpp b/applications/GeoMechanicsApplication/custom_processes/apply_c_phi_reduction_process.cpp
index f9b7b3e4304..8931010ff5a 100644
--- a/applications/GeoMechanicsApplication/custom_processes/apply_c_phi_reduction_process.cpp
+++ b/applications/GeoMechanicsApplication/custom_processes/apply_c_phi_reduction_process.cpp
@@ -28,6 +28,8 @@ void ApplyCPhiReductionProcess::ExecuteInitializeSolutionStep()
     mReductionFactor = mPreviousReductionFactor - mReductionIncrement;
     KRATOS_ERROR_IF(mReductionFactor <= 0.01)
         << "Reduction factor should not drop below 0.01, calculation stopped." << std::endl;
+    KRATOS_ERROR_IF(mReductionIncrement <= 0.001)
+    << "Reduction increment should not drop below 0.001, calculation stopped. Final safety factor = " << 1.0 / mPreviousReductionFactor << std::endl;
     KRATOS_INFO("ApplyCPhiReductionProces::ExecuteInitializeSolutionStep")
         << "Try a c-phi reduction factor " << mReductionFactor << " (safety factor "
         << 1. / mReductionFactor << ") Previous reduction = " << mPreviousReductionFactor
diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/custom_processes/test_apply_c_phi_reduction_process.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_processes/test_apply_c_phi_reduction_process.cpp
index d4cd526bd6b..405b0c48dc8 100644
--- a/applications/GeoMechanicsApplication/tests/cpp_tests/custom_processes/test_apply_c_phi_reduction_process.cpp
+++ b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_processes/test_apply_c_phi_reduction_process.cpp
@@ -185,6 +185,14 @@ KRATOS_TEST_CASE_IN_SUITE(CheckFailureEmptyModelPartApplyCPhiReductionProcess, K
                                       "ApplyCPhiReductionProces has no elements in modelpart dummy")
 }
 
+KRATOS_TEST_CASE_IN_SUITE(CheckReturnsZeroForValidModelPartApplyCPhiReductionProcess, KratosGeoMechanicsFastSuite)
+{
+    Model model;
+    auto& r_model_part = PrepareCPhiTestModelPart(model);
+    ApplyCPhiReductionProcess process{r_model_part, {}};
+    KRATOS_CHECK_EQUAL(process.Check(), 0);
+}
+
 KRATOS_TEST_CASE_IN_SUITE(CheckFailureNegativeReductionFactorApplyCPhiReductionProcess, KratosGeoMechanicsFastSuite)
 {
     Model model;
@@ -200,4 +208,25 @@ KRATOS_TEST_CASE_IN_SUITE(CheckFailureNegativeReductionFactorApplyCPhiReductionP
         "Reduction factor should not drop below 0.01, calculation stopped.");
 }
 
+KRATOS_TEST_CASE_IN_SUITE(CheckFailureTooSmallReductionIncrementApplyCPhiReductionProcess, KratosGeoMechanicsFastSuite)
+{
+    Model model;
+    auto& r_model_part = PrepareCPhiTestModelPart(model);
+    // Number of cycles should be larger than one to trigger the step halving of the
+    // reduction increment
+    r_model_part.GetProcessInfo().GetValue(NUMBER_OF_CYCLES) = 10;
+    ApplyCPhiReductionProcess process{r_model_part, {}};
+    for (size_t i = 0; i < 6; ++i) {
+        process.ExecuteInitializeSolutionStep();
+        process.ExecuteFinalizeSolutionStep();
+    }
+
+    // After halving the initial reduction increment (0.1) for the seventh time, it becomes
+    // 0.00078125, so it should throw an exception
+    // The final safety factor can be calculated as (1.0 / (1 - 0.1 * 0.5 - 0.1 * 0.5^2 ... + 0.1*0.5^6)) = 1.10919
+    KRATOS_EXPECT_EXCEPTION_IS_THROWN(
+        process.ExecuteInitializeSolutionStep(),
+        "Reduction increment should not drop below 0.001, calculation stopped. Final safety factor = 1.10919");
+}
+
 } // namespace Kratos::Testing