Skip to content

Solve Boozer chartmap cold inverse on the full torus #397

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
krystophny merged 3 commits into from
Jul 16, 2026
Merged

Solve Boozer chartmap cold inverse on the full torus #397

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e40e06c
Add failing regressions for Boozer chartmap cold inverse (#395)
krystophny Jul 16, 2026
a00952b
Solve Boozer chartmap cold inverse on the full torus
krystophny Jul 16, 2026
2a4e81a
Align cold accept threshold with solver floors, harden seam test
krystophny Jul 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
104 changes: 72 additions & 32 deletions src/coordinates/libneo_coordinates_chartmap.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -803,7 +803,14 @@ subroutine chartmap_invert_cart_warm(self, x, u_guess, u, geom_status)

u(1) = min(max(u(1), 0.0_dp), 1.0_dp)
u(2) = modulo(u(2), TWOPI)
u(3) = modulo(u(3), zeta_period)
! The rz chart is wedge-native (its forward map reduces zeta internally);
! the co-rotating cart chart is only 2pi-periodic in zeta, so wrapping into
! the wedge would move the returned point by a field-period rotation.
if (self%has_spl_rz) then
u(3) = modulo(u(3), zeta_period)
else
u(3) = modulo(u(3), TWOPI)
end if

rscale = chartmap_radial_scale(self, u)
geom_status = chartmap_classify_root(u(1), res_best, rscale)
Expand Down Expand Up @@ -1105,10 +1112,13 @@ subroutine chartmap_invert_cart(self, x_target, u, ierr)
real(dp) :: res_norm, best_res
real(dp) :: zeta_period
real(dp) :: phi
real(dp) :: phi_shift, ca, sa
real(dp) :: x_wedge(3), x_round(3)
real(dp) :: accept_tol
real(dp) :: zeta_seed(5)
logical :: trace
integer :: k
integer :: ierr_try
integer :: ierr_seed, ierr_try

ierr = chartmap_from_cyl_err_max_iter
best_res = huge(1.0_dp)
Expand All @@ -1117,25 +1127,44 @@ subroutine chartmap_invert_cart(self, x_target, u, ierr)
zeta_period = TWOPI/real(self%num_field_periods, dp)
if (.not. self%has_spl_rz) then
trace = chartmap_trace_once_enabled()
call chartmap_initial_guess_cartesian_3d(self, x_target, u_try)
call chartmap_refine_cart_seed(self, x_target, zeta_period, u_try, &
res_norm, ierr)
if (ierr == chartmap_from_cyl_ok) then
call chartmap_solve_cart(self, x_target, u_try(1), u_try(2), u_try(3), &
ierr_try, trace)
if (ierr_try == chartmap_from_cyl_ok) then
ierr = chartmap_from_cyl_ok
end if
else
call chartmap_solve_cart(self, x_target, u_try(1), u_try(2), u_try(3), &
ierr_try, trace)
if (ierr_try == chartmap_from_cyl_ok) ierr = chartmap_from_cyl_ok
end if
! The co-rotating cart chart is field-period symmetric,
! x(rho, theta, zeta + k*zeta_period) = Rot_z(k*zeta_period) x(rho,
! theta, zeta), but the lab-frame position itself is only 2pi-periodic
! in zeta, so a target may sit in any field-period sector. Rotate it
! back by the whole field periods of its geometric angle: the rotated
! target has its root near the stored wedge, where the seed scan
! resolves, and the solved zeta shifts back exactly by the symmetry.
phi_shift = zeta_period*floor(modulo(atan2(x_target(2), x_target(1)), &
TWOPI)/zeta_period)
ca = cos(phi_shift)
sa = sin(phi_shift)
x_wedge(1) = ca*x_target(1) + sa*x_target(2)
x_wedge(2) = -sa*x_target(1) + ca*x_target(2)
x_wedge(3) = x_target(3)
call chartmap_initial_guess_cartesian_3d(self, x_wedge, u_try)
call chartmap_refine_cart_seed(self, x_wedge, zeta_period, u_try, &
res_norm, ierr_seed)
call chartmap_solve_cart(self, x_wedge, u_try(1), u_try(2), u_try(3), &
ierr_try, trace)
! Accept on the final residual, not on the solver exit path: the Newton
! step-size exit fires once steps reach tol_step in chart units, which
! on a reactor-size chart is a Cartesian residual of order
! |dx/du|*tol_step, well above tol_res yet far below the documented
! inverse tolerance. A target the chart cannot represent keeps a
! residual orders of magnitude above the accept tolerance and stays
! rejected. The threshold carries the solver's own convergence floors
! (tol_newton and the scale-aware term) so a root the Newton
! legitimately declares converged is never rejected here.
call chartmap_eval_cart(self, u_try, x_round)
res_norm = sqrt(sum((x_wedge - x_round)**2))
accept_tol = max(chartmap_invert_accept_tol, self%tol_newton, &
100.0_dp*epsilon(1.0_dp)*sqrt(sum(x_wedge**2)))
if (res_norm < accept_tol) ierr = chartmap_from_cyl_ok
if (ierr /= chartmap_from_cyl_ok) return
u = u_try
u(1) = min(max(u(1), 0.0_dp), 1.0_dp)
u(2) = modulo(u(2), TWOPI)
u(3) = modulo(u(3), zeta_period)
u(3) = modulo(u(3) + phi_shift, TWOPI)
return
end if

Expand Down Expand Up @@ -1208,7 +1237,6 @@ subroutine chartmap_solve_cart(self, x_target, rho, theta, zeta, ierr, trace)
real(dp) :: res_norm, res_norm_new
real(dp) :: tol_res, tol_step
real(dp) :: step_norm
real(dp) :: zeta_period
integer :: iter
logical :: success

Expand All @@ -1219,10 +1247,15 @@ subroutine chartmap_solve_cart(self, x_target, rho, theta, zeta, ierr, trace)
! for any orbit whose particle reaches rho -> 0. (X, Y) carries the same
! information with a Jacobian regular through rho = 0 (the standard near-axis
! chart, cf. flux_pseudocartesian and Pfefferle et al. arXiv:1412.5464).
! zeta iterates unwrapped: the lab-frame position of the co-rotating cart
! chart is only 2pi-periodic in zeta, so wrapping the trial into the stored
! wedge would teleport it by a field-period rotation. Callers reduce the
! solved zeta to their convention. The residual tolerance carries a
! scale-aware floor so reactor-size charts (|x| ~ 1e3 cm) are not asked for
! a residual below floating-point resolution.
ierr = chartmap_from_cyl_ok
zeta_period = TWOPI/real(self%num_field_periods, dp)
zeta = modulo(zeta, zeta_period)
tol_res = max(1.0e-10_dp, self%tol_newton)
tol_res = max(1.0e-10_dp, self%tol_newton, &
100.0_dp*epsilon(1.0_dp)*sqrt(sum(x_target**2)))
tol_step = max(1.0e-10_dp, self%tol_newton)
xw = rho*cos(theta)
yw = rho*sin(theta)
Expand All @@ -1236,9 +1269,9 @@ subroutine chartmap_solve_cart(self, x_target, rho, theta, zeta, ierr, trace)
ierr = chartmap_from_cyl_err_max_iter
exit
end if
call chartmap_line_search_xy(self, x_target, xw, yw, zeta, &
zeta_period, delta, res_norm, xw_new, &
yw_new, zeta_new, res_norm_new, success)
call chartmap_line_search_xy(self, x_target, xw, yw, zeta, delta, &
res_norm, xw_new, yw_new, zeta_new, &
res_norm_new, success)
if (.not. success) then
ierr = chartmap_from_cyl_err_max_iter
exit
Expand Down Expand Up @@ -1279,9 +1312,14 @@ subroutine chartmap_newton_delta_xy(self, x_target, xw, yw, zeta, delta, &
call chartmap_eval_residual_and_partials_cart(self, x_target, rho, theta, &
zeta, residual, res_norm, &
dx_drho, dx_dtheta, dx_dzeta)
rsafe = max(rho, 1.0e-30_dp)
cth = xw/rsafe
sth = yw/rsafe
! Take the gauge direction from theta, not from xw/rho: at an exact axis
! seed (xw = yw = 0) the quotient degenerates to cth = sth = 0, which zeroes
! both pseudo-Cartesian columns and freezes the iteration on the axis.
! atan2(0, 0) = 0 gives a well-defined gauge there. The radial floor caps
! the 1/rho amplification of axis-level spline noise in dx/dtheta.
rsafe = max(rho, 1.0e-8_dp)
cth = cos(theta)
sth = sin(theta)
dx_dxw = cth*dx_drho - (sth/rsafe)*dx_dtheta
dx_dyw = sth*dx_drho + (cth/rsafe)*dx_dtheta

Expand All @@ -1303,13 +1341,15 @@ subroutine chartmap_newton_delta_xy(self, x_target, xw, yw, zeta, delta, &
end subroutine chartmap_newton_delta_xy

! Backtracking line search in (X, Y, zeta). rho = sqrt(X^2 + Y^2) >= 0 is
! automatic, so only the outer bound needs a guard.
subroutine chartmap_line_search_xy(self, x_target, xw, yw, zeta, zeta_period, &
delta, res_norm, xw_new, yw_new, zeta_new, &
! automatic, so only the outer bound needs a guard. zeta moves unwrapped: the
! lab-frame position is only 2pi-periodic in zeta, so a wedge wrap would
! teleport the trial by a field-period rotation and stall the search.
subroutine chartmap_line_search_xy(self, x_target, xw, yw, zeta, delta, &
res_norm, xw_new, yw_new, zeta_new, &
res_norm_new, success)
class(chartmap_coordinate_system_t), intent(in) :: self
real(dp), intent(in) :: x_target(3)
real(dp), intent(in) :: xw, yw, zeta, zeta_period
real(dp), intent(in) :: xw, yw, zeta
real(dp), intent(in) :: delta(3), res_norm
real(dp), intent(out) :: xw_new, yw_new, zeta_new, res_norm_new
logical, intent(out) :: success
Expand All @@ -1329,7 +1369,7 @@ subroutine chartmap_line_search_xy(self, x_target, xw, yw, zeta, zeta_period, &
cycle
end if
theta_t = atan2(yw_new, xw_new)
zeta_new = modulo(zeta + alpha*delta(3), zeta_period)
zeta_new = zeta + alpha*delta(3)
uvec = [rho_t, theta_t, zeta_new]
call chartmap_eval_cart(self, uvec, x_trial)
res_norm_new = sqrt(sum((x_target - x_trial)**2))
Expand Down
74 changes: 74 additions & 0 deletions test/source/test_boozer_chartmap_roundtrip.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,10 @@ program test_boozer_chartmap_roundtrip
use new_vmec_stuff_mod, only: nper
use boozer_sub, only: get_boozer_coordinates, splint_boozer_coord
use boozer_chartmap, only: export_boozer_chartmap, load_boozer_from_chartmap
use libneo_coordinates, only: coordinate_system_t, &
make_chartmap_coordinate_system, &
chartmap_from_cyl_ok
use cylindrical_cartesian, only: cart_to_cyl
implicit none

real(dp), parameter :: PI = 3.14159265358979_dp
Expand All @@ -30,6 +34,10 @@ program test_boozer_chartmap_roundtrip
real(dp) :: phi_period, rel_err, max_err_bmod
integer :: i

integer, parameter :: n_coord_chk = 3
class(coordinate_system_t), allocatable :: coord_cs
real(dp) :: u_ref_chk(3, n_coord_chk), u_back_before(3, n_coord_chk)

use_B_r = .false.
call get_boozer_coordinates(wout_file, radial_spline_order=5, &
angular_spline_order=5, grid_refinment=3)
Expand All @@ -47,6 +55,7 @@ program test_boozer_chartmap_roundtrip
end do

call export_boozer_chartmap(chartmap_file)
call check_coordinate_inverse_before_load(chartmap_file)
call load_boozer_from_chartmap(chartmap_file)

do i = 1, n_test
Expand Down Expand Up @@ -74,6 +83,8 @@ program test_boozer_chartmap_roundtrip
end if
print *, "Boozer chartmap export/import preserves |B| within tolerance"

call check_coordinate_inverse_after_load()

contains

!> Interior test points: s in (0.1, 0.8), full theta period, scattered phi.
Expand All @@ -91,4 +102,67 @@ subroutine init_test_points(phi_per)
end do
end subroutine init_test_points

!> Issue #395 probe sequence on the exported production chartmap: an exact
!> forward image at interior rho (including rho = 0.5 and toroidal angles
!> outside the stored wedge) must be recovered by from_cyl. The result is
!> stored so the after-load check can pin invariance against the
!> module-level Boozer spline construction.
subroutine check_coordinate_inverse_before_load(path)
character(len=*), intent(in) :: path

real(dp) :: x(3), xcyl(3), u_back(3), x_back(3)
integer :: j, ierr

u_ref_chk(:, 1) = [0.5_dp, 1.9_dp, 0.35_dp*TWOPI]
u_ref_chk(:, 2) = [0.5_dp, 4.0_dp, 0.80_dp*TWOPI]
u_ref_chk(:, 3) = [0.3_dp, 0.7_dp, 0.55_dp*TWOPI]

call make_chartmap_coordinate_system(coord_cs, path)

do j = 1, n_coord_chk
call coord_cs%evaluate_cart(u_ref_chk(:, j), x)
call cart_to_cyl(x, xcyl)
call coord_cs%from_cyl(xcyl, u_back, ierr)
if (ierr /= chartmap_from_cyl_ok) then
print *, "coordinate inverse failed, ierr=", ierr, &
" at u_ref=", u_ref_chk(:, j)
error stop "chartmap coordinate inverse failed on exported chart"
end if
call coord_cs%evaluate_cart(u_back, x_back)
if (sqrt(sum((x_back - x)**2)) > 1.0e-6_dp) then
print *, "roundtrip residual=", sqrt(sum((x_back - x)**2)), &
" at u_ref=", u_ref_chk(:, j)
error stop "chartmap coordinate roundtrip residual too large"
end if
if (abs(u_back(1) - u_ref_chk(1, j)) > 1.0e-6_dp) then
print *, "recovered rho=", u_back(1), " want=", u_ref_chk(1, j)
error stop "chartmap coordinate inverse lost the radial label"
end if
u_back_before(:, j) = u_back
end do
print *, "Coordinate inverse recovers exported-chart points before load"
end subroutine check_coordinate_inverse_before_load

!> Building the module-level Boozer field splines must not change the
!> coordinate inverse of the standalone chartmap coordinate system.
subroutine check_coordinate_inverse_after_load()
real(dp) :: x(3), xcyl(3), u_back(3)
integer :: j, ierr

do j = 1, n_coord_chk
call coord_cs%evaluate_cart(u_ref_chk(:, j), x)
call cart_to_cyl(x, xcyl)
call coord_cs%from_cyl(xcyl, u_back, ierr)
if (ierr /= chartmap_from_cyl_ok) then
error stop "coordinate inverse regressed after Boozer spline load"
end if
if (any(u_back /= u_back_before(:, j))) then
print *, "u before=", u_back_before(:, j)
print *, "u after =", u_back
error stop "Boozer spline load changed the coordinate inverse"
end if
end do
print *, "Boozer spline load leaves the coordinate inverse unchanged"
end subroutine check_coordinate_inverse_after_load

end program test_boozer_chartmap_roundtrip
Loading
Loading