juFrac2D/verify.jl at main · MatMechLab/juFrac2D · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
# ===========================================================================
#  PhaseFieldFracture — a readable, pure-Julia phase-field fracture solver
#
#  Author:     Yang Bai — Materials Mechanics Laboratory (MMLab)
#  Contact:    yangbai90@outlook.com
#  Copyright:  © 2026 Materials Mechanics Laboratory (MMLab). All rights reserved.
# ===========================================================================

# ===========================================================================
#  verify.jl — self-contained correctness checks (run:  julia verify.jl)
# ===========================================================================
#
#  This script proves, with only the Julia standard library, that the Miehe
#  spectral tension/compression split in `strainsplit.jl` is implemented
#  correctly.  It checks, at several strain states:
#
#    (A) the split adds up:   ψ⁺ + ψ⁻ = ψ ,   σ⁺ + σ⁻ = D ε          (exact)
#    (B) the derivatives are consistent:  σ = ∂ψ/∂ε ,  C = ∂σ/∂ε     (vs finite diff.)
#    (C) the physics is right: pure compression gives ψ⁺ = σ⁺ = C⁺ = 0
#        (no tensile driving, full compressive stiffness kept), tangents are SPD
#    (D) `split = :none` reproduces the classic no-split model
#    (E) end-to-end: the SAME notched specimen cracks in TENSION but not in
#        COMPRESSION with the split — while with :none it cracks in both.
#
#  Every check prints PASS/FAIL and the script exits non-zero if any fails.

include("PhaseFieldFracture.jl")
using .PhaseFieldFracture
using .PhaseFieldFracture: spectral_stress_tangent, driving_energy_spectral,
                           driving_energy, constitutive, strain_energy
using LinearAlgebra, Printf

const PF = PhaseFieldFracture
npass = 0; nfail = 0
function check(name, ok; detail = "")
    global npass, nfail
    ok ? (npass += 1) : (nfail += 1)
    @printf("  [%s] %-52s %s\n", ok ? "PASS" : "FAIL", name, detail)
end

# A material to test with (steel-like, plane strain); mat.λ, mat.μ, mat.D are ready to use.
mat = Material(E = 210_000.0, ν = 0.30, Gc = 2.7, ℓ = 0.04)
λ, μ, D = mat.λ, mat.μ, mat.D
σ_scale = mat.E * 2e-3                         # ~420, a physical stress magnitude for scaling

# Central finite differences (engineering strain ε = [εₓₓ, ε_yy, γₓy]).
basis = [Float64.(1:3 .== i) for i in 1:3]
fd_grad(f, ε; h = 1e-8) = [(f(ε + h*e) - f(ε - h*e)) / (2h) for e in basis]
fd_jac(vf, ε; h = 1e-8) = hcat([(vf(ε + h*e) - vf(ε - h*e)) / (2h) for e in basis]...)

# split-energy helpers (ψ⁺ from the code; ψ⁻ = ψ − ψ⁺ by the sum identity we also test)
ψpos(ε) = driving_energy_spectral(λ, μ, ε)
ψneg(ε) = strain_energy(mat, ε) - ψpos(ε)

println("\n(A) split identities  ψ⁺+ψ⁻=ψ ,  σ⁺+σ⁻=Dε   (should be ~1e-16) ------------")
states = [[1e-3,0,0], [-1e-3,0,0], [1e-3,2e-3,0], [-1e-3,-2e-3,0],
          [1e-3,-0.5e-3,0.8e-3], [0,0,1e-3], [1e-3,1e-3,0], [7e-4,-3e-4,5e-4]]
eA = 0.0
for ε in states
    σp, σn, Cp, Cn = spectral_stress_tangent(λ, μ, ε)
    eψ = abs(ψpos(ε) + ψneg(ε) - strain_energy(mat, ε)) / max(strain_energy(mat, ε), 1e-20)
    eσ = norm(σp + σn - D*ε) / σ_scale
    # C⁺+C⁻ = D holds wherever trε ≠ 0 (at trε=0 the volumetric tangent term vanishes)
    eC = abs(ε[1] + ε[2]) < 1e-14 ? 0.0 : norm(Cp + Cn - D) / norm(D)
    global eA = max(eA, eψ, eσ, eC)
end
check("ψ⁺+ψ⁻=ψ and σ⁺+σ⁻=Dε at 8 strain states", eA < 1e-10, detail = @sprintf("max err %.1e", eA))

println("\n(B) consistent derivatives  σ=∂ψ/∂ε ,  C=∂σ/∂ε   (vs finite differences) ----")
smooth = [[2e-3,1e-3,0.5e-3], [-2e-3,-1e-3,0.3e-3], [3e-3,-1e-3,1e-3], [-6e-4,9e-4,-4e-4]]
eB = 0.0
for ε in smooth
    σp, σn, Cp, Cn = spectral_stress_tangent(λ, μ, ε)
    e1 = norm(σp - fd_grad(ψpos, ε)) / σ_scale
    e2 = norm(σn - fd_grad(ψneg, ε)) / σ_scale
    e3 = norm(Cp - fd_jac(e->spectral_stress_tangent(λ,μ,e)[1], ε)) / norm(D)
    e4 = norm(Cn - fd_jac(e->spectral_stress_tangent(λ,μ,e)[2], ε)) / norm(D)
    global eB = max(eB, e1, e2, e3, e4)
