Coupled step: an outer iteration that converges, and a coupled step across the Ampère gate

examples/README.md

@@ -16,16 +16,29 @@ of `Project.toml`, which needs Julia 1.11 or later). It also adds `Plots`, which
 the plotting extension; the mp4 uses the FFMPEG that ships with it. Outputs go to
 `examples/output/<name>/`.
 
+### Verdicts
+
+A script with a prediction checks itself against it. It saves its figure as
+`<name>__PASS.png` or `<name>__FAIL.png`, removes the figure of the other verdict, and heads
+the figure with what it checked. The output folder therefore shows the verdict of the last
+run. The scenario scripts have no prediction and save their figures without a verdict.
+
+Known failures, and the work each belongs to:
+
+| figure | what fails | cause |
+|---|---|---|
+| `force_balance_control/position_control__FAIL` | the column is not held; it reaches the wall | the force and velocity model: the global J×B force is an explicit kick on an accumulated velocity. The coupled solve is not the cause. |
+
 ## Scenarios
 
-| script | field and geometry | physics on | simulated time |
-|---|---|---|---|
-| `townsend_avalanche.jl` | single-quadrupole null, box wall | atomic reactions, transport | 0.8 ms |
-| `selfE_avalanche.jl` | single-quadrupole null, box wall | as `townsend_avalanche.jl`, plus E∥ cancellation, mean E×B, turbulent E×B mixing | 4 ms |
-| `current_diffusion.jl` | pure toroidal field | a current filament with Ampère, with and without the inductive E; single-filament L/R reference with the electrons' kinetic inductance | 2 × 20 ms |
-| `force_balance_control.jl` | pure toroidal field | J×B hoop force; curved vertical-field PID position control | 2 × 2 ms |
-| `full_startup.jl` | single-quadrupole null, box wall | every module except the global J×B force | 10 ms |
-| `kstar_reference.jl` | KSTAR, time-varying external field | self-E model, Ampère off | 40 ms |
+| script | field and geometry | physics on | simulated time | check |
+|---|---|---|---|---|
+| `townsend_avalanche.jl` | single-quadrupole null, box wall | atomic reactions, transport | 0.8 ms | scenario, no verdict |
+| `selfE_avalanche.jl` | single-quadrupole null, box wall | as `townsend_avalanche.jl`, plus E∥ cancellation, mean E×B, turbulent E×B mixing | 4 ms | scenario, no verdict |
+| `current_diffusion.jl` | pure toroidal field | a current filament with Ampère, with and without the inductive E; single-filament L/R reference with the electrons' kinetic inductance | 2 × 20 ms | the current follows the (L + L_kin)/R circuit within 1 % of saturation |
+| `force_balance_control.jl` | pure toroidal field | J×B hoop force; curved vertical-field PID position control | 2 × 2 ms | the controller holds the current centroid within 5 cm of 1.5 m after 1 ms |
+| `full_startup.jl` | single-quadrupole null, box wall | every module except the global J×B force | 10 ms | scenario, no verdict |
+| `kstar_reference.jl` | KSTAR, time-varying external field | self-E model, Ampère off | 40 ms | scenario, no verdict |
 
 ## Coupled-step verification (`coupled_step/`)
 
@@ -48,6 +61,29 @@ flux, and that flux is computed from the simulated current.
 | `column_pushed_toward_loop.jl` | the column pushed at 200 m/s toward a superconducting loop outside the wall | I_c = −Φ_p/L_c; the loop pushes the column back |
 | `column_shifted_in_shell.jl` | the column moved up one cell inside a shell of 24 superconducting filaments | the shell's flux-conserving currents push it back down |
 
+### The coupled solve's iteration
+
+Within a step the coupled solve iterates on the boundary flux and the coil currents (see
+`coupled_step/common.jl`). These scripts compare the default solve with the same equations
+iterated to convergence by the relaxed iteration, step by step, each step against the converged
+current of that step (floored at 1e-3 of the run's peak). A case whose converged run did not
+converge fails. Each figure shows four
+things:
+- where the column and the conductors sit;
+- the plasma current of the two runs;
+- the first step's induced-field error against the number of iterations, for the relaxed
+  iteration (the default before Anderson mixing) and for the default;
+- that error over the grid.
+
+| script | setup | check |
+|---|---|---|
+| `picard_center_column.jl` | a dense, hot column in the middle of the default domain | the default stays within 1 % of the converged current at every step |
+| `picard_kstar_inboard_limited.jl` | the KSTAR grid and first wall; a dense column 4 cm from the inboard wall | as above |
+| `picard_tight_box.jl` | a column filling a box wall one cell inside the grid | as above |
+| `picard_filament_shell.jl` | a dense column inside a shell of 24 copper filaments | as above |
+| `picard_regime_map.jl` | the KSTAR grid and first wall; inboard-limited columns of radius 0.25–0.45 m, from 1e16 m⁻³ and 2 eV to 1e19 m⁻³ and 20 eV | every case within 1 % |
+| `gate_crossing_dense_column.jl` | dense columns at rest under the default 1 A gate | the run steps as the run with the gate at 0, within 1 % |
+
 Each script writes one figure to `examples/output/coupled_step/`:
 
 ```

examples/common.jl

@@ -198,3 +198,33 @@ function animate2D(runs::Vector{<:Pair}, fields::Vector{Symbol}; file, fps = 10)
     return file
 end
 animate2D(RP::RAPID, fields::Vector{Symbol}; kw...) = animate2D(["" => RP], fields; kw...)
+
+# ── verdicts ───────────────────────────────────────────────────────────────────────
+# An example that checks itself against its prediction saves its figure as <name>__PASS.png
+# or <name>__FAIL.png, and removes the figure of the other verdict (and an unlabelled one), so
+# the output folder shows the last run's verdict. The verdict and `detail` also head the
+# figure and are printed.
+function save_with_verdict(fig, dir::AbstractString, name::AbstractString, pass::Bool, detail::AbstractString)
+    verdict = pass ? "PASS" : "FAIL"
+    for old in ("$name.png", "$(name)__PASS.png", "$(name)__FAIL.png")
+        isfile(joinpath(dir, old)) && rm(joinpath(dir, old))
+    end
+    # the heading wrapped to the figure's width (about 9 px per character at this size)
+    width = fig[:size][1]
+    words, lines = split("$verdict: $detail"), [""]
+    for w in words
+        if isempty(lines[end]) || length(lines[end]) + length(w) + 1 <= width / 9
+            lines[end] = isempty(lines[end]) ? String(w) : lines[end] * " " * w
+        else
+            push!(lines, String(w))
+        end
+    end
+    plot!(
+        fig; plot_title = join(lines, "\n"), plot_titlefontsize = 11, plot_titlefontcolor = pass ? :darkgreen : :red3,
+        plot_titlevspan = min(0.05 + 0.03 * (length(lines) - 1), 0.2),
+    )
+    file = joinpath(dir, "$(name)__$(verdict).png")
+    savefig(fig, file)
+    println(verdict, "  ", name, ": ", detail)
+    return file
+end

examples/coupled_step/coil_driven_column.jl

@@ -63,5 +63,12 @@ fig = plot(
     plot_layout(RP; title = "column and driving coil"), p1, p2;
     layout = @layout([a{0.38w} grid(2, 1)]), size = (1200, 620), margin = 4Plots.mm,
 )
-savefig(fig, joinpath(out, "coil_driven_column.png"))
+# The column follows the two-circuit model, and the run under the default gate steps as the
+# run without it.
+model_gap = maximum(abs.(open.Ip .- last.(model))) / maximum(abs, last.(model))
+gate_gap = maximum(abs.(default.Ip .- open.Ip)) / maximum(abs, open.Ip)
+save_with_verdict(
+    fig, out, "coil_driven_column", model_gap < 0.05 && gate_gap < 0.02,
+    @sprintf("threshold 0 strays up to %.1f %% from the two-circuit model (passes under 5 %%); the 1 A gate run strays up to %.1f %% from threshold 0 (passes under 2 %%)", 100model_gap, 100gate_gap),
+)
 println("outputs in ", out)

examples/coupled_step/column_pushed_toward_loop.jl

@@ -73,5 +73,9 @@ fig = plot(
     l1, l2, p1, p2, p3;
     layout = @layout([grid(2, 1){0.34w} grid(3, 1)]), size = (1200, 860), margin = 4Plots.mm,
 )
-savefig(fig, joinpath(out, "column_pushed_toward_loop.png"))
+gap = maximum(abs.(rec.Ic .- rec.Ic_flux)) / maximum(abs, rec.Ic_flux)
+save_with_verdict(
+    fig, out, "column_pushed_toward_loop", rec.Rc[end] > 1.6 && gap < 0.02,
+    @sprintf("the loop current follows −Φ_p/L_c within %.2g %% (needs 2 %%) while the column moves to R = %.2f m", 100gap, rec.Rc[end]),
+)
 println("outputs in ", out)

examples/coupled_step/column_shifted_in_shell.jl

@@ -5,7 +5,9 @@
 #   ΔI = −M⁻¹ (Φ(J_shifted) − Φ(J_before)),
 # M the filaments' inductance matrix and Φ the plasma flux through each, and those currents
 # push the column back down: F_Z = −∫ Jϕ B_R dV < 0. This restoring force is how a
-# conducting wall holds a vertical displacement until its currents decay.
+# conducting wall holds a vertical displacement until its currents decay. The column itself
+# answers the shell's currents inductively at first, so the force builds up to that value over
+# about 0.1 ms; the coupled solve iterated to convergence gives the same build-up.
 #
 #   julia --project=examples examples/coupled_step/column_shifted_in_shell.jl
 
@@ -80,5 +82,9 @@ fig = plot(
     plot_layout(RP; title = "column and shell filaments"), p1, p2;
     layout = @layout([a{0.38w} grid(2, 1)]), size = (1200, 680), margin = 4Plots.mm,
 )
-savefig(fig, joinpath(out, "column_shifted_in_shell.png"))
+gap = abs(rec.F[end] - pred.F[]) / abs(pred.F[])
+save_with_verdict(
+    fig, out, "column_shifted_in_shell", pred.F[] < 0 && rec.F[1] < 0 && gap < 0.02,
+    @sprintf("the shell pushes the column back from the first step; %.1f ms after the shift its force is within %.2g %% of the flux-conserving shell's (needs 2 %%)", rec.t[end] * 1.0e3, 100gap),
+)
 println("outputs in ", out)

examples/coupled_step/common.jl

@@ -16,9 +16,10 @@ include(joinpath(@__DIR__, "..", "common.jl"))
 function column(
         name; E0 = 0.3, Te = 1.0, n0 = 1.0e16, cenR = 1.5, cenZ = 0.0, radius = 0.3,
         threshold = 0.0, dt = 5.0e-6, t_end = 1.0e-3, moving = false,
+        manual = pure_toroidal(E0), wall_R = Float64[], wall_Z = Float64[],
     )
     config = SimulationConfig{Float64}(;
-        device_Name = "manual", manual = pure_toroidal(E0), NR = 30, NZ = 50, R0B0 = 3.0,
+        device_Name = "manual", manual, wall_R, wall_Z, NR = 30, NZ = 50, R0B0 = 3.0,
         prefilled_gas_pressure = 0.0,          # vacuum: no neutrals
         dt, t_end_s = t_end, snap0D_Δt_s = 10dt, snap2D_Δt_s = t_end,
         Output_path = output_dir(name),
@@ -119,11 +120,153 @@ function plot_layout(RP; J = current_density(RP), title = "layout")
     p = heatmap(
         G.R1D, G.Z1D, permutedims(Jmax > 0 ? J ./ Jmax : J);
         c = :balance, clims = (-1, 1), aspect_ratio = :equal, colorbar_title = "J / max|J|",
-        xlims = (min(G.R1D[1], minimum(rc) - 0.1), max(G.R1D[end], maximum(rc) + 0.1)),
+        xlims = (min(G.R1D[1], minimum(rc; init = Inf) - 0.1), max(G.R1D[end], maximum(rc; init = -Inf) + 0.1)),
         ylims = extrema(G.Z1D), xlabel = "R (m)", ylabel = "Z (m)", title, titlefontsize = 10,
         framestyle = :box,
     )
     plot!(p, vcat(RP.wall.R, RP.wall.R[1]), vcat(RP.wall.Z, RP.wall.Z[1]); c = :gray40, lw = 1, label = "wall")
-    scatter!(p, rc, zc; c = :orange, ms = 4, msw = 0, label = "loops")
+    isempty(rc) || scatter!(p, rc, zc; c = :orange, ms = 4, msw = 0, label = "loops")
     return p
 end
+
+# ── the coupled solve's outer iteration ────────────────────────────────────────────────
+# Within a step the coupled solve iterates on the boundary flux and the coil currents: solve
+# u∥ and ψ inside the domain with the boundary flux held, then recompute the boundary flux
+# (Green's functions) and the coil currents (circuits) from the new plasma current, and repeat.
+# The default mixes the iterates by Anderson (memory 8, at most 20 block solves); the relaxed
+# iteration (anderson_m = 0, boundary flux weighted by w = 0.5) was the default before.
+# CONVERGED_PICARD iterates the same equations to convergence with the relaxed iteration
+# (w = 0.05), to 1e-10 of the step's field and of each coil's change with no floors: what the
+# default should reproduce. A run whose reference did not converge has no expected result.
+const CONVERGED_PICARD = (
+    tolerance = 1.0e-10, max_iter = 20_000, relaxation_w = 0.05, anderson_m = 0, E_floor = 0.0, I_floor = 0.0,
+)
+
+# The largest gap between a plasma current history and the expected one, step by step against
+# the expected current of that step, floored at 1e-3 of its peak where it passes through zero.
+# Exact agreement is no gap, also where the expected current is zero throughout.
+function current_gap(I, I_ref)
+    scale = max.(abs.(I_ref), 1.0e-3 * maximum(abs, I_ref))
+    return maximum(ifelse(e == 0, 0.0, e / s) for (e, s) in zip(abs.(I .- I_ref), scale))
+end
+
+# Geometries. A KSTAR-like domain: the grid of the KSTAR field files (R 1.2–2.4 m, Z ±1.2 m)
+# with the KSTAR first wall (KSTAR_First_Wall.dat), whose inboard side is 6 cm (1.5 cells at
+# 30×50) inside the grid. A tight box: the default box wall one cell inside a 1.2 × 1.6 m grid.
+const KSTAR_WALL_R = [1.26, 1.632, 1.992, 2.256, 2.256, 1.992, 1.632, 1.26, 1.26]   # closed
+const KSTAR_WALL_Z = [1.13, 1.056, 0.732, 0.456, -0.456, -0.732, -1.056, -1.13, 1.13]
+kstar_like(E0) = ManualSetup{Float64}(R = (1.2, 2.4), Z = (-1.2, 1.2), BR = 0.0, BZ = 0.0, Eϕ = E0)
+tight_box(E0) = ManualSetup{Float64}(R = (1.0, 2.2), Z = (-0.8, 0.8), BR = 0.0, BZ = 0.0, Eϕ = E0, wall_margin_cells = 1)
+
+# A shell of `nfil` copper filaments on a circle of radius `r` around (cenR, 0), each of the
+# square cross-section that tiles the circle: a passive conducting structure inside the grid.
+function filament_shell!(RP; cenR, r, nfil = 24)
+    side = 2π * r / nfil
+    for k in 1:nfil
+        θ = 2π * (k - 0.5) / nfil + 0.05
+        R, Z = cenR + r * cos(θ), r * sin(θ)
+        add_loop!(RP, R, Z; a = side / sqrt(π), R = 1.68e-8 * 2π * R / side^2, name = "shell_$k")
+    end
+    initialize_coil_system!(RP)
+    return RP
+end
+
+quiet(f) = redirect_stdout(() -> redirect_stderr(f, devnull), devnull)
+
+# The column of `make()` run for `nsteps` steps twice, with the default Picard and with
+# CONVERGED_PICARD: the plasma current and the induced field after each step, and the Picard
+# counters.
+function default_vs_converged(make; nsteps)
+    return map((nothing, CONVERGED_PICARD)) do picard
+        RP = make()
+        isnothing(picard) || (RP.flags.ampere_picard = PicardSettings{Float64}(; picard...))
+        RP.t_end_s = nsteps * RP.dt
+        I, E = Float64[], Matrix{Float64}[]
+        record(rp) = (push!(I, plasma_current(rp, current_density(rp))); push!(E, copy(rp.fields.Eϕ_self)))
+        quiet(() -> run_simulation!(RP; callback_after_step = record))
+        (; RP, I, E, stats = deepcopy(RP.diagnostics.ampere_picard))
+    end
+end
+
+# The first step of the column of `make()`, solved from the same state and stopped after
+# L = 1…Lmax block solves, against the converged step: the error of the induced field,
+# max|Eϕ_L − Eϕ*| / max|Eϕ*|, for the relaxed iteration (anderson_m = 0, w = 0.5) and for the
+# default solve. Both keep the solve's failure policy: a residual that grows a thousandfold restarts
+# from the best iterate with half the mixing, so the relaxed iteration no longer runs away.
+function picard_error_by_iteration(make; Lmax = 30)
+    RP = make()
+    pla, F, csys = RP.plasma, RP.fields, RP.coil_system
+    # what run_simulation! does before its first step
+    pla.ne[RP.G.nodes.on_out_wall_nids] .= 0.0
+    pla.ni[RP.G.nodes.on_out_wall_nids] .= 0.0
+    initialize_coupled_fields!(RP)
+    RAPID2D.update_transport_quantities!(RP)
+    prepare_timestep!(RP)
+    saved = (
+        u = copy(pla.ue_para), ψ = copy(F.ψ_self), E = copy(F.Eϕ_self), Ep = copy(F.Eϕ_self_prev),
+        I = csys.n_total > 0 ? copy(get_all_currents(csys)) : Float64[],
+        Φ = csys.n_total > 0 ? copy(csys.coils.ψ_pla) : Float64[], t = csys.time_s,
+    )
+    function trial(; kw...)
+        pla.ue_para .= saved.u; F.ψ_self .= saved.ψ; F.Eϕ_self .= saved.E; F.Eϕ_self_prev .= saved.Ep
+        if csys.n_total > 0
+            set_all_currents!(csys, copy(saved.I)); csys.coils.ψ_pla = copy(saved.Φ); csys.time_s = saved.t
+        end
+        quiet(() -> RAPID2D.solve_combined_momentum_Ampere_equations_with_coils!(RP; kw...))
+        return copy(F.Eϕ_self)
+    end
+    E_star = trial(; CONVERGED_PICARD...)
+    err(m, L) = maximum(abs, trial(; tolerance = 0.0, max_iter = L, anderson_m = m, E_floor = 0.0, I_floor = 0.0) .- E_star)
+    scale = maximum(abs, E_star)
+    return (relaxed = [err(0, L) for L in 1:Lmax] ./ scale, default = [err(RP.flags.ampere_picard.anderson_m, L) for L in 1:Lmax] ./ scale)
+end
+
+# One figure for a case of the coupled solve: where the column and the conductors sit, the
+# plasma current step by step with the default and the converged solve, the first step's
+# induced-field error against the number of block solves (relaxed iteration and default), and
+# that error over the grid after the default solve. Passes when the default stays within 1 % of
+# the converged current at every step (`current_gap`).
+function picard_case(make, name, title; nsteps = 20)
+    def, conv = default_vs_converged(make; nsteps)
+    errs = picard_error_by_iteration(make)
+    t = (1:nsteps) .* def.RP.dt .* 1.0e6
+    gap = current_gap(def.I, conv.I)
+    pass = conv.stats.nunconverged == 0 && isfinite(gap) && gap <= 1.0e-2
+
+    blowup = maximum(abs, def.I) > 100 * maximum(abs, conv.I)
+    p1 = plot!(plot_layout(conv.RP; title); colorbar = false, titlefontsize = 9)
+    p2 = plot(
+        t, blowup ? abs.(conv.I) : conv.I; c = :black, lw = 3, label = "converged (expected)",
+        xlabel = "t (µs)", ylabel = blowup ? "|I_p| (A)" : "I_p (A)", yscale = blowup ? :log10 : :identity,
+        title = "plasma current", legend = :topleft,
+    )
+    plot!(p2, t, blowup ? max.(abs.(def.I), 1.0e-3) : def.I; c = :red3, ls = :dash, lw = 2, m = :circle, ms = 3, label = "default solve")
+    p3 = plot(
+        1:length(errs.relaxed), max.(errs.relaxed, 1.0e-16); yscale = :log10, c = :gray50, lw = 2, m = :circle, ms = 2,
+        label = "relaxed iteration (w = 0.5, halved when it grows)", xlabel = "block solves in the first step", ylabel = "max|ΔEϕ| / max|Eϕ*|",
+        title = "first step: error by iteration", ylims = (1.0e-12, 1.0e6),
+    )
+    plot!(p3, 1:length(errs.default), max.(errs.default, 1.0e-16); c = :red3, lw = 2, m = :circle, ms = 3, label = "default solve")
+    vline!(p3, [def.RP.flags.ampere_picard.max_iter]; c = :gray, ls = :dash, label = "default limit")
+    hline!(p3, [1.0e-3]; c = :green, ls = :dot, label = "tolerance (1e-3)")
+    G = conv.RP.G
+    ΔE = (def.E[1] .- conv.E[1]) ./ maximum(abs, conv.E[1])
+    lim = max(maximum(abs, filter(isfinite, ΔE)), 1.0e-12)
+    p4 = heatmap(
+        G.R1D, G.Z1D, permutedims(ΔE); c = :balance, clims = (-lim, lim), aspect_ratio = :equal,
+        xlabel = "R (m)", ylabel = "Z (m)", title = "step 1: (Eϕ − Eϕ*) / max|Eϕ*|", titlefontsize = 10,
+        framestyle = :box, xlims = extrema(G.R1D), ylims = extrema(G.Z1D),
+    )
+    plot!(p4, vcat(conv.RP.wall.R, conv.RP.wall.R[1]), vcat(conv.RP.wall.Z, conv.RP.wall.Z[1]); c = :gray40, lw = 1, label = "")
+    fig = plot(
+        p1, p2, p3, p4; layout = (1, 4), size = (1800, 540), margin = 5Plots.mm, left_margin = 10Plots.mm,
+        top_margin = 10Plots.mm, bottom_margin = 12Plots.mm,
+    )
+    detail = @sprintf(
+        "at every step the default solve is within %.2g %% of that step's converged current (floored at 1e-3 of the peak; passes under 1 %%); first step %.3g A vs %.3g A; %.1f block solves/step, %d unconverged (converged run: %.0f/step, %d unconverged)",
+        100gap, def.I[1], conv.I[1], def.stats.niter / def.stats.nsolve, def.stats.nunconverged,
+        conv.stats.niter / conv.stats.nsolve, conv.stats.nunconverged
+    )
+    save_with_verdict(fig, output_dir("coupled_step"), name, pass, detail)
+    return (; def, conv, errs, gap, pass)
+end