diff --git a/bench/bench_legacy_vs_modular.py b/bench/bench_legacy_vs_modular.py
index b83da6085..435dd729c 100644
--- a/bench/bench_legacy_vs_modular.py
+++ b/bench/bench_legacy_vs_modular.py
@@ -157,11 +157,11 @@ def bench_family(key, results_dir, repeats):
     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,
+        sim._core.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,
+        sim._core.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
diff --git a/docs/cpp_api/continuum_mechanics_overview.md b/docs/cpp_api/continuum_mechanics_overview.md
index 5db1aea97..00141ccab 100644
--- a/docs/cpp_api/continuum_mechanics_overview.md
+++ b/docs/cpp_api/continuum_mechanics_overview.md
@@ -60,15 +60,12 @@ User Material subroutines (UMAT) for finite element analysis, organized by strai
 - Isotropic, Orthotropic, Transversely Isotropic
 
 **Plasticity:**
-- `plastic_isotropic_ccp` - J2 plasticity with isotropic hardening
-- `plastic_kin_iso_ccp` - Combined kinematic-isotropic hardening
-- `plastic_chaboche_ccp` - Chaboche multi-kinematic hardening model
-- `Hill_isoh` - Hill anisotropic plasticity with isotropic hardening
-- `Hill_isoh_Nfast` - Hill plasticity with multiple hardening mechanisms
-- `Hill_chaboche_ccp` - Hill plasticity with Chaboche hardening
-- `Ani_chaboche_ccp` - General anisotropic plasticity with Chaboche hardening
-- `DFA_chaboche_ccp` - Distortion-based anisotropy with Chaboche hardening
-- `Generic_chaboche_ccp` - Generic anisotropic plasticity framework
+- `plastic_isotropic_ccp` - J2 plasticity with isotropic hardening (kept as a readable reference implementation)
+- `plastic_chaboche_ccp` - Chaboche multi-kinematic hardening model (kept as a readable reference implementation)
+- All other plasticity variants (kinematic-isotropic `EPKCP`, Hill `EPHIL`/`EPTRI`/`EPHIN`,
+  Hill-Chaboche `EPHAC`, anisotropic Chaboche `EPANI`/`EPDFA`/`EPCHG`) are served since 2.0
+  by the **modular UMAT engine** through `umat_legacy_modular` name adapters — legacy names
+  and props ABI preserved. See the UMAT catalog (`docs/simulation/umat_catalog.rst`).
 
 **Viscoelasticity:**
 - `Zener_fast` - Single Zener (Standard Linear Solid) element
