Repository navigation
Perf/optimization v2 - #148
Merged
Merged
Conversation
benchmark/benchmark_batch.py times the batch Python APIs (umat, umat_T, objective_rate, Lt_convert, stress_convert, Log_strain, rotations) on N points per call, plus sequences mirroring one fedoo global iteration in small strain and updated Lagrangian. Each point is the minimum time per call over adaptive repetitions; the summary command fits t = a + b*N on the serial runs, separating the fixed cost of a call from the cost of a point. benchmark/dump_results.py dumps every output array of the 0D solver (the benchmark cases, the legacy/MODUL families and the solves of five gallery examples, run live) and of the batch umat over loading and unloading, and compares two dumps bit for bit or within a tolerance, failing with status 1 otherwise. Performance changes are expected to publish this comparison. bench_legacy_vs_modular.py now benchmarks only EPICP, EPCHA and PRONK by default: the ten other legacy names are routed to the modular engine by a props adapter, so they compared the engine with itself. --all still runs them, flagged as adapters. The README documents the tools and replaces the February 2026 estimates of the five engine hotspots by their measured state on 2.2.0.
The Gauss tables of the Eshelby integrals were process-wide statics of ellipsoid_multi, rewritten by every L_eff and at the start of every ellipsoidal solve while solver_run releases the GIL: concurrent solves read each other's tables (IndexError, or wrong results without bounds checks), and nested schemes with different (mp, np) used the inner ones. Each (mp, np) now has immutable shared tables (eshelby_quadrature_points) that the schemes read from their own props and pass to fillT explicitly. DE_Homogeneous_E dereferenced a null pointer (MIHEN crashed the solver). MISCN failed on the second block of any loading path: the solver clears the effective tangent before priming a block and the self-consistent scheme inverted it as its reference medium; a cleared tangent is now seeded again as at start. n_threads is now honoured under GCD and OpenMP (it only took effect for 1) and Lt_convert takes it too: the default of 4 now means 4 threads. Also drop the dead re-accumulation of dE_Periodic_Layer and fix the nmse_per_response documentation (the tests are pooled per column).
umat_modular rebuilt and reconfigured a ModularUMAT at every call (stiffness, its Cholesky positive-definiteness check, typed tensors), which was most of an elastic integration point. Each thread now keeps up to eight configured instances in a small LRU keyed on the props (bit for bit), nstatev and the kinematic form, so the phases of a multiphase material keep theirs while their props alternate. The state still comes from statev at every call: run() now unpacks it on the start increment of a reused instance too, every per-call cache of the mechanisms is rewritten before it is read (the damage compliance is keyed on the exact stiffness), and an instance whose call throws is dropped. The small-strain route also skips work that is the identity: PKII is tau when F1 is exactly I, etot is summed without the Voigt rotation when DR is exactly I, and the sub-phases are no longer copied in to_start/set_start. The memory sink reserves its record vectors from the expected output count instead of copying every arma object at each growth. The plasticity mechanism reads p and EP through indices resolved at registration rather than string lookups, and the EPICP kernel skips the Fischer-Burmeister solve on an admissible trial and the flow Hessian when no multiplier is active; both shortcuts return exactly what the full path returns. All solver and batch outputs are byte-identical to the previous version. An elastic ELISO increment drops from 3.05 to 0.95 us, a batch ELISO point from 1.70 to 0.39 us and an EPKCP plastic point from 13.0 to 8.4 us; the peak memory of a 1e6-increment run drops from 4.8 to 3.2 GiB.
Gauss now evaluates the calG integrand on plain arrays, with the azimuthal sin/cos hoisted out of the loops and a per-ring pass before the accumulation; every entry of G receives the same operations in the same order, so the tensors are bitwise unchanged (x8 on 50x50). The Eshelby tensors integrated by a thread are kept in a small exact-key cache that outlives the phase objects L_eff and the solver rebuild, fillT returns its last S, T_loc and T when called again with bitwise the same tangents and geometry, the global-frame rotation is skipped, and Eshelby/T_II reuse the shared quadrature tables. L_eff releases the GIL around the homogenization so that identifications can run phases in threads. All outputs are bitwise identical; the composite SC identification goes from 74 s to 9.7 s.
Replace the LAPACK eigen-solves of the per-point finite-strain kinematics by a fixed-size cyclic Jacobi solver for symmetric 3x3 matrices (eig_sym_3x3, Kopp 2008), shared by the polar decompositions, the logarithmic strains, the objective rates and the exact box/material tangent maps. The log_R kinematics now take R1 and the eigenbasis of B1 from a single decomposition of C1 and apply A^R:D by spectral congruence; a corate_kinematics overload returns the strain increment in the same pass and the solver and the bindings use it. The global Newton solves its 6x6/7x7 system with a stack LU (small_lu) instead of a LAPACK inverse, Fischer-Burmeister solves its one- and two-mechanism systems in closed form, the SMA laws invert their isotropic or cubic Reuss mixtures in closed form (inv_cubic), and the isotropic constitutive tensors, the Voigt rotation operators and linearop_eigsym avoid heap temporaries. The modular closest-point return mapping accumulates its condensation matrix in a fixed 6x6 factorized once per iteration with small_lu, with work buffers sized once. The batch umat binding copies its outputs and zeroes its tangent per point inside the parallel region, and stress_convert runs in parallel on the slice memory. fillT no longer rotates an isotropic reference medium, so all orientation families share one Eshelby tensor. Results change at round-off level only: 1e-13 typically, up to 1e-10 on the free stresses of nearly incompressible hyperelastic solves, and 2e-6 relative on the tiny (1e-7) heat terms of the thermomechanical SMA. log_R objective rates go from 4.9 to 0.7 us per point and one fedoo updated-Lagrangian iteration is about 2x faster.
Laws built on F (the hyperelastic potentials and Saint-Venant-Kirchhoff) close on their spatial tangent, of which dS/dE is a plain pull-back, yet they always returned the box tangent of the solver's rate, which the Green-Lagrange and Biot solver routes and sim.umat's material or spatial tangent output immediately mapped back. Such a law now declares it in its output convention, and a caller can ask for dS/dE through corate_type_material_tangent, a tangent form carried in the corate slot that the corate dispatchers accept, so no kernel signature changes. The solver keeps the box contract of Lt outside its Newton loop and rebuilds the box once per increment. Alongside, the isochoric invariant tangent builds only the terms of non-zero potential derivatives (bit identical), the Biot route computes one polar decomposition of F1 per iteration and the start Biot stress once, the modular Hencky block gets its box tangent in closed form from the Lie tangent on the eigenbasis of the elastic strain, and RU_decomposition, VR_decomposition, isochoric_pstretch and stress_convert accept batches of 3x3 tensors, run in parallel without the GIL. Results change at round-off level only (1e-12 on the hyperelastic solves, more on nearly incompressible stresses). sim.umat NEOHC with a spatial tangent output is 3.1x faster and the hyperelastic identification example goes from 9.8 to 6.0 s.
Expose three numerical choices of the constitutive integration, all off by default and bit identical to the previous behaviour when off. A min-NCP complementarity step (ncp_step) replaces the Fischer-Burmeister step of a single-surface cutting plane when asked for (EPICP props[6], EPKCP props[7], MODUL through ModularMaterial(ncp="min")): the first step is the full cutting-plane step, so a von Mises return with isotropic or Prager hardening takes three local iterations instead of six, and the solve falls back to Fischer-Burmeister as soon as its error grows. Several coupled surfaces keep Fischer-Burmeister and Armstrong-Frederick back stresses are refused, since the cutting plane integrates them along its iterates and the full steps would change the return beyond the tolerance. EPICP props[7] selects the closed-form consistent tangent of the radial return instead of the generic 6x6 inversion. The global Newton accepts a relative tolerance, precision_rel in solver_params and solve(), giving max(precision, precision_rel |tau|). Mean-field schemes get opt-in numerical options as trailing props slots, all 0 by default (bitwise the previous behaviour): MIMTN [mp, np, n_matrix, closed_form, rule], MISCN [mp, np, n_matrix, start, closed_form, rule, anderson, rtol]. closed_form uses the sphere and spheroid Eshelby tensors when the reference medium is isotropic (and integrates otherwise, or near a_r = 1 where the spheroid formulas cancel); rule = 1 integrates on a half sphere with a trapezoidal rule in phi, half the points of the Gauss-Legendre rule for an equal or better accuracy; anderson and rtol accelerate and stop the self-consistent fixed point of L_eff on a relative criterion. discretize_odf gains weights="fourier", exact for Mori-Tanaka with an isotropic matrix with a handful of phases (non-negative moment- matching weights when the projected ones are negative, Simpson with a warning when the density is too peaked). The composite self-consistent identification drops from 9.3 s to 0.74 s, an orientation average from 180 phases to 5 at 28 times less cost. Every option is off by default: all solver, batch and homogenization outputs are bitwise those of the previous commit. The plane-stress predictor of the plan is left out (it needs a tangent before the first evaluation, which the condensed API does not return), and the EPCHA tangent mode is only documented.
identification() gains method="least_squares": bounded trust-region reflective least squares on parameters mapped to [0, 1] (logarithmically for positive bounds), from n_starts starting points, with a finite-difference Jacobian. calc_residuals() returns the residual vector whose sum of squares is calc_cost() for the mse, nmse and nmse_per_response metrics. On the Chaboche and composite examples it needs 16 to 31 times fewer model evaluations than differential evolution and reaches an equal or lower cost. The default remains differential_evolution, unchanged. n_threads evaluates the cost in a thread pool: a generation of differential evolution (updating="deferred") or the Jacobian columns of the least squares, whose result then does not depend on the thread count. It pays off because solver.solve and L_eff release the GIL; the default stays serial. KeyTemplate keeps the key template files in memory and writes each working file once per evaluation; copy/apply_parameters and copy/apply_constants now copy, read and write each file once instead of once per key. The Chaboche, composite and Mooney-Rivlin examples use the least squares: the Chaboche one goes from 220 s to 3.6 s, the Mooney-Rivlin one from 10 s to 4 s, and the composite one from 67 s to 6 s on 2.2.0 (1.6 s with the faster homogenization). The Chaboche pre-cycle runs 64 increments per step with sub-incrementation allowed, and the Mooney-Rivlin example shows the linear least-squares shortcut.
Finite-element codes rotate their Gauss-point stresses, strains and tangents by the (3,3,N) rotation increments objective_rate returns. The only way to do so from Python was to build a scipy Rotation from the matrices, a serial conversion to quaternions of about 1.5 us per point, and to apply it with einsum: in an updated-Lagrangian fedoo iteration this was a fifth of the time once the kernels ran in parallel. simcoon.rotate_voigt(values, DR, kind, active, n_threads) rotates (6, N) stress or strain vectors, or (6, 6, N) stiffness, compliance and concentration tensors, by one rotation matrix per point or one for all, on the layout of umat and objective_rate. It exposes the batch_rotate kernels of tensor.hpp, which now take n_threads, and runs them without the GIL. The result is that of Rotation.from_matrix(DR).apply_<kind> to a few ulps (at most 3e-15 relative). Rotating strain and stress goes from 1.68 to 0.126 us per point in serial and 0.013 in parallel (N = 1e5), a 6x6 stiffness from 0.75 to 0.145 and 0.011. Existing outputs are bitwise unchanged. The batch benchmark gains the fedooUL2 cases, the iteration as fedoo runs it with rotate_voigt.
Move bench/bench_legacy_vs_modular.py into benchmark/, update usages and imports, and adjust related documentation and comments. README.md (benchmark/README.md) is expanded with a temporary-folder note, updated profiling hotspots, and a migration table for laws moved to the modular engine. dump_results.py updates references to the new benchmark path, removes the old bench sys.path insertion, and prefixes solver case names with "benchmark/". A comment in legacy_adapters.hpp was updated to reflect the new script path. These changes consolidate benchmarking scripts and refresh docs and paths accordingly.
The complementarity function of the local cutting plane was a trailing props slot (EPICP props[6], EPKCP props[7], the MODUL slot after the last mechanism, ModularMaterial(ncp=...)). It is a setting of the run, not a material parameter: it is now given by name, solve(..., ncp="min") and sim.umat(..., ncp="min"), "fischer_burmeister" being the default. Like tangent_mode it travels in solver_params::ncp and state_variables::ncp, and EPICP, the modular adapter (EPKCP) and MODUL take it as a trailing argument. check_ncp_serves refuses "min" for any other law, for a thermomechanical block and for a mean-field model; MODUL still refuses it with an Armstrong-Frederick back stress (ModularUMAT::set_ncp). The props of these laws are back to their material parameters only. The closed-form consistent tangent of the radial return (EPICP props[7]) is no longer an option: tangent_mode 2 always takes it. It is the operator of the generic assembly to round-off, without the 6x6 inverse: an EPICP batch of 1e4 points runs in 20 ms instead of 35 ms (14 ms with ncp="min"). EPICP outputs differ from the previous ones at round-off (at most 8e-15 relative on the solver and batch dumps); every other law is bitwise unchanged.
The UMAT catalog explains how the two ncp choices work: the Kuhn-Tucker problem, the Fischer-Burmeister and min equations, why Fischer-Burmeister takes two thirds of the first cutting-plane step, the fallback, why several surfaces keep Fischer-Burmeister, and what ncp="min" does case by case. Fischer_Burmeister_m gets its Doxygen (equation, Jacobian, stopping criterion).
ncp="min" fell back to Fischer-Burmeister whenever the local system had more than one row. In MODUL a yield surface composed with viscoelastic (Prony) or damage mechanisms has 1 + N rows, but those extra rows are integrated outside the Newton (closed-form branches, damage fixed point): always satisfied (Phi = -Y_crit), no multiplier, no cross term. The option then silently had no effect. ncp_step now takes the min step row by row whenever the local Jacobian is diagonal, which covers one surface and one surface next to such rows; coupled systems (several yield surfaces) keep Fischer-Burmeister. A single row runs the same arithmetic as before. The default is bitwise unchanged (solver and batch dumps). With ncp="min", a cyclic uniaxial MODUL solve drops by 17 % with one Prony branch next to the yield surface, 29 % with three, 23 % with a damage mechanism, for the Fischer-Burmeister return to the local tolerance.
tangent_mode mixed two choices: the operator handed to the global Newton and, with mode 3, the local integrator (closest point instead of cutting plane). The run now chooses them apart, by name, next to ncp: integrator="cutting_plane" (default) or "closest_point", and tangent_mode "none", "continuum" or "algorithmic", the consistent operator of the integrator that ran. tangent_mode=3 / "closest_point" is no longer accepted: it raises ValueError pointing to integrator="closest_point". kernel_tangent_mode folds the pair into the code the kernels read (unchanged: 3 is still the closest point inside), for solve, sim.umat and sim.umat_T. The closest point is offered with its algorithmic operator only (its continuum operator and the explicit form are not provided: ValueError). ncp="min" is refused under the closest point, which keeps Fischer-Burmeister. INTEGRATORS and NCP_METHODS join the solver maps; tangent_closest_point leaves the public solver exports (it stays the kernel code a PYEXT law receives). Tests, the MODUL closest-point example and the docs move to integrator=; the tests that sweep the kernel codes go through solver_harness. tangent_kwargs. Solver and batch dumps are bitwise unchanged.
Brings in #146 (fix/toward-2.3: FE transfer and exception types) and the #147 merge. Conflicts resolved by keeping both sides: the ncp / integrator helpers next to total_formulation_serves_T, the ncp argument of the modular adapter inside master's logic_error handling, master's condense check in sim.umat / sim.umat_T before the integrator / tangent fold (whose range check replaces master's), the two sides' new objective-rate tests, and one comment for the same FD step change in test_umat_condensed. The modular cache test expects the configuration error master now raises (std::logic_error). Solver and batch dumps match the branch bit for bit.
Remove the legacy copy_parameters/apply_parameters and copy_constants/apply_constants file-based helpers and switch to the in-memory KeyTemplate workflow. Updated docs and examples to demonstrate KeyTemplate.write usage (including json_keys_identification and hyperelastic/hyperelastic_umat_identification). Cleaned up parameter.py and constant.py by deleting the old copy/apply implementations and related shutil usage, and adjusted python-setup/simcoon/__init__.py exports. Tests were updated to use KeyTemplate and to reflect clearer test names/expectations. Also increased ninc from 64 to 100 in identification path JSON fixtures to use finer discretization.
Add support for constants that hold one value per test and substitute the appropriate value when writing/rendering templates. Changes: - Add CHANGELOG.md with release notes and Unreleased section. - Update docs (identification.rst) to document constants with one value per test and KeyTemplate.write(..., test=i). - Modify python-setup/simcoon/parameter.py: - Add numpy and warnings imports; track warned single-value constants. - Implement _consts_of_test to select per-test constant values (or raise IndexError if out of range) and warn once for constants with a single value. - Extend render and write to accept optional test index, and factor rendering into _render. - Add tests (test_identification_options.py) covering per-test substitution, out-of-range errors, and single-value warning behavior. These changes allow writing one working directory per test with test-specific constant values and provide clearer errors/warnings for misuse.
Add check_integrator_serves to refuse integrator="closest_point" for laws that do not implement a closest-point update and wire it into the solver and UMAT dispatch (C++ and Python). kernel_tangent_mode now rejects tangent_mode==3 (migration message), and solver checks thermomechanical runs when validating integrator/tangent choices. Update Python wrappers/docs/CHANGELOG to describe the new rule. Refactor hyperelastic Python bindings to unify single-point and batch handling for isochoric_invariants/isochoric_pstretch (per_tensor_point helper), and add tests for batch invariants, single-point J handling, and integrator refusal. Test binaries updated.
Introduce utilities and optimizations for continuum mechanics: add cubic_coefficients and is_isotropic to detect cubic/isotropic 6x6 Voigt tensors, and spectral_congruence and add_linearop_eigsym for efficient spectral operations. Refactor linearop_eigsym to use an in-place accumulator to avoid 6x6 temporaries and improve performance. Wire in mandel_factors-based isotropy check and include tensor.hpp where needed. Improve eshelby memoization/caching logic and remove fragile static Armadillo objects (avoid DLL unload issues on Windows) by hoisting comments and making tables/structures persistent safely. Miscellaneous cleanups: clarified docs/comments (including changelog), small performance and readability tweaks across Python wrappers and several UMAT, homogenization, micromechanics and objective-rate routines, and expose closed-form radial-return tangent description in headers.
Add an ANSYS USERMAT (UMANS) compatibility test and Fortran example, and wire it into CMake so the UMANS plugin is built for tests. Remove unused C++ UMAT kernels (damage_weibull, Prony_Nfast_Plastic) and add related test binaries. Update documentation (CHANGELOG, external/ and catalog) to document plugin small-strain limitation, test locations, and rheology naming; fix multiple API docs and Doxygen comments. Correct several API/implementation details and bugs: Eshelby penny shear components, kinematics/Cayley note, objective-rate API (added inverse prototype), material/variant ctor and build semantics, Chaboche state indexing, and various header throws/validation notes. Miscellaneous improvements and small refactors across Continuum_mechanics and Simulation sources and tests to align code, docs and tests with these changes.
Introduce Ansys USERMAT support and refactor external-plugin handling. - Add fea_components and FE<->simcoon conversion helpers (fea_transfer) to map component orders and convert vectors/matrices between Fortran/FE layout and simcoon Voigt order. - Implement umat_plugin_ans using the new transfer helpers; provide two plugin builds (normal and layout-inspecting) and update CMake to build plugins from external/ instead of testBin/ locations. - Consolidate external example sources under external/, remove duplicated testBin copies, and update docs and CHANGELOG describing the USERMAT behaviour and layout checks (TUMABA/TUMANS). - Add a generic plugin loader helper and switch the small-strain plugin dispatch to use dynamic plugin objects; mark external plugins (UMEXT/UMABA/UMANS) as small-strain-only and make the finite-strain dispatch refuse them with a clear error. - Prevent sim.umat from handling external plugins directly (Python wrapper now raises for external plugin names) and add a unit test to assert that behaviour. - Misc: cleanup copy/destructor boilerplate in step classes (defaulted destructor), reorganize fea_transfer/umat_smart implementations and update energy/work bookkeeping when using external USERMAT. These changes centralise FE/Fortran data mapping, enable Ansys USERMAT testing (including data-layout verification), and ensure external plugins are loaded and executed via the solver plugin API only.
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
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.
This pull request updates the benchmarking and non-regression tools for comparing legacy UMATs and their MODUL (modular UMAT) equivalents. The main improvements clarify which legacy UMATs still use their own kernels, introduce clearer command-line options to select which families to benchmark, and enhance documentation to explain the benchmarking process and current engine performance hotspots.
Enhancements to legacy vs. modular benchmarking:
EPICP,EPCHA, andPRONKstill run a true legacy kernel; all other legacy UMAT names are now routed to the modular engine via a properties adapter and are only benchmarked on request (flagged as "adapter" in results) (bench/bench_legacy_vs_modular.py, bench/bench_legacy_vs_modular.pyR3-R12).LEGACY_KERNELStuple to explicitly list the UMATs with real legacy kernels (bench/bench_legacy_vs_modular.py, bench/bench_legacy_vs_modular.pyR74-R77).Command-line interface improvements:
--allflag to include adapter families, and improved help messages to reflect these changes (bench/bench_legacy_vs_modular.py, bench/bench_legacy_vs_modular.pyL408-R445).bench/bench_legacy_vs_modular.py, bench/bench_legacy_vs_modular.pyL408-R445).Documentation and non-regression tools:
benchmark/README.mdto document all benchmarking and non-regression scripts, including new sections on batch API benchmarks, non-regression array dumps, and a measured update on previously identified engine hotspots (benchmark/README.md, [1] [2].These changes improve clarity for users running benchmarks and help maintainers focus on meaningful performance comparisons.