Add vector_newton: coupled Newton for mixed scalar/tensor local systems - #17
Open
petlenz wants to merge 15 commits into
Open
Add vector_newton: coupled Newton for mixed scalar/tensor local systems#17petlenz wants to merge 15 commits into
petlenz wants to merge 15 commits into
Conversation
…e_nonnegative (Refs #13)
…ar/tensor solves (Refs #15)
The FetchContent fallback for tmech already existed but could never fire: find_package(tmech QUIET) succeeds against ANY installed version, so a stale system install always won. tmech's package config exposes no version constraint that could express "new enough", so there is no way to make find_package reject one — inverting the priority is the only fix. Had the fallback ever fired it would also have misbehaved: it set BUILD_TESTS, which is numsim-core's option, not tmech's. tmech uses TMECH_BUILD_TESTS and TMECH_BUILD_EXAMPLES, both defaulting ON, so fetching would have built tmech's entire gtest suite and examples into this project's build tree. Adds NUMSIM_USE_SYSTEM_TMECH (default OFF) to opt back into an installed copy via find_package(tmech REQUIRED). For co-developing the two repositories, -DFETCHCONTENT_SOURCE_DIR_TMECH=/path/to/tmech points at a local checkout — the same escape hatch already documented for numsim-core. GIT_TAG master matches the existing numsim-core convention. Verified: a clean configure fetches tmech at acbe92f, builds, and leaks no tmech tests into ctest.
…ms (Refs #14, #12, #15) The vector analogue of backward_euler::solve — a dense pivoted Newton for R(x)=0 over a heterogeneous unknown set. Unblocks local return maps that are coupled systems: multi-surface plasticity, kinematic hardening with a tensor back-stress, viscoplastic couplings. - solvers/vector_newton.h: evaluate-first iteration, so a solve that converges on its last allowed update is not falsely reported failed. Residual inf-norm convergence, partialPivLu with a backward-error guard against a singular Jacobian, honest converged(), and no hidden projection. Covers both the direct solve() and the graph-driven update() mode (#12), plus solve_with_factorization() so a consistent tangent can reuse the converged LU without forming an inverse. - core/unknown_spec.h: the runtime unknown description, kept out of solvers/ so the JSON layer depends on neither the solver nor tmech. - JSON: unknown_spec and block_ref readers modelled on the existing sop converter. vector_newton registers as a single factory type — the layout is dispatched from the "unknowns" parameter rather than baked into the template, so there is no instantiation per unknown combination. Topology follows backward_euler: the solver produces the unknowns and consumes the residual and Jacobian pieces, so the function material reads the trial state back over Local edges. Jacobian blocks are wired dense by default, making an omitted block a wiring error rather than a silently-wrong Jacobian; "zero_blocks" makes sparsity opt-in. Rank-4 blocks must be minor-symmetric, since Mandel stores 36 of 81 components. Major symmetry is not required, so a non-associative algorithmic tangent packs correctly — there is a test on a minor-symmetric, major-asymmetric block.
petlenz
force-pushed
the
feature/vector-solver
branch
from
July 28, 2026 21:26
fcf8756 to
02b6243
Compare
Fetch tmech instead of finding it; fix the fallback's option names
…luation
Review fixes for the coupled Newton solver.
Name pre-check, before any property is created. This has to run first because
the property registry cannot report the problem afterwards: add_property()
silently returns the EXISTING property when a name is taken and downcasts it to
the requested type. A duplicate unknown name would therefore alias two unknowns
onto one storage slot while N still counted them separately, and a duplicate
whose kind differed would reinterpret a property<double> as a property<tensor2>
— undefined behaviour, with no diagnostic anywhere. check_names() now rejects
empty lists, empty names, duplicates, and 'zero_blocks' entries naming an
undeclared unknown (a typo there would silently drop a real block from the
Jacobian).
It also rejects colliding generated names rather than banning underscores:
unknowns {a, b_c, a_b, c} produce "jacobian_a_b_c" from both (a, b_c) and
(a_b, c), so the check requires every generated property name to come out
distinct while leaving "back_stress" perfectly usable.
gather_jacobian() no longer calls update_source(). The material computes
residual and Jacobian in one callback and gather_residual() has just run it at
this iterate, so re-triggering per block re-ran compute N^2 times per iteration
for no change in the result.
solve_with_factorization() now throws unless the last solve converged; it would
otherwise hand back a factorization from an earlier iterate, or from an earlier
call entirely, as though it were the converged tangent.
Tests: name validation, factorization guard, a mistyped producer property
failing to wire (previously claimed as a design benefit but never tested), and
the mixed scalar/tensor system at Dim 2 — the width-3 sym_tensor path had no
coverage at all. 77 tests pass, up from 67.
Documents that a failed solve leaves the raw iterate in the output properties.
That matches backward_euler::solve, which also returns its raw iterate on the
failure paths, but the values here sit in graph properties that downstream
materials read regardless of converged().
…obian block Regression from the previous commit, which narrowed update_source() to the residual inputs on the grounds that the material computes residual and Jacobian in one callback. update_source() fires the callback of the SPECIFIC property it is called on. A material that binds its compute to a Jacobian block therefore never got re-evaluated at all: the solver gathered default-constructed properties, saw R == 0, and reported convergence at the initial iterate. Not a slow or failed solve — a confident wrong answer, which is the exact failure class the rest of this solver is built to avoid. Restores the update over every input, but hoisted into refresh_sources() and called once per iteration rather than once per gather. Inputs whose source has no callback cost an empty std::function check, so in the normal case compute still runs exactly once per iteration. Adds the regression test that caught it: a material binding compute to jacobian_x_x must solve identically to one binding it to residual_x. Every existing test material happened to bind to a residual property, which is why the suite stayed green through the regression.
The Newton iteration tolerates an inexact Jacobian — the root is pinned by
R == 0, so a wrong J costs iterations, not accuracy. The linearization does
not. The consistent tangent comes from the implicit function theorem,
J * dx/deps = -dR/deps, J = dR/dx at the converged point
which inverts J as the actual derivative. Zeroing out a block that is not
actually zero therefore yields a silently wrong dx/deps, hence a wrong
dsigma/deps and the loss of quadratic convergence in the global Newton.
solve_with_factorization() hands back exactly that J, so it inherits the
constraint.
The previous test made this worse by endorsing the unsound usage: it declared
dR_x/dy zero for a system where that entry is 1, called it "a lie", and
concluded that an omitted block "changes the path, not the answer" — true of
the root, false of the tangent.
Replaces it with a genuinely decoupled system (dR_x/dy identically zero), where
the omission is exact, and asserts that the reused factorization reproduces the
true J^-1 column — so the test now pins the property that actually matters.
Documents the constraint at the block wiring and on
solve_with_factorization(), and notes there that a multi-column right-hand side
is accepted, so the whole N x 6 strain derivative solves in one call.
Member
Author
Newton convergence on the test systemsMeasured through the solver itself rather than a reimplementation: run with Nonlinear 2×2 —
|
| k | ‖R‖∞ | x | y | r/r_prev² |
|---|---|---|---|---|
| 0 | 5.0000e+00 | +0.0000000000 | +0.0000000000 | |
| 1 | 2.5000e+01 | +5.0000000000 | +3.0000000000 | 1.000 |
| 2 | 5.7113e+00 | +2.6101694915 | +1.8983050847 | 0.009 |
| 3 | 1.1831e+00 | +1.5224833506 | +1.8651055884 | 0.036 |
| 4 | 1.8142e-01 | +1.0965508740 | +1.9789946553 | 0.130 |
| 5 | 8.4280e-03 | +1.0047466324 | +1.9989122235 | 0.256 |
| 6 | 2.2411e-05 | +1.0000126379 | +1.9999971348 | 0.316 |
| 7 | 1.5971e-10 | +1.0000000001 | +2.0000000000 | 0.318 |
| 8 | 0.0000e+00 | +1.0000000000 | +2.0000000000 |
Asymptotically r/r_prev² settles at ≈ 0.318 — textbook quadratic convergence.
Note the first step: the residual rises 5 → 25. At (0, 0) the Jacobian is [[0,1],[1,0]] — nonsingular, but a poor linearization point, so the full Newton step overshoots to x = 5. It recovers, but this is precisely the behaviour that a line search would damp, and is worth keeping in mind against the "no damping" item in #21.
Mixed, 1 scalar + 1 symmetric tensor (N = 7) — g + tr(B) = 2, B = g·P
| k | ‖R‖∞ | g |
|---|---|---|
| 0 | 2.0000e+00 | +0.000000000000 |
| 1 | 2.2204e-16 | +0.645161290323 |
One update to machine precision — because this system is linear in (g, B). Worth stating plainly: the committed MixedScalarAndSymmetricTensorSystem test validates the plumbing (block assembly, Mandel packing, the N×N solve, unpacking) but does not exercise iterative convergence at all.
Mixed, genuinely nonlinear (N = 7) — B = g·P − 0.1·g²·I
| k | ‖R‖∞ | g | r/r_prev² |
|---|---|---|---|
| 0 | 2.0000e+00 | +0.000000000000 | |
| 1 | 4.1623e-02 | +0.645161290323 | 0.010 |
| 2 | 2.1186e-04 | +0.691189444210 | 0.122 |
| 3 | 5.6022e-09 | +0.691426133042 | 0.125 |
| 4 | 2.2204e-16 | +0.691426139301 |
r/r_prev² settles at ≈ 0.125. Quadratic convergence is preserved with a tensor unknown in the system, which is the substantive result here: it confirms the Mandel packing keeps the flat system equivalent to the tensor one, exactly as the isometry argument predicts. A Voigt packing would not preserve the ∞-norm this is measured in.
Takeaways
- Quadratic convergence holds for scalar-only and for mixed scalar/tensor systems alike.
- The committed N = 7 test is linear, so it converges in one step and never exercises the loop. Adding the nonlinear variant above as a test would close that gap cheaply.
- Undamped Newton can increase the residual substantially on the first step from a poor initial guess (5 → 25 here). Relevant to vector_newton: robustness and coverage follow-ups from review #21.
The existing mixed test is linear in (g, B), so Newton reaches machine
precision in a single update. It pins the plumbing — block assembly, Mandel
packing, the NxN solve, unpacking — but never runs the iteration for a system
containing a tensor unknown, so a regression in the loop would not be caught
there.
Adds a quadratic term, R_B = B - g*P + c*g^2*I, behind a 'c' parameter that
defaults to 0 so the existing linear tests are unchanged. Measured convergence
from the zero iterate:
2.00e+00 -> 4.16e-02 -> 2.12e-04 -> 5.60e-09 -> 2.22e-16
with r/r_prev^2 settling around 0.125 — quadratic, with a tensor unknown in the
system. That is the property the Mandel packing has to preserve: the
convergence test is an infinity-norm over the packed vector, so a
non-isometric packing would not keep it equivalent to the tensor system.
The test asserts the closed-form root of D*c*g^2 - (1 + tr P)*g + a = 0, the
matching B, and — so it cannot quietly go vacuous — that a single update does
NOT converge.
…Refs #19) Differentiating the converged residual R(x(eps), eps) = 0 gives J * dx/deps = -dR/deps, J = dR/dx at the converged point which is solved from the factorization already computed there — no inverse is formed, and the whole N x W system is one multi-column back-substitution. The material completes the chain rule itself: dsigma/deps = dsigma/deps|_x + (dsigma/dx) : (dx/deps) Opt-in through a 'strain_derivative' parameter, default off, since it costs a gather and a solve that a caller wanting only the root does not need. When on, the solver consumes d_residual_<name>_d_strain and produces d_<name>_d_strain. dR_I/deps has the shape of a Jacobian block whose column kind is the strain, so jacobian_block_layout already serializes it and no new packing was needed. What was missing is the scatter direction: dx/deps comes back as an N x W matrix and has to become per-unknown tensors of rank r_I + 2 again. block_scatter_base and strain_derivative_layout add that, mirroring the row-major/column-major transposition that jacobian_block_layout owns on the way out. Verified against a closed form rather than only a round trip. For R_g = g + tr(B) - a - tr(eps), R_B = B - g*P the solution is g = (a + tr eps)/(1 + trP), hence dg/deps = I/(1+trP) and dB/deps = (P (x) I)/(1+trP) — matched to 1e-12. The case is chosen so the scalar row of dR/deps is nonzero while the tensor block is zero, and so dB/deps is minor-symmetric but major-ASYMMETRIC, with an assertion that the reference really is asymmetric so the transpose check cannot go vacuous. Layout-level tests cover scatter inverting gather at non-zero offset and stride for both the 1 x W row and W x W block shapes. Also registers bool with the JSON reader registry, which had no reader for it.
Library: the strided placement in jacobian_block_layout::gather and strain_derivative_layout::scatter_block was hand-written index arithmetic — precisely where a row-major/column-major slip hides. It is now expressed as Map-to-Map assignments, so Eigen performs the storage-order conversion. The Rx1 column case drops its scratch buffer entirely: those destination slots are contiguous, so tmech packs straight into them. Tests: nested component loops replaced by tmech::almost_equal for tensor comparisons, and by Eigen Map expressions (dot, matrix-vector, isApprox) for flat buffers. The dimension-templated test tensor is a literal per dimension instead of an index-computed fill; only tr(P) enters the closed forms, which the new literals preserve. One trap found while doing this, worth recording. A single row of a column-major buffer must be viewed as a vector with INNER stride ld, not as a fixed-size Matrix<T,1,C> with outer stride: Eigen treats Matrix<T,1,C> as a vector at compile time, linearizes it with inner stride 1, and silently ignores the outer stride — writing the row contiguously instead of across the buffer. The first version of this refactor did exactly that, and the layout unit test did not catch it because it read the row back through the same faulty strided view, so both sides were wrong consistently while the 3D mixed solver tests failed. The placement tests now read through a plain dense Map of the whole buffer, with no custom stride and no fixed-size vector type, so they cannot be self-consistent with a broken writer.
A UMAT host inverts every ownership the framework assumes. Strain comes from STRAN/DSTRAN instead of a stepper, history from STATEV instead of history_property, and the host decides when an increment has converged — while re-calling the material on iterates it later discards. The layer added here makes a finalized material_context behave as a stateless material-point evaluator: raw pointers in, raw pointers out, nothing retained between calls. That statelessness is the design, not an optimisation. Reloading every history variable from STATEV at the top of each call is what makes a repeated call on an unconverged global iterate give the same answer, and it is why one context per thread suffices rather than one per integration point. It also makes the plane-stress solve a one-liner: each trial out-of-plane strain has to integrate from the same converged state at t_n, which is just another unpack. No existing material is modified. An earlier draft converted time to a plain property so the STATEV wrapper would not pick it up; that would have broken autocatalytic_reaction and curing_rate, which legitimately need the old/new pair. Host-supplied strain and time stay history properties — external sources write BOTH sides directly — and the wrapper excludes them by an explicit (material, property) list instead. Every exclusion must match a real history property or construction throws: an unmatched entry would leave host-owned state in STATEV and shift every following slot. Layers, each usable without the one above it: tensor_conversion host pointer <-> canonical 6-slot buffer <-> tmech external_state_source host-driven "strain"/"time", bind() sets old and new statev_map history <-> flat STATEV, with the exclusion list material_point_evaluator unpack -> bind -> update -> read -> pack plane_stress_evaluator outer eps_33 Newton + static condensation umat_interface registry, per-thread contexts, the Fortran shim Every element family runs the same Dim=3 material; 2D is handled by widening the host's NTENS slots into the canonical buffer, not by a second policy. A genuinely 2D material would need a separate implementation per model, and plane stress needs a local solve regardless. Traps found, all now covered by tests that fail without the fix: tmech already ships the Abaqus orderings (abq_std, abq_exp), so the conversion is reuse rather than hand-rolled index tables. Three conventions cross there and only their combination is right: strain carries engineering shear, stress does not, and the rank-4 tangent takes no scaling at all because the factor 2 is already absorbed by the minor-symmetry sum. The tangent is packed row-major by tmech but DDSDDE is a Fortran array, so narrow_matrix transposes — invisible for a major-symmetric tangent, which is why the tests use a deliberately asymmetric one and a second test guards that the fixture stays asymmetric. Plane stress is not a prefix of the 3D slot map: its components are 11, 22, 12, so host slot 2 is the shear, not the out-of-plane normal. Assuming a prefix still yields a plausible-looking 3-vector. STATEV layout is sorted by (owner, property) rather than following the graph's execution order, which comes from a topological sort seeded by unordered_map iteration and is not stable across builds. A layout that moved between builds would silently reinterpret every existing restart file. Failures are classified rather than uniformly retried. fatal_error zeroes the outputs and terminates through XIT; anything else asks for a smaller increment via PNEWDT, which is the right default for exceptions escaping the constitutive models themselves. Getting this backwards is worse than not classifying: an unknown CMNAME originally returned normally with DDSDDE untouched, so Abaqus solved on with whatever the buffer happened to contain, and a misconfigured model requested a cutback forever against a fault no increment size could fix. Under NLGEOM=YES Abaqus rotates STRESS and STRAN but never touches user state, so statev_map rotates tensor-valued history itself. Skipping it leaves plastic strain in a stale frame and the error accumulates silently over the rotation history. Note tmech::dot is first-order only; the rank-2 rotation uses operator*, which contracts the last index of the left with the first of the right. The call struct carries spans, not bare pointers. Array lengths are dictated by the element family, and passing a 4-element array with the default solid3d case reads six — a mistake made while writing these tests, which surfaced only as a wrong value. It is now rejected. Verification beyond unit tests: DDSDDE finite-differenced through the raw pointer boundary in both the elastic and plastic regimes, which is what catches a shear factor, a slot permutation or a transpose; statelessness shown by rebuilding a fresh context every step and matching a persistent reference bit-for-bit through 25 plastic steps; plane stress audited by replaying the recovered eps_33 through an independent 3D evaluator and confirming sigma_33 vanishes; the eight-thread concurrency test run under ThreadSanitizer. Not covered: model configuration is a builder function, so wiring CMNAME to a JSON config file is still unbuilt. SSE and SPD are reported only when the model names its plastic strain, since splitting work into stored and dissipated parts is not possible without the elastic strain; when unset the host's energies are left alone rather than filled with a guess. RPL, DDSDDT, PREDEF and the deformation gradients are accepted and ignored.
…STATEV convention Follow-up to c9ea9ed, from a four-way review (Abaqus contract conformance, C++ safety, numerics, mutation testing). The contract items all came back conformant — argument order, component orderings, column-major DDSDDE, TIME(2), PNEWDT/XIT, the Fortran hidden string length — and the numerics verified against independent computations: the condensed plane-stress tangent matches finite differences of the full outer solve to 2.9e-6, the engineering-shear packing to 2.8e-17, and DROT is worth 15% error over 20 plastic steps. What follows is everything that did not come back clean. Setup faults are now classified at the boundary, not at each throw site. thread_state_for() runs the registered builder — arbitrary graph construction this layer does not own — and the core library reports a mistyped source name or an unknown material type as a plain std::runtime_error. Those reached the generic handler and asked for a cutback, against a fault no increment size can fix, rebuilt and re-failed on every call from every thread. This was the third instance of the same bug (c9ea9ed fixed two); classifying per throw site cannot work when the throwing code belongs to another layer, so construction is now wrapped whole and every failure inside it is fatal. DDSDDE is zeroed on the cutback path too, not only the fatal one. PNEWDT is the minimum over all calls for the iteration, so Abaqus finishes the element loop and assembles BEFORE acting on the request — previously with whatever the buffer happened to contain. A zero tangent is a poor stiffness and that is the accepted trade: the increment is discarded either way, and predictably soft beats arbitrarily wrong. umat_dispatch can no longer terminate instead of reporting. Two allocations escaped its handlers: the CMNAME string was built above the try, and the fatal message was concatenated inside a catch — where an exception is not caught by the sibling catch(...) and so escaped a noexcept function *before* the fatal handler ran, losing exactly the diagnostic that path exists to print. SSE and SPD are both accumulated incrementally, as Abaqus's own materials do. SSE was absolute (1/2 sigma:(eps - eps_p)) while SPD was derived from its change, which is only self-consistent if the host's incoming SSE is the stored energy at t_n — and Abaqus initialises SSE to zero whatever the initial state. Under *INITIAL CONDITIONS, TYPE=STRESS the whole pre-existing stored energy was booked once as NEGATIVE plastic dissipation (measured: SPD = -1.35 on a step with no plastic flow) and the offset persisted in ALLPD for the analysis. Rank-2 state variables are stored in STATEV with engineering shear. The library round-trips either convention exactly, so this was invisible internally and wrong only where it matters: against STRAN, against Abaqus's PE/LE output, and against an initial plastic state seeded through *INITIAL CONDITIONS, TYPE=SOLUTION. describe() now names the convention. Note this assumes tensor history is strain-like, which holds for everything the framework carries; a back stress would want the unscaled form and a per-property decision. Also: CMNAME is folded to a common case on both registration and lookup, since Abaqus input is case-insensitive and delivers the name upper-cased; the fatal handler is atomic (TSan-confirmed race); and plane stress rejects a DROT that tilts out of the plane, which the scalar out-of-plane strain in STATEV cannot represent. Tests: 159 -> 176. Mutation testing found 10 defects that no test caught, all now closed and each re-verified by re-running its mutation. Three are worth recording because the tests looked convincing and were not: The energy balance was a tautology. SSE and SPD are updated as `sse += dSSE; spd += dW - dSSE`, so their sum telescopes to the accumulated work for ANY definition of dSSE — defining the elastic strain as the total strain passed it. The split is now pinned by an absolute SSE on an elastic step and by zero dissipation while elastic, in both the 3D and plane-stress paths. The objectivity test used u::rotate as its own oracle: it rotated the inputs and computed the expectation with the function under test, so R a R^T -> R^T a R cancelled and passed. Compounding it, a 90-degree rotation of a diagonal tensor cannot tell the two senses apart, and an invariant check passes for any orthogonal map. There is now a 30-degree fixture with a shear component whose expectation comes from the closed-form rotation algebra. The STATEV encoding was unconstrained: every test round-tripped through the same pack/unpack pair, so swapping the shear convention or permuting two slots passed them all. The raw slot values are now asserted directly. Smaller ones: the whole rank-4 branch of statev_map was unreachable from tests (no production material has rank-4 history — there is now a test-only material that does); TIME(1) and TIME(2) were indistinguishable because the fixture passed identical values in both slots and no model registered at the ABI used a time source; the warm-start test passed against an implementation that never warm-starts, since both runs were then cold and EXPECT_LE held; and the NTENS-vs-NDI/NSHR guard was dead code, though removing it would have built 6-element spans over a 4-element host array. One comment was actively wrong and is corrected: the finite-difference tangent check claimed to catch a row/column-major transpose. It does not — J2's consistent tangent is major-symmetric, so the transpose is invisible there. That contract rests entirely on the deliberately major-asymmetric fixture in test_umat_conversion, which makes it a single point of coverage worth keeping. Verified: 176/176 in Debug and Release (NDEBUG), 77 tests under ASan+UBSan, the ABI and threading tests under ThreadSanitizer.
The remaining findings from the four-way review — documentation, plus one performance item that turned out to be in the hottest path there is. The macro's "accepted but not used" note listed four arguments and ignored twelve. The omission that matters is PROPS/NPROPS: this shim takes material constants from the registered builder, not from *USER MATERIAL CONSTANTS, which is the single most surprising fact about the design and was written down nowhere. A deck that varies PROPS while naming the same model silently gets the builder's constants for all of them. The full list is now enumerated with the conditions under which ignoring each stops being safe — RPL and friends in coupled thermal-stress, thermal-electrical-structural and adiabatic procedures (the old note said only "coupled temperature-displacement", which is narrower than the contract); TEMP/DTEMP meaning temperature-dependent properties are unavailable, while thermal expansion still works because Abaqus reduces STRAN/DSTRAN to mechanical strain; JSTEP(4) meaning a linear-perturbation step gets a nonlinear update and an elastic-plastic Jacobian; and the deck-side transverse-shear-stiffness requirement for beams and shells. CMNAME handling no longer allocates. trim_fortran_name built an 80-character std::string on every call — past the small-string buffer, so one heap allocation plus a hash per integration point per global iteration. It now returns a view, folding happens into a stack buffer, and the registry uses heterogeneous lookup. The buffer is NUL-terminated so the same storage feeds both the lookup and the C formatting in the error paths, which also removes the last allocation from the fatal path. Three comments were narrower than the truth and are corrected. The row/column-major transpose is observable only under *USER MATERIAL, UNSYMM, since Abaqus otherwise uses just the symmetric part of DDSDDE — worth saying, because that is precisely the configuration a non-associative model needs. The NDI/NSHR note claimed plane strain "happens to pass eps_33 = 0"; generalized plane strain (CPEG) reports the same NDI/NSHR with a nonzero eps_33, and the code is correct for it because it consumes whatever slot 2 holds. And the element families that fall through to the fatal error — beams, trusses, axisymmetric shells, cohesive elements — are now named, since that fallthrough is the intended outcome rather than a gap. Two facts that were true but undocumented. The `false` template argument in the rank-4 conversions is INERT: tmech's abq_std honours _ShearStrain only on the rank-2 path, so the tangent's no-scaling property comes from tmech's internals rather than from the argument written there, and would silently acquire a factor 2 if tmech ever implemented it — the physical check in the tests is what actually pins it. And the minor-symmetry guard on the tangent is an assert, so absent under -DNDEBUG; the consequence is bounded and now stated precisely (STRESS is unaffected, only the Jacobian degrades), which is why statev_map's equivalent throws instead. SCD is documented as never written, and its test assertion made meaningful: it now passes a nonzero value in and requires it back unchanged, rather than comparing zero against a zero-initialised variable — which a do-nothing implementation satisfied. Also added: null guards on CMNAME and TIME, matching the ones the other pointer arguments already had, and a note on why thread_state's member order matters on use rather than on teardown. Verified: 176/176 in Debug and Release, 77 tests under ASan+UBSan, the ABI and threading tests under ThreadSanitizer, and -Wall -Wextra -Wconversion -Wsign-conversion -Wshadow -Wpedantic clean across the umat headers (the remaining warnings all originate in vendored tmech).
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.
Summary
Adds a coupled Newton solver for
R(x) = 0over a mixed scalar/tensor unknown set — the vector analogue ofbackward_euler::solve— plus the Mandel serialization layer it needs.Closes #13. Implements #14, #12, and the solver half of #15.
Commits
8e84f13backward_euler::solve()becomes a general scalar Newton; opt-insolve_nonnegative()7a1c7beunknown_layout: tensor ↔ flat serializationfcf8756vector_newton+ JSON configDesign
Per-property layouts behind a type-erased base. Every piece of the local system — each unknown, each residual entry, each Jacobian block — carries its own layout. The solver holds a vector of layouts for the residual and a matrix for the Jacobian, and knows only
layout_base. Consequences:input_property<T>::wire()'sdynamic_casttype-checks each one; a rank mismatch fails to wire rather than corrupting the system.Ntemplate parameter, so it registers as a single factory type; the layout is dispatched from the"unknowns"parameter rather than baked into the template.Symmetry must be declared, not inferred.
tmech::tensor<T,Dim,2>is symmetry-agnostic storage (eveninv()takes a symmetrysequenceat the call site), andproperty_traitscarries no shape metadata. Hence the explicit kind tags, with the runtime counterparts incore/unknown_spec.hso the JSON layer depends on neither the solver nor tmech.Block layouts are parameterized on the kind pair, not the value type. Blocks
(scalar, sym)and(sym, scalar)are bothtensor<T,3,2>but pack as a 1×6 row and a 6×1 column respectively — the value type cannot disambiguate.Mandel, not Voigt.
A:B == mandel(A)·mandel(B)and(D:X) == mandel4(D)·mandel2(X), so the flat LU solve and the ∞-norm convergence test are equivalent to the tensor system. Applying the same weights uniformly to residual, unknown, and both sides of every block is what makes this work.Iteration
Evaluate-first, so the residual after the final update is the one tested — a solve converging on its last allowed update is not falsely reported failed. Residual ∞-norm convergence,
partialPivLuplus a backward-error guard against a singular Jacobian (cheaper than rank-revealing pivoting every iteration, and finiteness alone would not catch it), honestconverged(), and no hidden projection — which is what #13 was about.JSON
{"type": "vector_newton", "name": "solver", "function": "kin", "unknowns": [{"name": "dgamma", "kind": "scalar"}, {"name": "backstress", "kind": "sym_tensor"}], "zero_blocks": [["dgamma", "backstress"]]}No
"dim"— the dimension is fixed by the Traits policy, which keeps the kind set small enough for an exhaustive switch. Blocks are wired dense by default so an omitted block is a wiring error, not a silently-wrong Jacobian; sparsity is opt-in.Testing
67/67 pass (was 49).
test_j2_plasticityandtest_drucker_pragerare the regression signal that de-clampingsolve()changed nothing for existing callers.Three things the tests corrected during development:
(C:N)⊗(M:C)shape — confirming it lands untransposed at non-zero offset and stride.backward_euler: the solver produces the unknowns and consumes residual/Jacobian, so the function material reads the trial state back over Local edges.Not included
The physics integration tests from #14/#15 — two-surface plasticity and J2 with tensor back-stress — are new constitutive models rather than tests, and are left as follow-ups.
solve_with_factorization()lands as a hook for #12's consistent tangent, but the tangent-block wiring waits for a consumer; it holds a stale factorization if the solve did not converge, so callers must checkconverged().