Skip to content

Perf/micro optim - #150

Open
chemiskyy wants to merge 11 commits into
masterfrom
perf/micro-optim
Open

chemiskyy wants to merge 11 commits into
masterfrom
perf/micro-optim

Conversation

@chemiskyy

Copy link
Copy Markdown
Member

This pull request updates the benchmarking and non-regression tooling for the solver project, clarifying the current state of legacy and modular UMAT kernels, and improving usability and documentation. The most important changes are:

Documentation and Benchmarking Workflow Updates:

  • The benchmark/README.md has been significantly expanded to clarify the purpose and workflow of each benchmarking and non-regression script, including detailed instructions and tables tracking the migration of laws to the modular engine. It now documents the state of previously identified C++ engine hotspots and their resolution status. [1] [2]

Legacy vs Modular Kernel Benchmarking Improvements:

  • The script bench/bench_legacy_vs_modular.py has been renamed to benchmark/bench_legacy_vs_modular.py and enhanced to:
    • Clearly distinguish between legacy kernels and modular engine adapters, both in code and in output tables.
    • Add a new --all option to benchmark all families, including those routed to the modular engine via adapters, and flag them accordingly. [1] [2] [3] [4]
    • Default to benchmarking only the true legacy kernels (EPICP, EPCHA, PRONK), improving the focus and clarity of the results.
    • Output an improved results table with a new "legacy route" column indicating whether a family uses a legacy kernel or an adapter.

These changes make benchmarking more transparent, reproducible, and easier to interpret, especially in the context of the ongoing migration to the modular engine.

The micro-optimisation commits of perf/micro-optim (01-08) rebased onto the
2.2 cleanup, with their comments trimmed to the house style:

- Trim temporaries from the shared tensor kernels
- Skip masked tangent products and heap temporaries
- Remove repeated work from mechanical kernels
- Build anisotropic criteria operators on the stack
- Evaluate SMA criteria once per local iteration
- Drop duplicate flow evaluations in thermal J2 laws
- Trim the modular engine's per-iteration overhead
- Rebuild finite-strain kernels on the stack

Results bit-identical to the base (solver and batch dumps). Serial us/pt
against the base: MODUL plastic 4.74 -> 3.06, EPKCP 5.76 -> 3.60, EPICP
1.93 -> 1.55, MODUL elastic 0.37 -> 0.31, NEOHC 2.76 -> 1.30, fedoo UL
MODUL 7.3 -> 5.1.
The accumulation pass of Gauss now holds the 36 sums of G as 18 pairs
of rows in two-lane values (a GCC/Clang vector type, a plain pair
elsewhere), so that they stay in registers and every operation serves
two entries. Each lane goes through the scalar expressions, contracted
as before, and the sum over the points of a ring stays sequential, so
G is bitwise unchanged with Apple clang (conda-forge clang 22 contracts the former
scalar loop differently: round-off, at most 2e-15 relative); the 50x50 integration drops from 36 to 25 us.
The contraction of G in Eshelby and T_II runs on plain arrays instead
of a bounds-checked Mat<int> table.
Mori-Tanaka keeps the last inverse of its normalisation tensor, which an
elastic medium hands back bitwise at every pass, and builds its sums in
fixed-size matrices; the periodic-layer scheme builds the Voigt
operators of each layer once per call instead of at each use (same
operators, same products, no identity shortcut); the loops
over the sub-phases take them by reference instead of copying them.
Results are bitwise unchanged; an elastic Mori-Tanaka increment is 1.5
times faster, a laminate with a plastic layer 1.2 times.
to_start and set_start of a phase dispatch through the virtual methods
instead of casting shared pointers, local2global and global2local cast
references (no atomic reference counts on every constitutive call), and
the state variables skip the rotation calls when the material frame is
the identity. The sub-phase outputs iterate by reference. Results are
bitwise unchanged.
Every attempt of an increment starts from the start state: set_start has
just accepted the previous increment, or assess_inc has rolled back a cut
one. The to_start that preceded each constitutive call of the Newton loops
therefore copied identical values on the first trial; it now runs from the
second trial on, and the single call of the fully strain-driven mechanical path no
longer rolls back (the thermomechanical one keeps its to_start). The invariant is documented where it is relied on and
asserted in debug builds. Lt_2_K copies the stress-controlled rows element by
element instead of through row subviews.

