diff --git a/include/numsim-materials/materials/drucker_prager_plasticity.h b/include/numsim-materials/materials/drucker_prager_plasticity.h index 819da4f..e0b6fcf 100644 --- a/include/numsim-materials/materials/drucker_prager_plasticity.h +++ b/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; } diff --git a/include/numsim-materials/materials/j2_kinematic_plasticity.h b/include/numsim-materials/materials/j2_kinematic_plasticity.h index cf87252..a3d605e 100644 --- a/include/numsim-materials/materials/j2_kinematic_plasticity.h +++ b/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 " diff --git a/include/numsim-materials/materials/j2_plasticity.h b/include/numsim-materials/materials/j2_plasticity.h index 7b8a7dd..21ed115 100644 --- a/include/numsim-materials/materials/j2_plasticity.h +++ b/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 " diff --git a/include/numsim-materials/materials/j2_rk_plasticity.h b/include/numsim-materials/materials/j2_rk_plasticity.h index 419bdad..00d7ec7 100644 --- a/include/numsim-materials/materials/j2_rk_plasticity.h +++ b/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)); } } diff --git a/include/numsim-materials/materials/neo_hooke.h b/include/numsim-materials/materials/neo_hooke.h index 58646e0..554a5c0 100644 --- a/include/numsim-materials/materials/neo_hooke.h +++ b/include/numsim-materials/materials/neo_hooke.h @@ -4,6 +4,7 @@ #include #include #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, 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)}; diff --git a/include/numsim-materials/mean_field/kernels.h b/include/numsim-materials/mean_field/kernels.h index ee11c1d..a7e9eaf 100644 --- a/include/numsim-materials/mean_field/kernels.h +++ b/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, 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)); } diff --git a/tests/test_drucker_prager.cpp b/tests/test_drucker_prager.cpp index e336262..09f80ad 100644 --- a/tests/test_drucker_prager.cpp +++ b/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 direction, T increment, T H_mod) { + dp_graph(std::vector direction, T increment, T H_mod, T tolerance = T{1e-10}) { param_type p; p.insert("name", "stepper"); p.insert("increment", increment); @@ -105,6 +106,7 @@ struct dp_graph { p.clear(); p.insert("name", "solver"); + p.insert("tolerance", tolerance); ctx.create>(p); p.clear(); @@ -718,3 +720,28 @@ TEST(DPBranch, StressIsContinuousAcrossTheBranchBoundary) { } } // namespace + +namespace { +template +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(); }); +} diff --git a/tests/test_finite_strain.cpp b/tests/test_finite_strain.cpp index 2a4e925..2318750 100644 --- a/tests/test_finite_strain.cpp +++ b/tests/test_finite_strain.cpp @@ -100,6 +100,16 @@ TEST(SaintVenantKirchhoff, TangentMatchesFiniteDifferences) { expect_tangent_matches_finite_differences>(ctx, src); } +// An inverted element is the classic reason for a smaller step. +TEST(NeoHooke, ANonPositiveDeterminantAsksForASmallerStep) { + ctx_type ctx; + auto& src = add_model>(ctx); + tensor2 F{tmech::eye()}; + 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>(ctx); diff --git a/tests/test_j2_kinematic_plasticity.cpp b/tests/test_j2_kinematic_plasticity.cpp index e57e648..26b7560 100644 --- a/tests/test_j2_kinematic_plasticity.cpp +++ b/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* src{}; nm::vector_newton* 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("name", "stepper"); src = &ctx.create>(p); @@ -83,7 +86,7 @@ struct rig { p.clear(); p.insert("name", "hardening"); p.insert("source", "kin"); - p.insert("K", H_iso); + p.insert("K", H_mod); ctx.create>(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 diff --git a/tests/test_j2_plasticity.cpp b/tests/test_j2_plasticity.cpp index cfab3bd..8424371 100644 --- a/tests/test_j2_plasticity.cpp +++ b/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; using param_type = policy::ParameterHandler; using tensor2 = tmech::tensor; +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("name", "stepper"); @@ -257,6 +259,7 @@ ctx_type* build_softening(ctx_type& ctx, T hardening) { p.clear(); p.insert("name", "solver"); + p.insert("tolerance", tolerance); ctx.create>(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) { diff --git a/tests/test_mean_field_kernels.cpp b/tests/test_mean_field_kernels.cpp index 087ca5f..08206c7 100644 --- a/tests/test_mean_field_kernels.cpp +++ b/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("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) { diff --git a/tests/test_rk_convergence.cpp b/tests/test_rk_convergence.cpp index 5b5a8ec..10a3a20 100644 --- a/tests/test_rk_convergence.cpp +++ b/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*>(failed_first.find("j2"))->template set_parameter("max_iter", 50); failed_first.update(); // the stepper is at its second increment now