PhaseFieldFracture

A small, readable, pure‑Julia solver for 2‑D brittle fracture using the phase‑field (AT2) method with a small‑strain, plane‑strain elastic bulk and a tension/compression split of the elastic energy (Miehe's spectral decomposition) so that only tensile strain drives — and is degraded by — the crack. It has zero external dependencies — it uses only the Julia standard library (LinearAlgebra, SparseArrays, Printf) — and writes its own ParaView (.vtu/.pvd) and CSV output by hand.

The reference problem is the classic single‑edge‑notched tension (SENT) test: a square plate with a horizontal pre‑crack is pulled apart until the crack runs across the ligament.

• Author: Yang Bai — Materials Mechanics Laboratory (MMLab)
• Contact: yangbai90@outlook.com
• Copyright: © 2026 Materials Mechanics Laboratory (MMLab). All rights reserved.

---

Table of contents

1. What this code does
2. Quick start
3. The physical model
4. The numerical method
5. Code architecture (file by file)
6. Public API reference
7. Usage & customization
8. Output files & post‑processing
9. The SENT benchmark & expected results
10. Parallel execution
11. Performance & scaling
12. Critical review: assumptions & limitations
13. Extending the code
14. Troubleshooting / FAQ
15. References
16. License & citation

---

1. What this code does

Given a rectangular plate, a material, and a loading program, the solver computes the coupled evolution of

• the displacement field u(x) (plane‑strain linear elasticity), and
• the damage / phase field φ(x) ∈ [0,1] (φ = 0 intact, φ = 1 fully broken),

as the imposed boundary displacement is increased in quasi‑static load steps. It tracks an irreversible history field H (the largest elastic energy density ever seen at each integration point) that drives — and never un‑drives — the damage, so cracks grow but never heal.

Per load step it reports a force–displacement point (reaction vs. imposed displacement) and periodically saves full nodal fields as a ParaView time series you can animate.

Everything is written for readability: each source file corresponds to exactly one mathematical ingredient of the model, and every function carries a docstring explaining the mathematics it implements.

---

2. Quick start

Requirements: a Julia install (developed and tested on Julia 1.12; any recent 1.x should work). No packages to add — the three modules it uses ship with Julia itself.

# from the directory that contains run.jl
julia run.jl

This runs the shipped SENT example (a 100×100 mesh, 50 load steps). It prints a per‑step log to the terminal and writes:

output/field_0000.vtu … field_0050.vtu     # nodal snapshots
output/timeseries.pvd                       # ParaView time series (open THIS)
HistoryOutput.txt                           # CSV force/damage history

Open output/timeseries.pvd in ParaView to watch the crack grow; plot HistoryOutput.txt (columns described below) to get the force–displacement curve.

Check correctness first (recommended for students):

julia verify.jl     # proves the tension/compression split is correct, in ~1 min

verify.jl checks the split identities and derivatives against finite differences and runs a tiny tension‑vs‑compression demonstration — a good place to start reading, and to confirm your install works.

Tip — faster first run. Julia compiles on first use. To try the model quickly, edit run.jl to a coarser mesh and fewer steps, e.g. rectangle_mesh(40, 40; …) and Settings(nsteps = 20, …).

---

3. The physical model

3.1 Energy functional

The model is the standard variational (Bourdin/Francfort–Marigo) brittle‑fracture regularization, with the tension/compression split of Miehe et al. (2010). The total energy of the body Ω is

E[u, φ] = ∫_Ω [ g(φ) ψ⁺(ε(u)) + ψ⁻(ε(u)) ] dΩ  +  Gc ∫_Ω ( φ²/(2ℓ) + (ℓ/2) |∇φ|² ) dΩ
          └──────── stored elastic energy ────────┘   └──── regularized fracture energy (AT2) ────┘

with

Symbol	Meaning	In code
ψ⁺, ψ⁻	tensile / compressive parts of the strain‑energy density (§3.3)	driving_energy, spectral_stress_tangent (strainsplit.jl)
ψ(ε) = ½ εᵀ D ε	total strain‑energy density, ψ = ψ⁺ + ψ⁻	strain_energy (elasticity.jl)
g(φ) = (1−φ)² + k	stiffness degradation function (acts on ψ⁺ only)	degrade (elasticity.jl)
D	plane‑strain elasticity matrix built from (E, ν)	Material (elasticity.jl)
Gc	critical energy release rate (fracture toughness)	Material.Gc
ℓ	phase‑field length scale (crack‑band width)	Material.ℓ
k	small residual stiffness so g(1) = k > 0 keeps Kᵤ regular	Material.k

Only the tensile energy ψ⁺ is multiplied by g(φ): a fully broken point loses its tensile stiffness but keeps its compressive stiffness, so crack faces can reopen but not interpenetrate, and compression alone cannot drive damage (§3.3). Setting split = :none in Material recovers the classic model in which the whole energy ψ drives and is degraded (ψ⁺ = ψ, ψ⁻ = 0).

The φ²/(2ℓ) + (ℓ/2)|∇φ|² term is the Ambrosio–Tortorelli “AT2” crack‑surface density: as ℓ → 0 its integral Γ‑converges to Gc times the crack length, so φ is a smeared, mesh‑resolvable stand‑in for a sharp crack of width ≈ ℓ.

3.2 Elasticity (plane strain)

Engineering strain ε = [εₓₓ, ε_yy, γₓy] with γₓy = ∂u/∂y + ∂v/∂x, and stress σ = D ε with

        E                ⎡ 1−ν    ν      0    ⎤
D = ───────────────── ·  ⎢  ν    1−ν     0    ⎥          (plane strain, 0 ≤ ν < ½)
    (1+ν)(1−2ν)          ⎣  0     0   (1−2ν)/2 ⎦

Equivalently, in Lamé form with λ = Eν/[(1+ν)(1−2ν)] and μ = E/[2(1+ν)], σ = λ (tr ε) I + 2μ ε and ψ = ½ λ (tr ε)² + μ (ε:ε). The split in §3.3 is built on this Lamé form. Only the tensile part of the stiffness is degraded by g(φ), so a fully broken point (φ = 1) keeps its full compressive stiffness and retains only the residual tensile stiffness k.

3.3 Tension–compression split (Miehe spectral decomposition)

The problem it fixes. If the whole energy ψ drives the crack and the whole stiffness is degraded, then (i) a crack can nucleate and grow under pure compression, and (ii) the faces of an open crack can freely interpenetrate, because the material has lost all stiffness even against closing load. Real brittle cracks only open and grow under tension.

The remedy (Miehe et al. 2010). Diagonalise the strain tensor into its principal strains εₐ and principal directions nₐ (its spectral decomposition),

ε = Σₐ εₐ (nₐ ⊗ nₐ) ,

and split it with the Macaulay brackets ⟨x⟩₊ = max(x,0), ⟨x⟩₋ = min(x,0) (note ⟨x⟩₊ + ⟨x⟩₋ = x), keeping only positive / only negative principal strains:

ε⁺ = Σₐ ⟨εₐ⟩₊ (nₐ ⊗ nₐ) ,     ε⁻ = Σₐ ⟨εₐ⟩₋ (nₐ ⊗ nₐ) .

The elastic energy then splits into a tensile and a compressive part (the exact definitions this code implements):

ψ⁺(ε) = ½ λ ⟨tr ε⟩₊²  +  μ Σₐ ⟨εₐ⟩₊² ,
ψ⁻(ε) = ½ λ ⟨tr ε⟩₋²  +  μ Σₐ ⟨εₐ⟩₋² ,        ψ⁺ + ψ⁻ = ψ  (since ⟨x⟩₊²+⟨x⟩₋²=x²).

Only ψ⁺ is degraded, ψ(ε,φ) = g(φ) ψ⁺(ε) + ψ⁻(ε), and only ψ⁺ drives the damage (§3.4). Differentiating the degraded energy gives the stress and its consistent tangent (needed because the split makes σ a nonlinear function of ε):

σ = g(φ) σ⁺ + σ⁻ ,     σ⁺ = λ⟨tr ε⟩₊ I + 2μ ε⁺ ,     σ⁻ = λ⟨tr ε⟩₋ I + 2μ ε⁻ ,
C = g(φ) C⁺ + C⁻ ,     C⁺ = λ H(tr ε) I⊗I + 2μ ℙ⁺ ,   C⁻ = λ H(−tr ε) I⊗I + 2μ (𝕀 − ℙ⁺) ,

where H is the Heaviside step and ℙ⁺ = ∂ε⁺/∂ε is the fourth‑order positive projection tensor of the strain (built in closed form from the eigenprojections; see strainsplit.jl for the exact formula and the derivation).

Plane strain. Out of plane ε_zz = 0, so 0 is automatically a third principal strain; since ⟨0⟩₊ = ⟨0⟩₋ = 0 it contributes nothing, and the whole split reduces to the in‑plane 2×2 strain tensor and its two principal strains — no special out‑of‑plane term is needed.

Where it lives. strainsplit.jl implements ψ⁺ (driving_energy_spectral) and (σ⁺, σ⁻, C⁺, C⁻) (spectral_stress_tangent). The functions are verified against finite differences (σ = ∂ψ/∂ε, C = ∂σ/∂ε) and against the exact identities ψ⁺+ψ⁻ = ψ, σ⁺+σ⁻ = Dε — run julia verify.jl.

3.4 Damage evolution & the history field

Minimizing E over φ at fixed u, and replacing the tensile driving energy ψ⁺ by the irreversible history field

H(x, t) = max_{s ≤ t}  ψ⁺(ε(u(x, s)))

gives a linear elliptic equation for the damage (this is what makes each sub‑solve cheap):

(2H + Gc/ℓ) φ  −  Gc ℓ Δφ  =  2H          in Ω

with φ = 1 prescribed on the pre‑crack and natural (zero‑flux) conditions elsewhere. The history field enters as both a reaction coefficient and a source.

Irreversibility (no healing) is enforced in two complementary ways:

1. H is a running maximum in time, so the driving force never decreases; and
2. after each damage solve the nodal field is projected, φ ← clamp(max(φ, φ_prev), 0, 1), so damage cannot fall below its value at the start of the load step.

3.5 Coupled system

The two stationarity conditions,

δE/δu = 0   →   div( σ ) = 0 ,   σ = g(φ) σ⁺(ε) + σ⁻(ε)   (degraded elasticity)
δE/δφ = 0   →   (2H + Gc/ℓ) φ − Gc ℓ Δφ = 2H             (AT2 damage)

are coupled through g(φ) (damage softens the tensile response) and H(ψ⁺(ε(u))) (tensile strain drives the damage). They are solved by alternate minimization (see §4.3). Note the mechanical equilibrium is now nonlinear in u — because σ⁺, σ⁻ depend nonlinearly on ε through the spectral split — so the displacement sub‑problem is solved by Newton's method with the consistent tangent C = g(φ) C⁺ + C⁻.

---

4. The numerical method

4.1 Finite element discretization

• Element: 4‑node bilinear quadrilateral (Q4), isoparametric.

• Reference shape functions on (ξ,η) ∈ [−1,1]²:

  N₁ = ¼(1−ξ)(1−η)   N₂ = ¼(1+ξ)(1−η)   N₃ = ¼(1+ξ)(1+η)   N₄ = ¼(1−ξ)(1+η)
  

• Quadrature: 2×2 Gauss (four points at ±1/√3, each weight 1), exact for the bilinear element. dΩ = detJ · weight.

• Isoparametric map & physical gradients: with the Jacobian J = ∂(x,y)/∂(ξ,η), the physical shape‑function gradients are ∇N = ∇_ref N · J⁻¹ (implemented as the right‑division ref∇N / J). A non‑positive detJ aborts the run (a node‑ordering guard).

• Strain–displacement operator B (3×8) such that ε = B uₑ (strain_matrix).

• Fields: displacement is Q4 vector‑valued (2 dofs/node, interleaved as (2n−1, 2n) = (uₓ, u_y)); damage is Q4 scalar (1 dof/node).

4.2 Assembly

Both stiffness matrices are assembled from COO triplets (I, J, V) and handed to sparse(I, J, V, …), which sums duplicates:

• Elasticity (assemble_mechanics): at the current (u, φ) it assembles both the tangent stiffness K = Σₑ ∫ Bᵀ C B dΩ and the internal force f_int = Σₑ ∫ Bᵀ σ dΩ, where σ = g(φ) σ⁺ + σ⁻ and C = g(φ) C⁺ + C⁻ come from the split (constitutive), with φ interpolated to each Gauss point φ(ξ,η) = N·φₑ. Each cell contributes an 8×8 block (64 entries). f_int is the Newton residual; its sum over the top‑edge dofs is the reaction force.

• Phase field (phasefield_system): on each cell,

  Kₑ[a,b] = ∫ ( (2H + Gc/ℓ) Nₐ N_b  +  Gc ℓ ∇Nₐ·∇N_b ) dΩ      (mass + diffusion)
  fₑ[a]   = ∫ ( 2H Nₐ ) dΩ                                       (source)
  

  a 4×4 block plus a 4‑vector source, with H taken per Gauss point.

Both Kᵤ and the phase‑field K are symmetric positive definite (the residual k and the Gc/ℓ reaction term keep them regular even at full damage).

4.3 Staggered (alternate‑minimization) solver

The coupled minimization is nonlinear only through the g(φ)ψ product. The code uses the robust, industry‑standard staggered scheme: freeze one field, solve the (now linear) problem for the other, and iterate. Per load step:

repeat until converged (or max_iter):
    (1) solve  div σ(ε(u), φ) = 0   with Dirichlet BCs      → displacement u
        by NEWTON:  repeat  K Δu = −f_int,  u ← u + Δu  until ‖f_int_free‖ small
    (2) H_trial ← max(H, ψ⁺(ε(u)))  at every Gauss point    → tensile driving force
    (3) solve  K(H_trial) φ = f(H_trial)                    → damage φ,
        then project: φ ← clamp(max(φ, φ_committed), 0, 1), φ = 1 on the notch
    convergence test:  max( ‖Δu‖/‖u‖ , ‖Δφ‖/‖φ‖ ) < tol
commit history:  H ← H_trial

Step (1) is a Newton loop because the spectral split makes the stress nonlinear in the strain; the potential g ψ⁺ + ψ⁻ is convex in ε, so Newton with the consistent tangent converges robustly (and in a single step when split = :none, where σ = g D ε is linear).

Loading is a linear displacement ramp: u_y^top = umax · step / nsteps. The bottom edge is clamped (uₓ = u_y = 0); the top edge is given the vertical displacement u_y with uₓ left free (Poisson contraction allowed).

4.4 Linear solver & boundary conditions

Each linear system is solved by solve_dirichlet(K, f, bc), which performs a static condensation into free/fixed dofs and solves the reduced system with Julia’s sparse backslash \ (a direct sparse LU factorization via SuiteSparse):

x_fixed = prescribed values
K_ff x_free = f_free − K_fc x_fixed        # the reduced system actually solved

---

5. Code architecture (file by file)

The module is deliberately split so each file is one ingredient of the model. PhaseFieldFracture.jl is the umbrella module that includes the rest in dependency order and exports the public names.

File	Lines¹	Responsibility	Key definitions
PhaseFieldFracture.jl	~50	Module: includes, exports, top‑level docs	module PhaseFieldFracture, export …
mesh.jl	~100	Structured Q4 mesh of a rectangle; boundary node sets; dof maps; pre‑crack nodes	Mesh, rectangle_mesh, crack_line_nodes, nnodes, ncells, xdof, ydof, grid_node
element.jl	~90	Q4 shape functions, isoparametric map, quadrature, B‑matrix	GAUSS, shape, element, strain_matrix, cell_dofs
strainsplit.jl	~200	Miehe spectral tension/compression split: ψ⁺, split stress σ±, consistent tangent C±	principal_strains, driving_energy_spectral, spectral_stress_tangent
elasticity.jl	~200	Plane‑strain material (+λ,μ,split), degradation g(φ), constitutive combine, tangent+residual assembly, Newton displacement solve	Material, degrade, driving_energy, constitutive, assemble_mechanics, displacement_bcs, solve_displacement
phasefield.jl	~115	AT2 system assembly, history update (driven by ψ⁺), damage solve with irreversibility	phasefield_system, update_history!, solve_phase
solver.jl	~175	Settings/Problem types, Dirichlet solve, staggered load‑stepping driver	Settings, Problem, rel_change, solve_dirichlet, run_simulation
output.jl	~210	Scalar monitors, CSV history, hand‑written VTU + PVD writers	reaction_top, crack_tip_x, max_free_phi, write_history_header, log_step, write_vtu, write_pvd
verify.jl	~110	Self‑contained correctness checks (split identities, FD derivatives, tension‑vs‑compression)	(script, not a module)
run.jl	~50	Driver script for the SENT example	(script, not a module)

¹ approximate.

Dependency flow: mesh → element → strainsplit → elasticity → phasefield → solver → output, with run.jl (and the standalone verify.jl) sitting on top. strainsplit.jl is pure math (no mesh/FE dependency) so its formulas can be unit‑tested in isolation; nothing depends on output.jl except the driver, so you can swap the I/O layer freely.

---

6. Public API reference

Everything below is exported by using .PhaseFieldFracture.

Types

Mesh                     # coords, cells, nx, ny, and boundary node vectors top/bottom/left/right
Material(; E, ν, Gc, ℓ, k = 1e-8, split = :spectral)   # builds D, λ, μ from (E, ν)
                         # split = :spectral (Miehe tension/compression split) | :none (no split)
Settings(; nsteps = 50, umax = 8e-3, max_iter = 50, tol = 1e-6, out_every = 5)
Problem(mesh, material, settings, crack_nodes, outdir, histfile)   # positional constructor

Settings and Material use keyword constructors; Problem is a plain struct built positionally.

Mesh construction

rectangle_mesh(nx, ny; x0 = -0.5, x1 = 0.5, y0 = -0.5, y1 = 0.5) -> Mesh
crack_line_nodes(mesh; crack_y = 0.0, tip_x = 0.0) -> Vector{Int}
nnodes(mesh) -> Int
ncells(mesh) -> Int

crack_line_nodes returns the nodes of a straight horizontal pre‑crack on y = crack_y from the left edge to x = tip_x; these are pinned to φ = 1. Use an even ny so that mesh nodes fall exactly on the crack line (the function errors out otherwise).

Running a simulation

run_simulation(problem::Problem) -> (u, φ, H)

Runs the full staggered load‑stepping loop and returns the final displacement vector u (length 2·nnodes), damage vector φ (length nnodes), and history field H (size 4 × ncells). Along the way it writes the CSV history, the .vtu snapshots, and the .pvd collection.

Post‑processing monitors

reaction_top(mesh, f_int) -> Float64     # Σ vertical internal force over top-edge nodes
crack_tip_x(mesh, φ; threshold = 0.95) -> Float64   # right-most x with φ ≥ threshold
max_free_phi(φ, crack_nodes) -> Float64  # largest damage OUTSIDE the pre-crack (growth monitor)

Here f_int = Σₑ ∫ Bᵀ σ dΩ is the assembled internal force returned by run_simulation's displacement solve; its vertical sum over the top edge is the reaction (it is no longer Kᵤ·u, since the split makes the material nonlinear).

Output writers

write_vtu(path, mesh, u, φ) -> path      # one ParaView UnstructuredGrid snapshot (ASCII XML)
write_pvd(path, entries)    -> path      # collection tying snapshots into a time series

---

7. Usage & customization

The driver run.jl is the template. Its four blocks are all you normally touch:

include("PhaseFieldFracture.jl")
using .PhaseFieldFracture

# 1) Geometry & mesh  (keep ny even so nodes lie on the crack line y = 0)
mesh  = rectangle_mesh(100, 100; 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)

# 2) Material  (consistent units: E in MPa, lengths in mm, Gc in N/mm)
material = Material(E = 210_000.0, ν = 0.30, Gc = 2.7, ℓ = 0.04, k = 1e-8)

# 3) Solver settings
settings = Settings(nsteps = 50, umax = 8e-3, max_iter = 50, tol = 1e-6, out_every = 5)

# 4) Assemble & run
problem = Problem(mesh, material, settings, crack, "output", "HistoryOutput.txt")
u, φ, H = run_simulation(problem)

Choosing parameters

• Mesh resolution vs. ℓ. The length scale must be resolved by the mesh: keep the element size h ≲ ℓ/2. The shipped run has h = 1/100 = 0.01 and ℓ = 0.04, i.e. ~4 elements across the crack band — good. If you refine ℓ, refine the mesh with it.
• Units. The code is unit‑agnostic; just be consistent. The example uses MPa / mm / (N/mm), giving E = 210 GPa, Gc = 2.7 kJ/m² (steel‑like).
• Load steps. More nsteps (smaller Δu) improves the staggered convergence near the peak load and resolves the softening branch better.
• out_every. Snapshots are saved every out_every steps (and always on the final step). Snapshot 0 is the regularized notch before loading.

Setting up a different problem

• Different geometry / loading: change the domain in rectangle_mesh and the BCs in displacement_bcs (elasticity.jl). Any Dirichlet program can be expressed as a Dict(dof => value) and fed through solve_dirichlet.
• Different notch: pass a different crack_y/tip_x to crack_line_nodes, or supply your own Vector{Int} of pre‑broken nodes to Problem.
• No pre‑crack: pass an empty Int[] (nucleation will then be governed by the stress field and ℓ).

---

8. Output files & post‑processing

CSV load history (HistoryOutput.txt)

One header line plus one row per load step (comma‑separated):

Column	Name	Meaning
1	step	load‑step index (0 = pre‑load snapshot)
2	imposed_top_uy	prescribed top displacement (the “time”)
3	reaction_y_top	vertical reaction on the top edge (the force)
4	crack_tip_x	right‑most x with φ ≥ 0.95 (crack‑tip position)
5	max_phi	maximum damage anywhere (usually 1 on the notch)
6	max_free_phi	maximum damage outside the pre‑crack (growth monitor)
7	staggered_iterations	staggered iterations used this step
8	converged	true/false — did the staggered loop meet tol?

Columns 2 and 3 are your force–displacement curve. A quick plot with pure Julia:

using DelimitedFiles, Printf
data = readdlm("HistoryOutput.txt", ','; comments = true, comment_char = '#')
uy, R = data[:, 2], data[:, 3]
# … feed (uy, R) to your favourite plotting package …

ParaView time series (.vtu + .pvd)

• Each field_XXXX.vtu is a VTK UnstructuredGrid (ASCII XML, written by hand) holding the mesh plus two nodal fields: phase_field (scalar φ) and displacement (3‑vector, w = 0). Cells are VTK_QUAD (type 9); connectivity is 0‑based per the VTK convention.
• timeseries.pvd is a VTK Collection that lists the snapshots and assigns each a “time” equal to the imposed top displacement, so ParaView’s animation slider scrubs through the loading history.

To view: open timeseries.pvd in ParaView, colour by phase_field, and press play. Apply Warp By Vector on displacement to see the deformed shape.

---

9. The SENT benchmark & expected results

The shipped example is the single‑edge‑notched tension test — the standard phase‑field verification case of Miehe et al. (2010), now with the same spectral split they use:

• Square plate [−0.5, 0.5]², horizontal pre‑crack from the left edge to the centre (x = 0, y = 0).
• Bottom edge clamped; top edge pulled vertically to umax = 8×10⁻³.
• The crack grows straight to the right and splits the specimen.

Expected physics. The reaction rises almost linearly, reaches a peak load, then drops sharply as the crack propagates unstably across the ligament — the hallmark brittle response. crack_tip_x jumps from 0 to the right edge when the crack runs; max_free_phi climbs from ≈ 0 toward 1.

Illustrative numbers from a coarse, fast smoke run (30×30 mesh, 8 steps — not the shipped 100×100/50‑step configuration, which is sharper) show the characteristic shape clearly:

step	u_y	reaction	crack_tip_x	max_free_phi	converged
1	1.0e‑3	191.7	0.00	0.43	✓
2	2.0e‑3	373.1	0.00	0.46	✓
3	3.0e‑3	532.9	0.00	0.50	✓
4	4.0e‑3	654.5 ← peak	0.00	0.59	✗ (hit max_iter)
5	5.0e‑3	557.3	0.00	0.95	✓
6	6.0e‑3	519.8	0.00	0.97	✗
7	7.0e‑3	16.5 ← crack ran	0.50	0.98	✓
8	8.0e‑3	11.0	0.50	0.99	✓

The load collapses by ~40× once the crack traverses. Note that a few steps near the critical load stall at max_iter — this is the well‑known slow convergence of the staggered scheme at the onset of unstable growth (see §12), not a bug; finer load steps mitigate it.

The table above is a legacy no‑split smoke run; under tension the spectral split gives an almost identical, slightly stiffer curve (a coarse 24×24/10 run peaks near ~695 with the split vs ~678 without). The split’s decisive effect is under compression / crack closure: verify.jl loads the same notch in tension and in compression and shows the crack runs only in tension (Δφ ≈ 0.65) while compression does not (Δφ ≈ 0.03); with split = :none compression spuriously cracks (Δφ ≈ 0.65) because then ψ(ε) = ψ(−ε).

---

10. Parallel execution

Honest status: this code is currently single‑threaded and serial. There is no Threads.@threads, @distributed, MPI, or GPU code anywhere in it. This section explains (a) what parallelism you already get for free, and (b) exactly how to add more, because the code is structured to make that easy.

10.1 What runs in parallel today (implicitly)

The only parallelism you get out of the box is inside the linear‑algebra kernels:

• Sparse direct solve (\). The dominant cost per iteration is the sparse LU factorization in SuiteSparse. Its dense sub‑kernels call multithreaded BLAS (OpenBLAS by default). You can size that pool:

  using LinearAlgebra
  BLAS.set_num_threads(8)          # or export OPENBLAS_NUM_THREADS=8 before launch

  This helps the numerical factorization but does not parallelize the symbolic work or the assembly loops.

Starting Julia with -t N (JULIA_NUM_THREADS) has no effect until you add threaded code — see below.

10.2 The embarrassingly parallel case: parameter sweeps

If you want to run many independent simulations (mesh‑refinement studies, ℓ/Gc sweeps, multiple notch positions), that is trivially parallel and needs no code change — just give each run its own outdir/histfile. Two easy ways:

Shell level (simplest):

# runs several independent parameterized drivers at once; each writes its own output
for L in 0.02 0.03 0.04 0.06; do
  julia sweep.jl $L &          # sweep.jl reads ℓ from ARGS and sets a unique outdir
done
wait

Julia Distributed (one controller, N workers):

using Distributed
addprocs(8)
@everywhere include("PhaseFieldFracture.jl")
@everywhere using .PhaseFieldFracture

ℓs = [0.02, 0.03, 0.04, 0.06]
pmap(ℓs) do ℓ
    mesh  = rectangle_mesh(100, 100)
    crack = crack_line_nodes(mesh)
    mat   = Material(E = 210_000.0, ν = 0.3, Gc = 2.7, ℓ = ℓ, k = 1e-8)
    set   = Settings(nsteps = 50, umax = 8e-3)
    out   = "output_l$(ℓ)"
    run_simulation(Problem(mesh, mat, set, crack, out, joinpath(out, "hist.txt")))
end

On this 96‑core machine that is the highest‑value parallelism available — near linear speedup with zero risk.

10.3 Parallelizing a single simulation (how to add it)

Two hotspots dominate a single run: matrix assembly and the linear solve.

(a) Thread the element assembly loops. The COO tangent assembly in assemble_mechanics is already structured for threading: element e writes only to its own disjoint block [64(e−1)+1 : 64e] of the (I, J, V) arrays, so there is no write conflict. The change is mechanical — index the triplet block by e instead of by a shared running counter, then annotate the loop:

function assemble_tangent_threaded(mesh, mat, φ, u)
    ndof, ncell = 2*nnodes(mesh), ncells(mesh)
    I = Vector{Int}(undef, 64*ncell)
    J = Vector{Int}(undef, 64*ncell)
    V = Vector{Float64}(undef, 64*ncell)
    Threads.@threads for e in 1:ncell           # each e owns a disjoint 64-entry block
        nodes = @view mesh.cells[e, :]
        X     = mesh.coords[nodes, :]
        φe    = @view φ[nodes]
        dofs  = cell_dofs(nodes)
        ue    = u[dofs]
        Kₑ = zeros(8, 8)
        for (ξ, η) in GAUSS
            N, ∇N, detJ = element(X, ξ, η)
            B = strain_matrix(∇N)
            g = degrade(mat, dot(N, φe))
            _, C = constitutive(mat, B * ue, g)         # degraded tangent from the split
            Kₑ .+= (B' * C * B) .* detJ
        end
        base = 64*(e - 1)
        for a in 1:8, b in 1:8
            base += 1
            I[base] = dofs[a]; J[base] = dofs[b]; V[base] = Kₑ[a, b]
        end
    end
    return sparse(I, J, V, ndof, ndof)
end

Launch with julia -t 8 run.jl (or -t auto). The caveat: the internal force f_int[dofs[a]] += fₑ[a] in assemble_mechanics (and the phase‑field source f[nodes[a]] += fₑ[a] in phasefield_system) is a genuine data race under threads because neighbouring cells share nodes. Fix it with a per‑thread buffer reduced at the end, atomics, or by scattering the vector through the same triplet mechanism as K. The stiffness triplets themselves are race‑free.

(b) Speed up / parallelize the solve. Both matrices are symmetric positive definite, but generic \ uses an LU. Switching to a Cholesky is a large, free win, and reusing the symbolic factorization avoids re‑analyzing the (fixed) sparsity pattern each iteration:

using SparseArrays, LinearAlgebra
F = cholesky(Symmetric(K_ff))     # CHOLMOD; ~2× faster than LU for SPD
x_free = F \ rhs

For very large meshes, replace the direct solve with a threaded iterative solver (preconditioned conjugate gradients + algebraic multigrid) — CG’s sparse mat‑vec products parallelize well and scale to problems a direct factor cannot hold.

10.4 Summary

Level	Effort	Payoff	Notes
Threaded BLAS in the solve	set one env var	modest	already available
Parameter sweeps (multi‑run)	none / shell / pmap	~linear	best value now
Threaded assembly (@threads)	small edit	good on many cores	mind the f‑vector race
SPD cholesky + cached symbolic	small edit	large per‑solve	matrices are SPD
Iterative CG + AMG solver	larger edit	enables huge meshes	for scaling out

---

11. Performance & scaling

The cost is dominated by the per‑iteration re‑assembly and re‑factorization of the two sparse systems. Each staggered iteration solves the mechanics by Newton (each Newton step reassembles the tangent K and residual f_int because the split makes them depend on u, and factorizes K), then reassembles and factorizes the phase‑field K (because H changed). Newton typically needs only 2–5 steps (the split’s stress is piecewise linear, so the "active set" of tensile/compressive directions settles quickly), but it does multiply the mechanical solve cost by that factor versus the old single linear solve. This is where the time and memory go.

• Assembly is O(ncells) and highly allocation‑heavy as written (fresh small matrices per Gauss point, for clarity over speed).
• Direct solve of a 2‑D problem with N dofs is roughly O(N^{1.5}) time and O(N log N) memory with a good fill‑reducing ordering (SuiteSparse provides one automatically).
• Total work ≈ nsteps × (avg staggered iters) × (assemble + factor + solve). The staggered iteration count spikes near the critical load (it hit the max_iter = 50 cap at the peak in the smoke run), so most of the time is spent around crack initiation.

Measured data point (coarse smoke test, 30×30 mesh = 961 nodes, 8 steps, on this machine): ≈ 20 s wall including ~30 % first‑call compilation, ~8 GiB total allocations. The shipped 100×100 / 50‑step run is substantially heavier (≈ 10× the cells, ≈ 6× the steps, and each solve is larger) — expect minutes, not seconds. Practical speedups, in order of value: (1) run sweeps in parallel (§10.2); (2) switch the solve to cholesky and cache the symbolic factorization (§10.3b); (3) thread the assembly (§10.3a); (4) reduce per‑Gauss‑point allocations (pre‑allocate B, Kₑ, reuse buffers).

---

12. Critical review: assumptions & limitations

This is a clean, correct reference implementation — its priorities are readability and self‑containment, not production performance or full physical generality. The important modeling choices to be aware of:

1. Tension–compression split — now implemented (Miehe spectral). Only the tensile energy ψ⁺ drives and is degraded, so the model no longer damages under pure compression and broken crack faces transmit compressive contact (strainsplit.jl, §3.3). Two caveats remain: (a) the spectral split is used; the alternative Amor et al. (2009) volumetric–deviatoric split is not provided (it is cheaper but less accurate in shear — easy to add as a second constitutive branch); and (b) no split fully suppresses shear‑driven crack growth under confined compression, so very high compressive/shear loads can still nucleate some tip damage — this is physical and a known limitation of the basic split, not a bug. Set split = :none to recover the original no‑split model for comparison.

2. AT2 has no elastic threshold. With the φ²/(2ℓ) term, damage begins to grow at any nonzero load (note max_free_phi ≈ 0.43 already at the tiny first step), so there is no truly linear‑elastic regime and the pre‑peak stiffness is slightly reduced. The AT1 model (a φ/… linear term with a [0,1] bound constraint; Pham/Marigo, Tanné et al.) introduces a genuine elastic phase and a sharper crack profile, at the cost of needing a bound‑constrained solve.

3. Staggered convergence is slow near instability. Alternate minimization is only linearly convergent and can stall exactly when the crack becomes unstable (seen as converged = false at the peak in the smoke run). It is very robust but can need many iterations. Alternatives: a monolithic Newton solve (faster but less robust, needs globalization), over‑relaxation / Anderson acceleration of the staggered map, or adaptive load stepping that shrinks Δu near the peak.

4. Irreversibility is enforced by projection, not a variational inequality. φ ← clamp(max(φ, φ_prev), 0, 1) plus the monotone history H is adequate and common for AT2, but it is not a rigorous KKT/bound‑constrained enforcement. For AT1 (or strict irreversibility) a proper bound‑constrained solver is preferable.

5. Direct solver, re‑assembled and re‑factorized every iteration. Fine and robust for these mesh sizes, but the single biggest performance lever (§§10–11). The matrices are SPD, so cholesky and a cached symbolic factorization are easy wins; iterative solvers are needed to scale to large 3‑D problems.

6. Fixed linear load ramp. No arc‑length / adaptive control, so snap‑back branches of the force–displacement curve cannot be traced and the steps at the instability are the hardest to converge.

7. Scope. 2‑D, plane strain, small strain, single isotropic material, Q4 elements, structured rectangular meshes, displacement‑controlled Dirichlet loading only. No body forces, tractions, thermal/multiphysics coupling, unstructured meshes, or higher‑order elements. These are deliberate simplifications, not oversights.

8. ASCII VTU output. Human‑readable and dependency‑free, but the files are large and slower to write/read than binary/appended‑base64 VTK for big meshes.

9. No automated tests and no Project.toml. For reproducibility you may want to pin the Julia version and add a small verification test (e.g. a patch‑test / known peak‑load check).

None of these are correctness bugs in what the code claims to do — the FEM kernels, assembly, BCs, and staggered logic were exercised end‑to‑end and produce the expected brittle SENT response. They are the boundaries of the model’s applicability and the natural roadmap for extension.

---

13. Extending the code

Because each file is one ingredient, most extensions are local:

Goal	Where to work
Volumetric–deviatoric (Amor) split	add a branch to constitutive (elasticity.jl) and driving_energy; return σ±, C± as in strainsplit.jl
Anisotropic / other degradation g(φ)	degrade (elasticity.jl) — feeds constitutive and phasefield_system
AT1 model	the source/reaction terms in phasefield_system; add a bound‑constrained solve in solve_phase
Different BCs / loading	displacement_bcs (elasticity.jl) and the ramp in run_simulation (solver.jl)
Body forces / tractions	add to the RHS built in solve_displacement
New geometry / notch	rectangle_mesh, crack_line_nodes (mesh.jl)
Faster / parallel solve	solve_dirichlet (solver.jl) — swap in cholesky, cache symbolic, thread assembly
Binary VTK output	write_vtu (output.jl)

---

14. Troubleshooting / FAQ

• “No crack nodes found — use an even ny …” crack_line_nodes needs a row of mesh nodes on y = crack_y. With the default domain that means an even ny.
• “Non‑positive Jacobian determinant …” A cell is degenerate or its nodes are mis‑ordered. rectangle_mesh always produces valid CCW cells, so this only appears if you build a custom mesh — check the counter‑clockwise node ordering.
• “plane strain requires 0 ≤ ν < 0.5” D is singular at ν = ½ (incompressible). Use ν < 0.5.
• A step prints “did not converge in N iterations.” Expected near the peak load; the run continues using the last iterate. Reduce Δu (raise nsteps) or raise max_iter for a cleaner curve.
• First run feels slow. That is Julia compiling. The second run in the same session is much faster; for repeated use, keep a session open or build a sysimage.
• ParaView shows nothing animating. Open the .pvd, not an individual .vtu; the .pvd is what carries the time series.

---

15. References

• G. A. Francfort, J.-J. Marigo (1998). Revisiting brittle fracture as an energy minimization problem. J. Mech. Phys. Solids 46(8), 1319–1342.
• B. Bourdin, G. A. Francfort, J.-J. Marigo (2000/2008). Numerical experiments in revisited brittle fracture / The variational approach to fracture.
• L. Ambrosio, V. M. Tortorelli (1990). Approximation of functionals depending on jumps by elliptic functionals via Γ‑convergence. (the AT1/AT2 regularizations)
• C. Miehe, M. Hofacker, F. Welschinger (2010). A phase field model for rate‑independent crack propagation: Robust algorithmic implementation based on operator splits. Comput. Methods Appl. Mech. Engrg. 199, 2765–2778. (the spectral split implemented here; the staggered scheme; the SENT benchmark)
• C. Miehe, F. Welschinger, M. Hofacker (2010). Thermodynamically consistent phase‑field models of fracture: Variational principles and multi‑field FE implementations. Int. J. Numer. Methods Engng 83, 1273–1311. (companion paper; the ψ⁺/ψ⁻ energy split used in strainsplit.jl)
• C. Miehe (1998). Comparison of two algorithms for the computation of fourth‑order isotropic tensor functions. Comput. Struct. 66, 37–43. (the spectral derivative giving the positive projection ℙ⁺ = consistent tangent)
• H. Amor, J.-J. Marigo, C. Maurini (2009). Regularized formulation of the variational brittle fracture with unilateral contact … JMPS 57(8). (volumetric–deviatoric split)
• K. Pham, H. Amor, J.-J. Marigo, C. Maurini (2011). Gradient damage models and their use to approximate brittle fracture. (AT1, elastic threshold)
• E. Tanné, T. Li, B. Bourdin, J.-J. Marigo, C. Maurini (2018). Crack nucleation in variational phase‑field models of brittle fracture. JMPS 110.

---

16. License & citation

© 2026 Materials Mechanics Laboratory (MMLab). All rights reserved.

Author: Yang Bai — yangbai90@outlook.com. If you use this code in academic work, please cite it as the PhaseFieldFracture pure‑Julia solver by Yang Bai (MMLab, 2026), and cite the primary references above for the underlying model.