Repository navigation
Feature/cpp return mapping - #131
Merged
Merged
Conversation
Introduce a closest-point projection (CPP) return-mapping implementation and wire its exact consistent tangent into the UMAT framework. Add new API and implementation: include/simcoon/Continuum_mechanics/Umat/return_mapping.hpp and src/Continuum_mechanics/Umat/return_mapping.cpp (closest_point_return_mapping, ReturnMappingResult, cpp_consistent_tangent, and caller hooks). Extend strain and plasticity mechanisms with ClosestPointIngredients and closest_point_ingredients callbacks and add backward_euler_factors to kinematic hardening classes to provide linearisation at the refreshed state. Update tangent_assembly, parameter/ state handling and ModularUMAT so tangent_mode 3 (tangent_closest_point) is accepted; when a kernel does not provide a CPP operator it falls back to the algorithmic operator. Add tests and Python bindings to accept and exercise tangent_mode 3, update CMakeLists to build the new source, and refresh docs to describe the new closest-point integrator and its tangent behavior. Several hardening implementations (Prager/Armstrong–Frederick/Chaboche) now provide backward-Euler factors used by the CPP solver.
Implement a closest-point (tangent_mode == 3) return-mapping branch and consistent tangent support for the modular UMAT. Key changes: - Modular UMAT: wire a new return_mapping_cpp() path selected when ndi==3, elasticity has constant stiffness and all multiplier-carrying mechanisms support closest-point; return_mapping() now returns bool and reports convergence; store closest-point result in cpp_result_ and use it in compute_tangent() to produce the exact Jacobian. Adjust run() to reject on closest-point non-convergence and update drift-guard logic. - New return_mapping_cpp() implementation: builds per-row closest-point ingredients, runs closest_point_return_mapping(), applies inert (non-multiplier) rows afterwards, commits a self-consistent stress and keeps the converged solve for the tangent. - Mechanism API changes (strain/plasticity/viscoelastic/damage): add carries_multipliers() and supports_closest_point() where appropriate; refresh_state/closest_point_ingredients contract updated to support the closest-point integrator and to throw on default implementations. Damage and viscoelastic mechanisms marked non-multiplier-carrying. - Plasticity: implement supports_closest_point(), rework refresh_state for CPP to derive state from total multipliers and couple n<->X via Newton. Add constitutive include. - Tangent assembly: compute_tangent() handles tangent_closest_point specially (uses cpp_consistent_tangent when available) and degrades to algorithmic operator when closest-point did not run. - Numerical solver: add Fischer_Burmeister_residual declaration to num_solve.hpp for line-search/ratings of trial points. - Python bindings/tests: update python module docs for tangent_mode constants and add extensive test_modular_closest_point.py exercising mode 3 behavior and comparisons with cutting-plane modes. These changes enable solving stress, multipliers and backward-Euler state simultaneously for mechanisms with flow Hessians, producing an exact discrete-update Jacobian and improved behavior for cases where the closest-point formulation applies.
Add examples/mechanical/MODUL_closest_point.py and update the examples README to document the closest-point (tangent_mode=3) demo. Implement several UMAT improvements to support the closest-point integrator and reduce work in the common small-strain / elastic cases: - modular_umat.hpp: document admissible-trial behaviour and add Y_crit parameter to the CPP entry. - modular_umat.cpp: evaluate trial constraints up-front and commit inert rows for an admissible trial; pass Y_crit to return_mapping_cpp; avoid rotating mechanisms when DR is identity; integrate an elastic-trial short-circuit and ensure constraint re-evaluation between cutting-plane iterations; skip Hessian assembly when no active multipliers. - internal_variable_collection.cpp: skip quaternion rotation when DR == I to avoid unnecessary work. - plasticity_mechanism.cpp: short-circuit updates and energy/work computation when plastic increments are zero. These changes improve correctness and performance for the closest-point projection, avoid redundant operations in the small-strain/elastic guard, and reduce unnecessary allocations/updates for zero updates.
Add support for a nonlinear (hyperelastic) elastic block in the closest-point return-mapping: introduce ReturnStateHooks::elastic_response and eps_el_tr so the solver can evaluate the elastic response and its tangent at every iterate. Change flow_state_coupling semantics to return strain-typed dLambda (the helper multiplies by the current elastic tangent), use a local Lc (the block tangent at the iterate) when assembling condensation and kappa, and remove the previous constant-stiffness short-circuit so mode-3 can dispatch a nonlinear block. Update docs/comments about the elastic block and robustness of the semi-smooth Newton. Wire the hooks in ModularUMAT when elasticity is non-constant and add a Python test (test_closest_point_with_hyperelastic_block_is_exact) to verify exactness of the composite tangent with a Neo-Hookean block.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closest-point projection (
tangent_mode = 3): thereturn_mappinghelper and the modular UMAT on itReplaces the cutting-plane local loop by a closest-point projection (CPP), law by law, so that the
Simo-Hughes operator becomes the exact Jacobian of the discrete update. Today it is exact for J2 +
isotropic hardening only; with kinematic hardening or a Hill/DFA/Drucker criterion the flow direction
rotates between iterates, the cutting-plane update is path dependent and the mode-2 operator is
1e-4..1e-2 off. Fix the integrator, never the tangent. Modes 0/1/2 are untouched (bitwise).
This branch carries the first two steps of the series (the helper, then the modular UMAT); the dedicated
kernels (EPICP, EPCHA, thermomechanical twins, SMA unified_T / unified_TR) follow on the same helper.
1. The helper and
tangent_mode = 3accepted everywhereinclude/simcoon/Continuum_mechanics/Umat/return_mapping.hpp+src/.../return_mapping.cpp:closest_point_return_mapping(sigma_tr, L, mechanisms, hooks, Y_crit, control)— the fully implicitsystem {σ − σ_tr + Σ Δs^j L Λ^j(σ, V̂) = 0 ; Δs ≥ 0, Φ ≤ 0, Δs·Φ = 0} with the internal state resolved
by backward Euler at every iterate, condensed through M = I + L Σ Δs^j dΛ^j/dσ to the N×N multiplier
system handed to the EXISTING
Fischer_Burmeister_m;cpp_consistent_tangent(result, L)=assemble_algorithmic_tangenton the converged ingredients (exact, since the CPP map IS the implicitupdate the algorithmic tangent linearises). Elastic guard returns the trial bit-identically.
Ported from
feature/algorithmic_tangent(35d0dbd) with the robustness policy reduced to what therules allow: Δs ≥ 0 projection, backtracking on the merit (
Fischer_Burmeister_residual, new innum_solve: the error loop ofFischer_Burmeister_mverbatim, without its N×N solve, + ‖R_σ‖/σ_ref),deactivation snap after 20 iterations; the Newton-step cap, the loose acceptance and the env-gated
fprintf are gone. A trial point whose inner state solve fails or is not finite is a rejected line-search
step (halved); a LAPACK failure inside
Fischer_Burmeister_mis caught; non-convergence →converged = false→ the caller's tnew_dt = 0.5 step cut, never an exception out of a UMAT.tangent_mode = 3accepted by the solver,sim.umat,sim.umat_T. A kernel not wired on the helper(EPICP, EPCHA, SMA,
_Tfor now) returns the mode-2 operator (compute_tangent_operator) — documentedin solver.rst / umat_catalog.rst / parameter.hpp.
Treturn_mapping(J2 closed-form radial return 1e-10; elastic guard bit-identical; inactive-rowmasking; Hill anisotropic Q-quadratic convergence; Hill exact consistent tangent vs FD + symmetric
discrete map),
Tcpp_ccp_robustness(6 cases × 7 step scales up to 100× yield: 0 failures for eitherintegrator, mean iterations CCP 6.07 / CPP 6.69 — identical to the July figures, so the dropped guards
cost nothing; CCP↔CPP stress difference O(Δε²) under step halving). Python:
test_tangent_modes_converge_to_same_state[closest_point],test_closest_point_mode_is_accepted_everywhere.2. The modular UMAT integrates with it
StrainMechanism:ClosestPointIngredients(total derivatives of a constraint row at the refreshedstate: Φ, ñ = dΦ/dσ, Λ, D = dΛ/dσ, dΛ/dΔs, K = dΦ/dΔs); static capabilities
carries_multipliers()(default true; viscoelastic and damage false) and
supports_closest_point()(Plasticity =has_flow_hessian());refresh_statebuilds the backward-Euler state AND its linearisation in one pass,closest_point_ingredients()returns it (the baserefresh_statethrows if a mechanism claims supportwithout implementing it). The partial slots the cutting-plane loop uses are untouched.
KinematicHardening::backward_euler_factors(Prager / Armstrong-Frederick / Chaboche): the backstress isaffine in the normal through the backward-Euler branch maps, X = X₀(Δp) + γ T n, with
γ = Δp Σ (2/3)C_i/(1+D_iΔp), β = Σ (2/3)C_i (n − D_i α_i)/(1+D_iΔp).
PlasticityMechanism::refresh_state: the n ↔ X coupling is a 6×6 Newton on A = I + γ H T (zero steps forvon Mises from the relaxed-backstress normal, a few for Hill/DFA/Ani) instead of the former fixed point
(non-contractive for stiff recall, ratio ≈ Δp ΣC/σ_eq); then D = A⁻¹H, ñ = n − γ Dᵀ T n,
K = −n·T A⁻¹β − R'(p), dn/dΔp = −A⁻¹ H T β (analytic; FD only in the tests). Limits in the class note
(γ = 0 → partial forms; Prager → ñ = n but D ≠ H, so the partial Hessian was inexact even there).
ModularUMAT::return_mapping_cpp: adapter over the helper — mechanisms' rows become its callbacks, inertrows (viscoelastic, damage) are evaluated once at the converged effective stress (strain equivalence),
an admissible trial returns before any callback is built, the stress is committed from the refreshed
state;
compute_tangentmode 3 =cpp_consistent_tangenton the converged solve (cpp_result_).Non-convergence → the existing reject block (tnew_dt = 0.5), no commit-at-maxiter, drift guard moot.
ElasticityModule::has_constant_stiffness()gates the branch (the helper freezes L).cutting-plane loop + mode-2 operator.
Measured
Mode-3 Lt vs central FD of the increment map (relative), mode 2 in parentheses:
J2+Voce 6e-11 (7e-11, already exact), J2+Prager 5e-11 (1.2e-4), J2+AF 9e-11 (2.4e-4),
J2+Chaboche-2 8e-11 (1.7e-4), Hill+Voce 6e-11 (1.6e-3), Hill+AF 8e-11 (1.1e-3).
Exact operator symmetric for isotropic/Prager (1e-16), NOT for AF/Chaboche (6e-6..5e-5 relative
asymmetry, as predicted: dynamic recovery is not generalised-standard). CCP↔CPP one-increment stress gap
O(h²) (ratios > 3 under halving); over a path both are first order (gap halves). Two von Mises rows under
a large reversed reload: every increment at full size with mode 3.
Speed (C++
ModularUMAT::run, Release, per increment, best of 3×20000; mode 2 → mode 3):Plastic increments pay the inner-consistent scheme — a backward-Euler state solve and a 6×6 condensation per iterate against the cutting-plane loop's scalar update: 2–4× is intrinsic (the plan's "+20–60 %" was wrong); the review pass took it from 2.8–5.9× to the figures above (single trial evaluation, elastic guard before any callback, FB residual without a solve, M⁻¹ products reused, fast inner solve). Measured and rejected: warm-starting the inner Newton from the previous normal (Hill −7 %, J2 +40 %: the relaxed-backstress start is exact for von Mises).
Side gain for modes 1/2 (bitwise-neutral, verified on the regression matrix): on an admissible trial state the cutting-plane loop no longer assembles the Jacobian and solves the FB system for a zero step, and
compute_tangentno longer runs the masked algorithmic assembly — elastic increments 2.2–3.0 us → 0.6–1.1 us.sim.solver.solve, 2×50 (2×10) increments: strain control 1.5×/1.9×/1.9× (1.4×/1.4×/1.5×) for J2+Voce/J2+AF/Hill+AF; stress control 1.1×/1.2×/1.5× (1.1×/1.3×/1.7×) — the exact operator saves global Newton iterations under stress control without paying for the per-point cost on a single material point. Modes 0/1/2: timing unchanged (≤ 1.5 %, noise), bitwise identical.Verification
Treturn_mapping5/5,Tcpp_ccp_robustness2/2,ModularUMATClosestPoint.NonConvergenceRequestsStepCut(maxiter = 1: tnew_dt = 0.5, statev / sigma untouched, Lt elastic; default budget converges); the AF
refresh_statetest now exercises the Newton; 33/33 suites.test_modular_closest_point.py(14 tests: exact tangent ×6 + symmetry classification,radial-return equivalence with mode 2, O(h²) gap, elastic increment bitwise across modes 1/2/3, composites
with viscoelasticity/damage exact at mode 3, Tresca degradation == mode 2 bitwise, two-mechanism reload
without step cut) + 2 in
test_solver_run.py; test_core 992 passed, 1 skipped.finite_strain_framework/scripts/regression_baseline.py --compare: 100/100 IDENTICAL (modes 0/1/2 untouched).