diff --git a/src/Equilibrium/Axisymmetric/Direct.jl b/src/Equilibrium/Axisymmetric/Direct.jl index e931091ac..1e15b4943 100644 --- a/src/Equilibrium/Axisymmetric/Direct.jl +++ b/src/Equilibrium/Axisymmetric/Direct.jl @@ -116,77 +116,219 @@ the Julia spline implementation. bf_out.bzz = -bf_out.psirz / r end +# Numerical controls of the magnetic-axis search. +const AXIS_MARCH_STEPS = 200 # midplane march steps allowed before giving up on a B_z sign change +const AXIS_NEWTON_ITERATIONS = 200 +const AXIS_SINGULAR_DET = 1e-20 # absolute floor on |det ∂(B_R,B_Z)/∂(R,Z)| in T²/m²; it does not scale with ψ₀ or machine size +const AXIS_STEP_RTOL = 1e-12 # Newton step test, relative to R +const AXIS_RESIDUAL_RTOL = 1e-10 # bound on |B_p| relative to B_scale, applied on the restart path only +const AXIS_SEED_CELL_TOL = 0.1 # the bisection sweep has settled when a pass moves the point by less than this many cells +const AXIS_SEED_PASSES = 50 +const AXIS_PSI_RATIO = (0.5, 2.0) # accepted band of ψ(axis)/ψ₀; every reader puts ψ = ψ₀ at the axis +const SEPARATRIX_STEPS = 2000 # the Fortran allows more Newton steps here than for the axis solve + +_axis_jacobian_det(bf::DirectBField) = bf.brr * bf.bzz - bf.brz * bf.bzr + +# Poloidal-field scale |ψ₀| / (R ΔR), about a sixth of the edge poloidal field. +_axis_field_scale(rp::DirectRunInput, r::Float64) = abs(rp.psio) / (r * (rp.rmax - rp.rmin)) + """ - direct_position!(raw_profile) + _midplane_axis_guess(rp, bfield, sq_in_deriv) -> (r, z, step) -Finds the key geometric locations of the equilibrium: the magnetic axis (O-point) -and the inboard/outboard separatrix crossings on the midplane. It also updates the -spline representing the poloidal flux `ψ(R,Z)` based on the new magnetic axis location. -This function performs the same overall function as the Fortran `direct_position` -subroutine with better iteration control and error handling. We have also added a -helper function for separatrix finding. +Initial guess for the axis: march outward along the grid's midplane `Z`, in steps of a twentieth of the grid width, until +`B_z ≥ 0`. +""" +function _midplane_axis_guess(rp::DirectRunInput, bfield::DirectBField, sq_in_deriv) + r = (rp.rmax + rp.rmin) / 2.0 + z = (rp.zmax + rp.zmin) / 2.0 + step = (rp.rmax - rp.rmin) / 20.0 + for _ in 1:AXIS_MARCH_STEPS + direct_get_bfield!(bfield, r, z, rp.psi_in, rp.sq_in, sq_in_deriv, rp.psio; derivs=1) + bfield.bz >= 0.0 && return r, z, step + r += step + end + error("Took too many iterations to get bz=0.") +end -## Arguments: +""" + _axis_newton(rp, bfield, sq_in_deriv, r, z, cap) -> (r, z, failure) - - `raw_profile`: A `DirectRunInput` object containing splines and parameters. +Newton iteration for the O-point `B_R = B_Z = 0`, each step capped at `cap`. `failure` is `nothing` once a step passes +`|δR|, |δZ| ≤ AXIS_STEP_RTOL × R`, else why the iteration stopped. A singular Jacobian stops it, where the Fortran divides +by the determinant unguarded. +""" +function _axis_newton(rp::DirectRunInput, bfield::DirectBField, sq_in_deriv, r::Float64, z::Float64, cap::Float64) + for _ in 1:AXIS_NEWTON_ITERATIONS + direct_get_bfield!(bfield, r, z, rp.psi_in, rp.sq_in, sq_in_deriv, rp.psio; derivs=2) + det = _axis_jacobian_det(bfield) + abs(det) < AXIS_SINGULAR_DET && return r, z, "stopped on a singular Jacobian, |det ∂(B_R,B_Z)/∂(R,Z)| < $AXIS_SINGULAR_DET" + # Δx = -J⁻¹ F + δr = (bfield.brz * bfield.bz - bfield.bzz * bfield.br) / det + δz = (bfield.bzr * bfield.br - bfield.brr * bfield.bz) / det + step = hypot(δr, δz) + step > cap && ((δr, δz) = (δr * cap / step, δz * cap / step)) + r += δr + z += δz + abs(δr) <= AXIS_STEP_RTOL * abs(r) && abs(δz) <= AXIS_STEP_RTOL * abs(r) && return r, z, nothing + end + return r, z, "stopped without convergence in $AXIS_NEWTON_ITERATIONS iterations" +end -## Returns: +""" + _axis_bisection_seed(rp, bfield, sq_in_deriv, r, z, halfwidth; passes) -> (r, z, failure) - - `ro`: R-coordinate of the magnetic axis [m]. - - `zo`: Z-coordinate of the magnetic axis [m]. - - `rs1`: R-coordinate of the inboard separatrix crossing [m]. - - `rs2`: R-coordinate of the outboard separatrix crossing [m]. - - `psi_in_new` : returns psi_in renormalized by * psio/psi(ro,zo) +Seed for the Newton restart: alternate 1-D bisections on `B_Z` along `R` and on `B_R` along `Z`, each bracketed by doubling +`halfwidth` up to the grid edges. The sweep is a nonlinear Gauss–Seidel iteration that contracts by ψ_RZ²/(ψ_RR ψ_ZZ) per +pass, so stopping once a pass moves the point by less than `AXIS_SEED_CELL_TOL` cells leaves it within about a cell of the +axis for contractions up to 0.9. `failure` names a missing sign change or a sweep that did not settle in `passes` passes. """ -function direct_position!(raw_profile::DirectRunInput) +function _axis_bisection_seed(rp::DirectRunInput, bfield::DirectBField, sq_in_deriv, r::Float64, z::Float64, halfwidth::Float64; + passes::Int=AXIS_SEED_PASSES) + _bz(rr, zz) = (direct_get_bfield!(bfield, rr, zz, rp.psi_in, rp.sq_in, sq_in_deriv, rp.psio; derivs=1); bfield.bz) + _br(rr, zz) = (direct_get_bfield!(bfield, rr, zz, rp.psi_in, rp.sq_in, sq_in_deriv, rp.psio; derivs=1); bfield.br) + hr = (rp.rmax - rp.rmin) / (length(rp.psi_in_xs) - 1) + hz = (rp.zmax - rp.zmin) / (length(rp.psi_in_ys) - 1) + # Nearest sign change of f around x, bisected to a hundredth of a cell h; `nothing` if the whole grid line has none. + function _bracketed_root(f, x, xmin, xmax, h) + hw = halfwidth + while true + lo, hi = max(x - hw, xmin), min(x + hw, xmax) + sign(f(lo)) * sign(f(hi)) <= 0 && return find_zero(f, (lo, hi), Roots.Bisection(); xatol=0.01 * h) + (lo == xmin && hi == xmax) && return nothing + hw *= 2 + end + end + for _ in 1:passes + r_new = _bracketed_root(rr -> _bz(rr, z), r, rp.rmin, rp.rmax, hr) + r_new === nothing && return r, z, "B_z has no sign change along Z = $z" + z_new = _bracketed_root(zz -> _br(r_new, zz), z, rp.zmin, rp.zmax, hz) + z_new === nothing && return r, z, "B_r has no sign change along R = $r_new" + moved = max(abs(r_new - r) / hr, abs(z_new - z) / hz) + r, z = r_new, z_new + moved < AXIS_SEED_CELL_TOL && return r, z, nothing + end + return r, z, "the bisection sweep did not settle in $passes passes" +end - bfield = DirectBField() - max_iterations = 200 - sq_in_deriv = deriv1(raw_profile.sq_in) +""" + _axis_failure(rp, bfield, sq_in_deriv, r, z; bound_residual) -> failure - # For an axis initial guess, we find zero crossing of Bz on the midplane - # Note the Fortran had this wrapped in a if ro == 0 block, which we omit here - # because I don't think it's ever used. If needed, it can be re-added. - r = (raw_profile.rmax + raw_profile.rmin) / 2.0 - z = (raw_profile.zmax + raw_profile.zmin) / 2.0 - dr = (raw_profile.rmax - raw_profile.rmin) / 20.0 - - for _ in 1:max_iterations - direct_get_bfield!(bfield, r, z, raw_profile.psi_in, raw_profile.sq_in, sq_in_deriv, raw_profile.psio; derivs=1) - if bfield.bz >= 0.0 - break - end - r += dr +Why `(r, z)` is not an acceptable magnetic axis, or `nothing`. It must lie on the ψ grid (a root off it is an artifact of +spline extrapolation), be an O-point (`det ∂(B_R,B_Z)/∂(R,Z) = det(Hess ψ)/R² > 0`; an X-point has `det < 0`) that is a +maximum of ψ (`ψ_RR < 0` and `ψ > 0`), and carry ψ within `AXIS_PSI_RATIO` of ψ₀, which rejects a secondary maximum +elsewhere in the grid. With `bound_residual`, `|B_p| ≤ AXIS_RESIDUAL_RTOL × B_scale` must hold as well. +""" +function _axis_failure(rp::DirectRunInput, bfield::DirectBField, sq_in_deriv, r::Float64, z::Float64; bound_residual::Bool) + (rp.rmin <= r <= rp.rmax && rp.zmin <= z <= rp.zmax) || + return "converged at (R, Z) = ($r, $z), outside the ψ grid R ∈ [$(rp.rmin), $(rp.rmax)], Z ∈ [$(rp.zmin), $(rp.zmax)]" + direct_get_bfield!(bfield, r, z, rp.psi_in, rp.sq_in, sq_in_deriv, rp.psio; derivs=2) + det = _axis_jacobian_det(bfield) + residual = hypot(bfield.br, bfield.bz) + det > 0 || return "converged at (R, Z) = ($r, $z), which is not an O-point (det ∂(B_R,B_Z)/∂(R,Z) = $det ≤ 0, |B_p| = $residual T)" + (bfield.psirr < 0 && bfield.psi > 0) || + return "converged at (R, Z) = ($r, $z), an O-point that is not the maximum of ψ (ψ = $(bfield.psi), ψ_RR = $(bfield.psirr))" + lo, hi = AXIS_PSI_RATIO + lo < bfield.psi / rp.psio < hi || + return "converged at (R, Z) = ($r, $z), a maximum of ψ that is not the main axis (ψ/ψ₀ = $(bfield.psi / rp.psio), outside ($lo, $hi))" + if bound_residual + bound = AXIS_RESIDUAL_RTOL * _axis_field_scale(rp, r) + residual <= bound || return "converged at (R, Z) = ($r, $z) with |B_p| = $residual T, above $AXIS_RESIDUAL_RTOL × B_scale = $bound T" end + return nothing +end - # If we never exited early, the loop failed to find bz = 0 - !(bfield.bz >= 0) && error("Took too many iterations to get bz=0.") +# Field diagnostics at an accepted axis, for the restart warning. +function _axis_diagnostics(rp::DirectRunInput, bfield::DirectBField, sq_in_deriv, r::Float64, z::Float64) + direct_get_bfield!(bfield, r, z, rp.psi_in, rp.sq_in, sq_in_deriv, rp.psio; derivs=2) + residual = hypot(bfield.br, bfield.bz) + return "at (R, Z) = ($r, $z) with |B_p| = $residual T ($(residual / _axis_field_scale(rp, r)) × B_scale) and " * + "det ∂(B_R,B_Z)/∂(R,Z) = $(_axis_jacobian_det(bfield)) T²/m²" +end - # Now, use Newton iteration to find the O-point (magnetic axis) where Br=0 and Bz=0 - dr, dz = 0.0, 0.0 - for _ in 1:max_iterations - direct_get_bfield!(bfield, r, z, raw_profile.psi_in, raw_profile.sq_in, sq_in_deriv, raw_profile.psio; derivs=2) - det = bfield.brr * bfield.bzz - bfield.brz * bfield.bzr - if abs(det) < 1e-20 - error("Jacobian matrix is singular near ($r, $z).") +""" + _find_magnetic_axis(rp, bfield, sq_in_deriv, r0, z0, march_step) -> (r, z) + +The magnetic axis, by the first attempt whose Newton result passes `_axis_failure`: plain Newton from the midplane guess +`(r0, z0)`, accepted on the step test alone so that smooth files keep the historical arithmetic, then Newton restarted from +a bisection seed with steps capped at one march step and the residual bound added. A restart warns with the field at the +axis; two failures error. +""" +function _find_magnetic_axis(rp::DirectRunInput, bfield::DirectBField, sq_in_deriv, r0::Float64, z0::Float64, march_step::Float64) + restart_cap = march_step # Newton step cap on the restart path + seed_halfwidth = march_step # initial bracket half-width of the bisection seed + attempts = ( + (label="Newton from the midplane guess", seed=() -> (r0, z0, nothing), cap=Inf, bound_residual=false), + (label="Newton restarted from a first-derivative bisection seed", + seed=() -> _axis_bisection_seed(rp, bfield, sq_in_deriv, r0, z0, seed_halfwidth), cap=restart_cap, bound_residual=true) + ) + failures = String[] + for a in attempts + rs, zs, why = a.seed() + if why !== nothing + push!(failures, "$(a.label): the seed failed, $why") + continue end - # Δx = -J⁻¹ F - dr = (bfield.brz * bfield.bz - bfield.bzz * bfield.br) / det - dz = (bfield.bzr * bfield.br - bfield.brr * bfield.bz) / det - r += dr - z += dz - if abs(dr) <= 1e-12 * abs(r) && abs(dz) <= 1e-12 * abs(r) - @info "Magnetic axis found at R = $(@sprintf("%.3f", r)), Z = $(@sprintf("%.3f", z))" - break + r, z, why = _axis_newton(rp, bfield, sq_in_deriv, rs, zs, a.cap) + why === nothing && (why = _axis_failure(rp, bfield, sq_in_deriv, r, z; bound_residual=a.bound_residual)) + if why === nothing + isempty(failures) || @warn "$(join(failures, "; ")); magnetic axis found by $(a.label) " * + "$(_axis_diagnostics(rp, bfield, sq_in_deriv, r, z)). Check ψ near the axis in the equilibrium file." + return r, z end + push!(failures, "$(a.label) $why") end + error("Failed to find magnetic axis: $(join(failures, "; ")).") +end - if !(abs(dr) <= 1e-12 * abs(r) && abs(dz) <= 1e-12 * abs(r)) - error("Failed to find magnetic axis after $max_iterations iterations.") +""" + _separatrix_crossing(rp, bfield, sq_in_deriv, ro, zo, r_edge, label) -> r + +Midplane separatrix crossing on the `r_edge` side of the axis: unbracketed Newton `δR = −ψ/ψ_R` toward ψ = 0 of the +renormalized spline, as in the Fortran, so that it extrapolates past the grid edge for limited equilibria whose box hugs +the LCFS. Each failed attempt moves the initial guess, a blend of `r_edge` and the axis, inward. +""" +function _separatrix_crossing(rp::DirectRunInput, bfield::DirectBField, sq_in_deriv, ro::Float64, zo::Float64, r_edge::Float64, label::String) + for restart in 0:5 + r = ((3 - 0.5 * restart) * r_edge + ro) / (4 - 0.5 * restart) + for _ in 1:SEPARATRIX_STEPS + direct_get_bfield!(bfield, r, zo, rp.psi_in, rp.sq_in, sq_in_deriv, rp.psio; derivs=1) + dr = -bfield.psi / bfield.psir + r += dr + if abs(dr) <= AXIS_STEP_RTOL * abs(r) + @info "$label separatrix found at R = $(@sprintf("%.3f", r))" + return r + end + end end + error("Took too many steps to find $label separatrix.") +end - ro = r - zo = z +""" + direct_position!(raw_profile) + +Finds the magnetic axis (O-point) and the inboard and outboard separatrix crossings on its midplane, and replaces the +`ψ(R, Z)` spline in `raw_profile` by one renormalized so that ψ at the axis equals ψ₀. The axis is a root of `∇ψ = 0` +that passes the Newton step test `|δR|, |δZ| ≤ 1e-12 R` and the checks of `_axis_failure`, found by plain Newton from a +midplane guess or, when that fails, by a restart from a bisection seed with capped steps (see `_find_magnetic_axis`). +The restart warns and a second failure errors. A g-file whose header signs disagree with its ψ grid errors here, where +the Fortran would have flipped the sign through the renormalization. + +## Arguments: + + - `raw_profile`: A `DirectRunInput` object containing splines and parameters; its `psi_in` spline is replaced. + +## Returns: + + - `ro`, `zo`: R and Z coordinates of the magnetic axis [m]. + - `rs1`, `rs2`: R-coordinates of the inboard and outboard separatrix crossings [m]. +""" +function direct_position!(raw_profile::DirectRunInput) + bfield = DirectBField() + sq_in_deriv = deriv1(raw_profile.sq_in) + + r0, z0, march_step = _midplane_axis_guess(raw_profile, bfield, sq_in_deriv) + ro, zo = _find_magnetic_axis(raw_profile, bfield, sq_in_deriv, r0, z0, march_step) + @info "Magnetic axis found at R = $(@sprintf("%.3f", ro)), Z = $(@sprintf("%.3f", zo))" # Renormalize psi based on the value at the magnetic axis direct_get_bfield!(bfield, ro, zo, raw_profile.psi_in, raw_profile.sq_in, sq_in_deriv, raw_profile.psio; derivs=0) @@ -197,32 +339,8 @@ function direct_position!(raw_profile::DirectRunInput) # Because DirectRunInput is a mutable struct, we can update the spline here raw_profile.psi_in = cubic_interp((x_coords, y_coords), new_psi_fs; extrap=ExtendExtrap()) - # Replicate Fortran direct.f direct_position: unbracketed Newton iteration dr=-ψ/ψ_r toward - # ψ=0 (the renormalized separatrix). Newton extrapolates past the grid edge, so it finds the - # separatrix even for limited equilibria whose grid box hugs the LCFS (ψ>0 to the edge, no - # sign change to bracket). Fortran allows more steps here than for the axis solve (2000 vs 200). - max_separatrix_steps = 2000 - function find_separatrix_crossing(r_edge, label) - # On each failed attempt, shift the initial guess (a weighted blend of r_edge and the axis - # ro) from near the grid edge toward the axis (Fortran direct.f restart counter ird = 0..5). - for restart in 0:5 - r = ((3 - 0.5 * restart) * r_edge + ro) / (4 - 0.5 * restart) - for _ in 1:max_separatrix_steps - direct_get_bfield!(bfield, r, zo, raw_profile.psi_in, raw_profile.sq_in, sq_in_deriv, raw_profile.psio; derivs=1) - dr = -bfield.psi / bfield.psir - r += dr - if abs(dr) <= 1e-12 * abs(r) - @info "$label separatrix found at R = $(@sprintf("%.3f", r))" - return r - end - end - end - error("Took too many steps to find $label separatrix.") - end - - # Find inboard (rs1) and outboard (rs2) separatrix positions - rs1 = find_separatrix_crossing(raw_profile.rmin, "inboard") - rs2 = find_separatrix_crossing(raw_profile.rmax, "outboard") + rs1 = _separatrix_crossing(raw_profile, bfield, sq_in_deriv, ro, zo, raw_profile.rmin, "inboard") + rs2 = _separatrix_crossing(raw_profile, bfield, sq_in_deriv, ro, zo, raw_profile.rmax, "outboard") return ro, zo, rs1, rs2 end diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index fc7b01980..6016d9f88 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -16,7 +16,12 @@ using TOML psilow=0.01, psihigh=0.994 ) - global plasma_eq_efit = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(efit_config) + # A healthy g-file takes the plain Newton path to the magnetic axis, without the bisection restart. + logger = Test.TestLogger(; min_level=Base.CoreLogging.Warn) + global plasma_eq_efit = Base.CoreLogging.with_logger(logger) do + GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(efit_config) + end + @test !any(occursin("bisection", log.message) for log in logger.logs) @test plasma_eq_efit isa GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium @test 6.5 < plasma_eq_efit.ro < 7.5 # Physical sanity check diff --git a/test/runtests_equil_axis_newton.jl b/test/runtests_equil_axis_newton.jl new file mode 100644 index 000000000..b144f9846 --- /dev/null +++ b/test/runtests_equil_axis_newton.jl @@ -0,0 +1,134 @@ +@testset "Direct equilibrium: magnetic-axis search robustness" begin + using GeneralizedPerturbedEquilibrium.Equilibrium + using GeneralizedPerturbedEquilibrium.Equilibrium: EquilibriumConfig, SolovevConfig, DirectRunInput, sol_run, direct_position! + using FastInterpolations + EQ = GeneralizedPerturbedEquilibrium.Equilibrium + + fallback = r"restarted from a first-derivative bisection seed" + # Relative Newton step tolerance of the axis search, reused as the bound on the residual |∇ψ|. + newton_tol = 1e-12 + + # A `DirectRunInput` carrying the cubic spline of `psi_fn(R, Z)` on the grid `rs × zs`, with `psio` as the declared axis flux. + function direct_input(psi_fn, rs, zs, psio) + nr = length(rs) - 1 + sol = SolovevConfig(; mr=nr, mz=nr) + cfg = EquilibriumConfig(; eq_type="sol") + psi_in = cubic_interp((rs, zs), [psi_fn(R, Z) for R in rs, Z in zs]; extrap=ExtendExtrap()) + return DirectRunInput(cfg, sol_run(cfg, sol).sq_in, psi_in, rs, zs, rs[1], rs[end], zs[1], zs[end], psio, 1, 1, nothing) + end + + # Up-down symmetric Solovev ψ whose axis sits off the grid's R-midpoint and off every node, plus a + # conical cusp -A·ρ·exp(-(ρ/w)⁴) at the axis: the cells there carry a curvature spike, and the + # cusp strength A is set so that a Newton step overshoots the axis by `ncell` grid cells. + # `r_in`, `r_out` set the grid's extent about r0 in units of a; with r_in > r_out the grid's R-midpoint lies inboard + # of the axis, so the midplane march takes a step and brackets the axis in R. `z_shift` moves the grid, and with it + # the march's midplane, in Z in units of a. + function cusp_solovev(; ncell, nr=96, r_in=1.4, r_out=1.6, z_shift=0.0) + (; elongation, a, r0, q0, b0fac) = SolovevConfig() + e = elongation + psio = e * r0 * b0fac * a^2 / (2 * q0 * r0) + psifac = psio / (a * r0)^2 + rs = collect(range(r0 - r_in * a, r0 + r_out * a; length=nr + 1)) + zh = collect(range(0.0, 1.5e * a; length=nr ÷ 2 + 1)) + zs = z_shift * a .+ vcat(-reverse(zh[2:end]), zh) + h = rs[2] - rs[1] + A = ncell * 2 * psifac * r0^2 * h + cusp(ρ) = A * ρ * exp(-(ρ / (0.5a))^4) + psi(R, Z) = psio - psifac * ((R * Z)^2 / e^2 + (R^2 - r0^2)^2 / 4) - cusp(hypot(R - r0, Z)) + return direct_input(psi, rs, zs, psio), r0, h + end + + # Quadratic ψ = ψ₀ − c (x² + 2b·x·y + y²/e²) about (r0, 0), exact on the cubic spline; `sgn = -1` turns the + # axis into a minimum of ψ. With b·e near 1 the Hessian is strongly tilted (ψ_RZ ≠ 0), where the Gauss–Seidel + # bisection sweep contracts by (b e)² per pass. `psio_declared` is the ψ₀ the input claims, against the maximum of 1. + function quadratic_psi(; b, sgn=1.0, nr=64, psio_declared=1.0) + (; elongation, a, r0) = SolovevConfig() + e = elongation + c = 1.0 / a^2 + rs = collect(range(r0 - 1.5a, r0 + 1.5a; length=nr + 1)) + zs = collect(range(-1.5e * a, 1.5e * a; length=nr + 1)) + psi(R, Z) = 1.0 - sgn * c * ((R - r0)^2 + 2b * (R - r0) * Z + Z^2 / e^2) + return direct_input(psi, rs, zs, psio_declared), r0, rs[2] - rs[1], e + end + + # The returned axis must be an O-point of the ψ spline with a vanishing gradient. + function test_axis_is_o_point(rp, ro, zo) + d(i, j) = rp.psi_in((ro, zo); deriv=DerivOp(i, j)) + @test d(2, 0) * d(0, 2) - d(1, 1)^2 > 0 + @test hypot(d(1, 0), d(0, 1)) < newton_tol * rp.psio / (rp.rmax - rp.rmin) + end + + # Newton on the cusp is chaotic: whether it converges, cycles or leaves the grid at a given cusp strength depends + # on the machine's rounding, so which runs take the restart is not deterministic; that some do over this many is. + # Whatever it does, the search must return the O-point at the cusp or fail loudly, and the smooth ψ (ncell = 0) + # must take the plain Newton path on every grid. + @testset "curvature spike at the axis: a verified axis or a loud failure" begin + grids = ((;), (; r_in=1.6, r_out=1.4), (; r_in=1.7, r_out=1.3, z_shift=0.7)) + # On the shifted grid the march's midplane lies off the axis's, where the Solovev B_z changes sign at + # R² = r0² - 2Z²/e², inboard of the axis: the march stops short and the axis lies outside its last step. + (; elongation) = SolovevConfig() + rp, r0, _ = cusp_solovev(; ncell=0.0, grids[3]...) + r_mid, z_mid, dr = (rp.rmin + rp.rmax) / 2, (rp.zmin + rp.zmax) / 2, (rp.rmax - rp.rmin) / 20 + @test r_mid < sqrt(r0^2 - 2 * z_mid^2 / elongation^2) < r_mid + dr < r0 + + n_fallback = 0 + for grid in grids, ncell in 0.0:0.25:3.0 + rp, r0, h = cusp_solovev(; ncell, grid...) + logger = Test.TestLogger(; min_level=Base.CoreLogging.Warn) + outcome = try + Base.CoreLogging.with_logger(() -> direct_position!(rp), logger) + catch err + err + end + if outcome isa Exception + @test ncell > 0 + @test occursin("Failed to find magnetic axis", sprint(showerror, outcome)) + else + ro, zo = outcome[1], outcome[2] + test_axis_is_o_point(rp, ro, zo) + @test hypot(ro - r0, zo) < h + fell_back = any(occursin(fallback, log.message) for log in logger.logs) + n_fallback += fell_back + if ncell == 0 + @test !fell_back + get(grid, :z_shift, 0.0) == 0.0 && @test abs(zo) <= eps(ro) + end + end + end + @test n_fallback > 0 + end + + @testset "Newton converging off the ψ grid is an error" begin + # A smooth ψ on a grid that lies wholly outboard of the axis: Newton follows the spline's extrapolation inboard. + rp, r0, _ = cusp_solovev(; ncell=0.0, r_in=-0.2, r_out=1.6) + @test rp.rmin > r0 + @test_throws r"outside the ψ grid" direct_position!(rp) + end + + @testset "axis checks accept the main maximum of ψ and reject a minimum or a secondary maximum" begin + rp, r0, _, _ = quadratic_psi(; b=0.0) + bf = EQ.DirectBField() + @test EQ._axis_failure(rp, bf, deriv1(rp.sq_in), r0, 0.0; bound_residual=true) === nothing + rp_min, _, _, _ = quadratic_psi(; b=0.0, sgn=-1.0) + msg = EQ._axis_failure(rp_min, bf, deriv1(rp_min.sq_in), r0, 0.0; bound_residual=false) + @test msg isa String && occursin("not the maximum of ψ", msg) + rp_low, _, _, _ = quadratic_psi(; b=0.0, psio_declared=10.0) + msg = EQ._axis_failure(rp_low, bf, deriv1(rp_low.sq_in), r0, 0.0; bound_residual=false) + @test msg isa String && occursin("not the main axis", msg) + end + + @testset "restart seed and capped Newton recover a tilted axis" begin + rp, r0, h, e = quadratic_psi(; b=0.9 / 1.6) + bf, sqd = EQ.DirectBField(), deriv1(rp.sq_in) + w = (rp.rmax - rp.rmin) / 20 + # Start ten cells away in both directions; a few passes at (b e)² = 0.81 leave it several cells off, so a + # sweep allowed one pass must report that it did not settle. + @test occursin("did not settle", EQ._axis_bisection_seed(rp, bf, sqd, r0 + 10h, 10h, w; passes=1)[3]) + rb, zb, failure = EQ._axis_bisection_seed(rp, bf, sqd, r0 + 10h, 10h, w) + @test failure === nothing + @test hypot(rb - r0, zb) < h + ro, zo, failure = EQ._axis_newton(rp, bf, sqd, rb, zb, w) + @test failure === nothing + @test abs(ro - r0) < 1e-10 && abs(zo) < 1e-10 + end +end diff --git a/test/suites.jl b/test/suites.jl index 32fdaa645..225eb2fd8 100644 --- a/test/suites.jl +++ b/test/suites.jl @@ -4,6 +4,7 @@ const TEST_SUITES = [ "runtests_fouriertransforms.jl" => "core", "runtests_vacuum.jl" => "core", "runtests_equil.jl" => "core", + "runtests_equil_axis_newton.jl" => "core", "runtests_sfl_sampling.jl" => "core", "runtests_grid_refinement.jl" => "api", "runtests_core_capped_knots.jl" => "core",