diff --git a/docs/simulation/umat_catalog.rst b/docs/simulation/umat_catalog.rst
index c6b570af4..bb25313d3 100644
--- a/docs/simulation/umat_catalog.rst
+++ b/docs/simulation/umat_catalog.rst
@@ -10,9 +10,12 @@ names — **the calling convention is identical for all of them** (same
 - **modular (native)**: the composable ``MODUL`` engine, configured from a
   props stream (see :mod:`simcoon.modular` for the Python builder).
 - **modular (adapter)**: a legacy name whose dedicated kernel was removed
-  after its equivalence with a ``MODUL`` configuration was proven
-  (bit-identical results, dedicated tests in ``test_modular.py``); a
-  translator maps the legacy props to the modular configuration at each call.
+  after its equivalence with a ``MODUL`` configuration was proven against the
+  retained reference kernels — bit-identical for the elastic, power-law and
+  Hill families; within 2e-3 relative for the Voce/Chaboche families, whose
+  legacy incremental ``Hp`` update differs from the modular closed form
+  (tests in ``test_modular.py`` and ``Treference_umats``); a translator maps
+  the legacy props to the modular configuration at each call.
 - **legacy (kept)**: a dedicated, self-contained implementation kept either
   for pedagogy (readable single-file reference of the CCP return mapping) or
   because no modular equivalent exists.
diff --git a/examples/umat_prototype_tensor.cpp b/examples/umat_prototype_tensor.cpp
deleted file mode 100644
index e15912a77..000000000
--- a/examples/umat_prototype_tensor.cpp
+++ /dev/null
@@ -1,243 +0,0 @@
-/*
- * PROTOTYPE: Elastic-plastic UMAT using tensor2 / tensor4
- * ========================================================
- *
- * Side-by-side comparison with the current raw arma::vec / arma::mat style.
- * NOT meant to compile or replace anything — just to show how the code reads.
- */
-
-#include <armadillo>
-#include <simcoon/Continuum_mechanics/Functions/tensor.hpp>
-#include <simcoon/Continuum_mechanics/Functions/constitutive.hpp>
-#include <simcoon/Continuum_mechanics/Functions/contimech.hpp>
-#include <simcoon/Simulation/Maths/rotation.hpp>
-
-using namespace arma;
-using namespace simcoon;
-
-
-// ============================================================================
-//  Isotropic elasticity — 3 lines of mechanics
-// ============================================================================
-//
-//  BEFORE (raw arma):
-//      L = L_iso(E, nu, "Enu");                         // what type? who knows
-//      vec Eel = Etot + DEtot - alpha*(T+DT-T_init);    // stress voigt? strain voigt?
-//      sigma = el_pred(L, Eel, ndi);                     // hope you got the convention right
-//      Lt = L;
-//
-//  AFTER (typed):
-
-void umat_elasticity_iso_typed(
-    const vec &Etot, const vec &DEtot,
-    vec &sigma_out, mat &Lt_out, mat &L_out,
-    const vec &props, double T_init)
-{
-    // --- The actual mechanics: 3 lines ---
-    tensor4 L(L_iso(props(0), props(1), "Enu"), Tensor4Type::stiffness);
-    tensor2 sigma = L * strain(Etot + DEtot - props(2)*(T_init)));
-    Lt_out = L_out = L.mat();
-
-    sigma_out = sigma.voigt();
-}
-
-
-// ============================================================================
-//  NEW STYLE — elastic-plastic with isotropic hardening (CCP)
-// ============================================================================
-
-void umat_plasticity_iso_CCP_typed(
-    const vec &Etot, const vec &DEtot,
-    vec &sigma_out, mat &Lt_out, mat &L_out,
-    vec &sigma_in,
-    const mat &DR,
-    const vec &props, vec &statev,
-    const double &T, const double &DT,
-    const bool &start, const int &solver_type)
-{
-    double E        = props(0);
-    double nu       = props(1);
-    double alpha_v  = props(2);
-    double sigmaY   = props(3);
-    double k        = props(4);
-    double m        = props(5);
-
-    double T_init = statev(0);
-    double p      = statev(1);
-
-    // --- Typed stiffness ---
-    tensor4 L(L_iso(E, nu, "Enu"), Tensor4Type::stiffness);
-
-    // --- Plastic strain: tagged as VoigtType::strain ---
-    // Can't accidentally rotate it with rotate_stress anymore.
-    tensor2 EP = tensor2::from_voigt(
-        vec::fixed<6>(statev.memptr() + 2), VoigtType::strain
-    );
-
-    // --- Rotation of internal variables ---
-    // EP.rotate() dispatches to apply_strain automatically.
-    // With raw code you'd write:  EP = rotate_strain(EP, DR);
-    //   ... and if you wrote rotate_stress(EP, DR) by mistake: silent bug.
-    Rotation dR = Rotation::from_matrix(DR);
-    EP = EP.rotate(dR);
-
-    if (start) {
-        T_init = T;
-        EP = tensor2::zeros(VoigtType::strain);
-        p = 0.0;
-    }
-
-    // --- CTE as a strain tensor ---
-    tensor2 alpha_tensor = tensor2::identity(VoigtType::strain) * alpha_v;
-
-    // --- Elastic predictor ---
-    tensor2 Etot_t = strain(Etot + DEtot);
-
-    tensor2 thermal = alpha_tensor * (T + DT - T_init);
-    tensor2 Eel = Etot_t - thermal - EP;
-
-    // --- sigma = L : Eel ---
-    // Output is automatically VoigtType::stress. No ambiguity.
-    tensor2 sigma = L * Eel;
-
-    // --- Yield function ---
-    double Hp = (p > 1e-12) ? k * std::pow(p, m) : 0.0;
-    double dHpdp = (p > 1e-12) ? m * k * std::pow(p, m - 1) : 0.0;
-
-    double phi = Mises(sigma) - Hp - sigmaY;
-
-    // --- Return mapping (CCP loop) ---
-    tensor2 sigma_start = sigma;
-    tensor2 EP_start = EP;
-
-    if (phi > 0.0) {
-        // Flow direction: dev(sigma) / ||dev(sigma)||
-        // This is a stress-like quantity — dev() returns stress Voigt
-        vec::fixed<6> n_v = dev(sigma);
-        double norm_n = arma::norm(n_v);
-        if (norm_n > 1e-14) n_v /= norm_n;
-
-        // Flow direction as stress tensor
-        tensor2 N = tensor2::from_voigt(n_v, VoigtType::stress);
-
-        // kappa = L : N  → gives stress (stiffness contracts stress-normal → stress)
-        tensor2 kappa = L.contract(
-            tensor2::from_voigt(n_v, VoigtType::strain)  // n as strain for L:n
-        );
-
-        // Newton iterations on the plastic multiplier dp
-        double dp = 0.0;
-        for (int iter = 0; iter < 25; ++iter) {
-            Hp = (p + dp > 1e-12) ? k * std::pow(p + dp, m) : 0.0;
-            dHpdp = (p + dp > 1e-12) ? m * k * std::pow(p + dp, m - 1) : 0.0;
-
-            phi = Mises(sigma) - Hp - sigmaY;
-            if (std::abs(phi) < 1e-10) break;
-
-            double denom = arma::dot(vec(dev(sigma)) / Mises(sigma), vec(kappa.voigt())) + dHpdp;
-            double ddp = phi / denom;
-            dp += ddp;
-
-            // Update plastic strain and stress
-            // EP += ddp * N  — both are strain-tagged, arithmetic preserves tag
-            EP = EP + tensor2::from_voigt(vec::fixed<6>(n_v * ddp), VoigtType::strain);
-
-            Eel = Etot_t - thermal - EP;
-            sigma = L * Eel;
-        }
-        p += dp;
-    }
-
-    // --- Tangent modulus ---
-    if (solver_type == 0 || solver_type == 2) {
-        if (phi > -1e-10 && Mises(sigma) > 1e-10) {
-            // Elastoplastic tangent
-            vec::fixed<6> n_v = dev(sigma);
-            double norm_n = arma::norm(n_v);
-            if (norm_n > 1e-14) n_v /= norm_n;
-
-            // kappa = L : n (as strain)
-            vec::fixed<6> kappa_v = L.mat() * n_v;
-
-            // Denominator: n : L : n + dH/dp
-            double denom = arma::dot(n_v, kappa_v) + dHpdp;
-
-            // Lt = L - (kappa * kappa^T) / denom
-            // Type is preserved: stiffness - stiffness = stiffness
-            tensor4 correction(
-                mat::fixed<6,6>(kappa_v * kappa_v.t() / denom),
-                Tensor4Type::stiffness
-            );
-            tensor4 Lt = L - correction;
-            Lt_out = Lt.mat();
-        } else {
-            Lt_out = L.mat();
-        }
-    } else if (solver_type == 1) {
-        sigma_in = -L.mat() * EP.voigt();
-    }
-
-    // --- Export ---
-    sigma_out = sigma.voigt();
-    L_out = L.mat();
-
-    statev(0) = T_init;
-    statev(1) = p;
-    vec::fixed<6> ep_v = EP.voigt();
-    for (int i = 0; i < 6; ++i) statev(2 + i) = ep_v(i);
-}
-
-
-// ============================================================================
-//  Neo-Hookean (finite strain) — 5 lines of mechanics
-// ============================================================================
-//
-//  BEFORE: 15 lines of FTensor boilerplate + manual /J everywhere
-//  AFTER:  push_forward(F) includes the 1/J metric by default
-
-void umat_neohookean_typed(
-    const mat::fixed<3,3> &F,
-    vec &sigma_out, mat &c_out,
-    const vec &props)
-{
-    double mu  = props(0) / (2.0 * (1.0 + props(1)));
-    double lam = props(0) * props(1) / ((1.0 + props(1)) * (1.0 - 2.0 * props(1)));
-    double J   = arma::det(F);
-
-    // --- 5 lines: stress + tangent ---
-    tensor2 I = tensor2::identity(VoigtType::stress);
-    tensor2 b = stress(F * F.t());
-    tensor2 sigma = (mu * (b - I) + (lam * std::log(J)) * I) * (1.0 / J);
-
-    tensor4 C_ref = lam * tensor4::volumetric() + 2.0 * mu * tensor4::identity();
-    c_out = C_ref.push_forward(F).mat();       // metric=true: includes 1/J
-
-    sigma_out = sigma.voigt();
-}
-
-
-// ============================================================================
-//  Rotated anisotropic elasticity — rotation dispatch is automatic
-// ============================================================================
-//
-//  BEFORE: you must remember QS for stiffness, QE for compliance
-//  AFTER:  .rotate(R) picks the right rule from the type tag
-
-void umat_elastic_aniso_typed(
-    const vec &Etot, const vec &DEtot,
-    vec &sigma_out, mat &Lt_out, mat &L_out,
-    const vec &props, const double &theta_deg)
-{
-    tensor4 L_mat(L_isotrans(props(0), props(1), props(2), props(3), props(4), 1),
-                  Tensor4Type::stiffness);
-
-    Rotation R = Rotation::from_axis_angle(theta_deg * datum::pi / 180.0, 3);
-    tensor4 L_glob = L_mat.rotate(R);   // dispatches to QS * L * QS^T
-
-    // compliance rotation uses QE automatically — no manual dispatch
-    tensor4 M_glob = L_mat.inverse().rotate(R);  // QE * M * QE^T
-
-    tensor2 sigma = L_glob * strain(Etot + DEtot));
-    sigma_out = sigma.voigt();
-    Lt_out = L_out = L_glob.mat();
-}
diff --git a/examples/umat_tutorial/umat_tutorial_J2.cpp b/examples/umat_tutorial/umat_tutorial_J2.cpp
index d3a0fc18f..cb4c5182e 100644
--- a/examples/umat_tutorial/umat_tutorial_J2.cpp
+++ b/examples/umat_tutorial/umat_tutorial_J2.cpp
@@ -72,14 +72,14 @@ void umat_tutorial_J2(const std::string & /*umat_name*/, const vec &Etot,
     if (start) {
         T_init = T;
         p = 0.;
-        EP = tensor2::zeros(VoigtType::strain);
+        EP = tensor2::zeros(Tensor2Type::strain);
         sigma = zeros(6);
         Wm = 0.; Wm_r = 0.; Wm_ir = 0.; Wm_d = 0.;
     }
 
     // Objectivity: the stored plastic strain co-rotates with the material
     // frame. The typed rotate() picks the strain rotation kernel (factor-2
-    // shear convention) from the tensor's own VoigtType.
+    // shear convention) from the tensor's own Tensor2Type.
     EP = EP.rotate(Rotation::from_matrix(DR));
 
     const vec sigma_start = sigma;
@@ -100,7 +100,7 @@ void umat_tutorial_J2(const std::string & /*umat_name*/, const vec &Etot,
     const double f_trial = Mises(sig_trial) - (sigma_Y + H * p);
 
     double Dp = 0.;
-    tensor2 Lam = tensor2::zeros(VoigtType::strain);
+    tensor2 Lam = tensor2::zeros(Tensor2Type::strain);
 
     if (f_trial > 0.) {
         // -------------------------------------------------------------- 6.
diff --git a/examples/umats/external/umat_plugin_ext_EPICP.cpp b/examples/umats/external/umat_plugin_ext_EPICP.cpp
index cf7251379..b58c9a36b 100644
--- a/examples/umats/external/umat_plugin_ext_EPICP.cpp
+++ b/examples/umats/external/umat_plugin_ext_EPICP.cpp
@@ -19,7 +19,7 @@
  *   ``L - dyadic(kappa, kappa/Bhat)``      L - (kappa (x) kappa) / Bhat
  *
  * Overhead: the typed objects are fixed-size (stack, no heap). Construction uses
- * the ``VoigtType`` *enum* (no per-call string parsing), and the only Voigt->tensor
+ * the ``Tensor2Type`` *enum* (no per-call string parsing), and the only Voigt->tensor
  * conversion sits *before* the Newton loop — the hot loop never touches a string
  * tag or rebuilds a tensor from a raw vector.
  *
@@ -101,15 +101,15 @@ class LIB_EXPORT umat_epicp_tensor : public umat_plugin_ext_api {
         L_out = L.mat();
 
         // ---- Wrap incoming state as tensors (the only Voigt->tensor conversions) ----
-        tensor2 sig       = tensor2::from_voigt(sigma, VoigtType::stress);
-        tensor2 EP        = tensor2::from_voigt(vec(statev.subvec(2, 7)), VoigtType::strain);
+        tensor2 sig       = tensor2::from_voigt(sigma, Tensor2Type::stress);
+        tensor2 EP        = tensor2::from_voigt(vec(statev.subvec(2, 7)), Tensor2Type::strain);
         tensor2 sig_start = sig;
         tensor2 EP_start  = EP;
 
         // ---- Trial mechanical strain (built once, reused every iteration) ----
         vec alpha_v = {alpha, alpha, alpha, 0., 0., 0.};
         tensor2 eps_tot = tensor2::from_voigt(
-            vec(Etot + DEtot - alpha_v * (T + DT - T_init)), VoigtType::strain);
+            vec(Etot + DEtot - alpha_v * (T + DT - T_init)), Tensor2Type::strain);
 
         // ---- Hardening at the start of the step ----
         double Hp = 0., dHpdp = 0.;
@@ -123,7 +123,7 @@ class LIB_EXPORT umat_epicp_tensor : public umat_plugin_ext_api {
         sig = L.contract(eps_tot - EP);
 
         // ---- Return mapping (CCP) — scalars stay primitive, tensors stay typed ----
-        tensor2 N(VoigtType::strain), kappa(VoigtType::stress);
+        tensor2 N(Tensor2Type::strain), kappa(Tensor2Type::stress);
         double error = 1.;
         for (int iter = 0;
              iter < simcoon::maxiter_umat && error > simcoon::precision_umat;
@@ -165,7 +165,7 @@ class LIB_EXPORT umat_epicp_tensor : public umat_plugin_ext_api {
         // ---- Work (typed double contractions  sigma : d_eps) ----
         tensor2 sig_sum = sig_start + sig;
         tensor2 DEP     = EP - EP_start;
-        tensor2 DEtot_t = tensor2::from_voigt(DEtot, VoigtType::strain);
+        tensor2 DEtot_t = tensor2::from_voigt(DEtot, Tensor2Type::strain);
         double  Dp      = p - p_old;
         double  A_p     = -Hp;
 
diff --git a/examples/umats/external/umat_plugin_ext_NEOHC.cpp b/examples/umats/external/umat_plugin_ext_NEOHC.cpp
index d9ae53b64..572f17628 100644
--- a/examples/umats/external/umat_plugin_ext_NEOHC.cpp
+++ b/examples/umats/external/umat_plugin_ext_NEOHC.cpp
@@ -114,8 +114,8 @@ class LIB_EXPORT umat_neohookean_tensor : public umat_plugin_ext_api {
         mat invC_mat = inv(C_mat);
 
         // ---- PKII stress using Tensor2 (enum tags: no string parsing in the UMAT hot path) ----
-        tensor2 invC_t = tensor2::from_voigt(t2v_stress(invC_mat), VoigtType::stress);
-        tensor2 I_t    = tensor2::identity(VoigtType::stress);
+        tensor2 invC_t = tensor2::from_voigt(t2v_stress(invC_mat), Tensor2Type::stress);
+        tensor2 I_t    = tensor2::identity(Tensor2Type::stress);
         tensor2 S = mu * (I_t - invC_t) + lambda * lnJ * invC_t;
 
         // ---- Cauchy stress via push-forward ----
diff --git a/include/simcoon/Continuum_mechanics/Functions/tensor.hpp b/include/simcoon/Continuum_mechanics/Functions/tensor.hpp
index fd4468b2c..69e151f8e 100644
--- a/include/simcoon/Continuum_mechanics/Functions/tensor.hpp
+++ b/include/simcoon/Continuum_mechanics/Functions/tensor.hpp
@@ -22,6 +22,12 @@
 #include <optional>
 #include <Fastor/Fastor.h>
 
+/**
+ * @file tensor.hpp
+ * @brief Typed 2nd- and 4th-order tensor API: type-tagged Voigt/Mandel conversions,
+ * convention-aware contractions and rotations (Voigt order [11,22,33,12,13,23], units MPa).
+ */
+
 namespace simcoon {
 
 // Forward declarations
diff --git a/include/simcoon/Simulation/Solver/step_meca.hpp b/include/simcoon/Simulation/Solver/step_meca.hpp
index 83142e7e0..d156dc653 100755
--- a/include/simcoon/Simulation/Solver/step_meca.hpp
+++ b/include/simcoon/Simulation/Solver/step_meca.hpp
@@ -68,6 +68,14 @@ class step_meca : public step
     
     using step::generate;
     virtual void generate(const double&, const arma::vec &, const arma::vec&, const double&);
+    /**
+     * @brief Generate the increment tables for a fully kinematic step (control types 5/6).
+     * @param[in] mTime current (absolute) time at the start of the step
+     * @param[in] mF current 3x3 deformation gradient \f$ \mathbf{F} \f$ (the step anchor)
+     * @param[in] mT current temperature
+     * Modes 1/2 interpolate geodesically between mF and BC_meca via the real matrix
+     * logarithm/exponential; mode 3 reads the 9-component F table (absolute time column).
+     */
     virtual void generate_kin(const double&, const arma::mat &, const double &);
     virtual void assess_inc(const double &, double &, const double &, phase_characteristics &, double &, const double &, const arma::mat &, const int &);
     
diff --git a/include/simcoon/Simulation/Solver/step_thermomeca.hpp b/include/simcoon/Simulation/Solver/step_thermomeca.hpp
index 5c466977e..0810b0b04 100755
--- a/include/simcoon/Simulation/Solver/step_thermomeca.hpp
+++ b/include/simcoon/Simulation/Solver/step_thermomeca.hpp
@@ -67,6 +67,13 @@ class step_thermomeca : public step
     
     using step::generate;
     virtual void generate(const double&, const arma::vec&, const arma::vec&, const double&);
+    /**
+     * @brief Generate the increment tables for a fully kinematic thermomechanical step
+     * (control types 5/6); same geodesic-F conventions as step_meca::generate_kin.
+     * @param[in] mTime current (absolute) time at the start of the step
+     * @param[in] mF current 3x3 deformation gradient \f$ \mathbf{F} \f$ (the step anchor)
+     * @param[in] mT current temperature
+     */
     virtual void generate_kin(const double&, const arma::mat&m, const double &);    
     virtual void assess_inc(const double &, double &, const double &, phase_characteristics &, double &, const double &, const arma::mat &, const int &);
     
diff --git a/pyproject.toml b/pyproject.toml
index b3833db13..9a85a774b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -8,7 +8,7 @@ build-backend = "scikit_build_core.build"
 
 [project]
 name = "simcoon"
-version = "1.14.0"
+version = "2.0.0b1"
 description = "Simulation in Mechanics and Materials: Interactive Tools - A library for the simulation of multiphysics systems and heterogeneous materials"
 readme = "README.md"
 license = {text = "GPL-3.0-or-later"}
diff --git a/python-setup/simcoon/solver/files.py b/python-setup/simcoon/solver/files.py
index 36862fd57..e2f7f3f87 100644
--- a/python-setup/simcoon/solver/files.py
+++ b/python-setup/simcoon/solver/files.py
@@ -199,9 +199,37 @@ def from_file(path_data: str = "data", pathfile: str = "path.txt") -> Tuple[List
         parse_step = _parse_meca_step if btype == 1 else _parse_thermomeca_step
         steps = [parse_step(tk, control_type, path_data) for _ in range(nstep)]
         blocks.append(Block(steps=steps, control_type=control_type, ncycle=ncycle))
+    _reanchor_tabular_times(blocks)
     return blocks, T_init
 
 
+def _reanchor_tabular_times(blocks) -> None:
+    """Shift restarting mode-3 time axes so from_file matches the C++ file path.
+
+    The legacy FILE convention lets each mode-3 table restart its time axis at
+    0 (first row = anchor at the current state); the C++ file reader shifts
+    such an axis to continue from the running simulation time. In-memory
+    ``tabular`` arrays are bound to the strict absolute-time contract instead,
+    so a parsed legacy table must be re-anchored here — otherwise a file set
+    that runs fine through the file-driven solver is rejected by solve().
+    Tables whose time column already continues absolutely are left untouched
+    (shift only when the first row is EARLIER than the running time, the same
+    guard as the C++ reader).
+    """
+    t_run = 0.0
+    for b in blocks:
+        for _ in range(b.ncycle):
+            for s in b.steps:
+                if s.mode in ("tabular", 3) and s.tabular is not None:
+                    t0 = float(s.tabular[0, 0])
+                    if t0 < t_run - 1e-12:
+                        s.tabular = s.tabular.copy()
+                        s.tabular[:, 0] += t_run - t0
+                    t_run = float(s.tabular[-1, 0])
+                else:
+                    t_run += float(s.time)
+
+
 def material_from_file(path_data: str = "data", materialfile: str = "material.dat") -> dict:
     """Parse a legacy material definition file into solve() keyword arguments.
 
diff --git a/simcoon-python-builder/include/simcoon/python_wrappers/Libraries/Solver/solver_run.hpp b/simcoon-python-builder/include/simcoon/python_wrappers/Libraries/Solver/solver_run.hpp
index fd5b7156e..d340740fa 100644
--- a/simcoon-python-builder/include/simcoon/python_wrappers/Libraries/Solver/solver_run.hpp
+++ b/simcoon-python-builder/include/simcoon/python_wrappers/Libraries/Solver/solver_run.hpp
@@ -8,7 +8,27 @@ namespace py = pybind11;
 
 namespace simpy {
 
-//In-memory solver: takes loading blocks as a list of dicts, returns results as a dict of numpy arrays
+/**
+ * @brief In-memory solver entry point for the Python bindings (no file I/O).
+ *
+ * Runs the same incremental solver loop as the file-based binding, but takes the
+ * loading programme as Python data and returns the full converged history in memory.
+ *
+ * @param[in] blocks_py list of block dicts (control_type, ncycle, steps[] with
+ *            cBC_meca/BC_meca/mode/times/tabular... — see simcoon.solver.Block.to_dict)
+ * @param[in] T_init initial temperature
+ * @param[in] umat_name 5-char UMAT name (dispatch key)
+ * @param[in] props_py material properties (units MPa)
+ * @param[in] nstatev number of internal state variables
+ * @param[in] psi_rve,theta_rve,phi_rve material orientation Euler angles
+ * @param[in] solver_type solver scheme (0 = Newton)
+ * @param[in] corate_type objective rate (0 Jaumann, 1 Green-Naghdi, 2 log, 3 log_R, 4 Truesdell, 5 log_F)
+ * @param[in] params_py numeric controls (div/mul_tnew_dt, miniter, maxiter, inforce,
+ *            precision, lambda_solver, tangent_mode)
+ * @param[in] record_tangent whether to record the per-increment tangent history
+ * @return dict of numpy arrays keyed "Strain", "Stress", "F", "R", "Statev", "Wm", "Temp",
+ *         "Time", ... plus "status" (0 = completed) — consumed by simcoon.solver.SolverResults
+ */
 py::dict solver_run(const py::list &blocks_py, const double &T_init,
                     const std::string &umat_name, const py::array_t<double> &props_py,
                     const int &nstatev, const double &psi_rve, const double &theta_rve, const double &phi_rve,
diff --git a/simcoon-python-builder/src/python_wrappers/python_module.cpp b/simcoon-python-builder/src/python_wrappers/python_module.cpp
index a940b49b4..5857394eb 100755
--- a/simcoon-python-builder/src/python_wrappers/python_module.cpp
+++ b/simcoon-python-builder/src/python_wrappers/python_module.cpp
@@ -230,6 +230,12 @@ PYBIND11_MODULE(_core, m)
     m.def("stress_convert", &stress_convert, "sigma"_a, "F"_a, "converter_key"_a, "J"_a = 0., "copy"_a = true, simcoon_docs::stress_convert);
 
     // umat
+    // NOTE (single/direct callers): kernels may request a step cut instead of
+    // integrating a too-large increment (e.g. the modular engine rejects a
+    // non-finite or runaway return-mapping result and leaves statev untouched
+    // with an elastic Lt). The simcoon solver honors this via tnew_dt < 1 and
+    // retries with a smaller increment; couplers driving umat directly must
+    // subdivide the increment themselves in that situation.
     m.def("umat", &launch_umat, "umat_name"_a, "etot"_a, "Detot"_a, "F0"_a, "F1"_a, "sigma"_a, "DR"_a, "props"_a, "statev"_a, "time"_a, "dtime"_a, "Wm"_a, "temp"_a = pybind11::none(), "ndi"_a = 3, "n_threads"_a = 4, "tangent_mode"_a = simcoon::tangent_default);
     m.def("umat_T", &launch_umat_T, "umat_name"_a, "etot"_a, "Detot"_a, "sigma"_a, "DR"_a, "props"_a, "statev"_a, "time"_a, "dtime"_a, "Wm"_a, "Wt"_a, "T"_a, "DT"_a, "ndi"_a = 3, "n_threads"_a = 4, "tangent_mode"_a = simcoon::tangent_default);
 
diff --git a/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.cpp b/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.cpp
index 98aea36b4..8aafd436d 100755
--- a/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.cpp
+++ b/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.cpp
@@ -634,8 +634,16 @@ void umat_sma_unified_T(const string &umat_name, const vec &Etot, const vec &DEt
     Bhat(1,1) = sum(dPhiRdsigma%kappa_j[1]) - K(1,1);
 
     const std::vector<vec> dPhidsigma_l = { dPhiFdsigma, dPhiRdsigma };
+    // tangent_algorithmic is clamped to continuum for the SMA kernels: the
+    // finite-difference transformation-flow Hessian makes the local system
+    // near-singular on the transformation plateau (rcond ~ 1e-17; whether the
+    // backend flags it is LAPACK-dependent — see the thermomechanical twin,
+    // which crashed CI through this path). Re-enable together with the exact
+    // CPP Hessian rework.
+    const int tangent_mode_eff = (tangent_mode == tangent_algorithmic)
+        ? tangent_continuum : tangent_mode;
     const ContinuumTangent ct = compute_tangent_operator(
-        tangent_mode, Bhat, kappa_j, dPhidsigma_l, Ds_j, L,
+        tangent_mode_eff, Bhat, kappa_j, dPhidsigma_l, Ds_j, L,
         [&]() -> std::vector<mat> {  // lazy: evaluated only in algorithmic mode
             // Simo-Hughes algorithmic tangent (closest-point). The forward transformation strain-flow
             // Lambda_eps^F = Hcur(sigma)*dDrucker(sigma) + DM*sigma + Dalpha_T couples to stress through
diff --git a/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_TR.cpp b/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_TR.cpp
index 479b09b44..bab164dc2 100644
--- a/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_TR.cpp
+++ b/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_TR.cpp
@@ -657,8 +657,16 @@ void umat_sma_unified_TR(const string &umat_name, const vec &Etot, const vec &DE
     }
 
     const std::vector<vec> dPhidsigma_l = { dPhiFdsigma, dPhiRdsigma, dPhiReodsigma };
+    // tangent_algorithmic is clamped to continuum for the SMA kernels: the
+    // finite-difference transformation-flow Hessian makes the local system
+    // near-singular on the transformation plateau (rcond ~ 1e-17; whether the
+    // backend flags it is LAPACK-dependent — see the thermomechanical twin,
+    // which crashed CI through this path). Re-enable together with the exact
+    // CPP Hessian rework.
+    const int tangent_mode_eff = (tangent_mode == tangent_algorithmic)
+        ? tangent_continuum : tangent_mode;
     const ContinuumTangent ct = compute_tangent_operator(
-        tangent_mode, Bhat, kappa_j, dPhidsigma_l, Ds_j, L,
+        tangent_mode_eff, Bhat, kappa_j, dPhidsigma_l, Ds_j, L,
         [&]() -> std::vector<mat> {  // lazy: evaluated only in algorithmic mode
             // Simo-Hughes algorithmic tangent (closest-point), 3-mechanism. dLambda/dsigma by central
             // FD of the transformation flows: forward Hcur(sigma)*dDrucker(sigma); reorientation
diff --git a/src/Simulation/Solver/step.cpp b/src/Simulation/Solver/step.cpp
index 593412dc2..09df10736 100755
--- a/src/Simulation/Solver/step.cpp
+++ b/src/Simulation/Solver/step.cpp
@@ -151,7 +151,10 @@ int step::mode3_ninc() const
     ifstream pathinc(file, ios::in);
     if(!pathinc)
     {
-        cout << "Error: cannot open the file " << file << "\n Please check if the file is correct and if you have added the extension\n";
+        // Throw instead of falling through: getline on a never-opened stream
+        // fails without ever setting eofbit, so the counting loop below would
+        // spin forever.
+        throw simcoon::exception_solver("step " + std::to_string(number) + " (mode 3): cannot open the increment file '" + file + "' — check the file name (including extension) and the data path");
     }
     while (!pathinc.eof())
     {
diff --git a/src/Simulation/Solver/step_thermomeca.cpp b/src/Simulation/Solver/step_thermomeca.cpp
index f54256324..8f793e1ec 100755
--- a/src/Simulation/Solver/step_thermomeca.cpp
+++ b/src/Simulation/Solver/step_thermomeca.cpp
@@ -15,7 +15,7 @@
  
  */
 
-///@file step_thermomeca.hpp
+///@file step_thermomeca.cpp
 ///@brief object that defines a thermomechanical step
 ///@version 1.0
 
diff --git a/test/Libraries/Umat/Treference_umats.cpp b/test/Libraries/Umat/Treference_umats.cpp
index 56714d3cd..49e7c28d0 100644
--- a/test/Libraries/Umat/Treference_umats.cpp
+++ b/test/Libraries/Umat/Treference_umats.cpp
@@ -59,9 +59,13 @@ using UmatFn = std::function<void(
 
 // Drive one kernel through a cyclic confined-uniaxial strain path
 // (E11: 0 -> +amp -> -amp -> 0, other components held at zero) and return the
-// sigma11 history. tangent_continuum on both sides for determinism.
+// sigma11 history. Continuum tangent on both sides — NOTE the reference
+// kernels are verbatim PRE-2.0 copies, so their continuum mode is the OLD
+// literal 0, not tangent_continuum (=1, which they would read as algorithmic).
+// The mode only selects the Lt output in these kernels (stress is unaffected),
+// but the intent is a continuum-vs-continuum comparison.
 vec drive(const UmatFn& umat, const std::string& name, const vec& props,
-          int nstatev, double amp, int n_quarter) {
+          int nstatev, double amp, int n_quarter, int tangent_mode) {
     vec Etot = zeros(6);
     vec sigma = zeros(6);
     vec statev = zeros(nstatev);
@@ -85,7 +89,7 @@ vec drive(const UmatFn& umat, const std::string& name, const vec& props,
         umat(name, Etot, DEtot, sigma, Lt, L, DR,
              static_cast<int>(props.n_elem), props, nstatev, statev,
              T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d,
-             3, 3, start, tnew_dt, tangent_continuum);
+             3, 3, start, tnew_dt, tangent_mode);
         Etot += DEtot;
         Time += DTime;
         start = false;
@@ -97,8 +101,10 @@ vec drive(const UmatFn& umat, const std::string& name, const vec& props,
 void expect_equiv(const UmatFn& reference, const std::string& name,
                   const vec& props, int nstatev, double rel_tol,
                   double amp = 0.02, int n_quarter = 25) {
-    const vec ref = drive(reference, name, props, nstatev, amp, n_quarter);
-    const vec mod = drive(umat_legacy_modular, name, props, nstatev, amp, n_quarter);
+    const vec ref = drive(reference, name, props, nstatev, amp, n_quarter,
+                          /*pre-2.0 continuum literal*/ 0);
+    const vec mod = drive(umat_legacy_modular, name, props, nstatev, amp,
+                          n_quarter, tangent_continuum);
     const double peak = std::max(std::abs(ref.max()), std::abs(ref.min()));
     ASSERT_GT(peak, 0.);
     const double dev = abs(ref - mod).max() / peak;
diff --git a/test/Libraries/Umat/reference_kernels/Ani_chaboche_ccp.cpp b/test/Libraries/Umat/reference_kernels/Ani_chaboche_ccp.cpp
index 0642d9afa..a4d65071f 100644
--- a/test/Libraries/Umat/reference_kernels/Ani_chaboche_ccp.cpp
+++ b/test/Libraries/Umat/reference_kernels/Ani_chaboche_ccp.cpp
@@ -19,7 +19,7 @@
  
  */
 
-///@file plastic_kin_iso_ccp.cpp
+///@file Ani_chaboche_ccp.cpp
 ///@brief User subroutine for elastic-plastic materials in 1D-2D-3D case
 ///@brief This subroutines uses a convex cutting plane algorithm
 ///@brief Linear Kinematical hardening coupled with a power-law hardenig is considered
diff --git a/test/Libraries/Umat/reference_kernels/DFA_chaboche_ccp.cpp b/test/Libraries/Umat/reference_kernels/DFA_chaboche_ccp.cpp
index 638ab7e8e..1ad8c52b2 100644
--- a/test/Libraries/Umat/reference_kernels/DFA_chaboche_ccp.cpp
+++ b/test/Libraries/Umat/reference_kernels/DFA_chaboche_ccp.cpp
@@ -19,7 +19,7 @@
  
  */
 
-///@file plastic_kin_iso_ccp.cpp
+///@file DFA_chaboche_ccp.cpp
 ///@brief User subroutine for elastic-plastic materials in 1D-2D-3D case
 ///@brief This subroutines uses a convex cutting plane algorithm
 ///@brief Linear Kinematical hardening coupled with a power-law hardenig is considered
diff --git a/test/Libraries/Umat/reference_kernels/Generic_chaboche_ccp.cpp b/test/Libraries/Umat/reference_kernels/Generic_chaboche_ccp.cpp
index fff9c486e..6590be2d8 100644
--- a/test/Libraries/Umat/reference_kernels/Generic_chaboche_ccp.cpp
+++ b/test/Libraries/Umat/reference_kernels/Generic_chaboche_ccp.cpp
@@ -19,7 +19,7 @@
  
  */
 
-///@file plastic_kin_iso_ccp.cpp
+///@file Generic_chaboche_ccp.cpp
 ///@brief User subroutine for elastic-plastic materials in 1D-2D-3D case
 ///@brief This subroutines uses a convex cutting plane algorithm
 ///@brief Linear Kinematical hardening coupled with a power-law hardenig is considered
diff --git a/test/Libraries/Umat/reference_kernels/Hill_chaboche_ccp.cpp b/test/Libraries/Umat/reference_kernels/Hill_chaboche_ccp.cpp
index c7775b23f..a1f888948 100644
--- a/test/Libraries/Umat/reference_kernels/Hill_chaboche_ccp.cpp
+++ b/test/Libraries/Umat/reference_kernels/Hill_chaboche_ccp.cpp
@@ -19,7 +19,7 @@
  
  */
 
-///@file plastic_kin_iso_ccp.cpp
+///@file Hill_chaboche_ccp.cpp
 ///@brief User subroutine for elastic-plastic materials in 1D-2D-3D case
 ///@brief This subroutines uses a convex cutting plane algorithm
 ///@brief Linear Kinematical hardening coupled with a power-law hardenig is considered
diff --git a/test/Libraries/Umat/reference_kernels/Hill_isoh.cpp b/test/Libraries/Umat/reference_kernels/Hill_isoh.cpp
index f1dc92e7d..667d377ef 100644
--- a/test/Libraries/Umat/reference_kernels/Hill_isoh.cpp
+++ b/test/Libraries/Umat/reference_kernels/Hill_isoh.cpp
@@ -19,7 +19,7 @@
  
  */
 
-///@file Iso_Nfast.cpp
+///@file Hill_isoh.cpp
 ///@brief User subroutine for plasticity model in 3D case
 ///@author Chemisky, Chatzigeorgiou
 ///@version 1.0
diff --git a/test/Libraries/Umat/reference_kernels/Hill_isoh_Nfast.cpp b/test/Libraries/Umat/reference_kernels/Hill_isoh_Nfast.cpp
index 359ffda17..f6f08c231 100644
--- a/test/Libraries/Umat/reference_kernels/Hill_isoh_Nfast.cpp
+++ b/test/Libraries/Umat/reference_kernels/Hill_isoh_Nfast.cpp
@@ -19,7 +19,7 @@
  
  */
 
-///@file Iso_Nfast.cpp
+///@file Hill_isoh_Nfast.cpp
 ///@brief User subroutine for Multiple plasticity model in 3D case
 ///@author Chemisky, Chatzigeorgiou
 ///@version 1.0
diff --git a/test/Libraries/Umat/reference_kernels/elastic_orthotropic.cpp b/test/Libraries/Umat/reference_kernels/elastic_orthotropic.cpp
index 8e83347b9..029d4ced7 100644
--- a/test/Libraries/Umat/reference_kernels/elastic_orthotropic.cpp
+++ b/test/Libraries/Umat/reference_kernels/elastic_orthotropic.cpp
@@ -19,7 +19,7 @@
  
  */
 
-///@file Elastic_orthotropic.cpp
+///@file elastic_orthotropic.cpp
 ///@brief User subroutine for ortothropic elastic materials in 3D case
 
 #include <iostream>
diff --git a/test/Libraries/Umat/reference_kernels/elastic_transverse_isotropic.cpp b/test/Libraries/Umat/reference_kernels/elastic_transverse_isotropic.cpp
index 5ba719fd1..b3e306ec3 100644
--- a/test/Libraries/Umat/reference_kernels/elastic_transverse_isotropic.cpp
+++ b/test/Libraries/Umat/reference_kernels/elastic_transverse_isotropic.cpp
@@ -19,7 +19,7 @@
  
  */
 
-///@file elastic_transversely_isotropic.cpp
+///@file elastic_transverse_isotropic.cpp
 ///@brief User subroutine for transversely isotropic elastic materials in 3D case
 ///@version 1.0