Edge dW scan: unbracketed Newton in free_compute_wv_spline throws ConvergenceFailed

Open
  • Is something wrong?
#648 0 comments 0 reactions 0 assignees View on GitHub

Maintainers usually reply within 3 days

Nobody has claimed this yet.

Assessment

Difficulty
3/5
Estimated time
1-2 days
Newbie friendliness
74/100
Issue type
Bug
Clarity
Clearly specified
Activity status
Active
Tech stack
julia

Research direction

Read free_compute_wv_spline in src/ForceFreeStates/Free.jl and follow the call path through _handle_edge_dW_scan!, scan_edge_dW!, and findmax_dW_edge! in EulerLagrange.jl. Check how q_spline behaves over [dW_edge_scan_start, psilim], then reproduce the n = 1 edge scan if possible. Done means the scan completes when an unbracketed Newton solve would fail and existing converged cases retain the same roots.

Written by the indexing model from the issue text.

Description

Summary

free_compute_wv_spline (src/ForceFreeStates/Free.jl) inverts q(ψ) at each edge-scan node using an unbracketed Roots.Newton() started from a ψ-linear guess. On one equilibrium in a parameter scan, Newton did not converge. The whole run then aborted with Roots.ConvergenceFailed("Algorithm failed to converge") during the n = 1 edge δW scan.

The same code is still on develop, at about line 222.

Where

qi   = qedge + (intr.qlim - qedge) * ((i - 1) / npsi)
psii = ctrl.dW_edge_scan_start + (intr.psilim - ctrl.dW_edge_scan_start) * ((i - 1) / npsi)
psi_array[i] = find_zero((psi -> profiles.q_spline(psi) - qi, psi -> profiles.q_deriv(psi)),
                         psii, Roots.Newton())

Call path: _handle_edge_dW_scan! → scan_edge_dW! → findmax_dW_edge! (EulerLagrange.jl:1140) → free_compute_wv_spline (Free.jl:232).

Context

  • Build: 728a7e7d5 (develop).
  • Settings: n = 1, no wall, psihigh = 0.998, mpsi = 256, grid_type = "core_edge_packed", dW_edge_scan_start = 0.97, and psilim set from rational_offset_fraction = 0.2.
  • The equilibrium is an H-mode with q95 ≈ 4, where q rises steeply over the last few percent in ψ. The g-file's FPOL is single precision, and GPEC already warns about that.
  • The failure came from one point of a 96-point scan. Its neighbors, with slightly lower and slightly higher pedestal-region current, ran normally, so the failure looks like Newton's sensitivity to its starting guess rather than a bad equilibrium.
  • I can't attach the equilibrium publicly, but I can share it privately.

Likely cause

The target q values are spaced evenly in q, while the starting guesses are spaced evenly in ψ. Where q(ψ) is strongly convex near the edge, the guess sits far from the root. A Newton step can then overshoot past psilim, into the spline's ExtendExtrap region, or oscillate.

Suggested fix

Over [dW_edge_scan_start, psilim], q is monotonic and the root is bracketed by construction. A bracketed solve would always converge:

psi_array[i] = find_zero(psi -> profiles.q_spline(psi) - qi,
                         (ctrl.dW_edge_scan_start, intr.psilim), Roots.Brent())

Alternatively, keep Newton and fall back to the bracketed solve when it fails. Wherever Newton converges today, either option returns the same root, so no existing result changes.

Dominant language
Julia
Stars
11
Forks
7
Avg merge
3d 22h
Merged PRs (30d)
70

Getting set up

  • No Dockerfile or Docker Compose file
  • Has a pull request template
  • No contributing guide

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

More from OpenFUSIONToolkit/GPEC

All issues in OpenFUSIONToolkit/GPEC

Similar issues

More Julia issues

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.