Skip to content

Coupled step: a direct solve as the reference, and fixes after #28 - #29

Merged
mgyoo86 merged 9 commits into
masterfrom
feature/coupled-step-direct-reference
Oct 6, 2026
Merged

mgyoo86 merged 9 commits into
masterfrom
feature/coupled-step-direct-reference

Conversation

@mgyoo86

@mgyoo86 mgyoo86 commented Oct 5, 2026 •

Copy link
Copy Markdown
Member

Summary

A follow-up to #28. The coupled solve gains an opt-in direct method. It solves a step's equations at once, without iterating, and the default iterative solve is now checked against it. With this PR:

  • A reference that does not depend on the iteration converging. flags.ampere_picard.method = DirectOuterSolve() assembles the outer map's linear part and solves it in one Newton step. It replaces the earlier references: the relaxed iteration run to $10^{-10}$, which took hundreds of block solves per step, and the uncalled alternative solver. Its cost suits tests and coarse grids, so the default stays AndersonOuterSolve().
  • The outer solve's method is a policy type, as IonTransportPolicy is: AndersonOuterSolve(; memory, relaxation_w), the default, and DirectOuterSolve(), the reference. Anderson's parameters move into it from PicardSettings.
  • The uncalled alternative solver is removed, with its stopping test and the old block assembly.
  • A solve that throws leaves u∥, the fields and the coils as it found them.
  • Copilot's review comments on Coupled step: an outer iteration that converges, and a coupled step across the Ampère gate #28 are addressed (below).

Each item has a test.

The direct solve

The outer unknowns are the edge flux and the coil currents, $x = (\psi_b, I_c)$. With the step's coefficients held at $t^n$, every equation of the step is linear in its unknowns, so the map the iteration solves is affine, $g(x) = T x + c$. Here $T$ is the response of $x$ to itself through the plasma current, and $c$ collects the forcing. The step solves

$$ \left( \mathbb{1} - T \right) x = c $$

  • $T$ column by column (coupled_map_jacobian). Column $j$ is what a unit $x_j$ gives back. That takes one back-substitution with the step's block LU, which is already factorized, and then the Green's functions and the circuits.
  • One Newton step from the first iterate lands on the fixed point, unless an eigenvalue of $T$ is exactly 1. Where the iteration diverges or stalls makes no difference to it. The unknowns are weighted into one unit, as in Anderson mixing; in raw webers and amperes the system is badly conditioned. Later evaluations only remove rounding (NewtonStepper).
  • What it checks: it shows whether the iteration found the step's solution. It does not check the discretization; the examples with analytic predictions do that.
  • Cost: $N_b + N_c$ back-substitutions per step ($N_b$: edge nodes, $N_c$: coils).

anderson_solve! becomes fixed_point_solve!. It drives either AndersonMixer or NewtonStepper under the same acceptance policy; outer_stepper builds one from the policy.

Review comments on #28

  • A failed solve changed the state before it threw. The predictor's last field, the friction and an unset coil memory are now written on acceptance only. Results are unchanged bit for bit.
  • The gate snapshot allocates on every step below the gate. Left as it is: it is a few percent of that step's allocation and negligible in its time.
  • Two example headers described the old default and a jump at the gate that the re-solve had removed. Both are corrected.

Tests and examples

  • fixed_point_test.jl: Newton's step on an affine map that the relaxed iteration diverges on, with unknowns over twelve decades; a step that stalls; residuals that are not finite; the driver's handling of a stalled evaluation.
  • ampere_picard_test.jl:
    • the default solve against the direct one in the hard cases (tight box, filament shell, mixed coils);
    • the direct solve takes each hard step in two evaluations;
    • the relaxed iteration, Anderson mixing and the direct solve reach the same fixed point;
    • a solve that throws leaves the state unchanged.
  • coil_plasma_flux_test.jl: with a powered coil and a resistive loop, the direct solve takes the step the iteration converges to.
  • circuit_equations_test.jl: the coils' source on the grid is linear in their currents, however small.
  • examples/coupled_step/picard_*: the default solve against the direct solve.

