From e40e06cc51dff85c7b47b93736d06b7f55ba290c Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 14:58:25 +0200 Subject: [PATCH 1/3] Add failing regressions for Boozer chartmap cold inverse (#395) --- .../source/test_boozer_chartmap_roundtrip.f90 | 74 ++++++ .../test_chartmap_boozer_coordinates.f90 | 220 ++++++++++++++++++ test/source/test_chartmap_invert_cart.f90 | 37 +++ 3 files changed, 331 insertions(+) diff --git a/test/source/test_boozer_chartmap_roundtrip.f90 b/test/source/test_boozer_chartmap_roundtrip.f90 index aa7d4ba6..67d8df96 100644 --- a/test/source/test_boozer_chartmap_roundtrip.f90 +++ b/test/source/test_boozer_chartmap_roundtrip.f90 @@ -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 @@ -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) @@ -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 @@ -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. @@ -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 diff --git a/test/source/test_chartmap_boozer_coordinates.f90 b/test/source/test_chartmap_boozer_coordinates.f90 index 0f71d753..d8789327 100644 --- a/test/source/test_chartmap_boozer_coordinates.f90 +++ b/test/source/test_chartmap_boozer_coordinates.f90 @@ -8,6 +8,7 @@ program test_chartmap_boozer_coordinates implicit none character(len=*), parameter :: filename = "chartmap_boozer.nc" + character(len=*), parameter :: filename_reactor = "chartmap_boozer_reactor.nc" integer :: ierr, nerrors character(len=2048) :: message @@ -22,6 +23,9 @@ program test_chartmap_boozer_coordinates call run_roundtrip_checks(filename, nerrors) + call write_reactor_boozer_file(filename_reactor) + call run_full_torus_cold_inverse_checks(filename_reactor, nerrors) + if (nerrors > 0) then print *, "FAILED: ", nerrors, " error(s) detected in Boozer chartmap tests" error stop 1 @@ -135,6 +139,222 @@ pure real(dp) function angular_diff(a, b, period) angular_diff = abs(modulo(a - b + 0.5_dp*period, period) - 0.5_dp*period) end function angular_diff + ! Regression for issue #395: the cold inverse of a rotating-frame Boozer chart + ! must recover interior points whose toroidal angle lies anywhere on the full + ! torus, not only inside the stored field-period wedge. The lab-frame position + ! of the co-rotating chart is 2pi-periodic in zeta, so zeta is compared modulo + ! TWOPI and the recovered u must reproduce the Cartesian target. + subroutine run_full_torus_cold_inverse_checks(path, nerrors) + character(len=*), intent(in) :: path + integer, intent(inout) :: nerrors + + class(coordinate_system_t), allocatable :: cs + + call make_chartmap_coordinate_system(cs, path) + + select type (ccs => cs) + type is (chartmap_coordinate_system_t) + call check_cold_interior_points(ccs, nerrors) + call check_cold_outside_target(ccs, nerrors) + call check_cold_near_axis_target(ccs, nerrors) + class default + print *, " FAIL: reactor chartmap did not load as chartmap type" + nerrors = nerrors + 1 + end select + end subroutine run_full_torus_cold_inverse_checks + + subroutine check_cold_interior_points(ccs, nerrors) + type(chartmap_coordinate_system_t), intent(in) :: ccs + integer, intent(inout) :: nerrors + + real(dp) :: u(3), u_back(3), x(3), xcyl(3), x_back(3) + real(dp) :: du_theta, du_zeta + integer :: ir, it, iz, ierr_local, npoints, nfail + real(dp), parameter :: rho_vals(4) = [0.12_dp, 0.45_dp, 0.5_dp, 0.82_dp] + real(dp), parameter :: theta_vals(5) = [0.0_dp, 1.4_dp, 3.0_dp, 4.9_dp, & + 6.2_dp] + real(dp), parameter :: zeta_fracs(6) = [0.0_dp, 0.13_dp, 0.29_dp, 0.48_dp, & + 0.77_dp, 0.999_dp] + real(dp), parameter :: tol_u = 1.0e-8_dp + real(dp), parameter :: tol_x = 1.0e-8_dp + + npoints = 0 + nfail = 0 + do iz = 1, size(zeta_fracs) + do it = 1, size(theta_vals) + do ir = 1, size(rho_vals) + u = [rho_vals(ir), theta_vals(it), zeta_fracs(iz)*TWOPI] + ! Every third point enters many field periods outside the + ! stored wedge; the inverse must recover zeta modulo TWOPI. + if (mod(npoints, 3) == 0) u(3) = u(3) + 3.0_dp*TWOPI + npoints = npoints + 1 + + call ccs%evaluate_cart(u, x) + xcyl(1) = sqrt(x(1)**2 + x(2)**2) + xcyl(2) = atan2(x(2), x(1)) + xcyl(3) = x(3) + + call ccs%from_cyl(xcyl, u_back, ierr_local) + if (ierr_local /= chartmap_from_cyl_ok) then + print *, " FAIL: cold from_cyl ierr=", ierr_local, & + " at u=", u + nfail = nfail + 1 + cycle + end if + + du_theta = angular_diff(u_back(2), u(2), TWOPI) + du_zeta = angular_diff(u_back(3), u(3), TWOPI) + call ccs%evaluate_cart(u_back, x_back) + if (abs(u_back(1) - u(1)) > tol_u .or. du_theta > tol_u .or. & + du_zeta > tol_u .or. & + sqrt(sum((x_back - x)**2)) > tol_x) then + print *, " FAIL: cold from_cyl roundtrip mismatch" + print *, " u =", u + print *, " u_back=", u_back + print *, " |dx| =", sqrt(sum((x_back - x)**2)) + nfail = nfail + 1 + end if + end do + end do + end do + + nerrors = nerrors + nfail + if (nfail == 0) then + print *, " PASS: ", npoints, " cold full-torus roundtrips recovered" + end if + end subroutine check_cold_interior_points + + ! A target genuinely outside the mapped volume must return a nonzero status + ! instead of masquerading as an interior root. + subroutine check_cold_outside_target(ccs, nerrors) + type(chartmap_coordinate_system_t), intent(in) :: ccs + integer, intent(inout) :: nerrors + + real(dp) :: u_edge(3), x_edge(3), x(3), u_back(3) + real(dp) :: e_cov(3, 3) + integer :: ierr_local + + u_edge = [1.0_dp, 0.9_dp, 0.7_dp*TWOPI] + call ccs%evaluate_cart(u_edge, x_edge) + call ccs%covariant_basis(u_edge, e_cov) + x = x_edge + 0.5_dp*e_cov(:, 1) + + call ccs%from_cart(x, u_back, ierr_local) + if (ierr_local == chartmap_from_cyl_ok) then + print *, " FAIL: past-edge target accepted as interior root, u=", u_back + nerrors = nerrors + 1 + else + print *, " PASS: past-edge target returns nonzero status" + end if + end subroutine check_cold_outside_target + + ! A near-axis target below the chart resolution must either locate with a tiny + ! radius that reproduces the target, or fail explicitly; it must never return + ! success on an unrelated radial branch. + subroutine check_cold_near_axis_target(ccs, nerrors) + type(chartmap_coordinate_system_t), intent(in) :: ccs + integer, intent(inout) :: nerrors + + real(dp) :: u(3), u_back(3), x(3), x_back(3) + integer :: ierr_local + + u = [1.0e-7_dp, 2.3_dp, 0.6_dp*TWOPI] + call ccs%evaluate_cart(u, x) + call ccs%from_cart(x, u_back, ierr_local) + if (ierr_local == chartmap_from_cyl_ok) then + call ccs%evaluate_cart(u_back, x_back) + if (u_back(1) > 1.0e-4_dp .or. sqrt(sum((x_back - x)**2)) > 1.0e-6_dp) then + print *, " FAIL: near-axis target accepted on wrong branch" + print *, " u_back=", u_back + nerrors = nerrors + 1 + return + end if + end if + print *, " PASS: near-axis target not accepted on an unrelated branch" + end subroutine check_cold_near_axis_target + + ! Reactor-scale endpoint-excluded rotating-frame Boozer chart, nfp = 5, with a + ! field-period-symmetric Boozer toroidal shift, sized so the inverse exercises + ! the production probe/refine/solve path at reactor dimensions. + subroutine write_reactor_boozer_file(path) + character(len=*), intent(in) :: path + + integer :: ncid + integer :: dim_rho, dim_theta, dim_zeta + integer :: var_rho, var_theta, var_zeta, var_x, var_y, var_z, var_nfp + integer, parameter :: nrho = 25, ntheta = 48, nzeta = 32, nfp = 5 + real(dp), parameter :: r0 = 1000.0_dp + real(dp), parameter :: a = 250.0_dp + real(dp) :: rho(nrho), theta(ntheta), zeta(nzeta) + real(dp) :: x(nrho, ntheta, nzeta), y(nrho, ntheta, nzeta), & + z(nrho, ntheta, nzeta) + real(dp) :: rho_val, theta_val, zeta_val + real(dp) :: phi_geom, rmaj, zpos, period + integer :: ir, it, iz + + period = TWOPI/real(nfp, dp) + + do ir = 1, nrho + rho(ir) = real(ir - 1, dp)/real(nrho - 1, dp) + end do + do it = 1, ntheta + theta(it) = TWOPI*real(it - 1, dp)/real(ntheta, dp) + end do + do iz = 1, nzeta + zeta(iz) = period*real(iz - 1, dp)/real(nzeta, dp) + end do + + do iz = 1, nzeta + zeta_val = zeta(iz) + do it = 1, ntheta + theta_val = theta(it) + do ir = 1, nrho + rho_val = rho(ir) + phi_geom = zeta_val + 0.12_dp*rho_val*sin(theta_val) + & + 0.03_dp*sin(real(nfp, dp)*zeta_val) + rmaj = r0 + a*rho_val*cos(theta_val) + & + 0.05_dp*a*rho_val*cos(theta_val - real(nfp, dp)*zeta_val) + zpos = a*rho_val*sin(theta_val) + & + 0.03_dp*a*rho_val*sin(real(nfp, dp)*zeta_val) + + x(ir, it, iz) = rmaj*cos(phi_geom) + y(ir, it, iz) = rmaj*sin(phi_geom) + z(ir, it, iz) = zpos + end do + end do + end do + + call nc_check(nf90_create(trim(path), NF90_NETCDF4, ncid)) + call nc_check(nf90_put_att(ncid, NF90_GLOBAL, "zeta_convention", "boozer")) + call nc_check(nf90_put_att(ncid, NF90_GLOBAL, "rho_convention", "rho_tor")) + call nc_check(nf90_def_dim(ncid, "rho", nrho, dim_rho)) + call nc_check(nf90_def_dim(ncid, "theta", ntheta, dim_theta)) + call nc_check(nf90_def_dim(ncid, "zeta", nzeta, dim_zeta)) + call nc_check(nf90_def_var(ncid, "rho", NF90_DOUBLE, [dim_rho], var_rho)) + call nc_check(nf90_def_var(ncid, "theta", NF90_DOUBLE, [dim_theta], var_theta)) + call nc_check(nf90_def_var(ncid, "zeta", NF90_DOUBLE, [dim_zeta], var_zeta)) + call nc_check(nf90_def_var(ncid, "x", NF90_DOUBLE, & + [dim_rho, dim_theta, dim_zeta], var_x)) + call nc_check(nf90_def_var(ncid, "y", NF90_DOUBLE, & + [dim_rho, dim_theta, dim_zeta], var_y)) + call nc_check(nf90_def_var(ncid, "z", NF90_DOUBLE, & + [dim_rho, dim_theta, dim_zeta], var_z)) + call nc_check(nf90_def_var(ncid, "num_field_periods", NF90_INT, var_nfp)) + call nc_check(nf90_put_att(ncid, var_x, "units", "cm")) + call nc_check(nf90_put_att(ncid, var_y, "units", "cm")) + call nc_check(nf90_put_att(ncid, var_z, "units", "cm")) + call nc_check(nf90_enddef(ncid)) + + call nc_check(nf90_put_var(ncid, var_rho, rho)) + call nc_check(nf90_put_var(ncid, var_theta, theta)) + call nc_check(nf90_put_var(ncid, var_zeta, zeta)) + call nc_check(nf90_put_var(ncid, var_x, x)) + call nc_check(nf90_put_var(ncid, var_y, y)) + call nc_check(nf90_put_var(ncid, var_z, z)) + call nc_check(nf90_put_var(ncid, var_nfp, nfp)) + call nc_check(nf90_close(ncid)) + end subroutine write_reactor_boozer_file + subroutine write_boozer_file(path) character(len=*), intent(in) :: path diff --git a/test/source/test_chartmap_invert_cart.f90 b/test/source/test_chartmap_invert_cart.f90 index b3d00ade..a78cf1b3 100644 --- a/test/source/test_chartmap_invert_cart.f90 +++ b/test/source/test_chartmap_invert_cart.f90 @@ -68,6 +68,7 @@ subroutine run_invert_checks(path, nerrors) zeta_period = TWOPI/real(ccs%num_field_periods, dp) call check_interior_roundtrip(ccs, zeta_period, nerrors) + call check_full_torus_roundtrip(ccs, nerrors) call check_seam_roundtrip(ccs, zeta_period, nerrors) call check_near_axis(ccs, zeta_period, nerrors) call check_axis(ccs, zeta_period, nerrors) @@ -118,6 +119,42 @@ subroutine check_interior_roundtrip(ccs, zeta_period, nerrors) if (nfail == 0) print *, " PASS: interior roundtrip recovers u0, LOCATED" end subroutine check_interior_roundtrip + ! Full-torus points (issue #395): the lab-frame position of the co-rotating + ! Boozer chart is 2pi-periodic in zeta, not field-period-periodic, so the + ! inverse must recover zeta modulo TWOPI and reproduce the Cartesian target + ! for angles in every field-period sector. + subroutine check_full_torus_roundtrip(ccs, nerrors) + type(chartmap_coordinate_system_t), intent(in), target :: ccs + integer, intent(inout) :: nerrors + real(dp) :: u0(3), x(3), u(3), u_guess(3), x_back(3) + real(dp) :: dth, dze + integer :: k, status, nfail + real(dp), parameter :: zeta_fracs(4) = [0.15_dp, 0.4_dp, 0.65_dp, 0.9_dp] + real(dp), parameter :: tol = 1.0e-6_dp + + nfail = 0 + do k = 1, size(zeta_fracs) + u0 = [0.52_dp, 2.6_dp, zeta_fracs(k)*TWOPI] + call ccs%evaluate_cart(u0, x) + u_guess = [u0(1) + 0.07_dp, u0(2) + 0.15_dp, u0(3) - 0.05_dp] + call ccs%invert_cart(x, u_guess, u, status) + + dth = angular_diff(u(2), u0(2), TWOPI) + dze = angular_diff(u(3), u0(3), TWOPI) + call ccs%evaluate_cart(u, x_back) + if (status /= CHARTMAP_LOCATED .or. abs(u(1) - u0(1)) > tol .or. & + dth > tol .or. dze > tol .or. & + sqrt(sum((x_back - x)**2)) > tol) then + print *, " FAIL: full-torus roundtrip u0=", u0 + print *, " u=", u, " status=", status + print *, " |dx|=", sqrt(sum((x_back - x)**2)) + nfail = nfail + 1 + end if + end do + nerrors = nerrors + nfail + if (nfail == 0) print *, " PASS: full-torus roundtrip recovers zeta mod 2pi" + end subroutine check_full_torus_roundtrip + ! Seam points: zeta within a small fraction of a field-period boundary, where a ! global Boozer phi guess would be a full period stale. invert_cart reseeds the ! toroidal angle and must still locate. From a00952ba8cc81cee5afd5213cc2af2244c3d01d2 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 15:13:38 +0200 Subject: [PATCH 2/3] Solve Boozer chartmap cold inverse on the full torus The cold Cartesian inverse of the co-rotating cart chart constrained the logical toroidal angle to the stored field-period wedge (seed scan, line search wrap, and final reduction), but the lab-frame position of that chart is only 2pi-periodic in zeta, so any target outside the first field-period sector could never converge and an exact forward-map roundtrip failed with err_max_iter. Rotate the target back by the whole field periods of its geometric angle (an exact chart symmetry), iterate zeta unwrapped, and shift the solved zeta back, returning it modulo 2pi. Accept the cold root on its final Cartesian residual against the documented accept tolerance instead of the solver exit path, whose step-size exit sits above the reachable residual on reactor-size charts. Give the Newton residual tolerance a scale-aware floor, fix the axis-degenerate gauge in the pseudo-Cartesian step, and stop wedge-wrapping the warm inverse result for co-rotating charts. Fixes #395 --- .../libneo_coordinates_chartmap.f90 | 99 +++++++++++++------ .../test_chartmap_boozer_coordinates.f90 | 9 +- 2 files changed, 73 insertions(+), 35 deletions(-) diff --git a/src/coordinates/libneo_coordinates_chartmap.f90 b/src/coordinates/libneo_coordinates_chartmap.f90 index f224160d..8a83bc80 100644 --- a/src/coordinates/libneo_coordinates_chartmap.f90 +++ b/src/coordinates/libneo_coordinates_chartmap.f90 @@ -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) @@ -1105,10 +1112,12 @@ 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) :: 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) @@ -1117,25 +1126,40 @@ 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. + call chartmap_eval_cart(self, u_try, x_round) + res_norm = sqrt(sum((x_wedge - x_round)**2)) + if (res_norm < chartmap_invert_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 @@ -1208,7 +1232,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 @@ -1219,10 +1242,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) @@ -1236,9 +1264,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 @@ -1279,9 +1307,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 @@ -1303,13 +1336,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 @@ -1329,7 +1364,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)) diff --git a/test/source/test_chartmap_boozer_coordinates.f90 b/test/source/test_chartmap_boozer_coordinates.f90 index d8789327..e3b5f295 100644 --- a/test/source/test_chartmap_boozer_coordinates.f90 +++ b/test/source/test_chartmap_boozer_coordinates.f90 @@ -175,8 +175,11 @@ subroutine check_cold_interior_points(ccs, nerrors) 6.2_dp] real(dp), parameter :: zeta_fracs(6) = [0.0_dp, 0.13_dp, 0.29_dp, 0.48_dp, & 0.77_dp, 0.999_dp] - real(dp), parameter :: tol_u = 1.0e-8_dp - real(dp), parameter :: tol_x = 1.0e-8_dp + ! The documented inverse tolerance: chartmap_invert_accept_tol on the + ! Cartesian residual, with the coordinate error bounded by residual over + ! the smallest local basis-vector scale. + real(dp), parameter :: tol_u = 1.0e-6_dp + real(dp), parameter :: tol_x = 1.0e-6_dp npoints = 0 nfail = 0 @@ -263,7 +266,7 @@ subroutine check_cold_near_axis_target(ccs, nerrors) call ccs%from_cart(x, u_back, ierr_local) if (ierr_local == chartmap_from_cyl_ok) then call ccs%evaluate_cart(u_back, x_back) - if (u_back(1) > 1.0e-4_dp .or. sqrt(sum((x_back - x)**2)) > 1.0e-6_dp) then + if (u_back(1) > 1.0e-4_dp .or. sqrt(sum((x_back - x)**2)) > 1.0e-4_dp) then print *, " FAIL: near-axis target accepted on wrong branch" print *, " u_back=", u_back nerrors = nerrors + 1 From 2a4e81a040e3c700de2cfa8913044d7dbe0be55b Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 15:25:26 +0200 Subject: [PATCH 3/3] Align cold accept threshold with solver floors, harden seam test Review findings: the fixed cold accept tolerance could sit below the Newton's own convergence target for coarse tol_newton or very large charts, rejecting a legitimately converged root; carry the same floors in the acceptance. The warm seam test compared zeta modulo one field period and never re-evaluated the Cartesian point, masking exactly the wrong-wedge regression it guards; compare modulo 2pi and assert the reconstructed position. --- src/coordinates/libneo_coordinates_chartmap.f90 | 9 +++++++-- test/source/test_chartmap_invert_cart.f90 | 10 +++++++--- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/src/coordinates/libneo_coordinates_chartmap.f90 b/src/coordinates/libneo_coordinates_chartmap.f90 index 8a83bc80..19ad72d2 100644 --- a/src/coordinates/libneo_coordinates_chartmap.f90 +++ b/src/coordinates/libneo_coordinates_chartmap.f90 @@ -1114,6 +1114,7 @@ subroutine chartmap_invert_cart(self, x_target, u, ierr) 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 @@ -1151,10 +1152,14 @@ subroutine chartmap_invert_cart(self, x_target, u, ierr) ! |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. + ! 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)) - if (res_norm < chartmap_invert_accept_tol) ierr = chartmap_from_cyl_ok + 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) diff --git a/test/source/test_chartmap_invert_cart.f90 b/test/source/test_chartmap_invert_cart.f90 index a78cf1b3..20b931e9 100644 --- a/test/source/test_chartmap_invert_cart.f90 +++ b/test/source/test_chartmap_invert_cart.f90 @@ -162,7 +162,7 @@ subroutine check_seam_roundtrip(ccs, zeta_period, nerrors) type(chartmap_coordinate_system_t), intent(in), target :: ccs real(dp), intent(in) :: zeta_period integer, intent(inout) :: nerrors - real(dp) :: u0(3), x(3), u(3), u_guess(3) + real(dp) :: u0(3), x(3), u(3), u_guess(3), x_back(3) real(dp) :: dth, dze integer :: k, status, nfail real(dp), parameter :: seam_fracs(3) = [0.0_dp, 0.999_dp, 0.5_dp] @@ -176,9 +176,13 @@ subroutine check_seam_roundtrip(ccs, zeta_period, nerrors) u_guess = [u0(1), u0(2), u0(3) + zeta_period] call ccs%invert_cart(x, u_guess, u, status) dth = angular_diff(u(2), u0(2), TWOPI) - dze = angular_diff(u(3), u0(3), zeta_period) + ! zeta compares modulo 2pi: a result off by a field period would pass + ! a wedge-modular comparison yet reconstruct a rotated position. + dze = angular_diff(u(3), u0(3), TWOPI) + call ccs%evaluate_cart(u, x_back) if (status /= CHARTMAP_LOCATED .or. abs(u(1) - u0(1)) > tol .or. & - dth > tol .or. dze > tol) then + dth > tol .or. dze > tol .or. & + sqrt(sum((x_back - x)**2)) > tol) then print *, " FAIL: seam roundtrip u0=", u0, " u=", u, " status=", status nfail = nfail + 1 end if