Add J2 plasticity analytical tangent

Also remove FPE trapping as it causes more
problems due to how FPEs are used in some
Julia packages.

Author: lxmota

src/Norma.jl

@@ -21,9 +21,6 @@ end
 
 function run(sim::Simulation)
     start_time = time()
-    if get(sim.params, "enable FPE", false) == true
-        enable_fpe_traps()
-    end
     evolve(sim)
     elapsed_time = time() - start_time
     norma_log(0, :done, "Simulation Complete")

src/constitutive.jl

@@ -567,7 +567,8 @@ function _j2_stress(
 end
 
 # Consistent algorithmic tangent AA_{iJkL} = ∂P_{iJ}/∂F_{kL} via forward FD.
-function _j2_tangent(
+# 9 extra stress evaluations; used for benchmarking against the analytical tangent.
+function _j2_tangent_fd(
     material::J2Plasticity, F::SMatrix{3,3,Float64,9}, state_old::Vector{Float64},
     P0::SMatrix{3,3,Float64,9}
 )
@@ -589,9 +590,151 @@ function _j2_tangent(
     return SArray{Tuple{3,3,3,3},Float64,4,81}(AA)
 end
 
+# ---------------------------------------------------------------------------
+# Analytical tangent AA_{iJkL} = ∂P_{iJ}/∂F_{kL}.
+#
+# P = Fᵉ_new⁻ᵀ M_new Fᵖ_new⁻ᵀ,  Fᵉ_new = F Fᵖ_new⁻¹.
+#
+# Approximation: Fᵖ_new is treated as FROZEN when differentiating P with
+# respect to F.  The consequence:
+#   • Elastic step  (Fᵖ_new = Fᵖ_old, truly constant): result is EXACT.
+#   • Plastic step  (Fᵖ_new depends on F): the term ∂Fᵖ_new/∂F is omitted,
+#     introducing an error O(Δεᵖ) relative to the full consistent tangent.
+#     In practice this error is small (~0.1–1 % of the tangent norm) and
+#     Newton convergence remains rapid.
+#
+# ∂P/∂F = geometric + material contributions:
+#
+#   Geometric:  –(F⁻ᵀ)_{iL} P_{kJ}
+#               from ∂Fᵉ_new⁻ᵀ/∂F at fixed Fᵖ_new and M_new.
+#
+#   Material:   2 Σ_{ABCD} (Fᵉ_new⁻ᵀ)_{iA} ℂ_{ABCD} (Fᵖ_new⁻ᵀ)_{BJ}
+#                          (Fᵉ_tr)_{kC} (Fᵖ_old⁻¹)_{DL}
+#               from ∂M_new/∂Cᵉ_tr · ∂Cᵉ_tr/∂F; the factor of 2 comes
+#               from symmetry of ℂ in (C,D) and ∂Cᵉ_tr/∂F.
+#               Note: Cᵉ_tr = Fᵉ_trᵀ Fᵉ_tr uses Fᵖ_OLD (fixed), so this
+#               chain is exact.
+#
+# ℂ_{ABCD} = ∂M_new/∂Cᵉ_tr assembled in the spectral basis {nᵢ} of Cᵉ_tr:
+#
+#   Material part:  Σᵢⱼ [ĉᵢⱼ/(2cⱼ)] (nᵢ⊗nᵢ)_{AB} (nⱼ⊗nⱼ)_{CD}
+#   Geometric part: Σᵢ≠ⱼ θᵢⱼ (nᵢ⊗nⱼ)_sym_{AB} (nᵢ⊗nⱼ)_sym_{CD}
+#
+#   cᵢ = eigenvalues of Cᵉ_tr,  nᵢ = eigenvectors,
+#   θᵢⱼ = (mᵢ_new − mⱼ_new)/(cᵢ − cⱼ)  [L'Hôpital limit when cᵢ ≈ cⱼ].
+#
+# Principal-space algorithmic moduli ĉᵢⱼ = ∂mᵢ_new/∂εⱼ_tr:
+#   Elastic:  ĉᵢⱼ = λ + 2μ δᵢⱼ
+#   Plastic:  ĉᵢⱼ = (λ + μβ) + 2μ(1 – 3β/2) δᵢⱼ + (2μβ – 4μ²/(3μ+H)) Nᵢ Nⱼ
+#             β = 2μ Δεᵖ / σvm_tr,  Nᵢ = 1.5 mdev_i_tr / σvm_tr.
+# ---------------------------------------------------------------------------
+function _j2_tangent_analytical(
+    material::J2Plasticity,
+    F::SMatrix{3,3,Float64,9},
+    state_old::Vector{Float64},
+    P::SMatrix{3,3,Float64,9},
+    state_new::Vector{Float64},
+)
+    λ = material.λ
+    μ = material.μ
+    H = material.H
+
+    # Recover kinematics
+    Fp_old = reshape(state_old[1:9], 3, 3)
+    eqps_old = state_old[10]
+    Fp_new = reshape(state_new[1:9], 3, 3)
+    Δεᵖ = state_new[10] - eqps_old
+
+    Fm = Matrix{Float64}(F)
+    Fp_old_inv = inv(Fp_old)
+    Fp_new_inv = inv(Fp_new)
+    Fe_tr = Fm * Fp_old_inv
+    Fe_new = Fm * Fp_new_inv
+    Fe_new_inv_T = inv(Fe_new)'
+    F_inv_T = inv(Fm)'
+
+    # Spectral decomposition of trial Ce
+    Ce_tr = Symmetric(Fe_tr' * Fe_tr)
+    eig = eigen(Ce_tr)
+    c = eig.values      # principal squared stretches
+    Q = eig.vectors     # columns are eigenvectors
+
+    # Trial principal Hencky strains and Mandel stresses
+    ε_tr = 0.5 .* log.(c)
+    tr_ε = sum(ε_tr)
+    m_tr = [λ * tr_ε + 2μ * ε_tr[i] for i in 1:3]
+    m_dev_tr = m_tr .- sum(m_tr) / 3
+    σvm_tr = sqrt(1.5) * norm(m_dev_tr)
+
+    # Principal-space algorithmic moduli ĉᵢⱼ and updated principal stresses
+    f_tr = σvm_tr - (material.σy + H * eqps_old)
+    if f_tr ≤ 0.0
+        # Elastic step
+        m_new = m_tr
+        ĉ = [λ + 2μ * Float64(i == j) for i in 1:3, j in 1:3]
+    else
+        # Plastic step
+        N_pr = 1.5 .* m_dev_tr ./ σvm_tr   # principal values of flow direction
+        m_new = m_tr .- 2μ * Δεᵖ .* N_pr
+        β = 2μ * Δεᵖ / σvm_tr
+        A = λ + μ * β
+        B = 2μ * (1.0 - 1.5β)
+        Ccoef = 2μ * β - 4μ^2 / (3μ + H)
+        ĉ = [A + B * Float64(i == j) + Ccoef * N_pr[i] * N_pr[j] for i in 1:3, j in 1:3]
+    end
+
+    # Assemble ℂ_{ABCD} = ∂M_new/∂Ce_tr in the principal basis
+    CC = zeros(3, 3, 3, 3)
+
+    # Material part: Σᵢⱼ [ĉᵢⱼ/(2cⱼ)] (nᵢ⊗nᵢ)_{AB} (nⱼ⊗nⱼ)_{CD}
+    for α in 1:3, β_idx in 1:3
+        coeff = ĉ[α, β_idx] / (2c[β_idx])
+        nα = Q[:, α]
+        nβ = Q[:, β_idx]
+        for A in 1:3, B in 1:3, C in 1:3, D in 1:3
+            CC[A, B, C, D] += coeff * nα[A] * nα[B] * nβ[C] * nβ[D]
+        end
+    end
+
+    # Geometric part: Σᵢ≠ⱼ θᵢⱼ (nᵢ⊗nⱼ)_sym_{AB} (nᵢ⊗nⱼ)_sym_{CD}
+    tol = 1e-10
+    for α in 1:3, β_idx in 1:3
+        α == β_idx && continue
+        nα = Q[:, α]
+        nβ = Q[:, β_idx]
+        Δc = c[α] - c[β_idx]
+        if abs(Δc) > tol * (c[α] + c[β_idx])
+            θ = (m_new[α] - m_new[β_idx]) / Δc
+        else
+            # L'Hôpital: lim_{cβ→cα} (mα−mβ)/(cα−cβ) = (ĉ_{αα}−ĉ_{βα})/(2cα)
+            θ = (ĉ[α, α] - ĉ[β_idx, α]) / (2c[α])
+        end
+        for A in 1:3, B in 1:3, C in 1:3, D in 1:3
+            sym_AB = 0.5 * (nα[A] * nβ[B] + nα[B] * nβ[A])
+            sym_CD = 0.5 * (nα[C] * nβ[D] + nα[D] * nβ[C])
+            CC[A, B, C, D] += θ * sym_AB * sym_CD
+        end
+    end
+
+    # Assemble AA_{iJkL} = –(F⁻ᵀ)_{iL} P_{kJ}
+    #   + 2 Σ_{ABCD} (Fe_new⁻ᵀ)_{iA} CC_{ABCD} (Fp_new⁻ᵀ)_{BJ} (Fe_tr)_{kC} (Fp_old⁻¹)_{DL}
+    # Factor of 2 arises from symmetry of CC in (C,D) and ∂Ce_tr/∂F.
+    Fp_new_inv_T = Fp_new_inv'
+    AA = MArray{Tuple{3,3,3,3},Float64}(undef)
+    for i in 1:3, J in 1:3, k in 1:3, L in 1:3
+        val = -F_inv_T[i, L] * P[k, J]
+        for A in 1:3, B in 1:3, C in 1:3, D in 1:3
+            val += 2 * Fe_new_inv_T[i, A] * CC[A, B, C, D] * Fp_new_inv_T[B, J] *
+                   Fe_tr[k, C] * Fp_old_inv[D, L]
+        end
+        AA[i, J, k, L] = val
+    end
+    return SArray{Tuple{3,3,3,3},Float64,4,81}(AA)
+end
+
 function constitutive(material::J2Plasticity, F::SMatrix{3,3,Float64,9}, state_old::Vector{Float64})
     W, P, state_new = _j2_stress(material, F, state_old)
-    AA = _j2_tangent(material, F, state_old, P)
+    AA = _j2_tangent_analytical(material, F, state_old, P, state_new)
     return W, P, AA, state_new
 end
 

src/utils.jl

@@ -149,31 +149,6 @@ function configure_logger()
     return nothing
 end
 
-# Constants for Linux/glibc FPE trap masks
-const FE_INVALID = 1    # 0x01
-const FE_DIVBYZERO = 4  # 0x04
-const FE_OVERFLOW = 8   # 0x08
-
-"""
-    enable_fpe_traps()
-
-Attempts to enable floating-point exceptions (FPE traps) on supported platforms.
-Catches invalid operations, divide-by-zero, and overflow. Only supported on Linux x86_64.
-Warns that parser behavior may break due to this setting.
-"""
-function enable_fpe_traps()
-    if Sys.islinux() && Sys.ARCH == :x86_64
-        norma_log(0, :warning, "Enabling FPE traps can break Julia's parser if done too early")
-        norma_log(0, :warning, "(e.g., before Meta.parse()).")
-        mask = FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW
-        ccall((:feenableexcept, "libm.so.6"), Cuint, (Cuint,), mask)
-        norma_log(0, :info, "Floating-point exceptions enabled (invalid, div-by-zero, overflow)")
-    else
-        norma_log(0, :warning, "FPE trap support not available on this platform: $(Sys.KERNEL) $(Sys.ARCH)")
-    end
-    return nothing
-end
-
 function format_time(seconds::Float64)::String
     days = floor(Int, seconds / 86400)  # 86400 seconds in a day
     seconds %= 86400

test/benchmark-j2-tangent.jl

@@ -0,0 +1,110 @@
+# Norma: Copyright 2025 National Technology & Engineering Solutions of
+# Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS,
+# the U.S. Government retains certain rights in this software. This software
+# is released under the BSD license detailed in the file license.txt in the
+# top-level Norma.jl directory.
+
+# Benchmark: finite-difference vs analytical consistent tangent for J2Plasticity.
+#
+# For each of two representative states (elastic step, plastic step):
+#   • Verifies that the analytical and FD tangents agree to within FD truncation error.
+#   • Times both implementations over N_REPS repetitions and reports the speedup.
+#
+# Run standalone:
+#   cd test && julia --project=.. benchmark-j2-tangent.jl
+
+using LinearAlgebra
+using Printf
+using StaticArrays
+using Test
+
+include("../src/Norma.jl")
+
+const N_REPS = 500
+
+# ── helpers ────────────────────────────────────────────────────────────────────
+
+function time_reps(f, args...; n=N_REPS)
+    f(args...)            # warm-up / JIT
+    t = @elapsed for _ in 1:n
+        f(args...)
+    end
+    return t / n          # seconds per call
+end
+
+function max_rel_err(A::SArray{Tuple{3,3,3,3}}, B::SArray{Tuple{3,3,3,3}})
+    return norm(A - B) / max(norm(A), norm(B), 1.0)
+end
+
+# ── material ───────────────────────────────────────────────────────────────────
+
+params = Norma.Parameters(
+    "model"             => "j2 plasticity",
+    "elastic modulus"   => 200.0e9,
+    "Poisson's ratio"   => 0.3,
+    "density"           => 7800.0,
+    "yield stress"      => 250.0e6,
+    "hardening modulus" => 20.0e9,
+)
+mat = Norma.J2Plasticity(params)
+
+# ── deformation states ─────────────────────────────────────────────────────────
+
+# State 1 – elastic: small uniaxial stretch (well below yield, σy = 250 MPa)
+F_el = SMatrix{3,3,Float64,9}(
+    1.0,  0.0,  0.0,
+    0.0,  1.0,  0.0,
+    0.0,  0.0,  1.001,
+)
+s_el = Norma.initial_state(mat)   # Fp = I, eqps = 0
+
+# State 2 – plastic: larger strain, pre-loaded plastic state
+F_pre = SMatrix{3,3,Float64,9}(
+    0.98,  0.02,  0.0,
+    0.0,   0.97,  0.01,
+    0.0,   0.0,   1.05,
+)
+_, _, s_pre = Norma._j2_stress(mat, F_pre, s_el)
+
+F_pl = SMatrix{3,3,Float64,9}(
+    0.97,  0.03,  0.0,
+    0.0,   0.96,  0.02,
+    0.0,   0.0,   1.08,
+)
+
+# ── benchmark ──────────────────────────────────────────────────────────────────
+
+println("\n── J2 Plasticity Tangent Benchmark  (N_REPS = $N_REPS) ──\n")
+@printf("  %-16s  %12s  %12s  %8s  %12s\n",
+        "case", "FD (μs)", "Analytic (μs)", "Speedup", "Rel error")
+println("  ", "─"^68)
+
+@testset "J2 Tangent: analytical vs FD" begin
+    for (label, F, s_old) in [
+            ("elastic step", F_el, s_el),
+            ("plastic step", F_pl, s_pre),
+        ]
+
+        W, P, s_new = Norma._j2_stress(mat, F, s_old)
+
+        AA_fd = Norma._j2_tangent_fd(mat, F, s_old, P)
+        AA_an = Norma._j2_tangent_analytical(mat, F, s_old, P, s_new)
+
+        t_fd = time_reps(Norma._j2_tangent_fd, mat, F, s_old, P)
+        t_an = time_reps(Norma._j2_tangent_analytical, mat, F, s_old, P, s_new)
+
+        speedup = t_fd / t_an
+        err     = max_rel_err(AA_fd, AA_an)
+
+        @printf("  %-16s  %12.2f  %12.2f  %7.1fx  %12.2e\n",
+                label, t_fd*1e6, t_an*1e6, speedup, err)
+
+        # The analytical tangent uses the frozen-Fᵖ approximation (Fᵖ_new treated
+        # as constant when differentiating P).  For elastic steps this is exact
+        # (no plastic flow, so frozen-Fᵖ = true tangent); for plastic steps the
+        # error is O(Δεᵖ) relative to the full consistent tangent.
+        tol = label == "elastic step" ? 1e-4 : 1e-2
+        @test err < tol
+    end
+end
+println()

test/constitutive.jl

@@ -215,3 +215,34 @@ end
     σvm = sqrt(1.5) * norm(σdev)
     @test isapprox(σvm, σy; rtol=1e-3)
 end
+
+@testset "J2Plasticity FD Tangent" begin
+    # Verify _j2_tangent_fd agrees with the analytical tangent for both
+    # elastic and plastic steps. This keeps coverage on the FD helper,
+    # which is retained as a reference implementation for future models.
+    params = Norma.Parameters(
+        "elastic modulus" => 200.0e9, "Poisson's ratio" => 0.3,
+        "density" => 7800.0, "yield stress" => 250.0e6, "hardening modulus" => 20.0e9,
+    )
+    mat = Norma.J2Plasticity(params)
+    state0 = Norma.initial_state(mat)
+
+    # Elastic step
+    F_el = @SMatrix [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.001]
+    W, P, state_new = Norma._j2_stress(mat, F_el, state0)
+    AA_fd = Norma._j2_tangent_fd(mat, F_el, state0, P)
+    AA_an = Norma._j2_tangent_analytical(mat, F_el, state0, P, state_new)
+    @test size(AA_fd) == (3, 3, 3, 3)
+    @test norm(AA_fd - AA_an) / max(norm(AA_fd), norm(AA_an), 1.0) < 1e-4
+
+    # Plastic step
+    F_pre = @SMatrix [0.98 0.02 0.0; 0.0 0.97 0.01; 0.0 0.0 1.05]
+    _, _, s_pre = Norma._j2_stress(mat, F_pre, state0)
+    F_pl = @SMatrix [0.97 0.03 0.0; 0.0 0.96 0.02; 0.0 0.0 1.08]
+    W2, P2, state_new2 = Norma._j2_stress(mat, F_pl, s_pre)
+    AA_fd2 = Norma._j2_tangent_fd(mat, F_pl, s_pre, P2)
+    AA_an2 = Norma._j2_tangent_analytical(mat, F_pl, s_pre, P2, state_new2)
+    @test size(AA_fd2) == (3, 3, 3, 3)
+    # Frozen-Fᵖ approximation introduces O(Δεᵖ) error on plastic steps
+    @test norm(AA_fd2 - AA_an2) / max(norm(AA_fd2), norm(AA_an2), 1.0) < 1e-2
+end

test/runtests.jl

@@ -66,7 +66,7 @@ const indexed_test_files = [
     (50, "smoothing.jl"),
     (51, "nnopinf-schwarz-overlap-cuboid-hex8.jl"),
     (52, "single-static-solid-j2-plasticity.jl"),
-    (53, "utils.jl"), # Must go last due to FPE traps
+    (53, "utils.jl"),
 ]
 
 # Extract test file names

test/utils.jl

@@ -63,14 +63,4 @@ using Logging
         end
     end
 
-    @testset "Enable Fpe Traps" begin
-        # This is platform-specific and side-effect prone
-        # So we test only that it runs without error
-        try
-            Norma.enable_fpe_traps()
-            @test true
-        catch e
-            @test false
-        end
-    end
 end