Behaviour changes

  • Runs with the default settings: unchanged, apart from the deposition below.
  • Removed:
    • solve_coupled_momentum_Ampere_equations_with_coils!. It eliminated u∥, took a full LU every iteration, and added the in-grid coils' source without $A_u$.
    • picard_step_converged.
    • The export combine_Au_and_ΔGS_sparse_matrices, which CoupledBlock replaced.
  • PicardSettings: anderson_m and relaxation_w become method = AndersonOuterSolve(; memory, relaxation_w), or DirectOuterSolve(). A relaxation weight above 1 (over-relaxation) is now allowed; it must be finite and positive.
  • Coil needs a positive self-inductance. Zero, or leaving it out (it defaulted to zero), is refused when the coil is made. A toroidal loop always has some, and without it the coils' inductance matrix is indefinite; the direct solve's weights also divided by it.
  • Coil currents on the grid: a current below machine epsilon (in amperes) used to be skipped when it was deposited. It is now deposited, so the coils' source is linear in their currents.

Remaining (unchanged from #28)

Time levels at the step's start and end, plasma motion reaching the coils one step late, an older LU as a preconditioner, and the nonlinear coupling.

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.
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).
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_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.
…embly

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.
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.
…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.
@codecov

codecov Bot commented Oct 5, 2026 •

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 98.82353% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 93.64%. Comparing base (0bbc207) to head (c52ddab).
⚠️ Report is 1 commits behind head on master.

Files with missing lines Patch % Lines
src/numerics/fixed_point.jl 95.83% 1 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master      #29      +/-   ##
==========================================
- Coverage   93.69%   93.64%   -0.05%     
==========================================
  Files          47       50       +3     
  Lines        5075     5130      +55     
==========================================
+ Hits         4755     4804      +49     
- Misses        320      326       +6     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

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.

Copilot AI left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copilot review overview

🟡 Changes recommended

Direct solving fails for valid zero-inductance coils, and policy conversion can bypass its finite-positive invariant.

Review effort: Balanced
Findings: 2 Medium severity

Open (2)
What changed in this PR

Adds a direct reference solver for the coupled Ampère step while retaining Anderson mixing as the default.

Changes:

  • Introduces outer-solve policies and a shared fixed-point driver.
  • Adds direct affine-map solving with safer state acceptance.
  • Updates coil deposition, tests, and examples.
File Description
src/​types.jl Defines outer-solve policies and settings.
src/​RAPID2D.jl Includes fixed-point numerics.
src/​physics/​physics.jl Implements direct coupled solving.
src/​numerics/​fixed_point.jl Adds the shared solver and Newton stepper.
src/​numerics/​anderson.jl Adapts Anderson mixing to the shared driver.
src/​coils/​circuit_equations.jl Makes coil deposition linear and adds non-mutating memory access.
test/​unit/​physics/​ampere_picard_test.jl Tests policies, direct solving, and rollback.
test/​unit/​numerics/​fixed_point_test.jl Tests Newton and fixed-point behavior.
test/​unit/​numerics/​anderson_test.jl Updates shared-driver tests.
test/​unit/​coils/​coil_plasma_flux_test.jl Compares iterative and direct coil steps.
test/​unit/​coils/​circuit_equations_test.jl Tests small-current linearity.
examples/​README.md Documents direct-reference comparisons.
examples/​coupled_step/​common.jl Uses the direct reference in examples.
examples/​coupled_step/​picard_regime_map.jl Updates regime-map comparisons.
examples/​coupled_step/​picard_tight_box.jl Updates reference description.
examples/​coupled_step/​picard_kstar_inboard_limited.jl Updates reference description.
examples/​coupled_step/​picard_filament_shell.jl Updates reference description.
examples/​coupled_step/​picard_center_column.jl Documents current solver behavior.
examples/​coupled_step/​coil_driven_column.jl Documents gate-crossing behavior.

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

Comment thread src/physics/physics.jl
Comment thread src/types.jl
… 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.
@mgyoo86
mgyoo86 merged commit 72751ec into master Oct 6, 2026
5 checks passed
@mgyoo86
mgyoo86 deleted the feature/coupled-step-direct-reference branch October 6, 2026 17:10
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants