Add Drucker-Prager plasticity with apex return - #10
Conversation
- small_strain_plasticity<Traits, YieldFunction>: generic return mapping - Consistent tangent via implicit function theorem (generic, not per yield fn) - Internal Newton via solver.solve(eval) — solver is an external material - j2_yield_function: stateless policy with 6 static methods - linear_isotropic_hardening, exponential_isotropic_hardening - material_ref<T>: lazy material references resolved at finalize() - add_material_ref<T>(name) in material_interface - wire_materials() called before wire_inputs() in finalize() - backward_euler: dual-mode (graph-driven update + direct solve call) - 6 new tests (elastic, yielding, yield surface, deviatoric, hardening, tangent) - All 23 tests pass
- material_ref: debug assert on null dereference - material_interface: defaulted destructor, batch missing-material errors - object_store: string_view in find() - property_engine: dump() to stderr - json_parameter_converter: warnings to stderr - material_context: consistent include guard - test_j2: tighten yield surface tolerance from 10 MPa to 1 MPa - test_property_graph: add missing-parameter and circular-dependency error tests - 25/25 tests pass
- butcher_tableau: runtime data + 7 factory functions
- explicit_rk_integrator: no solver needed
- dirk_integrator: Newton per diagonal stage
- implicit_rk_integrator: coupled Newton (Gaussian elimination)
- curing_rate: pure rate function for autocatalytic curing (RK-compatible)
- autocatalytic_reaction: added rate/rate_derivative outputs
- 12 new tests: convergence order (exponential decay) + curing simulation
- Forward Euler order 1, RK4 order 4, implicit midpoint order 2,
Crank-Nicolson order 2, Gauss-Legendre order 4
- Curing converges with RK4, implicit midpoint, and Gauss-Legendre
- All 37 tests pass
…leau) - rk_integrator: single class handles explicit/DIRK/fully implicit (dispatches at construction based on tableau structure) - small_strain_plasticity: optional tableau parameter for multi-stage return mapping. Without tableau → solver.solve() (classical). With tableau → multi-stage RK stages. - Delete: explicit_rk_integrator, dirk_integrator, implicit_rk_integrator, rk_plasticity, plasticity_integrator, j2_constitutive_law - solver_source now optional (not needed with tableau) - 40/40 tests pass
- plasticity_utils.h: free functions compute_trial() and compute_tangent() - small_strain_plasticity: lean single-stage, no tableau overhead - rk_plasticity: multi-stage with tableau, uses same utils - No dead member storage in the simple class - 40/40 tests pass
…ial return is exact)
petlenz
left a comment
petlenz
left a comment
There was a problem hiding this comment.
Review from the numsim-codegen side
Reviewing with the codegen hat on: numsim-codegen's NumSimMaterialTarget generates constitutive materials that target exactly this runtime's solver/property contract (<function>::rate/rate_derivative and residual/jacobian over Local edges, driven by update_source()). So this solver layer is shared infrastructure and changes ripple into generated code. Overall this is careful, well-structured plasticity — findings are mostly about the generic solver picking up domain-specific behavior, plus one tangent-consistency question.
Done well
- Scalar return mapping in full-tensor representation. The local unknown is the scalar multiplier
Δλ; the tensor update is closed-form (ε_p += Δλ·N). That's why a scalar solver suffices and there is no Voigt/Mandel flattening. (So the Mandel + vector-solver gaps I filed — #11/#12 — are for a different class, coupled tensor-valued local systems, not this. This PR correctly sidesteps them.) - Convergence-failure handling is robust in
small_strain_plasticity::compute()— checksconverged(), falls back to apex, and throws if both fail rather than using a bad iterate. Better rigor than codegen's current in-function Newton (numsim-codegen#85). - The conservative apex pre-check (
dλ_max = F_trial/G_eff) to skip a doomed smooth Newton is a nice optimization.
Findings
HIGH — backward_euler is a generic name with domain-specific behavior baked in. (solvers/backward_euler.h)
update()ends withm_delta = std::abs(m_delta);— "curing degree can only increase" (:80). A curing-specific monotonicity assumption hardcoded into a generic solver.solve()clamps every result withstd::max(x, value_type{0})— "negative plastic-multiplier unphysical" (:93,94,98). A plasticity-specific assumption, also in the generic solver.
So "backward_euler" is really two domain-specific solvers wearing a generic name. This blocks codegen reuse: a generated material whose state can legitimately decrease or go negative would be silently corrupted by abs()/max(·,0). Suggest hoisting the sign/clamp policy out of the solver (a clamp_nonnegative parameter defaulting off, or keep the clamp in the plasticity caller).
MEDIUM — update() and solve() disagree on convergence reporting. solve() sets m_converged; update() never does — on max_iter exhaustion it silently proceeds with the last iterate. The graph-driven path (curing, and any codegen rate material wired through it) has no non-convergence signal. Mirror converged() into update().
MEDIUM — Drucker-Prager algorithmic tangent uses the trial flow direction. materials/small_strain_plasticity.h:168-169. For DP the converged flow direction differs from trial, so N_trial gives a tangent that is not the consistent algorithmic tangent — degrading host-FE Newton from quadratic toward linear near yielding. The material advertises a tangent output; is the trial-direction tangent an intentional documented approximation, or should DP assemble at the converged direction?
LOW — magic constants in the solver. The m_delta = 5e-12 seed (:59), the 1e-30 singular-Jacobian threshold (:94), and the 5-iteration NaN-damping cap should be named/parameterized.
Cross-repo note
The residual/jacobian (and rate/rate_derivative) Local-edge contract itself is clean and is what codegen emits against. The one thing currently preventing a generated material from reusing backward_euler as-is is the HIGH finding (the baked-in sign policy). Hoist that out and this solver becomes directly consumable by codegen output — the goal of the graph-coupled architecture.
| auto x = x0; | ||
| for (int i = 0; i < m_max_iter; ++i) { | ||
| auto [r, dr] = eval(x); | ||
| if (std::abs(r) < m_tol) { m_converged = true; return std::max(x, value_type{0}); } |
There was a problem hiding this comment.
HIGH (codegen-relevant): std::max(x, 0) bakes a plasticity assumption (Δλ ≥ 0) into a generically-named solver. A codegen-generated material solving a general residual whose root is legitimately negative would be silently clamped to 0. Suggest moving the clamp into the plasticity caller, or gating it behind a clamp_nonnegative parameter (default off).
| m_stress = tmech::dcontract(C_e, m_strain.get() - m_eps_p.new_value()); | ||
|
|
||
| // Tangent at trial state (return mapping uses N_trial). | ||
| // For J2, trial = converged. For DP, they differ. |
There was a problem hiding this comment.
The trial-direction tangent is not the consistent algorithmic tangent for DP (N differs trial vs converged), which degrades host-FE Newton convergence near yielding. Intentional documented approximation, or should the DP tangent be assembled at the converged flow direction? The material advertises a tangent output, so worth pinning down.
The pull_request trigger filtered on branches: [main], so a PR targeting a feature branch matched nothing. Work here lands through stacks -- #27..#31 all target a feature branch -- and not one of them has ever been built or tested by CI. The last run of any kind was main, three weeks ago. Dropping the filter runs the job for every pull request whatever its base. The push trigger keeps its main filter, so branch pushes add no load: a stacked branch is covered by its own PR.
EIGEN_BUILD_TESTING is only honoured by Eigen after 3.4.0; the pinned tag gates on BUILD_TESTING, so the existing guard did nothing. Any build that FETCHES Eigen -- every build on a machine without it installed, i.e. CI -- got Eigen's whole suite in its ctest run. Measured here: 987 tests, 835 of them failing, against 46 of our own. Both variables are set so a bumped tag stays covered. Our tests are unaffected: they register through enable_testing(), not BUILD_TESTING. Reproduced with -DCMAKE_DISABLE_FIND_PACKAGE_Eigen3=ON, which is what a clean runner does. 46/46 after.
The previous fix forced BUILD_TESTING OFF to stop a fetched Eigen registering
its ~780 tests into our ctest run. That worked, but BUILD_TESTING is a GLOBAL
variable, and it only left our own tests standing because they register through
a bare enable_testing(). Anyone switching this project to include(CTest) would
have silently dropped the entire suite.
SOURCE_SUBDIR names a directory with no CMakeLists.txt, so MakeAvailable
populates Eigen and stops -- add_subdirectory is never called and none of
Eigen's CMake runs. Eigen is header-only, so the source dir is the include
path and nothing is lost. It is the treatment tmech and nlohmann already get
here: take the headers, leave the build system alone.
Verified with -DCMAKE_DISABLE_FIND_PACKAGE_Eigen3=ON:
- 46/46, ours alone
- 46/46 again with -DBUILD_TESTING=ON, so the coupling is gone rather than
merely satisfied
- _deps/eigen-build contains 0 files
A raw ${eigen_SOURCE_DIR} in the INTERFACE_INCLUDE_DIRECTORIES of an EXPORTED
target is rejected at generate time:
Target "numsim-materials" INTERFACE_INCLUDE_DIRECTORIES property contains
path: .../build/_deps/eigen-src which is prefixed in the build directory.
BUILD_INTERFACE scopes it to the build tree, which is all it can describe -- a
consumer of the INSTALLED package supplies its own Eigen, as it already does
for tmech.
Missed locally because the check piped configure to /dev/null and relied on
&&, and CMake still writes usable build files after a generate error: the
build and all 46 tests ran green on top of a failed configure.
The project called enable_testing() directly, so BUILD_TESTING -- the switch consumers expect to reach -- did not exist. include(CTest) declares it and calls enable_testing() itself. Gated on PROJECT_IS_TOP_LEVEL: embedded in a superproject, that project owns the dashboard targets, and its BUILD_TESTING choice should reach us rather than be re-declared. Embedded without one, BUILD_TESTING defaults ON so behaviour is unchanged for existing consumers. Both switches are kept and either turns tests off: BUILD_TESTING is how a superproject silences every subproject at once, NUMSIM_BUILD_TESTS only ours. This was NOT adoptable before the previous commit. Suppressing a fetched Eigen's ~780 tests meant forcing BUILD_TESTING OFF globally, which under include(CTest) would have silenced our own suite as well. Eigen's CMake no longer runs at all, so the name is free. Nothing else gates on it -- numsim-core uses BUILD_TESTS, tmech TMECH_BUILD_TESTS, nlohmann JSON_BuildTests, all forced off by name. Verified: default 46/46; BUILD_TESTING=OFF and NUMSIM_BUILD_TESTS=OFF each skip the suite; BUILD_TESTING=ON gives 46 again rather than Eigen's.
petlenz
left a comment
petlenz
left a comment
There was a problem hiding this comment.
Critical review — probed rather than read. One high finding, plus a pre-existing packaging defect that shows up here because this is the earliest open PR.
| severity | finding |
|---|---|
| high | the apex return is reached by no test in the suite |
| medium | the INSTALLED package cannot be consumed (pre-existing on main) |
The installed package is unusable
Not introduced by this PR — I reproduced it on main — but every PR ships install(EXPORT ...) and CI never runs install, so nothing catches it:
find_package(numsim-materials REQUIRED)
CMake Error: ... set numsim-materials_FOUND to FALSE ...
The following imported targets are referenced, but are missing:
numsim-core::numsim-core
numsim-core is FetchContent'd into our build and never exported, so the generated numsim-materialsTargets.cmake names a target no consumer can resolve. The install succeeds; the result is just not usable. A find_dependency(numsim-core) in numsim-materialsConfig.cmake.in plus exporting or requiring it would close it — out of scope here, but worth an issue since it silently affects every release.
While checking that, I confirmed the normal developer path still works after today's dependency churn: with Eigen and nlohmann found rather than fetched, configure/build/test are clean at 216/216 on the tip of the stack. I had been testing only the fetched path.
| /// Check if the standard return overshoots the DP cone apex. | ||
| /// When G_shear*Δλ ≥ q_trial, the deviatoric correction flips direction. | ||
| /// G_shear is the plain shear modulus (not G_eff). | ||
| bool needs_apex_return(T G_shear, T dlambda, T sqrt_j2) const { |
There was a problem hiding this comment.
Finding (high): nothing in the test suite ever reaches the apex branch.
I instrumented this predicate with a counter and ran every test binary in the suite:
NO test binary in the suite ever reaches the apex branch
APEX_HITS=0, including test_drucker_prager itself. So needs_apex_return, apex_modified_sig_eq, apex_effective_modulus, apex_plastic_strain and apex_tangent — roughly fifty lines, and the subject of this PR's most recent commit — are executed by nothing.
That matters more here than it usually would, because the code itself says the apex is the hard part:
The return map is nonsmooth at the apex, so this is not a unique classical derivative. [...] the apex branch is rate-indifferent and has no stiffness.
A branch that is documented as nonsmooth, non-unique and stiffness-free is exactly the one whose behaviour you want pinned by a test, and DPConvergence.TangentErrorIsBounded cannot pin it: it drives a single path at 0.02 increments that stays on the smooth cone the whole way, so the tangent it validates is the cone tangent.
tests/debug_apex.cpp exists but is an add_executable, not add_numsim_test — CI compiles it and never runs it, and it asserts nothing. That is the shape that makes this easy to miss: there IS apex code in tests/, it just is not a test.
What would close it: a load path in high triaxial tension that provably enters the branch (assert on a state the apex produces, e.g. q == 0 within tolerance), plus the tangent checker run on that path so the apex tangent is compared against a numerical derivative rather than asserted by construction. If the honest answer is that the apex tangent cannot be FD-verified because the map is nonsmooth there, then a test that documents the chosen convention is still worth more than none.
Nothing in the suite reached it. Instrumenting needs_apex_return() with a counter and running every test binary gave APEX_HITS=0: apex_modified_sig_eq, apex_effective_modulus, apex_plastic_strain and apex_tangent were executed by no test. tests/debug_apex.cpp is an add_executable, not add_numsim_test, so CI compiled it and never ran it. Every existing path is uniaxial and stays on the smooth cone, which is why DPConvergence.TangentErrorIsBounded cannot cover the apex: it validates the cone tangent. The apex sits on the hydrostatic axis, and tensor_component_stepper moves one component at a time, so the tests carry a small hydrostatic driver. Two tests: hydrostatic tension drives the deviatoric stress to zero, which the smooth branch cannot do (there the return is proportional to a nonzero s_trial); and the resulting state is admissible -- pressure capped by the hardening-shifted cone tip k/eta, plastic volume change positive for beta > 0. Same counter after the change: APEX_HITS=37. cmake: the installed package could not be consumed. find_package() failed with 'the following imported targets are referenced, but are missing: numsim-core::numsim-core' -- the exported set names every target linked INTERFACE and the generated Config re-found none of them. It now re-finds whatever was linked at build time; a FETCHED dependency is header-only, reaches the consumer through BUILD_INTERFACE, and is correctly absent from the list. Pre-existing on main, and invisible because CI never runs install.
|
The Root-caused and filed upstream as NumSim-Stack/numsim-core#18 — The two defects are independent; fixing this one exposed that one. Consuming an installed numsim-materials needs both. |
1 participant
Summary
Adds Drucker-Prager (DP) plasticity with non-associative flow and apex return
to the existing J2 framework. The same
small_strain_plasticitytemplate nowhandles both J2 and DP via a yield-function policy. Also fixes a class of
tangent-consistency bugs and adds extensive documentation.
What's new
Drucker-Prager yield function
drucker_prager_yield_function<T, Dim>policy with frictionη, dilatancyβ, bulk modulusK_bulkη ≠ β) — yield normalM = dF/dσand flownormal
N = dG/dσare distinctdrucker_prager_plasticity<Traits>mirrorsj2_plasticity<Traits>Apex return
√J₂below zero (pressure-dominatedloading)
Newton iterations
r = η·p_trial - K·η·β·Δκ - k - H = 0K·H' / (K·η·β + H') · I⊗Ihas_apex_returnconcept — J2 path is unaffectedTangent consistency fixes
mismatch between
modified_equivalent_stressandyield_normal√J₂ + η·I₁paired withs/(2q) + η/3·I(factor-of-3 mismatchin the volumetric coupling)
√J₂ + η·pconvention everywhereuses
N_trial)~10⁻⁴ to ~10⁻¹⁰ as a side effect
Refactoring
plasticity_utils.h:evaluate_at_state,compute_trial,compute_tangent,make_IIdevsmall_strain_plasticity::compute()is now a flat dispatch: elastic /apex / smooth — each path is a named private member function
solve_scalar_return(phi, G_eff, kappa_n)is shared between smooth andapex Newton (matches the unified residual structure)
α → κ(avoids clash with DP frictioncoefficient; matches the document notation)
α → ηin codeDocumentation
docs/small_strain_plasticity.md: full derivation of yield functions,return-mapping algorithm, consistent tangent (smooth and apex), policy
interface, consistency requirements, and common failure modes
Tests and tooling
test_drucker_prager.cpp: 6 tests covering elastic-before-yield, yielding,volumetric plastic strain, pressure sensitivity, consistent tangent, and
tangent convergence (max error 5.2e-10)
scripts/plot_plasticity.py+tests/plot_data.cpp: stress-strain andtangent-accuracy plots for three loading paths (ε₁₁, ε₂₂, ε₁₂)
tests/debug_apex.cpp: standalone diagnostic for the apex-overshootcondition
Bug fixes (also applied to existing code)
rk_plasticity.h:jacobian(m_G, ...)should bejacobian(effective_modulus(m_G), ...)— pre-existing scaling bugbackward_euler::solve(): now reports convergence viaconverged()andclamps result to
≥ 0(no backward plastic flow)small_strain_plasticity::compute(): throws on non-convergence (was silentwrong-result bug); falls back to apex if smooth Newton fails and apex
support exists
apex_tangent: scale-relative zero-denominator threshold(
eps · (|Kηβ| + |H'|)) instead ofnumeric_limits::min()flow_normal_stress_derivativenow takessig_devdirectly — avoidscancellation error from reconstructing
sfrom non-associativeNevaluate_at_statethreshold scaled toσ_0— prevents subnormal-stressdivisions in downstream
1/J₂termsTest plan
fixes)
bug fixed)