end
check("σ⁺,σ⁻,C⁺,C⁻ match finite differences", eB < 1e-5, detail = @sprintf("max err %.1e", eB))

println("\n(C) physics of the split  ---------------------------------------------------")
εc = [-2e-3, -1e-3, 0.3e-3]                     # fully compressive state (both εₐ<0, trε<0)
σp, σn, Cp, Cn = spectral_stress_tangent(λ, μ, εc)
check("compression ⇒ ψ⁺ = 0 (no tensile driving)", ψpos(εc) == 0.0)
check("compression ⇒ σ⁺ = 0 and C⁺ = 0", norm(σp) == 0.0 && norm(Cp) == 0.0)
check("compression ⇒ full compressive stiffness kept (σ⁻=Dε, C⁻=D)",
      norm(σn - D*εc)/σ_scale < 1e-12 && norm(Cn - D)/norm(D) < 1e-12)
εt = [2e-3, 1e-3, 0.5e-3]                        # fully tensile state
σp, σn, Cp, Cn = spectral_stress_tangent(λ, μ, εt)
# in pure tension nothing lives in the compressive part: ψ⁺=ψ, σ⁺=Dε, σ⁻≈0
check("tension ⇒ ψ⁺ = ψ and σ⁺ = Dε (compressive part empty)",
      abs(ψpos(εt) - strain_energy(mat, εt)) / strain_energy(mat, εt) < 1e-12 &&
      norm(σp - D*εt) / σ_scale < 1e-12 && norm(σn) / σ_scale < 1e-12)
εm = [3e-3, -1e-3, 1e-3]                          # mixed state
_, _, Cp, Cn = spectral_stress_tangent(λ, μ, εm)
check("tangents C⁺, C⁻ are positive semi-definite",
      minimum(eigvals(Symmetric(Cp))) ≥ -1e-6 && minimum(eigvals(Symmetric(Cn))) ≥ -1e-6)

println("\n(D) split = :none reproduces the classic model  -----------------------------")
matn = Material(E = 210_000.0, ν = 0.30, Gc = 2.7, ℓ = 0.04, split = :none)
eD = 0.0
for ε in smooth
    σ, C = constitutive(matn, ε, 0.7)            # degraded stress/tangent with g = 0.7
    global eD = max(eD, norm(σ - 0.7*(D*ε))/σ_scale, norm(C - 0.7*D)/norm(D),
                    abs(driving_energy(matn, ε) - strain_energy(matn, ε)))
end
check(":none gives σ = gDε, C = gD, ψ⁺ = ½εᵀDε", eD < 1e-12, detail = @sprintf("max err %.1e", eD))

println("\n(E) end-to-end: the same notch cracks in TENSION, not in COMPRESSION  --------")
function growth(; split, umax)
    mesh  = rectangle_mesh(20, 20; x0=-0.5, x1=0.5, y0=-0.5, y1=0.5)
    crack = crack_line_nodes(mesh; crack_y=0.0, tip_x=0.0)
    m     = Material(E=210_000.0, ν=0.30, Gc=2.7, ℓ=0.05, k=1e-8, split=split)
    set   = Settings(nsteps=6, umax=umax, max_iter=60, tol=1e-6, out_every=10_000)
    prob  = Problem(mesh, m, set, crack, tempname(), tempname())
    # zero-load notch halo baseline (AT2 damage around the pinned crack, before loading)
    H0    = zeros(4, ncells(mesh)); φ0 = zeros(nnodes(mesh)); for n in crack; φ0[n]=1.0; end
    halo  = max_free_phi(PF.solve_phase(mesh, m, H0, φ0, crack), crack)
    local φ
    redirect_stdout(devnull) do; _, φ, _ = run_simulation(prob); end
    return max_free_phi(φ, crack) - halo          # damage GROWTH above the static halo
end
gt  = growth(split=:spectral, umax=+8e-3)
gc  = growth(split=:spectral, umax=-8e-3)
gcn = growth(split=:none,     umax=-8e-3)
@printf("    spectral tension     Δφ = %.3f   (want > 0.4, crack runs)\n", gt)
@printf("    spectral compression Δφ = %.3f   (want < 0.1, crack does NOT run)\n", gc)
@printf("    none     compression Δφ = %.3f   (want > 0.4, spurious cracking)\n", gcn)
check("split: tension cracks but compression does not", gt > 0.4 && gc < 0.1)
check("no-split: compression spuriously cracks (the bug the split fixes)", gcn > 0.4)

@printf("\n==================  %d passed, %d failed  ==================\n", npass, nfail)
nfail == 0 ? println("ALL CHECKS PASS ✓") : (println("SOME CHECKS FAILED ✗"); exit(1))