solver.jl

# ===========================================================================
#  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.
# ===========================================================================

# ===========================================================================
#  Staggered solver: monotonic loading + alternate minimization
# ===========================================================================
#
#  The coupled problem "find u and φ that minimize the energy" is nonlinear
#  because of the g(φ)ψ coupling.  We solve it the standard robust way:
#  hold one field fixed and solve the (now linear) problem for the other, and
#  alternate.  Per load step:
#
#        ┌─ solve displacement u   with φ fixed        (Newton; nonlinear once
#        │                                               the tension/compression
#        │                                               split is active)
#        │  update history H(u)                         (irreversible driving force)
#        └─ solve damage    φ   with H fixed            (linear AT2 phase field)
#
#  repeated until neither field changes, then the history is committed and the
#  imposed displacement is increased for the next step.

# ---------------------------------------------------------------------------
#  Settings and problem definition
# ---------------------------------------------------------------------------

"Numerical settings for the load stepping and the staggered iteration."
Base.@kwdef struct Settings
    nsteps::Int    = 50        # number of load steps
    umax::Float64  = 8e-3      # final imposed top displacement u_y (reached at the last step)
    max_iter::Int  = 50        # cap on staggered iterations within one load step
    tol::Float64   = 1e-6      # convergence tolerance on the relative change of u and φ
    out_every::Int = 5         # save a .vtu snapshot every `out_every` steps
end

"Everything needed to run one simulation: geometry, material, settings, notch, output paths."
struct Problem
    mesh::Mesh                    # the finite element mesh
    material::Material            # elastic + fracture properties
    settings::Settings           # solver controls
    crack_nodes::Vector{Int}     # nodes of the prescribed initial crack (held at φ = 1)
    outdir::String               # folder for the .vtu / .pvd files
    histfile::String             # path of the CSV load-history file
end

# ---------------------------------------------------------------------------
#  Small numerical helpers
# ---------------------------------------------------------------------------

"Relative change between two successive iterates — the staggered convergence measure."
rel_change(new, old) = norm(new - old) / max(norm(new), 1.0)

"""
    solve_dirichlet(K, f, bc) -> x

Solve the linear system `K x = f` subject to Dirichlet (prescribed-value)
constraints `bc`, a dictionary `dof => value`.

Splitting the dofs into *free* (unknown) and *fixed* (prescribed):

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

Here `K[free, fixed] * x[fixed]` moves the known columns to the right-hand side.
"""
function solve_dirichlet(K, f, bc::Dict{Int,Float64})
    ndof = length(f)
    isfixed = falses(ndof)                   # boolean mask: is this dof prescribed?
    for d in keys(bc)
        isfixed[d] = true
    end
    fixed = findall(isfixed)                 # indices of prescribed dofs
    free  = findall(.!isfixed)               # indices of free (unknown) dofs

    x = zeros(ndof)
    for (d, v) in bc                         # write the prescribed values into x
        x[d] = v
    end
    if !isempty(free)                        # solve the reduced system for the free dofs
        x[free] = K[free, free] \ (f[free] - K[free, fixed] * x[fixed])
    end
    return x
end

# ---------------------------------------------------------------------------
#  The main driver
# ---------------------------------------------------------------------------

"""
    run_simulation(problem) -> (u, φ, H)

Run the full staggered simulation.  Returns the final displacement `u`, damage
`φ`, and history field `H`.  Along the way it writes:

  * `problem.histfile`                    — CSV load history (one row per step),
  * `problem.outdir/field_XXXX.vtu`       — nodal fields at saved steps,
  * `problem.outdir/timeseries.pvd`       — ParaView collection over those steps.
"""
function run_simulation(problem::Problem)
    mesh, mat, set = problem.mesh, problem.material, problem.settings
    mkpath(problem.outdir)
    write_history_header(problem.histfile)

    # State vectors.  `u` holds the two displacement dofs per node; `φ` holds one
    # damage value per node.  The pre-crack nodes start fully broken.
    u = zeros(2 * nnodes(mesh))
    φ = zeros(nnodes(mesh))
    for n in problem.crack_nodes
        φ[n] = 1.0
    end

    # History field H: the maximum elastic energy density seen so far at every
    # Gauss point.  It is the irreversible driving force for damage.  `H` is the
    # committed value (carried between load steps); `H_trial` is a scratch copy
    # recomputed inside the staggered loop from the current displacement.
    H       = zeros(length(GAUSS), ncells(mesh))
    H_trial = similar(H)

    # Regularize the notch: with zero history, one phase-field solve spreads the
    # sharp φ = 1 crack line into the smooth AT2 damage profile of width ~ℓ.
    φ = solve_phase(mesh, mat, H, φ, problem.crack_nodes)

    # Snapshot 0 — the cracked body before any load is applied.  `pvd` accumulates
    # (time, filename) pairs for the ParaView collection written at the very end.
    f_int = zeros(2 * nnodes(mesh))                          # internal force (= 0 at u = 0)
    write_vtu(joinpath(problem.outdir, "field_0000.vtu"), mesh, u, φ)
    pvd = Tuple{Float64,String}[(0.0, "field_0000.vtu")]
    log_step(problem, 0, 0.0, u, φ, f_int, 0, true)

    for step in 1:set.nsteps
        top_uy      = set.umax * step / set.nsteps    # ramp the imposed displacement linearly
        φ_committed = copy(φ)                          # damage floor for this step (irreversibility)
        converged   = false
        iters       = set.max_iter

        # ---- staggered (alternate-minimization) loop ----
        for it in 1:set.max_iter
            u_prev, φ_prev = copy(u), copy(φ)          # previous iterates, for convergence test

            u, f_int = solve_displacement(mesh, mat, φ, top_uy, u)                 # (1) elasticity (Newton)
            update_history!(H_trial, H, mesh, mat, u)                              # (2) H ← max(H, ψ⁺)
            φ = solve_phase(mesh, mat, H_trial, φ_committed, problem.crack_nodes)  # (3) damage

            # Stop once neither displacement nor damage changes appreciably.
            if max(rel_change(u, u_prev), rel_change(φ, φ_prev)) < set.tol
                converged = true
                iters = it
                break
            end
        end
        converged || @printf("  (step %d did not converge in %d iterations)\n", step, set.max_iter)

        # Commit the history field once the two fields have settled at this load.
        update_history!(H_trial, H, mesh, mat, u)
        H .= H_trial

        # Save a snapshot on the chosen cadence (and always on the final step).
        if step % set.out_every == 0 || step == set.nsteps
            fname = @sprintf("field_%04d.vtu", step)
            write_vtu(joinpath(problem.outdir, fname), mesh, u, φ)
            push!(pvd, (top_uy, fname))                # use the imposed load as the "time"
        end
        log_step(problem, step, top_uy, u, φ, f_int, iters, converged)

        @printf("step %3d/%d | u_y = %.4e | max φ (off-crack) = %.4f | iters = %d\n",
                step, set.nsteps, top_uy, max_free_phi(φ, problem.crack_nodes), iters)
    end

    # Tie all snapshots into one ParaView time series.
    write_pvd(joinpath(problem.outdir, "timeseries.pvd"), pvd)

    return u, φ, H
end