Repository navigation
Coupled step: an outer iteration that converges, and a coupled step across the Ampère gate - #28
Conversation
solve_timestep! now hands RP.flags.ampere_picard (tolerance, max_iter, relaxation_w) to the coupled solve, so a run can choose them. The defaults are the solver's own, so a default run is unchanged.
…d solve An example with a prediction now checks itself. It saves its figure as <name>__PASS.png or <name>__FAIL.png (save_with_verdict) and heads it with what it checked. The coupled_step examples, current_diffusion and force_balance_control carry such checks. coil_driven_column also checks the whole trajectory of the 1 A gate run, not only its end. New coupled_step scripts compare the default coupled solve with the same equations iterated to convergence: - a dense column in the middle of the domain; - a column against the inboard wall of the KSTAR grid and first wall; - a column filling a box one cell inside the grid; - a column inside a shell of copper filaments; - a map over column size and density in the KSTAR-like domain. gate_crossing_dense_column runs dense columns at rest under the 1 A gate against the gate at 0. The README lists each check and the known failures with their causes.
calculate_psi_by_green_function computed dpsi/dR only through m, at fixed Rd*Rs. But psi ~ sqrt(Rd Rs/m) f(m), so each R-derivative also carries psi/(2R). Without it the R-derivatives were about 20 % low in the cases measured. The Z-derivatives were right. Nothing in the time step uses these tables. The uncalled displacement-term function (the A.1b candidate) does; its docstring no longer warns about them. A new test compares all four derivatives with central differences, on array-shaped destinations and sources with currents other than one.
A small solver-independent helper for the coupled solve's outer iteration. anderson_step!(A, x, f) takes an iterate and its residual f = g(x) - x and returns the next iterate, a type-II Anderson step with memory m. The step uses a diagonal mixing B and diagonal least-squares weights W. With m = 0 it is the relaxed iteration x + B f. On an affine map it works as GMRES, so it converges where the relaxed iteration diverges. The mixer keeps the iterate with the smallest weighted residual. A residual that is not finite, or more than 1e3 times that smallest one, restarts from it: the history is dropped and B is halved. With the restarts used up it stops at the best iterate. With no finite iterate at all it reports failure. Tests use a 20-unknown affine map whose relaxed iteration diverges, and cover m = 0, restarts, non-finite residuals, dependent differences and exhausted restarts.
The combined u-psi-circuit solve iterates on x = (boundary flux, coil currents). Its relaxed iteration failed in three ways: - it diverged for a large or dense column near the grid edge, where one eigenvalue of the map's linear part falls below -3; - it crawled with a shell of filaments inside the grid (eigenvalues up to 0.85); - it started from a poor first iterate (u from A_u alone), so ten iterations left the first steps far off: on the KSTAR grid and first wall, the first two steps came out at 63 A and -46 A instead of 5.5 A and 11 A. Now: - the first iterate is a block solve with the predicted boundary flux; the separate A_u factorization is gone; - Anderson mixing (AndersonMixer, memory anderson_m = 8) drives the residual f = g(x) - x to zero. anderson_m = 0 is the old relaxed iteration; - the circuits' forcing is taken once per step; - the stopping test is on that residual (coupled_residual_converged): the boundary flux and the in-grid coils' source in induced-field units, and each coil current against its own change; - the accepted evaluation gives u and psi, and the coil currents of its J with psi_pla(J) as the coils' memory, so the circuit balance closes on solves cut short too; - an unconverged solve accepts its best finite evaluation. With none finite, it throws. flags.ampere_picard gains anderson_m, and max_iter now defaults to 20. The uncalled alternative solver keeps its own iteration. Tests compare the default solve with the same equations iterated to convergence by the relaxed iteration (w = 0.05) for three cases: - a column filling a box one cell inside the grid; - a dense column inside a filament shell; - a column with u-advection and coils of 1 mH and microhenries. Further tests check that the coils remember the flux of the accepted current, that the circuit balance holds on solves cut short, and that the fixed point does not depend on the mixing. The coupled_step examples now plot the relaxed and the default iteration side by side.
Below the gate the plasma is not a source of induction. A step whose u-update carries the current across the gate therefore accelerated the electrons without their self-inductance, by about L_p/L_kin. A dense column at rest (pre-ionized, or formed at I = 0) jumped to hundreds of amperes, or kiloamperes, in one step under the default 1 A gate, instead of a few amperes, and stayed there. Now such a step is solved again by the coupled solve, from the state it started at, when that solve can run (Ampere, the inductive E and u-evolution on). The below-gate trial's writes are put back first: u, J, the induced field and its parallel projection, and the coils' currents, memory and clock. The rest it writes is either reassigned by the coupled solve before use or a cache rebuilt from the same state. One re-solve per step; its result stands even if it ends below the gate. A current decaying through the gate, or crossing it by density growth as in an avalanche, takes the same path as before. Regression tests check that a dense column at rest under the 1 A gate steps exactly (==) as with the gate at 0. They also check: - the same with a coil whose voltage changes in time, a resistive loop and stale circuit matrices: coil currents, memory and clock; - a crossing taken with advance_timestep!, the coils' memory unset; - a run split just before the crossing.
…ilure left
After the Anderson and gate commits, every coupled_step example passes,
including the Picard cases and the gate crossing. The README's
known-failures table keeps only force_balance_control, whose position
control fails for a force-model reason. The Picard figures plot the
relaxed iteration with its restart rule ("halved when it grows") next
to the default solve.
…s settings The outer iteration's control flow moves into anderson_solve!, which accepts a converged evaluation only if the mixer kept it (:best or :ok). The stopping test requires a finite residual and a finite step field (Inf <= Inf passed before), and the solver marks an evaluation invalid when u_par or psi is not finite where the residual does not see it. The QR cutoff of the Anderson least squares is max(1e-10, n eps) of the largest pivot: 1e-10 alone kept round-off pivots in single precision. flags.ampere_picard gains E_floor and I_floor. The reference solves in the tests run with both at zero to 1e-10 and must converge; coils are compared one at a time against their own peak. New tests: the driver on synthetic maps (rejected evaluation, best-evaluation fallback at max_iter and on exhaustion, failed first evaluation), the solver's fallback against the solve cut at max_iter = 1, and each part of the stopping test on its own.
…entry one The coils remember a flux the column at rest no longer makes. The trial below the gate moves the memory to the present current's flux, so a crossing that did not restore it would start the coupled solve from the wrong change. The gated step must still equal the step with the gate at 0. Dropping the restore fails this test and no other.
…verged reference current_gap compares each step's plasma current with the converged current of that step (floor 1e-3 of its peak), not with the run's peak, so a start-up transient cannot hide behind a later large current. The converged reference runs to 1e-10 with no floors, and a case whose reference did not converge fails. The regime map counts a gap that is not finite as a failure.
The stopping test returns whether the residuals and the step field are finite, and the evaluation is valid only if they are, as well as u_par and psi. Before, a finite but huge psi could overflow -dpsi/(R dt) and still be kept as the best evaluation and accepted on fallback. Such a first evaluation now makes the solve throw. New test: a vacuum in a purely poloidal field entering with psi = 1e305 inside (it accepted 80 infinite field entries before). Another test makes the run's coil floor decide convergence, which the floor test without coils could not. Docstrings call Anderson a truncated relative of GMRES rather than GMRES.
Exact agreement is no gap, also where the expected current is zero throughout (0/0 failed before); any other gap against a zero reference fails. Captions and the README state that each step is compared with that step's converged current, floored at 1e-3 of the peak.
flags.ampere_picard was a six-field NamedTuple: changing one field needed merge(...), and nothing checked the values until the mixer threw mid-run. It is now a mutable struct checked on construction and on assignment, as ImplicitWeights is, so RP.flags.ampere_picard.max_iter = 40 works and max_iter = 0 fails where it is written. solve_timestep! passes NamedTuple(settings)... to the coupled solve. Also: the PicardStats docstring describes residuals and both ways to stop short, and the coupled solve's comments lose their stale step numbers.
…p does The snapshot the gate crossing restores is a hand-written list. This guard compares every number in the plasma, fields, transport and coils after the re-solved crossing with the step taken with the gate at 0, so a write the trial below the gate leaves behind fails here, by path, wherever it lands. An added write to Ti in the trial and a dropped coil-memory restore both fail it.
The block was rebuilt every step as OP.II + spdiagm(...) + dt*sparse(A_adv) and sparse(I, J, V). That arithmetic drops stored zeros, and the upwind side a face does not use is a stored zero on the wall pattern, so the pattern moved whenever a flow turned and the LU redid its symbolic analysis (200 of 351 coupled steps in full_startup at 30x50, on master as well). A_u is now written in place on the wall pattern (set_identity!, add_diagonal!, add_scaled!), as the ne and u_par solves do, and CoupledBlock fixes the 2N block's pattern on the first coupled solve; update_coupled_block! writes values only. One symbolic analysis per run; the coupled solve is about a quarter cheaper at 30x50. The block equals the old assembly; combine_Au_and_DeltaGS_sparse_matrices stays, uncalled.
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #28 +/- ##
===========================================
- Coverage 93.69% 79.68% -14.02%
===========================================
Files 47 49 +2
Lines 5075 5207 +132
===========================================
- Hits 4755 4149 -606
- Misses 320 1058 +738 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
There was a problem hiding this comment.
Copilot review overview
🟡 Changes recommended
The failed-solve path leaves persistent state mutated, and the below-gate snapshot introduces significant per-step allocations.
Review effort: Balanced
Findings: 2
Open (4)
What changed in this PR
Improves the coupled momentum–Ampère–circuit step’s convergence, failure handling, and gate-crossing behavior.
Changes:
- Adds Anderson mixing, residual-based convergence, and best-evaluation fallback.
- Reuses a fixed sparse block pattern and re-solves gate-crossing steps.
- Expands verification examples and derivative/regression tests.
| File | Description |
|---|---|
src/RAPID2D.jl |
Loads new numerical modules. |
src/types.jl |
Adds Picard settings and block storage. |
src/initialization.jl |
Initializes the momentum block operator. |
src/workflows.jl |
Re-solves current-gate crossings. |
src/physics/physics.jl |
Implements the revised coupled solve. |
src/numerics/anderson.jl |
Adds Anderson mixing and fallback. |
src/numerics/coupled_block.jl |
Adds fixed-pattern block assembly. |
src/utils/green_function.jl |
Corrects radial derivatives. |
src/diagnostics/types.jl |
Updates solve statistics documentation. |
test/unit/numerics/anderson_test.jl |
Tests mixing and failure behavior. |
test/unit/physics/ampere_picard_test.jl |
Tests convergence and settings. |
test/unit/utils/green_function_test.jl |
Checks derivatives numerically. |
test/unit/coils/coil_plasma_flux_test.jl |
Checks circuit balance on fallback. |
test/regression/coupled_step_test.jl |
Tests gate-crossing state restoration. |
examples/common.jl |
Adds self-checking figure verdicts. |
examples/README.md |
Documents verdicts and new examples. |
examples/current_diffusion.jl |
Adds an analytic verdict. |
examples/force_balance_control.jl |
Adds a position-control verdict. |
examples/coupled_step/common.jl |
Adds convergence comparison utilities. |
examples/coupled_step/two_coils.jl |
Adds circuit-recursion verification. |
examples/coupled_step/density_growth_dt.jl |
Adds timestep-convergence verification. |
examples/coupled_step/density_doubling.jl |
Adds current-jump verification. |
examples/coupled_step/column_shifted_in_shell.jl |
Adds shell-force verification. |
examples/coupled_step/column_pushed_toward_loop.jl |
Adds loop-flux verification. |
examples/coupled_step/coil_driven_column.jl |
Adds model and gate checks. |
examples/coupled_step/gate_crossing_dense_column.jl |
Demonstrates corrected gate crossing. |
examples/coupled_step/picard_center_column.jl |
Adds central-column convergence case. |
examples/coupled_step/picard_filament_shell.jl |
Adds conducting-shell case. |
examples/coupled_step/picard_kstar_inboard_limited.jl |
Adds KSTAR boundary case. |
examples/coupled_step/picard_regime_map.jl |
Maps convergence across regimes. |
examples/coupled_step/picard_tight_box.jl |
Adds tight-boundary case. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
* Coupled solve: a solve that throws leaves the state as it found it The predictor's last field (Eϕ_self_prev), the friction's part at t^n (Rue_ei) and an unset coil's plasma flux memory were written before the outer iteration. They are now written when an evaluation is accepted; the coil memory the circuits start from comes from coil_plasma_flux_memory, which writes nothing. Results are unchanged bit for bit. * Examples: headers describe the gate crossing and the Anderson default coil_driven_column said the 1 A run jumps before the coupled solve takes over, which the gate re-solve and the script's own check contradict. picard_center_column described the old relaxed default (w = 0.5, 10 solves). * Numerics: a fixed-point solve that drives any stepper, and Newton's step anderson_solve! becomes fixed_point_solve!, which calls fixed_point_step! on its stepper: AndersonMixer, or the new NewtonStepper, x + (1 - T)^-1 f with a factorization of W (1 - T) W^-1. On an affine map Newton's step lands on the fixed point at once; a residual no smaller than the best stalls, which ends the solve (accepted if it meets the stopping test). * Coupled solve: method = :direct, the step solved without iterating coupled_map_jacobian assembles the outer map's linear part T column by column (N_b + N_c back-substitutions on the step's block LU), and Newton's step on W (1 - T) W^-1 lands on the fixed point from the first iterate; later evaluations only remove rounding. flags.ampere_picard.method selects it (:anderson by default). The coupled-solve tests now compare the default with this direct solve instead of a relaxed iteration run to 1e-10. * Coupled solve: remove the uncalled alternative solver and the old assembly solve_coupled_momentum_Ampere_equations_with_coils! (u-par eliminated, a full LU per iteration, the relaxed iteration, and an in-grid coil source without A_u), its stopping test picard_step_converged, and the export combine_Au_and_DeltaGS_sparse_matrices, which CoupledBlock replaced. Their tests now check the direct solve against the iteration, and the fixed-pattern block against the block matrix sparse algebra assembles. * Examples: the picard cases check the default against the direct solve DIRECT_PICARD (method = :direct, checked to 1e-10 with no floors) replaces the relaxed iteration run to 1e-10 as the expected step. All five cases pass with the same gaps as before; the reference takes 2 block solves per step instead of a few hundred. * Review fixes: a linear coil source, Newton's best iterate, and exact docs - distribute_coil_currents_to_Jphi! skipped currents below eps (in amperes), so the in-grid coils' source was not linear in their currents, which the direct solve assumes; it now skips exact zeros only. - NewtonStepper keeps its best iterate and returns it on :exhausted, as the fixed_point_step! contract says. - The coupled solve's docstring: a throw leaves u-par, the fields and the coils unchanged (the step's caches may be rebuilt), and the direct solve needs no eigenvalue of T equal to 1. * Coupled solve: the outer solve's method is a policy type OuterSolvePolicy, as IonTransportPolicy: AndersonOuterSolve(; memory, relaxation_w), the default, and DirectOuterSolve(), the reference. Anderson's parameters move from PicardSettings into its policy, which checks them (memory >= 0, a finite positive weight; above 1 over-relaxes). outer_stepper resolves the policy, so the solve reads the default first and the reference code (outer_stepper for DirectOuterSolve, coupled_map_jacobian) sits below it. * Coil needs a positive self-inductance; the Anderson weight is checked as stored A Coil with zero self-inductance (also the old default when it was left out) is refused when it is made: a toroidal loop always has some, and without it the coils' inductance matrix is indefinite and the direct solve's weights divide by it. The check runs once per coil, not per step. AndersonOuterSolve checks relaxation_w after converting it to Float64.


Summary
The coupled step advances u∥, Ampère's free boundary and the coil circuits together, closing them with an outer iteration every step. That iteration could diverge or stall without notice, and the step that carried the current across the Ampère gate left out the plasma's self-inductance. With this PR:
<name>__PASS.pngor<name>__FAIL.png.Also: the Green-function R-derivatives gain the term from their$\sqrt{R_d R_s}$ factor (not used in the step yet), and the coupled solve's settings are a run flag,
flags.ampere_picard(aPicardSettings, checked when set).Each item has a test. What remains is listed at the end.
The outer iteration
With the step's coefficients (density, collision rates, field direction) held at$t^n$ , the coupled step solves four sets of equations together at $t^{n+1}$ :
The first two couple each cell only to its neighbours, so$u_\parallel$ and $\psi$ form one sparse block ($x = (\psi_b, I_c)$ :
CoupledBlock, on a pattern fixed for the run), factorized once per step. The last two couple every cell to every edge node and coil. They are closed by an outer iteration on those few global unknowns,At frozen coefficients$g$ is affine, $g(x) = T x + c$ ($T$ : its linear part, $c$ : the rest). The relaxed iteration $x \leftarrow x + w f$ ($w$ : the relaxation weight) multiplies each eigenmode of $T$ , with eigenvalue $\mu$ , by $1 - w (1 - \mu)$ . It failed in three ways:
No single$w$ handles the first two together, and a failed solve was accepted after one warning.
Now:
AndersonMixer, memory 8) fits the latest residualanderson_m = 0is the relaxed iteration.tolerancetimes the field the step induces. Each coil's part is measured against that coil's change over the step.anderson_solve!, apart from the physics.The step that crosses the gate
Below the gate the plasma's self-inductance$L_p$ is left out. That holds while it is small next to the electrons' inertia, the kinetic inductance $L_\mathrm{kin} = 2\pi R m_e / (n_e e^2 A)$ ($A$ : the column's cross-section), as in an avalanche. A column that is already dense when the drive comes on has $L_p \gg L_\mathrm{kin}$ : its first step overshot by about $L_p / L_\mathrm{kin}$ .
Now, when the coupled solve can run, a step below the gate whose u∥ update crosses the gate is restored to its starting state and solved once by the coupled solve. That step is then exactly the step taken with the gate at 0. A current that decays through the gate, or crosses it by density growth, takes the old path.
Tests and examples
anderson_test.jl: the mixer andanderson_solve!on synthetic maps: divergent maps, restarts, evaluations that are not finite or are rejected, fallback to the best, single precision.ampere_picard_test.jl: the default solve against the same equations iterated to convergence, in a tight box, in a filament shell, and with advection and coils of mixed inductance. It compares the current, the induced field, each coil, u∥ and ψ. The stopping test, the fallback and the settings' checks are also tested on their own, and the block keeps one pattern while the flow turns.coupled_step_test.jl: a dense column crossing the 1 A gate steps exactly (==) as with the gate at 0. The same holds with coils, a remembered coil flux,advance_timestep!and a split run. After the crossing, every number in the plasma, fields, transport and coils equals the gate-0 step's, which guards the restored snapshot.green_function_test.jl,coil_plasma_flux_test.jl: the derivatives against finite differences, and the circuits' flux balance on solves cut short.examples/coupled_step/:picard_*puts the default solve against the converged one:gate_crossing_dense_columnshows the crossing.Behaviour changes
full_startupis unchanged). Where it did not, they become the converged step's. Solves take fewer iterations, and the coupled solve is cheaper, since its LU reuses the symbolic analysis.max_iterdefaults to 20 (was 10).Remaining (next PRs)
solve_coupled_momentum_Ampere_equations_with_coils!(not called) adds the in-grid coils' source without