diff --git a/src/coordinates/libneo_coordinates_chartmap.f90 b/src/coordinates/libneo_coordinates_chartmap.f90 index f224160d..19ad72d2 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,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) @@ -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 @@ -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 @@ -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) @@ -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 @@ -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 @@ -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 @@ -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)) 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..e3b5f295 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,225 @@ 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] + ! 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 + 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-4_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..20b931e9 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. @@ -125,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] @@ -139,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