Perf/optimization v3

CHANGELOG.md

@@ -27,6 +27,10 @@ Known gaps and planned work, not yet in a release. An item moves to
 - The external plugins (UMEXT, UMABA, UMANS) serve the small-strain route only: a finite-strain
   entry in the plugin API (law fed the logarithmic strain, returning tau, as the internal
   kernels) would let them run under NLGEOM.
+- `tangent_assembly` takes per-mechanism `std::vector<arma::vec>` / `std::vector<arma::mat>`
+  (κ, ∂Φ/∂σ, ∂Λ/∂σ): every MODUL tangent copies the 6×6 flow Hessians to the heap, the
+  largest remaining allocation source in MODUL plasticity. Fixed-size or view-based arguments
+  need an API change of `assemble_*_tangent`.
 - Migration of the laws to the modular engine (MODUL), state in
   `benchmark/README.md`. Routed through the adapter (code 201): ELISO, ELIST,
   ELORT, EPKCP, EPHIL, EPTRI, EPHAC, EPANI, EPDFA, EPCHG, EPHIN. To do:

include/simcoon/Continuum_mechanics/Functions/fastor_bridge.hpp

@@ -98,8 +98,8 @@ inline arma::mat::fixed<3,3> fastor2_to_arma(const Fastor::Tensor<double,3,3> &T
  * or the type-free Mandel pair mandel_to_fastor4 / fastor4_to_mandel.
  */
 inline Fastor::Tensor<double,3,3,3,3> voigt_to_fastor4(const arma::mat::fixed<6,6> &L) {
+    // The four minor-symmetric writes below cover all 81 entries: no zero fill.
     Fastor::Tensor<double,3,3,3,3> C;
-    C.zeros();
 
     for (int i = 0; i < 3; ++i) {
         for (int j = i; j < 3; ++j) {
@@ -210,9 +210,9 @@ inline Fastor::Tensor<double,3,3,3,3> push_forward_4(
     const Fastor::Tensor<double,3,3,3,3> &C,
     const Fastor::Tensor<double,3,3> &F)
 {
+    // Every entry of each pass is written exactly once: no zero fill.
     // Pass 1: T1(i,J,M,N) = sum_L F(i,L) * C(L,J,M,N)
     Fastor::Tensor<double,3,3,3,3> T1;
-    T1.zeros();
     for (int i = 0; i < 3; ++i)
     for (int J = 0; J < 3; ++J)
     for (int M = 0; M < 3; ++M)
@@ -225,7 +225,6 @@ inline Fastor::Tensor<double,3,3,3,3> push_forward_4(
 
     // Pass 2: T2(i,s,M,N) = sum_J F(s,J) * T1(i,J,M,N)
     Fastor::Tensor<double,3,3,3,3> T2;
-    T2.zeros();
     for (int i = 0; i < 3; ++i)
     for (int s = 0; s < 3; ++s)
     for (int M = 0; M < 3; ++M)
@@ -238,7 +237,6 @@ inline Fastor::Tensor<double,3,3,3,3> push_forward_4(
 
     // Pass 3: T3(i,s,r,N) = sum_M F(r,M) * T2(i,s,M,N)
     Fastor::Tensor<double,3,3,3,3> T3;
-    T3.zeros();
     for (int i = 0; i < 3; ++i)
     for (int s = 0; s < 3; ++s)
     for (int r = 0; r < 3; ++r)
@@ -251,7 +249,6 @@ inline Fastor::Tensor<double,3,3,3,3> push_forward_4(
 
     // Pass 4: result(i,s,r,p) = sum_N F(p,N) * T3(i,s,r,N)
     Fastor::Tensor<double,3,3,3,3> result;
-    result.zeros();
     for (int i = 0; i < 3; ++i)
     for (int s = 0; s < 3; ++s)
     for (int r = 0; r < 3; ++r)

include/simcoon/Continuum_mechanics/Umat/Modular/damage_mechanism.hpp

@@ -107,6 +107,16 @@ class DamageMechanism final : public StrainMechanism {
     mutable bool M_cached_valid_;       ///< Whether the cached compliances are valid
     mutable arma::mat L_cached_;        ///< Stiffness the cached compliance was inverted from
     std::vector<tensor2> kappa_cache_{tensor2::zeros(Tensor2Type::strain)};   ///< zero: no multiplier carried
+    size_t D_index_ = 0;        ///< ivc_ index of "D", resolved by register_variables
+    size_t Y_max_index_ = 0;    ///< ivc_ index of "Y_max", resolved by register_variables
+
+    // last driving force and the effective stress it was evaluated at, for the cached compliance (reset with it)
+    mutable arma::vec::fixed<6> Y_sigma_;
+    mutable double Y_sigma_value_ = 0.0;
+    mutable bool Y_sigma_valid_ = false;
+
+    /// compute_driving_force(sigma_eff, M_cached_), memoised on the last stress
+    double driving_force(const arma::vec& sigma_eff) const;
 
 public:
     /**

include/simcoon/Continuum_mechanics/Umat/Modular/hardening.hpp

@@ -231,7 +231,8 @@ class KinematicHardening {
 
     /**
      * @brief Register internal variables for this hardening model
-     * @param ivc The owning mechanism's internal variable collection
+     * @param ivc The owning mechanism's internal variable collection; every other member taking a
+     *        collection must get this same one (read by the indices resolved here, debug-checked).
      */
     virtual void register_variables(InternalVariableCollection& ivc) = 0;
 
@@ -410,6 +411,7 @@ class PragerHardening final : public KinematicHardening {
 private:
     double C_;  ///< Kinematic hardening modulus
     std::string a_key_{"a"};  ///< IVC key for the back-strain α
+    size_t a_index_ = 0;      ///< its index in the collection, resolved by register_variables
 public:
     PragerHardening() : C_(0.0) {}
     void configure(const arma::vec& props, int& offset) override;
@@ -434,6 +436,7 @@ class ArmstrongFrederickHardening final : public KinematicHardening {
     double C_;  ///< Hardening parameter
     double D_;  ///< Dynamic recovery parameter
     std::string a_key_{"a"};  ///< IVC key for the back-strain α
+    size_t a_index_ = 0;      ///< its index in the collection, resolved by register_variables
 public:
     ArmstrongFrederickHardening() : C_(0.0), D_(0.0) {}
     void configure(const arma::vec& props, int& offset) override;
@@ -461,6 +464,7 @@ class ChabocheHardening final : public KinematicHardening {
     arma::vec C_;                    ///< Hardening parameters
     arma::vec D_;                    ///< Dynamic recovery parameters
     std::vector<std::string> a_keys_; ///< IVC keys for back-strains α_i ("a_0", "a_1", ...)
+    std::vector<size_t> a_index_;     ///< their indices in the collection, resolved by register_variables
 public:
     explicit ChabocheHardening(int N = 1) : N_(N), C_(arma::zeros(N)), D_(arma::zeros(N)) {}
     void configure(const arma::vec& props, int& offset) override;

include/simcoon/Continuum_mechanics/Umat/Modular/plasticity_mechanism.hpp

@@ -114,6 +114,9 @@ class PlasticityMechanism final : public StrainMechanism {
     mutable double H_total_{0.0};               ///< Hardening modulus (iso + kin)
     std::vector<ClosestPointIngredients> cpp_cache_{ClosestPointIngredients{}};  ///< built by refresh_state, read by closest_point_ingredients
 
+    /// Total backstress (engineering stress Voigt), summed in the X_branches_ buffer it overwrites (no allocation)
+    [[nodiscard]] arma::vec total_backstress_voigt() const;
+
 public:
     /**
      * @brief Constructor

include/simcoon/Continuum_mechanics/Umat/Modular/viscoelastic_mechanism.hpp

@@ -69,6 +69,8 @@ class ViscoelasticMechanism final : public StrainMechanism {
     // inside the FB hot loop).
     std::vector<std::string> ev_key_;   ///< "EV_<i>" per branch
     std::vector<std::string> v_key_;    ///< "v_<i>" per branch
+    std::vector<size_t> ev_index_;      ///< ivc_ index of "EV_<i>", resolved by register_variables
+    std::vector<size_t> v_index_;       ///< ivc_ index of "v_<i>", resolved by register_variables
 
     // Per-iteration caches (CCP: frozen flow direction from previous iteration)
     mutable std::vector<tensor2> kappa_t_;        ///< Zero fluxes: the rows carry no multiplier

include/simcoon/Continuum_mechanics/Umat/Modular/yield_criterion.hpp

@@ -55,6 +55,8 @@ class YieldCriterion {
 private:
     YieldType type_;
     arma::vec params_;      ///< Criterion-specific parameters
+    arma::mat P_;           ///< HILL / DFA / ANISOTROPIC: quadratic-form matrix of params_
+                            ///< (P_Hill / P_DFA / P_Ani), built once by configure_*
     bool configured_;
 
     /**
@@ -63,6 +65,18 @@ class YieldCriterion {
      */
     std::string type_string() const;
 
+    /**
+     * @brief equivalent_stress(sigma) and flow_direction(sigma) in one call (the quadratic criteria share P sigma).
+     */
+    void equivalent_stress_and_flow(const arma::vec& sigma, double& sigma_eq, arma::vec& flow) const;
+
+    /**
+     * @brief The same plus flow_hessian(sigma) (criteria with a flow Hessian only, else its exception).
+     */
+    void equivalent_stress_flow_hessian(const arma::vec& sigma, double& sigma_eq, arma::vec& flow,
+                                        arma::mat& hessian) const;
+    friend class PlasticityMechanism;
+
 public:
     /**
      * @brief Default constructor (von Mises)

include/simcoon/Simulation/Maths/small_dense.hpp

@@ -0,0 +1,85 @@
+/* This file is part of simcoon.
+
+ simcoon is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ simcoon is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with simcoon.  If not, see <http://www.gnu.org/licenses/>.
+
+ */
+
+///@file small_dense.hpp
+///@brief Fixed-order 6x6 products for the per-point kernels
+///@version 1.0
+
+#pragma once
+#include <armadillo>
+
+namespace simcoon{
+
+/** @addtogroup maths
+ *  @{
+ */
+
+/**
+ * @brief \f$ \mathbf{y} = \mathbf{A}\,\mathbf{x} \f$ for a column-major 6x6 @p A and a 6-vector @p x.
+ *
+ * Armadillo hands every product above 4x4 to BLAS; at this size the call costs more than the
+ * 36 multiply-adds, and the summation order depends on the BLAS backend. This loop has a fixed
+ * order (column by column), so the result is the same on every platform.
+ * @p y may alias @p x.
+ */
+inline void mul6(const double *A, const double *x, double *y) {
+    double r[6];
+    for (int i = 0; i < 6; ++i) r[i] = A[i]*x[0];
+    for (int k = 1; k < 6; ++k)
+        for (int i = 0; i < 6; ++i) r[i] += A[i + 6*k]*x[k];
+    for (int i = 0; i < 6; ++i) y[i] = r[i];
+}
+
+/// \f$ \mathbf{y} = \mathbf{A}\,\mathbf{x} \f$ for a 6x6 @p A and a 6-vector @p x (see the pointer form).
+inline void mul6(const arma::mat &A, const arma::vec &x, arma::vec &y) {
+    if (&y != &x) y.set_size(6);
+    mul6(A.memptr(), x.memptr(), y.memptr());
+}
+
+/// \f$ \mathbf{C} = \mathbf{A}\,\mathbf{B} \f$ for 6x6 operands, fixed summation order (see mul6).
+inline arma::mat::fixed<6,6> mul66(const arma::mat &A, const arma::mat &B) {
+    const double *a = A.memptr();
+    const double *b = B.memptr();
+    arma::mat::fixed<6,6> C;
+    double *c = C.memptr();
+    for (int j = 0; j < 6; ++j)
+        for (int i = 0; i < 6; ++i) {
+            double s = 0.0;
+            for (int k = 0; k < 6; ++k) s += a[i + 6*k]*b[k + 6*j];
+            c[i + 6*j] = s;
+        }
+    return C;
+}
+
+/// \f$ \mathbf{C} = \mathbf{A}\,\mathbf{B}^T \f$ for 6x6 operands, without forming the transpose.
+inline arma::mat::fixed<6,6> mul66_ABt(const arma::mat &A, const arma::mat &B) {
+    const double *a = A.memptr();
+    const double *b = B.memptr();
+    arma::mat::fixed<6,6> C;
+    double *c = C.memptr();
+    for (int j = 0; j < 6; ++j)
+        for (int i = 0; i < 6; ++i) {
+            double s = 0.0;
+            for (int k = 0; k < 6; ++k) s += a[i + 6*k]*b[j + 6*k];
+            c[i + 6*j] = s;
+        }
+    return C;
+}
+
+/** @} */ // end of maths group
+
+} //namespace simcoon

python-setup/simcoon/rotation.py

@@ -7,6 +7,7 @@
 methods (``apply_stress``, ``apply_stiffness``, etc.).
 """
 
+import math
 import warnings
 from typing import Dict, Sequence, Union
 
@@ -81,7 +82,16 @@ def from_mrp(cls, mrp):
 
     @classmethod
     def identity(cls, num=None):
-        return cls._ensure_cls(super().identity(num))
+        r = super().identity(num)
+        if type(r) is not cls:
+            # exact unit quaternions: nothing to renormalize through from_quat (the default
+            # orientation of every phase, built with each phase object)
+            r.__class__ = cls
+        if num is None:
+            # a single rotation cannot be modified in place: euler_angles reads its angles
+            # without decomposing it
+            r._exact_identity = True
+        return r
 
     @classmethod
     def random(cls, num=None, random_state=None):
@@ -501,17 +511,23 @@ def as_rotation(value: Orientation) -> Rotation:
         return value
     if isinstance(value, ScipyRotation):
         return Rotation.from_scipy(value)
+    return Rotation.from_euler(EULER_SEQ, _euler_triplet(value), degrees=True)
+
+
+def _euler_triplet(value):
+    """The Euler angles of a ``{"psi", "theta", "phi"}`` dict (missing angles are 0) or of a
+    3-sequence, in degrees, checked: unknown keys, a count other than 3 or a non-numeric value
+    raise ``ValueError``."""
     if isinstance(value, dict):
         unknown = set(value) - set(_ANGLES)
         if unknown:
             raise ValueError(f"orientation: unknown keys {sorted(unknown)}; expected {_ANGLES}")
-        angles = [float(value.get(k, 0.0)) for k in _ANGLES]
-    else:
-        angles = np.asarray(value, dtype=float).ravel()
-        if angles.size != 3:
-            raise ValueError(f"orientation: 3 Euler angles (psi, theta, phi) in degrees "
-                             f"expected, got {angles.size} values")
-    return Rotation.from_euler(EULER_SEQ, angles, degrees=True)
+        return [float(value.get(k, 0.0)) for k in _ANGLES]
+    angles = np.asarray(value, dtype=float).ravel()
+    if angles.size != 3:
+        raise ValueError(f"orientation: 3 Euler angles (psi, theta, phi) in degrees "
+                         f"expected, got {angles.size} values")
+    return angles
 
 
 def euler_angles(rotation: Orientation) -> Dict[str, float]:
@@ -525,11 +541,15 @@ def euler_angles(rotation: Orientation) -> Dict[str, float]:
     if not isinstance(rotation, ScipyRotation) and rotation is not None:
         # angles given as angles are written as given: no detour through a quaternion
         # (float noise, and a gimbal-locked triplet rewritten) for a no-op
-        as_rotation(rotation)   # validates the dict keys / the 3 values
-        values = ([float(rotation.get(k, 0.0)) for k in _ANGLES] if isinstance(rotation, dict)
-                  else [float(a) for a in np.asarray(rotation, dtype=float).ravel()])
+        values = [float(a) for a in _euler_triplet(rotation)]   # validates the dict keys / the 3 values
+        if not all(math.isfinite(a) for a in values):
+            as_rotation(rotation)   # the error of a rotation built from non-finite angles
         return dict(zip(_ANGLES, values))
     rot = as_rotation(rotation)
+    if getattr(rot, "_exact_identity", False):
+        # Rotation.identity() (the default orientation of every phase): the angles scipy
+        # returns for it, +0, without the decomposition
+        return {'psi': 0.0, 'theta': 0.0, 'phi': 0.0}
     with warnings.catch_warnings():
         warnings.simplefilter('ignore')   # scipy's "Gimbal lock detected" - see above
         psi, theta, phi = rot.as_euler(EULER_SEQ, degrees=True)

simcoon-python-builder/src/python_wrappers/Libraries/Continuum_mechanics/pyumat.cpp

@@ -51,6 +51,26 @@ struct step_cut_tag {};
 PYBIND11_CONSTINIT static py::gil_safe_call_once_and_store<py::object> g_step_cut;
 static py::object *g_fn = nullptr;   // touched only with the GIL held
 
+// The keyword names of the call, interned once (same lifetime rule as g_step_cut). The call
+// passes them with the values in one vector (PyObject_Vectorcall): building a kwargs dict
+// key by key, as a pybind11 keyword call does, was a third of the bridge.
+constexpr const char *k_kwnames[] = {"Etot", "DEtot", "sigma", "DR", "props", "statev", "T", "DT",
+                                     "Time", "DTime", "Wm", "ndi", "nshr", "start", "tangent_mode"};
+constexpr std::size_t k_nkw = sizeof(k_kwnames) / sizeof(k_kwnames[0]);
+PYBIND11_CONSTINIT static py::gil_safe_call_once_and_store<py::tuple> g_kwnames;
+
+const py::tuple &kwnames() {
+    return g_kwnames.call_once_and_store_result([]() {
+        py::tuple names(k_nkw);
+        for (std::size_t k = 0; k < k_nkw; k++) {
+            PyObject *s = PyUnicode_InternFromString(k_kwnames[k]);
+            if (!s) throw py::error_already_set();
+            names[k] = py::reinterpret_steal<py::str>(s);
+        }
+        return names;
+    }).get_stored();
+}
+
 using farr = py::array_t<double, py::array::f_style | py::array::forcecast>;
 
 // (rows, cols) column-major copy — cheaper than simpy::arma_to_numpy::mat_to_arr(arma::mat(m)), which
@@ -126,20 +146,28 @@ void python_umat_bridge(const std::string &umat_name, const arma::vec &Etot, con
     if (!g_fn || g_fn->is_none()) {
         throw simcoon::exception_solver("PYEXT: no Python UMAT registered");
     }
+    // own reference for the call: a law that unregisters or replaces itself must stay alive
+    const py::object fn = *g_fn;
 
     py::array_t<double> Wm4(4);
     {
         double *w = Wm4.mutable_data();
         w[0] = Wm; w[1] = Wm_r; w[2] = Wm_ir; w[3] = Wm_d;
     }
+    // keyword arguments, in the order of k_kwnames
+    const py::object values[k_nkw] = {np1d(Etot), np1d(DEtot), np1d(sigma), np2d(DR), np1d(props),
+                                      np1d(statev), py::float_(T), py::float_(DT), py::float_(Time),
+                                      py::float_(DTime), Wm4, py::int_(ndi), py::int_(nshr),
+                                      py::bool_(start), py::int_(tangent_mode)};
     py::object ret;
     try {
-        ret = (*g_fn)("Etot"_a = np1d(Etot), "DEtot"_a = np1d(DEtot), "sigma"_a = np1d(sigma),
-                      "DR"_a = np2d(DR), "props"_a = np1d(props),
-                      "statev"_a = np1d(statev),
-                      "T"_a = T, "DT"_a = DT, "Time"_a = Time, "DTime"_a = DTime,
-                      "Wm"_a = Wm4, "ndi"_a = ndi, "nshr"_a = nshr, "start"_a = start,
-                      "tangent_mode"_a = tangent_mode);
+        PyObject *argv[k_nkw + 1];   // argv[0]: the slot PY_VECTORCALL_ARGUMENTS_OFFSET lends the callee
+        argv[0] = nullptr;
+        for (std::size_t k = 0; k < k_nkw; k++) argv[k + 1] = values[k].ptr();
+        PyObject *res = PyObject_Vectorcall(fn.ptr(), argv + 1, PY_VECTORCALL_ARGUMENTS_OFFSET,
+                                            kwnames().ptr());
+        if (!res) throw py::error_already_set();
+        ret = py::reinterpret_steal<py::object>(res);
     } catch (py::error_already_set &e) {
         if (e.matches(g_step_cut.get_stored())) {
             // Requested step cut: leave sigma/statev/Wm untouched (the solver restores the