Results are bitwise unchanged. An elastic increment is about 7 percent
faster, a logarithmic-control one about 2.5 percent.
solver_run recorded through solver_memory_sink, which builds one arma
object per field and per record (a heap block for every 6x6 tangent),
then copied them row by row into the arrays returned to Python. The
binding now records into its own sink: one column-major (N, n) buffer per
field, sized for the increments of a complete run and filled in place, so
that the buffer itself becomes the NumPy array; it grows if need be and is
cut to the records of an aborted run. Keys, shapes, memory order and
values of the result are unchanged; solver_memory_sink stays as it was for
the C++ API. A new test pins the recorded history to the contract of
solver_memory_sink (one record per increment, sigma = tau / det(F1),
fields per state type, layout).

An elastic increment is about 30 percent faster, a plastic EPICP or a
logarithmic-control one about 10 percent, and the peak memory of a
1e6-increment run drops from 3.2 to 1.0 GiB.
The PYEXT bridge called the registered law with fifteen keyword
arguments through pybind11, which builds a kwargs dict key by key at every
call: a third of the bridge. It now uses PyObject_Vectorcall with a tuple
of the keyword names interned once. The law receives the same values in
the same order, and the bridge now holds its own reference to the law for
the duration of the call, which may unregister or replace it.

The bridge and solver overhead per call drops from 4.4 to 1.9 us; an
increment of a trivial Python law from 6.8 to 4.1 us.
The Python side passes Fortran-ordered cubes so that the binding can
borrow them, but the batch entries copied them, and copied and transposed
the (N, 6) Voigt arrays. Read-only cubes are now borrowed, and a
C-ordered (N, 6) array is read as the column-major (6, N) matrix it is;
other layouts are still copied. Outputs are unchanged.

Batch Tensor2 rotations, push-forwards and pull-backs are 26 to 34
percent faster, Tensor4 ones 8 to 11 percent.
euler_angles, called by every solve and for every phase handed to L_eff
or solve, validated Euler angles by building a scipy rotation it then
discarded, and decomposed the identity of every default phase
orientation. Angles are now checked directly, and Rotation.identity()
switches the class of scipy's exact unit quaternions instead of a
renormalizing round trip, and marks a single identity so that its angles
are returned without decomposition. Results and errors are unchanged.

euler_angles drops from 15 to 1.2 us for given angles, a phase object
builds in 3.6 instead of 15 us, a two-phase Mori-Tanaka L_eff takes 17
instead of 37 us, and every solve call saves about 14 us.
Armadillo hands every product above 4x4 to BLAS/LAPACK, whose call costs
more than the arithmetic at this size and whose summation order depends on
the backend:

- assemble_algorithmic_tangent: M^-1 [L | kappa] through small_lu<6>
  instead of a LAPACK inverse; L * sum(Ds dLambda/dsigma) through mul66
- assemble_tangent: the P / Lt products as fixed-order loops
- ElasticityModule::evaluate and PlasticityMechanism::compute_constraints:
  the 6x6 matrix-vector products through mul6
- small_dense.hpp gathers mul6 and the mul66 / mul66_ABt of tensor.cpp
- compute_tangent reserves its per-mechanism vectors

Not bit-identical: Lt within 1.5e-14 relative in the batch dumps, solver
paths within 3e-14 of each array's scale (Q of thermomechanical SMA_T 4e-11,
a small difference of terms of order 30). Serial us/pt: MODUL plastic
3.06 -> 2.16, EPKCP 3.60 -> 2.70, MODUL elastic 0.31 -> 0.28.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant