Feature/modular default

.github/workflows/build.yml

@@ -54,14 +54,10 @@ jobs:
       - name: Editable install
         run: |
           uv venv
-          # scikit-build-core pinned <1: the 1.x editable install changed where
-          # the compiled extension lands, which breaks the Windows DLL-vendoring
-          # step below (green on 0.12.2, DLL load failure on 1.0.2). Lift the pin
-          # together with a rework of that step. --no-build-isolation so the
-          # pinned backend in the venv is the one that actually builds the
-          # editable wheel (isolation would re-resolve pyproject's
-          # scikit-build-core>=0.11 to 1.x).
-          uv pip install "scikit-build-core<1" pybind11 numpy
+          # --no-build-isolation so the backend in the venv is the one that
+          # actually builds the editable wheel (and so later direct cmake
+          # rebuilds share its build tree).
+          uv pip install scikit-build-core pybind11 numpy
           uv pip install -e .[dev] --no-build-isolation
 
       - name: Vendor conda BLAS/LAPACK runtime next to the extension (Windows)
@@ -70,13 +66,18 @@ jobs:
           # Trigger the lazy scikit-build-core editable build so the extension and its
           # install dir exist (the import itself fails on the missing BLAS DLL -- expected).
           uv run python -c "import simcoon" 2>&1 | Out-Null
-          # simcoon/__init__.py adds its own package dir to the DLL search
-          # (os.add_dll_directory). Copy conda's reference BLAS/LAPACK + Fortran runtime
-          # there so _core resolves its transitive deps under the uv venv: unlike the conda
-          # Python used by ci.yml, the venv does not see C:\Miniconda\Library\bin, and
-          # Python 3.8+ ignores PATH for extension-module dependencies.
-          $dst = ".venv\Lib\site-packages\simcoon"
-          Copy-Item "$env:CONDA\Library\bin\*.dll" $dst -Force -ErrorAction SilentlyContinue
+          # Copy conda's reference BLAS/LAPACK + Fortran runtime NEXT TO _core.pyd:
+          # Windows resolves an extension module's dependencies from the pyd's own
+          # directory (PATH is ignored since Python 3.8, and the uv venv does not see
+          # C:\Miniconda\Library\bin like ci.yml's conda Python does). The pyd's
+          # location depends on the scikit-build-core editable layout (site-packages
+          # on 0.x, redirected to <build-dir>\install\platlib\simcoon on 1.x), so
+          # locate it instead of hardcoding either.
+          $pyds = Get-ChildItem -Path build, .venv\Lib\site-packages -Recurse -Filter "_core*.pyd" -ErrorAction SilentlyContinue
+          if (-not $pyds) { Write-Error "editable _core.pyd not found"; exit 1 }
+          foreach ($pyd in $pyds) {
+            Copy-Item "$env:CONDA\Library\bin\*.dll" $pyd.DirectoryName -Force -ErrorAction SilentlyContinue
+          }
           exit 0
 
       - name: Verify cmake rebuild works (Unix)

CMakeLists.txt

@@ -370,13 +370,6 @@ if(WIN32 AND BUILD_PYTHON_BINDINGS)
   endif()
 endif()
 
-# Compile public executables (with main functions)
-set(EXECUTABLES solver L_eff Elastic_props ODF PDF)
-foreach(exe_name ${EXECUTABLES})
-  add_executable(${exe_name} software/${exe_name}.cpp)
-  target_link_libraries(${exe_name} PRIVATE simcoon ${CMAKE_DL_LIBS})
-endforeach()
-
 ##Testing
 if(SIMCOON_BUILD_TESTS)
   # Test plugin libraries

bench/bench_legacy_vs_modular.py

@@ -0,0 +1,202 @@
+"""Benchmark legacy UMATs against their MODUL (modular UMAT) twin configurations.
+
+ZENER/ZENNK are intentionally absent: Zener_fast/Zener_Nfast are generalized
+KELVIN chains (branch driving force sigma - L_i EV_i), a different rheological
+model from the modular Prony (generalized Maxwell) — no twin exists.
+
+Methodology: INTERLEAVED A/B — each repetition runs the legacy kernel then the
+modular twin back-to-back, so load spikes hit both arms equally; the reported
+estimator is the MINIMUM over repetitions (least-noise wall time; naive medians
+proved 5x sensitive to background load). A max-deviation column doubles as an
+equivalence smoke check.
+
+Usage:
+    python bench/bench_legacy_vs_modular.py [--families ELISO,EPCHA,...]
+                                            [--repeats 8]
+
+Run from anywhere; paths resolve relative to the repo. Requires the simcoon
+python package (editable install). Results land in a temp dir and are removed.
+"""
+import argparse
+import os
+import statistics
+import tempfile
+import time
+from pathlib import Path
+
+import numpy as np
+
+import simcoon as sim
+from simcoon.modular import (
+    ModularMaterial, IsotropicElasticity, CubicElasticity,
+    TransverseIsotropicElasticity, OrthotropicElasticity,
+    Plasticity, VoceHardening, PowerLawHardening,
+    HillYield, DFAYield, AnisotropicYield,
+    PragerHardening, ChabocheHardening, Viscoelasticity,
+)
+
+REPO = Path(__file__).resolve().parents[1]
+EXAMPLES = REPO / "examples" / "mechanical"
+
+VE_TERMS = [(1500., 0.35, 3000., 1200.),
+            (800., 0.35, 30000., 12000.),
+            (400., 0.35, 300000., 120000.)]
+
+
+def _mat(elasticity, mechanisms=()):
+    return ModularMaterial(elasticity=elasticity, mechanisms=list(mechanisms))
+
+
+# family -> (legacy_name, legacy_props, legacy_nstatev, path_file, ModularMaterial)
+FAMILIES = {
+    "ELISO": ("ELISO", [210000., 0.3, 1.2e-5], 1, "MODUL_path.txt",
+              _mat(IsotropicElasticity(C1=210000., C2=0.3, alpha=1.2e-5))),
+    "ELIST": ("ELIST", [3, 230000., 15000., 0.02, 0.4, 50000., 0., 0.], 1,
+              "MODUL_path.txt",
+              _mat(TransverseIsotropicElasticity(EL=230000., ET=15000.,
+                                                 nuTL=0.02, nuTT=0.4,
+                                                 GLT=50000., axis=3))),
+    "ELORT": ("ELORT", [70000., 30000., 15000., 0.3, 0.3, 0.3,
+                        8000., 6000., 5000., 1e-5, 2e-5, 3e-5], 1,
+              "MODUL_path.txt",
+              _mat(OrthotropicElasticity(C1=70000., C2=30000., C3=15000.,
+                                         C4=0.3, C5=0.3, C6=0.3,
+                                         C7=8000., C8=6000., C9=5000.,
+                                         alpha1=1e-5, alpha2=2e-5, alpha3=3e-5))),
+    "EPICP": ("EPICP", [210000., 0.3, 0., 300., 1000., 0.3], 8, "MODUL_path.txt",
+              _mat(IsotropicElasticity(C1=210000., C2=0.3),
+                   [Plasticity(sigma_Y=300.,
+                               isotropic_hardening=PowerLawHardening(k=1000., m=0.3))])),
+    "EPKCP": ("EPKCP", [210000., 0.3, 0., 300., 1000., 1.0, 20000.], 33,
+              "MODUL_path.txt",
+              _mat(IsotropicElasticity(C1=210000., C2=0.3),
+                   [Plasticity(sigma_Y=300.,
+                               isotropic_hardening=PowerLawHardening(k=1000., m=1.0),
+                               kinematic_hardening=PragerHardening(C=1.5 * 20000.))])),
+    "EPCHA": ("EPCHA", [210000., 0.3, 0., 300., 200., 20.,
+                        30000., 172., 19500., 301.], 33, "MODUL_path.txt",
+              _mat(IsotropicElasticity(C1=210000., C2=0.3),
+                   [Plasticity(sigma_Y=300.,
+                               isotropic_hardening=VoceHardening(Q=200., b=20.),
+                               kinematic_hardening=ChabocheHardening(
+                                   terms=((30000., 172.), (19500., 301.))))])),
+    "EPHIL": ("EPHIL", [210000., 0.3, 0., 300., 5000., 1.0,
+                        0.5, 0.4, 0.6, 1.5, 1.5, 1.5], 33, "MODUL_path.txt",
+              _mat(IsotropicElasticity(C1=210000., C2=0.3),
+                   [Plasticity(sigma_Y=300.,
+                               yield_criterion=HillYield(F=0.5, G=0.4, H=0.6,
+                                                         L=1.5, M=1.5, N=1.5),
+                               isotropic_hardening=PowerLawHardening(k=5000., m=1.0))])),
+    "EPHAC": ("EPHAC", [210000., 0.3, 85000., 0., 300., 200., 20.,
+                        30000., 172., 19500., 301.,
+                        0.5, 0.4, 0.6, 1.5, 1.5, 1.5], 33, "MODUL_path.txt",
+              _mat(CubicElasticity(C1=210000., C2=0.3, C3=85000.),
+                   [Plasticity(sigma_Y=300.,
+                               yield_criterion=HillYield(F=0.5, G=0.4, H=0.6,
+                                                         L=1.5, M=1.5, N=1.5),
+                               isotropic_hardening=VoceHardening(Q=200., b=20.),
+                               kinematic_hardening=ChabocheHardening(
+                                   terms=((30000., 172.), (19500., 301.))))])),
+    "EPANI": ("EPANI", [210000., 0.3, 85000., 0., 300., 200., 20.,
+                        30000., 172., 19500., 301.,
+                        1.2, 1.1, 1.1, -0.6, -0.6, -0.5, 1.6, 1.5, 1.4], 33,
+              "MODUL_path.txt",
+              _mat(CubicElasticity(C1=210000., C2=0.3, C3=85000.),
+                   [Plasticity(sigma_Y=300.,
+                               yield_criterion=AnisotropicYield(
+                                   P11=1.2, P22=1.1, P33=1.1, P12=-0.6,
+                                   P13=-0.6, P23=-0.5, P44=1.6, P55=1.5,
+                                   P66=1.4),
+                               isotropic_hardening=VoceHardening(Q=200., b=20.),
+                               kinematic_hardening=ChabocheHardening(
+                                   terms=((30000., 172.), (19500., 301.))))])),
+    "EPDFA": ("EPDFA", [210000., 0.3, 85000., 0., 300., 200., 20.,
+                        30000., 172., 19500., 301.,
+                        0.5, 0.4, 0.6, 1.5, 1.5, 1.5, 0.1], 33, "MODUL_path.txt",
+              _mat(CubicElasticity(C1=210000., C2=0.3, C3=85000.),
+                   [Plasticity(sigma_Y=300.,
+                               yield_criterion=DFAYield(F=0.5, G=0.4, H=0.6,
+                                                        L=1.5, M=1.5, N=1.5,
+                                                        K=0.1),
+                               isotropic_hardening=VoceHardening(Q=200., b=20.),
+                               kinematic_hardening=ChabocheHardening(
+                                   terms=((30000., 172.), (19500., 301.))))])),
+    "EPCHG": ("EPCHG", [210000., 0.3, 85000., 0., 300., 2, 2, 0,
+                        150., 15., 50., 40., 30000., 172., 19500., 301.], 33,
+              "MODUL_path.txt",
+              _mat(CubicElasticity(C1=210000., C2=0.3, C3=85000.),
+                   [Plasticity(sigma_Y=300.,
+                               # legacy N-term iso couples through a single Hp
+                               # -> ONE effective Voce (b_eff=sum b_i,
+                               # Q_eff=sum(b_i Q_i)/sum(b_i))
+                               isotropic_hardening=VoceHardening(
+                                   Q=(15. * 150. + 40. * 50.) / 55., b=55.),
+                               kinematic_hardening=ChabocheHardening(
+                                   terms=((30000., 172.), (19500., 301.))))])),
+    # EPHIN legacy is defective for N >= 2 (NaN even for identical/inactive
+    # second surface) — bench the provable N = 1 case.
+    "EPHIN": ("EPHIN", [210000., 0.3, 0., 1,
+                        300., 3000., 1.0, 0.5, 0.4, 0.6, 1.5, 1.5, 1.5], 33,
+              "MODUL_path.txt",
+              _mat(IsotropicElasticity(C1=210000., C2=0.3),
+                   [Plasticity(sigma_Y=300.,
+                               yield_criterion=HillYield(F=0.5, G=0.4, H=0.6,
+                                                         L=1.5, M=1.5, N=1.5),
+                               isotropic_hardening=PowerLawHardening(k=3000., m=1.0))])),
+    "PRONK": ("PRONK", [3000., 0.35, 0., 3] + [x for t in VE_TERMS for x in t],
+              7 + 7 * 3, "PRONK_path.txt",
+              _mat(IsotropicElasticity(C1=3000., C2=0.35),
+                   [Viscoelasticity(terms=VE_TERMS)])),
+}
+
+
+def bench_family(key, results_dir, repeats):
+    name, props, nstatev, path_file, mat = FAMILIES[key]
+    props = np.asarray(props, dtype=float)
+    t_leg, t_mod = [], []
+    out_l, out_m = f"b_{key}_l.txt", f"b_{key}_m.txt"
+    for _ in range(repeats):
+        t0 = time.perf_counter()
+        sim.solver(name, props, nstatev, 0., 0., 0., 0, 1,
+                   "../data", results_dir, path_file, out_l)
+        t_leg.append(time.perf_counter() - t0)
+        t0 = time.perf_counter()
+        sim.solver("MODUL", mat.props, mat.nstatev, 0., 0., 0., 0, 1,
+                   "../data", results_dir, path_file, out_m)
+        t_mod.append(time.perf_counter() - t0)
+    # discard the first (warmup) pair; min = least-noise estimator
+    tl, tm = min(t_leg[1:]), min(t_mod[1:])
+    res = Path(EXAMPLES) / results_dir
+    a = np.loadtxt(res / out_l.replace(".txt", "_global-0.txt"))[:, 14]
+    b = np.loadtxt(res / out_m.replace(".txt", "_global-0.txt"))[:, 14]
+    dev = np.max(np.abs(a - b)) / max(np.max(np.abs(a)), 1e-12)
+    return tl * 1e3, tm * 1e3, dev
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--families", default=",".join(FAMILIES),
+                    help="comma-separated subset of: " + ",".join(FAMILIES))
+    ap.add_argument("--repeats", type=int, default=8)
+    args = ap.parse_args()
+
+    os.chdir(EXAMPLES)
+    with tempfile.TemporaryDirectory() as tmp:
+        link = EXAMPLES / f"_bench_{os.getpid()}"
+        link.symlink_to(tmp, target_is_directory=True)
+        try:
+            print(f"| family | legacy [ms] | MODUL [ms] | ratio | max dev |")
+            print(f"|---|---|---|---|---|")
+            for key in args.families.split(","):
+                key = key.strip()
+                if key not in FAMILIES:
+                    print(f"| {key} | (unknown family) | | | |")
+                    continue
+                tl, tm, dev = bench_family(key, link.name, args.repeats)
+                print(f"| {key} | {tl:.1f} | {tm:.1f} | {tm/tl:.2f} | {dev:.2e} |")
+        finally:
+            link.unlink()
+
+
+if __name__ == "__main__":
+    main()

docs/simulation/abaqus.rst

@@ -281,40 +281,56 @@ The following constitutive models are available through ``select_umat_M()``:
 
    * - Code
      - Model
-     - Properties
+     - Properties (in order)
      - State Variables
    * - ELISO
      - Isotropic elasticity
      - E, ν, α
-     - 1 (start flag)
+     - 1
    * - ELIST
      - Transversely isotropic elasticity
-     - E₁, E₂, ν₁₂, ν₂₃, G₁₂, α₁, α₂
+     - axis, EL, ET, νTL, νTT, GLT, αL, αT
      - 1
    * - ELORT
      - Orthotropic elasticity
      - E₁, E₂, E₃, ν₁₂, ν₁₃, ν₂₃, G₁₂, G₁₃, G₂₃, α₁, α₂, α₃
      - 1
    * - EPICP
-     - Isotropic plasticity (isotropic hardening)
-     - E, ν, α, σ_y, H
-     - 8 (p, Hp, ...)
+     - Von Mises plasticity, power-law isotropic hardening
+     - E, ν, α, σ_Y, k, m
+     - 8 (T_init, p, EP)
    * - EPKCP
-     - Kinematic + isotropic hardening
-     - E, ν, α, σ_y, H, C, γ
-     - 14
+     - Von Mises, power-law isotropic + Prager kinematic
+     - E, ν, α, σ_Y, k, m, kX
+     - 14 (T_init, p, EP, a)
    * - EPCHA
-     - Chaboche cyclic plasticity
-     - E, ν, α, σ_y, Q, b, C₁, γ₁, ...
-     - 14+
-   * - EPHIL
-     - Hill anisotropic plasticity (iso hardening)
-     - E, ν, α, σ_y, H, F, G, H, L, M, N
-     - 8
+     - Von Mises + Voce + 2× Armstrong-Frederick
+     - E, ν, α, σ_Y, Q, b, C₁, D₁, C₂, D₂
+     - 33
+   * - EPHIL / EPTRI
+     - Hill yield + power-law isotropic hardening
+     - E, ν, α, σ_Y, k, m, F, G, H, L, M, N
+     - 8 (T_init, p, EP)
    * - EPHAC
-     - Hill + Chaboche
-     - E, ν, α, σ_y, Q, b, C₁, γ₁, F, G, H, L, M, N
-     - 14+
+     - Cubic elasticity + Hill + Voce + 2× AF
+     - E, ν, G, α, σ_Y, Q, b, C₁, D₁, C₂, D₂, F, G, H, L, M, N
+     - 33
+   * - EPANI
+     - Cubic elasticity + anisotropic yield + Voce + 2× AF
+     - E, ν, G, α, σ_Y, Q, b, C₁, D₁, C₂, D₂, P₁₁..P₆₆ (9)
+     - 33
+   * - EPDFA
+     - Cubic elasticity + DFA yield + Voce + 2× AF
+     - E, ν, G, α, σ_Y, Q, b, C₁, D₁, C₂, D₂, F, G, H, L, M, N, K
+     - 33
+   * - EPCHG
+     - Generic Chaboche (selectable yield, N iso/kin terms)
+     - E, ν, G, α, σ_Y, N_iso, N_kin, criteria, (Q,b)×N, (C,D)×N, crit. params
+     - 33
+   * - EPHIN
+     - N Hill yield surfaces
+     - E, ν, α, N, per surface: σ_Y, k, m, F, G, H, L, M, N
+     - 1 + 7N
    * - SMAUT
      - SMA unified model
      - See SMA documentation
@@ -324,21 +340,29 @@ The following constitutive models are available through ``select_umat_M()``:
      - See SMA documentation
      - 24
    * - LLDM0
-     - Lemaitre-Chaboche damage
-     - E, ν, α, σ_y, H, S, s, D_c
+     - Lemaitre-Ladeveze-Dufailly damage
+     - See header documentation
      - 9
    * - ZENER
-     - Zener viscoelastic (single branch)
-     - E₀, E₁, η
-     - 7
+     - Kelvin viscoelastic (single branch)
+     - E₀, ν₀, α, E₁, ν₁, ηB₁, ηS₁
+     - 14
    * - ZENNK
-     - Zener viscoelastic (N branches)
-     - E₀, E₁, η₁, E₂, η₂, ...
-     - 1+6N
+     - Kelvin viscoelastic (N branches)
+     - E₀, ν₀, α, N, per branch: Eᵢ, νᵢ, ηBᵢ, ηSᵢ
+     - 7 + 7N
    * - PRONK
-     - Prony series viscoelastic
-     - G₀, K₀, g₁, τ₁, k₁, τ'₁, ...
-     - 1+12N
+     - Prony series viscoelastic (generalized Maxwell)
+     - E₀, ν₀, α, N, per branch: Eᵢ, νᵢ, ηBᵢ, ηSᵢ
+     - 7 + 7N
+   * - MODUL
+     - Composable modular UMAT
+     - self-describing stream (see :mod:`simcoon.modular`)
+     - model-dependent
+
+Several of these names are served by the modular engine through
+props-translating adapters — identical usage and results; see
+:doc:`umat_catalog` for per-name status and state-variable layout notes.
 
 .. list-table:: Micromechanics Models
    :header-rows: 1