Solver and return-map failures throw cutback_error

include/numsim-materials/materials/drucker_prager_plasticity.h

@@ -13,6 +13,7 @@
 #include "numsim-materials/materials/drucker_prager_yield_function.h"
 #include "numsim-materials/materials/plasticity_utils.h"
 #include "numsim-materials/solvers/local_newton.h"
+#include "numsim-materials/umat/errors.h"
 
 namespace numsim::materials {
 
@@ -104,7 +105,7 @@ class drucker_prager_plasticity final
     // the material, not about roundoff. Clamping it leaves the stress outside
     // the cone with no plastic flow and no error.
     if (smooth.converged && smooth.x < value_type{0})
-      throw std::runtime_error(
+      throw umat::cutback_error(
           "drucker_prager_plasticity: the return map produced a negative "
           "plastic multiplier, which means the hardening modulus is at or "
           "below -(G + K*eta*beta). The softening branch is unstable and this "
@@ -132,7 +133,7 @@ class drucker_prager_plasticity final
         m_yf.needs_apex_return(m_G, dlambda, ts.eval.sig_eq)) {
       // The apex is the remaining branch, and its own Newton can fail too.
       if (!do_apex_return(ts.eval.sig, C_e, kappa_n))
-        throw std::runtime_error(
+        throw umat::cutback_error(
             "drucker_prager_plasticity: both smooth and apex Newton failed");
       return;
     }

include/numsim-materials/materials/j2_kinematic_plasticity.h

@@ -10,6 +10,7 @@
 #include "numsim-materials/core/parameter_helpers.h"
 #include "numsim-materials/materials/plasticity_utils.h"
 #include "numsim-materials/solvers/vector_newton.h"
+#include "numsim-materials/umat/errors.h"
 
 namespace numsim::materials {
 
@@ -189,16 +190,16 @@ class j2_kinematic_plasticity final
 
     solver.solve();
     if (!solver.converged())
-      throw std::runtime_error("j2_kinematic_plasticity '" +
-                               std::string(base::name()) +
-                               "': the coupled return map did not converge");
+      throw umat::cutback_error("j2_kinematic_plasticity '" +
+                                std::string(base::name()) +
+                                "': the coupled return map did not converge");
 
     const auto dl = m_dl.get();
     // Same reasoning as j2_plasticity: the return map is only entered with
     // F_tr > 0, so a negative multiplier means the hardening is at or below
     // -(3G + C) and the local problem has no admissible solution.
     if (dl < value_type{0})
-      throw std::runtime_error(
+      throw umat::cutback_error(
           "j2_kinematic_plasticity '" + std::string(base::name()) +
           "': the return map produced a negative plastic multiplier: the "
           "hardening softens faster than 3G + C can absorb (for gamma = 0 that "

include/numsim-materials/materials/j2_plasticity.h

@@ -11,6 +11,7 @@
 #include "numsim-materials/core/material_ref.h"
 #include "numsim-materials/materials/plasticity_utils.h"
 #include "numsim-materials/solvers/local_newton.h"
+#include "numsim-materials/umat/errors.h"
 
 namespace numsim::materials {
 
@@ -148,7 +149,7 @@ class j2_plasticity final
     };
     const auto sol = m_solver.get().solve(eval);
     if (!sol.converged)
-      throw std::runtime_error(
+      throw umat::cutback_error(
           "j2_plasticity: return-mapping Newton failed to converge");
     // A negative multiplier is not a rounding artefact to be clamped away.
     // The return map is only entered with F_trial > 0, and the root is
@@ -162,7 +163,7 @@ class j2_plasticity final
     // Moderate softening, -G_eff < H' < 0, still gives a positive root and is
     // untouched.
     if (sol.x < value_type{0})
-      throw std::runtime_error(
+      throw umat::cutback_error(
           "j2_plasticity: the return map produced a negative plastic "
           "multiplier, which means the hardening modulus is at or below "
           "-3G. The softening branch is unstable and this local return map "

include/numsim-materials/materials/j2_rk_plasticity.h

@@ -10,6 +10,7 @@
 #include "numsim-materials/materials/yield_functions.h"
 #include "numsim-materials/materials/plasticity_utils.h"
 #include "numsim-materials/solvers/butcher_tableau.h"
+#include "numsim-materials/umat/errors.h"
 
 namespace numsim::materials {
 
@@ -187,7 +188,7 @@ class j2_rk_plasticity final
         // m_N_stage[i] is written on the converged path only; an unconverged stage would otherwise reuse the
         // previous call's flow direction. Throw, as the other return maps do.
         if (!stage_converged)
-          throw std::runtime_error(std::format("j2_rk_plasticity ('{}'): implicit stage {} of {} failed to converge in {} iterations (tolerance {:g})",
+          throw umat::cutback_error(std::format("j2_rk_plasticity ('{}'): implicit stage {} of {} failed to converge in {} iterations (tolerance {:g})",
                                                this->name(), i, s, m_max_iter, m_tol));
       }
     }

include/numsim-materials/materials/neo_hooke.h

@@ -4,6 +4,7 @@
 #include <cmath>
 #include <tmech/tmech.h>
 #include "numsim-materials/core/material_base.h"
+#include "numsim-materials/umat/errors.h"
 
 namespace numsim::materials {
 
@@ -46,7 +47,7 @@ class neo_hooke : public material_base<neo_hooke<Traits>, Traits> {
     const auto& Finv{m_F_inv.get()};
     const value_type J{m_J.get()};
     if (!(J > value_type(0)))
-      throw std::domain_error("neo_hooke: det F must be positive");
+      throw umat::cutback_error("neo_hooke: det F must be positive");
     const value_type lambda{m_K - value_type(2) / value_type(3) * m_G};
     const value_type lnJ{std::log(J)};
     const auto FinvT{tmech::trans(Finv)};

include/numsim-materials/mean_field/kernels.h

@@ -18,6 +18,7 @@
 #include "numsim-materials/core/tensor_identities.h"
 #include "numsim-materials/mean_field/kernel.h"
 #include "numsim-materials/mean_field/shapes.h"
+#include "numsim-materials/umat/errors.h"
 
 namespace numsim::materials::mean_field {
 
@@ -217,7 +218,7 @@ class self_consistent final : public kernel_base<self_consistent<Traits>, Traits
       }
       C = (1 - m_relaxation) * C + m_relaxation * C_new;
     }
-    throw std::runtime_error(std::format("self_consistent: no convergence in {} iterations (residual {:g}, tolerance {:g})",
+    throw umat::cutback_error(std::format("self_consistent: no convergence in {} iterations (residual {:g}, tolerance {:g})",
                                          m_max_iterations, residual, m_tolerance));
   }
 

tests/test_drucker_prager.cpp

@@ -13,6 +13,7 @@
 #include "numsim-materials/solvers/local_newton.h"
 #include "numsim-materials/solvers/butcher_tableau.h"
 #include "numsim-materials/postprocessing/numerical_diff_checker.h"
+#include "numsim-materials/umat/errors.h"
 
 namespace {
 
@@ -96,7 +97,7 @@ struct dp_graph {
   static constexpr T eta{0.1};
   static constexpr T beta{0.05};
 
-  dp_graph(std::vector<double> direction, T increment, T H_mod) {
+  dp_graph(std::vector<double> direction, T increment, T H_mod, T tolerance = T{1e-10}) {
     param_type p;
     p.insert<std::string>("name", "stepper");
     p.insert<T>("increment", increment);
@@ -105,6 +106,7 @@ struct dp_graph {
 
     p.clear();
     p.insert<std::string>("name", "solver");
+    p.insert<T>("tolerance", tolerance);
     ctx.create<numsim::materials::local_newton<policy>>(p);
 
     p.clear();
@@ -718,3 +720,28 @@ TEST(DPBranch, StressIsContinuousAcrossTheBranchBoundary) {
 }
 
 } // namespace
+
+namespace {
+template <typename Fn>
+void expect_cutback(Fn&& fn) {
+  EXPECT_THROW(
+      {
+        for (int i = 0; i < 20; ++i) {
+          fn();
+        }
+      },
+      numsim::materials::umat::cutback_error);
+}
+}  // namespace
+
+// Softening past -(G + K eta beta): no admissible multiplier, but a smaller step is worth a retry.
+TEST(DruckerPragerErrors, ANegativeMultiplierAsksForASmallerStep) {
+  dp_graph g({1.0, 0.0, 0.0, 0.0}, T{0.05}, T{-300.0});
+  expect_cutback([&] { g.ctx.update(); g.ctx.commit(); });
+}
+
+// A zero tolerance fails both the smooth and the apex Newton.
+TEST(DruckerPragerErrors, WhenBothNewtonsFailTheStepIsRetryable) {
+  dp_graph g({1.0, 0.0, 0.0, 0.0}, T{0.05}, T{500.0}, T{0});
+  expect_cutback([&] { g.ctx.update(); g.ctx.commit(); });
+}

tests/test_finite_strain.cpp

@@ -100,6 +100,16 @@ TEST(SaintVenantKirchhoff, TangentMatchesFiniteDifferences) {
   expect_tangent_matches_finite_differences<nm::saint_venant_kirchhoff<policy>>(ctx, src);
 }
 
+// An inverted element is the classic reason for a smaller step.
+TEST(NeoHooke, ANonPositiveDeterminantAsksForASmallerStep) {
+  ctx_type ctx;
+  auto& src = add_model<nm::neo_hooke<policy>>(ctx);
+  tensor2 F{tmech::eye<T, 3, 2>()};
+  F(0, 0) = -1;
+  src.bind(F, F);
+  EXPECT_THROW(ctx.update(), u::cutback_error);
+}
+
 TEST(NeoHooke, TangentMatchesFiniteDifferences) {
   ctx_type ctx;
   auto& src = add_model<nm::neo_hooke<policy>>(ctx);

tests/test_j2_kinematic_plasticity.cpp

@@ -14,6 +14,7 @@
 #include "numsim-materials/materials/j2_kinematic_plasticity.h"
 #include "numsim-materials/materials/linear_isotropic_hardening.h"
 #include "numsim-materials/solvers/vector_newton.h"
+#include "numsim-materials/umat/errors.h"
 #include "numsim-materials/umat/external_state_source.h"
 
 namespace {
@@ -58,12 +59,14 @@ t2 strain_b() {
 
 struct rig {
   static constexpr T K = 166.67, G = 76.92, sigma_0 = 1.0, H_iso = 5.0, C = 30.0;
+  T H_mod;
 
   ctx_type ctx;
   nm::external_strain_source<policy>* src{};
   nm::vector_newton<policy>* solver{};
 
-  explicit rig(T gamma, int max_iter = 50, bool graph_driven = false) {
+  explicit rig(T gamma, int max_iter = 50, bool graph_driven = false, T hardening = H_iso)
+      : H_mod(hardening) {
     param_type p;
     p.insert<std::string>("name", "stepper");
     src = &ctx.create<nm::external_strain_source<policy>>(p);
@@ -83,7 +86,7 @@ struct rig {
     p.clear();
     p.insert<std::string>("name", "hardening");
     p.insert<std::string>("source", "kin");
-    p.insert<T>("K", H_iso);
+    p.insert<T>("K", H_mod);
     ctx.create<nm::linear_isotropic_hardening<policy>>(p);
 
     p.clear();
@@ -254,6 +257,28 @@ TEST(J2Kinematic, RejectsAGraphDrivenSolver) {
   expect_throw_saying([&] { r.step(t2{0.1 * strain_a()}); }, "graph_driven = false");
 }
 
+/// A failed return map is retryable at a smaller step; a misconfigured solver is not.
+TEST(J2Kinematic, ANonConvergedReturnMapAsksForASmallerStep) {
+  rig af(40.0, /*max_iter=*/1);
+  EXPECT_THROW(af.step(strain_a()), nm::umat::cutback_error);
+}
+
+TEST(J2Kinematic, ANegativeMultiplierAsksForASmallerStep) {
+  rig soft(0.0, 50, false, /*H_iso=*/-1000.0);
+  EXPECT_THROW(soft.step(strain_a()), nm::umat::cutback_error);
+}
+
+TEST(J2Kinematic, AGraphDrivenSolverIsAConfigurationErrorNotACutback) {
+  rig r(40.0, 50, /*graph_driven=*/true);
+  try {
+    r.step(t2{0.1 * strain_a()});
+    FAIL() << "a graph-driven solver was accepted";
+  } catch (const nm::umat::cutback_error& e) {
+    FAIL() << "retrying a smaller step cannot fix the configuration: " << e.what();
+  } catch (const std::runtime_error&) {
+  }
+}
+
 /// Jacobian blocks against central differences of the residual, as directional
 /// derivatives along symmetric directions: J applied to a direction must equal
 /// the derivative of R along it. Directional, because the Mandel space the

tests/test_j2_plasticity.cpp

@@ -12,6 +12,7 @@
 #include "numsim-materials/solvers/local_newton.h"
 #include "numsim-materials/solvers/butcher_tableau.h"
 #include "numsim-materials/postprocessing/numerical_diff_checker.h"
+#include "numsim-materials/umat/errors.h"
 
 namespace {
 
@@ -20,6 +21,7 @@ using T = policy::value_type;
 using ctx_type = numsim::materials::material_context<policy>;
 using param_type = policy::ParameterHandler;
 using tensor2 = tmech::tensor<T, 3, 2>;
+using numsim::materials::umat::cutback_error;
 
 class J2PlasticityTest : public ::testing::Test {
 protected:
@@ -247,7 +249,7 @@ TEST_F(J2TangentTest, ConsistentTangentAllSteps) {
 /// Moderate softening, -3G < H' < 0, still has a positive root and must keep
 /// working.
 namespace {
-ctx_type* build_softening(ctx_type& ctx, T hardening) {
+ctx_type* build_softening(ctx_type& ctx, T hardening, T tolerance = T{1e-10}) {
   constexpr T G{76.92};
   param_type p;
   p.insert<std::string>("name", "stepper");
@@ -257,6 +259,7 @@ ctx_type* build_softening(ctx_type& ctx, T hardening) {
 
   p.clear();
   p.insert<std::string>("name", "solver");
+  p.insert<T>("tolerance", tolerance);
   ctx.create<numsim::materials::local_newton<policy>>(p);
 
   p.clear();
@@ -289,7 +292,21 @@ TEST(J2Softening, UnstableSofteningIsRejectedRatherThanClamped) {
           ctx.commit();
         }
       },
-      std::runtime_error);
+      cutback_error);
+}
+
+// A zero tolerance cannot be met, so the return map fails to converge.
+TEST(J2Softening, ANonConvergedReturnMapAsksForASmallerStep) {
+  ctx_type ctx;
+  build_softening(ctx, T{1000.0}, T{0});
+  EXPECT_THROW(
+      {
+        for (int i = 0; i < 10; ++i) {
+          ctx.update();
+          ctx.commit();
+        }
+      },
+      cutback_error);
 }
 
 TEST(J2Softening, ModerateSofteningStillYields) {

tests/test_mean_field_kernels.cpp

@@ -16,6 +16,7 @@
 #include "numsim-materials/io/json_parameter_converter.h"
 #include "numsim-materials/io/object_readers.h"
 #include "numsim-materials/mean_field/kernel_registry.h"
+#include "numsim-materials/umat/errors.h"
 
 namespace {
 
@@ -140,7 +141,7 @@ TEST(MeanFieldKernels, SelfConsistentThrowsWhenItDoesNotConverge) {
   PH p{self_consistent_parameters()};
   p.insert<int>("max_iterations", 2);
   auto k{make("self_consistent", p)};
-  EXPECT_THROW(concentrations(*k, ph), std::runtime_error);
+  EXPECT_THROW(concentrations(*k, ph), nm::umat::cutback_error);
 }
 
 TEST(MeanFieldKernels, SelfConsistentNeedsAShapePerPhase) {

tests/test_rk_convergence.cpp

@@ -9,6 +9,7 @@
 #include "numsim-materials/materials/linear_isotropic_hardening.h"
 #include "numsim-materials/materials/tensor_component_stepper.h"
 #include "numsim-materials/solvers/rk_integrator.h"
+#include "numsim-materials/umat/errors.h"
 
 namespace {
 
@@ -204,7 +205,7 @@ TEST(RkConvergence, AJ2RkPlasticStepConvergesWithEnoughIterations) {
 
 // the stage used to keep the previous call's flow direction and report nothing; the other return maps throw
 TEST(RkConvergence, AJ2RkStageThatRunsOutOfIterationsThrows) {
-  EXPECT_THROW(run_j2_rk(1, "implicit_euler"), std::runtime_error);
+  EXPECT_THROW(run_j2_rk(1, "implicit_euler"), nm::umat::cutback_error);
 }
 
 TEST(RkConvergence, TheJ2RkFailureNamesTheMaterialAndTheStage) {
@@ -224,7 +225,7 @@ TEST(RkConvergence, TheJ2RkFailureNamesTheMaterialAndTheStage) {
 TEST(RkConvergence, AFailedStageDoesNotInfluenceTheNextCall) {
   ctx_type failed_first, fresh;
   build_j2_rk(failed_first, 1, "implicit_euler");
-  EXPECT_THROW(failed_first.update(), std::runtime_error);
+  EXPECT_THROW(failed_first.update(), nm::umat::cutback_error);
   failed_first.revert(); // the host rejects the failed step (#187)
   dynamic_cast<nm::j2_rk_plasticity<policy>*>(failed_first.find("j2"))->template set_parameter<int>("max_iter", 50);
   failed_first.update(); // the stepper is at its second increment now