From 0640623edfbe7fbb8ee89a0ccf7b05609aacb4c4 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 28 Mar 2026 08:46:09 +0100 Subject: [PATCH] Support real-Cartesian Boozer chartmaps --- src/coordinates/libneo_coordinates_base.f90 | 7 +- .../libneo_coordinates_chartmap.f90 | 169 +++++++++++++- .../libneo_coordinates_chartmap_impl.inc | 169 +++++++++++++- .../libneo_coordinates_mapping.f90 | 17 +- .../libneo_coordinates_validator.f90 | 4 +- test/CMakeLists.txt | 9 + .../test_chartmap_boozer_coordinates.f90 | 216 ++++++++++++++++++ test/source/test_chartmap_validator.f90 | 87 +++++++ 8 files changed, 658 insertions(+), 20 deletions(-) create mode 100644 test/source/test_chartmap_boozer_coordinates.f90 diff --git a/src/coordinates/libneo_coordinates_base.f90 b/src/coordinates/libneo_coordinates_base.f90 index 9d6bb1ea..0a5f9f62 100644 --- a/src/coordinates/libneo_coordinates_base.f90 +++ b/src/coordinates/libneo_coordinates_base.f90 @@ -25,7 +25,7 @@ module libneo_coordinates_base " z(zeta,theta,rho)"//new_line('a')// & "- Optional: num_field_periods (integer >= 1)"//new_line('a')// & "- Required: zeta_convention (global attribute)"//new_line('a')// & - " allowed: cyl, vmec"//new_line('a')// & + " allowed: cyl, vmec, boozer"//new_line('a')// & "- Ranges:"//new_line('a')// & " rho in [0,1]"//new_line('a')// & " theta in [0,2pi)"//new_line('a')// & @@ -34,7 +34,10 @@ module libneo_coordinates_base "- Storage order:"//new_line('a')// & " file dims (zeta,theta,rho)"//new_line('a')// & "Fortran reads x(rho,theta,zeta)"//new_line('a')// & - "- Units: x,y,z in cm" + "- Units: x,y,z in cm"//new_line('a')// & + "- Contract:"//new_line('a')// & + " cyl/vmec require atan2(y,x)=zeta"//new_line('a')// & + " boozer stores true Cartesian geometry; zeta is a chart parameter" integer, parameter :: chartmap_from_cyl_ok = 0 integer, parameter :: chartmap_from_cyl_err_max_iter = 1 diff --git a/src/coordinates/libneo_coordinates_chartmap.f90 b/src/coordinates/libneo_coordinates_chartmap.f90 index d329f74e..7902b18f 100644 --- a/src/coordinates/libneo_coordinates_chartmap.f90 +++ b/src/coordinates/libneo_coordinates_chartmap.f90 @@ -440,9 +440,10 @@ subroutine initialize_chartmap(ccs, filename) ccs%rho_convention = rho_convention zeta_period = TWOPI/real(num_field_periods, dp) - if (zeta_convention /= CYL .and. zeta_convention /= VMEC) then + if (zeta_convention /= CYL .and. zeta_convention /= VMEC .and. & + zeta_convention /= BOOZER) then print *, "initialize_chartmap: unsupported zeta_convention" - print *, " must be cyl or vmec" + print *, " must be cyl, vmec, or boozer" error stop end if @@ -485,11 +486,9 @@ subroutine initialize_chartmap(ccs, filename) call chartmap_extend_theta(nrho, ntheta, nzeta, theta, x, y, z, theta_spl, & x_th, y_th, z_th) - allocate (zeta_spl(nzeta)) - zeta_spl = zeta - call move_alloc(x_th, x_spl) - call move_alloc(y_th, y_spl) - call move_alloc(z_th, z_spl) + periodic_cart(3) = .true. + call chartmap_extend_zeta(nrho, ntheta + 1, nzeta, zeta, zeta_period, & + x_th, y_th, z_th, zeta_spl, x_spl, y_spl, z_spl) x_min = [rho(1), theta_spl(1), zeta_spl(1)] x_max = [rho(nrho), theta_spl(size(theta_spl)), zeta_spl(size(zeta_spl))] @@ -642,6 +641,30 @@ subroutine chartmap_invert_cart(self, x_target, u, ierr) u = 0.0_dp 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 + 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) + return + end if + phi = atan2(x_target(2), x_target(1)) zeta_seed = [modulo(phi, zeta_period), 0.0_dp, 0.25_dp*zeta_period, & 0.50_dp*zeta_period, 0.75_dp*zeta_period] @@ -968,6 +991,11 @@ subroutine chartmap_initial_guess(self, x_target, zeta, rho, theta, ierr) ierr = chartmap_from_cyl_ok + if (.not. self%has_spl_rz) then + call chartmap_initial_guess_cartesian(self, x_target, zeta, rho, theta) + return + end if + uvec = [0.0_dp, 0.0_dp, zeta] call chartmap_eval_cart(self, uvec, vals) axis_x = vals @@ -993,6 +1021,133 @@ subroutine chartmap_initial_guess(self, x_target, zeta, rho, theta, ierr) rho = min(max(rho, 0.0_dp), 1.0_dp) end subroutine chartmap_initial_guess + subroutine chartmap_initial_guess_cartesian(self, x_target, zeta, rho, theta) + class(chartmap_coordinate_system_t), intent(in) :: self + real(dp), intent(in) :: x_target(3) + real(dp), intent(in) :: zeta + real(dp), intent(out) :: rho + real(dp), intent(out) :: theta + + real(dp) :: uvec(3) + real(dp) :: x_try(3) + real(dp) :: best_dist2, dist2 + integer :: nrho_probe, ntheta_probe + integer :: ir, it + + nrho_probe = min(max(self%nrho, 5), 9) + ntheta_probe = min(max(self%ntheta, 16), 32) + + best_dist2 = huge(1.0_dp) + rho = 0.0_dp + theta = 0.0_dp + + do it = 1, ntheta_probe + uvec(2) = TWOPI*real(it - 1, dp)/real(ntheta_probe, dp) + uvec(3) = zeta + do ir = 1, nrho_probe + uvec(1) = real(ir - 1, dp)/real(nrho_probe - 1, dp) + call chartmap_eval_cart(self, uvec, x_try) + dist2 = sum((x_target - x_try)**2) + if (dist2 < best_dist2) then + best_dist2 = dist2 + rho = uvec(1) + theta = uvec(2) + end if + end do + end do + end subroutine chartmap_initial_guess_cartesian + + subroutine chartmap_initial_guess_cartesian_3d(self, x_target, u) + class(chartmap_coordinate_system_t), intent(in) :: self + real(dp), intent(in) :: x_target(3) + real(dp), intent(out) :: u(3) + + real(dp) :: uvec(3) + real(dp) :: x_try(3) + real(dp) :: best_dist2, dist2 + real(dp) :: zeta_period + integer :: nrho_probe, ntheta_probe, nzeta_probe + integer :: ir, it, iz + + nrho_probe = min(max(self%nrho, 9), 17) + ntheta_probe = min(max(self%ntheta, 32), 64) + nzeta_probe = min(max(self%nzeta, 24), 48) + zeta_period = TWOPI/real(self%num_field_periods, dp) + + best_dist2 = huge(1.0_dp) + u = 0.0_dp + + do iz = 1, nzeta_probe + uvec(3) = zeta_period*real(iz - 1, dp)/real(nzeta_probe, dp) + do it = 1, ntheta_probe + uvec(2) = TWOPI*real(it - 1, dp)/real(ntheta_probe, dp) + do ir = 1, nrho_probe + uvec(1) = real(ir - 1, dp)/real(nrho_probe - 1, dp) + call chartmap_eval_cart(self, uvec, x_try) + dist2 = sum((x_target - x_try)**2) + if (dist2 < best_dist2) then + best_dist2 = dist2 + u = uvec + end if + end do + end do + end do + end subroutine chartmap_initial_guess_cartesian_3d + + subroutine chartmap_refine_cart_seed(self, x_target, zeta_period, u, res_norm, ierr) + class(chartmap_coordinate_system_t), intent(in) :: self + real(dp), intent(in) :: x_target(3) + real(dp), intent(in) :: zeta_period + real(dp), intent(inout) :: u(3) + real(dp), intent(out) :: res_norm + integer, intent(out) :: ierr + + real(dp) :: trial(3), best_u(3), x_trial(3) + real(dp) :: step_rho, step_theta, step_zeta + real(dp) :: best_dist2, dist2 + integer :: iter, ir, it, iz + + ierr = chartmap_from_cyl_ok + step_rho = 1.0_dp/16.0_dp + step_theta = TWOPI/64.0_dp + step_zeta = zeta_period/48.0_dp + + call chartmap_eval_cart(self, u, x_trial) + best_dist2 = sum((x_target - x_trial)**2) + + do iter = 1, 64 + best_u = u + do iz = -1, 1 + trial(3) = modulo(u(3) + real(iz, dp)*step_zeta, zeta_period) + do it = -1, 1 + trial(2) = modulo(u(2) + real(it, dp)*step_theta, TWOPI) + do ir = -1, 1 + trial(1) = min(max(u(1) + real(ir, dp)*step_rho, 0.0_dp), 1.0_dp) + call chartmap_eval_cart(self, trial, x_trial) + dist2 = sum((x_target - x_trial)**2) + if (dist2 < best_dist2) then + best_dist2 = dist2 + best_u = trial + end if + end do + end do + end do + + if (any(best_u /= u)) then + u = best_u + else + step_rho = 0.5_dp*step_rho + step_theta = 0.5_dp*step_theta + step_zeta = 0.5_dp*step_zeta + if (max(step_rho, max(step_theta, step_zeta)) < 1.0e-12_dp) exit + end if + end do + + call chartmap_eval_cart(self, u, x_trial) + res_norm = sqrt(sum((x_target - x_trial)**2)) + if (res_norm > 1.0e-4_dp) ierr = chartmap_from_cyl_err_max_iter + end subroutine chartmap_refine_cart_seed + subroutine chartmap_newton_delta(self, x_target, rho, theta, zeta, delta, & res_norm, & ierr) diff --git a/src/coordinates/libneo_coordinates_chartmap_impl.inc b/src/coordinates/libneo_coordinates_chartmap_impl.inc index e2d8353e..ff446221 100644 --- a/src/coordinates/libneo_coordinates_chartmap_impl.inc +++ b/src/coordinates/libneo_coordinates_chartmap_impl.inc @@ -403,9 +403,10 @@ subroutine initialize_chartmap(ccs, filename) ccs%rho_convention = rho_convention zeta_period = TWOPI/real(num_field_periods, dp) - if (zeta_convention /= CYL .and. zeta_convention /= VMEC) then + if (zeta_convention /= CYL .and. zeta_convention /= VMEC .and. & + zeta_convention /= BOOZER) then print *, "initialize_chartmap: unsupported zeta_convention" - print *, " must be cyl or vmec" + print *, " must be cyl, vmec, or boozer" error stop end if @@ -448,11 +449,9 @@ subroutine initialize_chartmap(ccs, filename) call chartmap_extend_theta(nrho, ntheta, nzeta, theta, x, y, z, theta_spl, & x_th, y_th, z_th) - allocate (zeta_spl(nzeta)) - zeta_spl = zeta - call move_alloc(x_th, x_spl) - call move_alloc(y_th, y_spl) - call move_alloc(z_th, z_spl) + periodic_cart(3) = .true. + call chartmap_extend_zeta(nrho, ntheta + 1, nzeta, zeta, zeta_period, & + x_th, y_th, z_th, zeta_spl, x_spl, y_spl, z_spl) x_min = [rho(1), theta_spl(1), zeta_spl(1)] x_max = [rho(nrho), theta_spl(size(theta_spl)), zeta_spl(size(zeta_spl))] @@ -608,6 +607,30 @@ subroutine chartmap_invert_cart(self, x_target, u, ierr) u = 0.0_dp 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 + 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) + return + end if + phi = atan2(x_target(2), x_target(1)) zeta_seed = [modulo(phi, zeta_period), 0.0_dp, 0.25_dp*zeta_period, & 0.50_dp*zeta_period, 0.75_dp*zeta_period] @@ -936,6 +959,11 @@ subroutine chartmap_initial_guess(self, x_target, zeta, rho, theta, ierr) ierr = chartmap_from_cyl_ok + if (.not. self%has_spl_rz) then + call chartmap_initial_guess_cartesian(self, x_target, zeta, rho, theta) + return + end if + uvec = [0.0_dp, 0.0_dp, zeta] call chartmap_eval_cart(self, uvec, vals) axis_x = vals @@ -961,6 +989,133 @@ subroutine chartmap_initial_guess(self, x_target, zeta, rho, theta, ierr) rho = min(max(rho, 0.0_dp), 1.0_dp) end subroutine chartmap_initial_guess +subroutine chartmap_initial_guess_cartesian(self, x_target, zeta, rho, theta) + class(chartmap_coordinate_system_t), intent(in) :: self + real(dp), intent(in) :: x_target(3) + real(dp), intent(in) :: zeta + real(dp), intent(out) :: rho + real(dp), intent(out) :: theta + + real(dp) :: uvec(3) + real(dp) :: x_try(3) + real(dp) :: best_dist2, dist2 + integer :: nrho_probe, ntheta_probe + integer :: ir, it + + nrho_probe = min(max(self%nrho, 5), 9) + ntheta_probe = min(max(self%ntheta, 16), 32) + + best_dist2 = huge(1.0_dp) + rho = 0.0_dp + theta = 0.0_dp + + do it = 1, ntheta_probe + uvec(2) = TWOPI*real(it - 1, dp)/real(ntheta_probe, dp) + uvec(3) = zeta + do ir = 1, nrho_probe + uvec(1) = real(ir - 1, dp)/real(nrho_probe - 1, dp) + call chartmap_eval_cart(self, uvec, x_try) + dist2 = sum((x_target - x_try)**2) + if (dist2 < best_dist2) then + best_dist2 = dist2 + rho = uvec(1) + theta = uvec(2) + end if + end do + end do +end subroutine chartmap_initial_guess_cartesian + +subroutine chartmap_initial_guess_cartesian_3d(self, x_target, u) + class(chartmap_coordinate_system_t), intent(in) :: self + real(dp), intent(in) :: x_target(3) + real(dp), intent(out) :: u(3) + + real(dp) :: uvec(3) + real(dp) :: x_try(3) + real(dp) :: best_dist2, dist2 + real(dp) :: zeta_period + integer :: nrho_probe, ntheta_probe, nzeta_probe + integer :: ir, it, iz + + nrho_probe = min(max(self%nrho, 9), 17) + ntheta_probe = min(max(self%ntheta, 32), 64) + nzeta_probe = min(max(self%nzeta, 24), 48) + zeta_period = TWOPI/real(self%num_field_periods, dp) + + best_dist2 = huge(1.0_dp) + u = 0.0_dp + + do iz = 1, nzeta_probe + uvec(3) = zeta_period*real(iz - 1, dp)/real(nzeta_probe, dp) + do it = 1, ntheta_probe + uvec(2) = TWOPI*real(it - 1, dp)/real(ntheta_probe, dp) + do ir = 1, nrho_probe + uvec(1) = real(ir - 1, dp)/real(nrho_probe - 1, dp) + call chartmap_eval_cart(self, uvec, x_try) + dist2 = sum((x_target - x_try)**2) + if (dist2 < best_dist2) then + best_dist2 = dist2 + u = uvec + end if + end do + end do + end do +end subroutine chartmap_initial_guess_cartesian_3d + +subroutine chartmap_refine_cart_seed(self, x_target, zeta_period, u, res_norm, ierr) + class(chartmap_coordinate_system_t), intent(in) :: self + real(dp), intent(in) :: x_target(3) + real(dp), intent(in) :: zeta_period + real(dp), intent(inout) :: u(3) + real(dp), intent(out) :: res_norm + integer, intent(out) :: ierr + + real(dp) :: trial(3), best_u(3), x_trial(3) + real(dp) :: step_rho, step_theta, step_zeta + real(dp) :: best_dist2, dist2 + integer :: iter, ir, it, iz + + ierr = chartmap_from_cyl_ok + step_rho = 1.0_dp/16.0_dp + step_theta = TWOPI/64.0_dp + step_zeta = zeta_period/48.0_dp + + call chartmap_eval_cart(self, u, x_trial) + best_dist2 = sum((x_target - x_trial)**2) + + do iter = 1, 64 + best_u = u + do iz = -1, 1 + trial(3) = modulo(u(3) + real(iz, dp)*step_zeta, zeta_period) + do it = -1, 1 + trial(2) = modulo(u(2) + real(it, dp)*step_theta, TWOPI) + do ir = -1, 1 + trial(1) = min(max(u(1) + real(ir, dp)*step_rho, 0.0_dp), 1.0_dp) + call chartmap_eval_cart(self, trial, x_trial) + dist2 = sum((x_target - x_trial)**2) + if (dist2 < best_dist2) then + best_dist2 = dist2 + best_u = trial + end if + end do + end do + end do + + if (any(best_u /= u)) then + u = best_u + else + step_rho = 0.5_dp*step_rho + step_theta = 0.5_dp*step_theta + step_zeta = 0.5_dp*step_zeta + if (max(step_rho, max(step_theta, step_zeta)) < 1.0e-12_dp) exit + end if + end do + + call chartmap_eval_cart(self, u, x_trial) + res_norm = sqrt(sum((x_target - x_trial)**2)) + if (res_norm > 1.0e-4_dp) ierr = chartmap_from_cyl_err_max_iter +end subroutine chartmap_refine_cart_seed + subroutine chartmap_newton_delta(self, x_target, rho, theta, zeta, delta, & res_norm, & ierr) diff --git a/src/coordinates/libneo_coordinates_mapping.f90 b/src/coordinates/libneo_coordinates_mapping.f90 index cd10f6df..045c2334 100644 --- a/src/coordinates/libneo_coordinates_mapping.f90 +++ b/src/coordinates/libneo_coordinates_mapping.f90 @@ -19,6 +19,7 @@ subroutine map_vmec_u_to_chartmap_u(vmec_cs, ccs, u_vmec, u_chart, ierr, message character(len=*), intent(out) :: message real(dp) :: xcyl(3) + real(dp) :: xcart(3) ierr = 0 message = "" @@ -35,8 +36,20 @@ subroutine map_vmec_u_to_chartmap_u(vmec_cs, ccs, u_vmec, u_chart, ierr, message return end if - ierr = 1 - message = "unsupported zeta_convention for vmec mapping" + call cyl_to_cart(xcyl, xcart) + call ccs%from_cart(xcart, u_chart, ierr) + if (ierr /= 0) then + write (message, '(a,i0)') "chartmap from_cart failed ierr=", ierr + end if end subroutine map_vmec_u_to_chartmap_u + pure subroutine cyl_to_cart(xcyl, xcart) + real(dp), intent(in) :: xcyl(3) + real(dp), intent(out) :: xcart(3) + + xcart(1) = xcyl(1)*cos(xcyl(2)) + xcart(2) = xcyl(1)*sin(xcyl(2)) + xcart(3) = xcyl(3) + end subroutine cyl_to_cart + end module libneo_coordinates_mapping diff --git a/src/coordinates/libneo_coordinates_validator.f90 b/src/coordinates/libneo_coordinates_validator.f90 index 980cfe3a..2bc47d4a 100644 --- a/src/coordinates/libneo_coordinates_validator.f90 +++ b/src/coordinates/libneo_coordinates_validator.f90 @@ -237,11 +237,11 @@ subroutine check_optional_zeta_convention(ncid, value_out, ierr, message) end if select case (trim(value)) - case ("cyl", "vmec") + case ("cyl", "vmec", "boozer") value_out = trim(value) case default ierr = 5 - message = "unsupported zeta_convention (must be cyl or vmec)" + message = "unsupported zeta_convention (must be cyl, vmec, or boozer)" end select end subroutine check_optional_zeta_convention diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f3dcf7cf..7f2abecb 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -185,6 +185,9 @@ target_link_libraries(test_chartmap_basis_debug.x neo ${MAGFIE_LIB}) add_executable(test_chartmap_validator.x source/test_chartmap_validator.f90) target_link_libraries(test_chartmap_validator.x neo) +add_executable(test_chartmap_boozer_coordinates.x source/test_chartmap_boozer_coordinates.f90) +target_link_libraries(test_chartmap_boozer_coordinates.x neo ${MAGFIE_LIB}) + # TODO: Update test for new psi_1..psi_7 interface # add_executable(test_analytical_gs_circular.x source/test_analytical_gs_circular.f90) # target_link_libraries(test_analytical_gs_circular.x ${MAGFIE_LIB} ${COMMON_LIBS}) @@ -415,6 +418,12 @@ set_tests_properties(test_chartmap_validator PROPERTIES FAIL_REGULAR_EXPRESSION "STOP" WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +add_test(NAME test_chartmap_boozer_coordinates + COMMAND test_chartmap_boozer_coordinates.x) +set_tests_properties(test_chartmap_boozer_coordinates PROPERTIES + FAIL_REGULAR_EXPRESSION "STOP" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + add_library(chartmap_test_utils STATIC source/chartmap_test_utils.f90) target_link_libraries(chartmap_test_utils PUBLIC neo) diff --git a/test/source/test_chartmap_boozer_coordinates.f90 b/test/source/test_chartmap_boozer_coordinates.f90 new file mode 100644 index 00000000..0f71d753 --- /dev/null +++ b/test/source/test_chartmap_boozer_coordinates.f90 @@ -0,0 +1,216 @@ +program test_chartmap_boozer_coordinates + use, intrinsic :: iso_fortran_env, only: dp => real64 + use libneo_coordinates, only: coordinate_system_t, chartmap_coordinate_system_t, & + make_chartmap_coordinate_system, validate_chartmap_file, & + chartmap_from_cyl_ok, BOOZER + use math_constants, only: TWOPI + use netcdf + implicit none + + character(len=*), parameter :: filename = "chartmap_boozer.nc" + integer :: ierr, nerrors + character(len=2048) :: message + + nerrors = 0 + + call write_boozer_file(filename) + call validate_chartmap_file(filename, ierr, message) + if (ierr /= 0) then + print *, " FAIL: validator rejected Boozer chartmap: ", trim(message) + error stop 1 + end if + + call run_roundtrip_checks(filename, nerrors) + + if (nerrors > 0) then + print *, "FAILED: ", nerrors, " error(s) detected in Boozer chartmap tests" + error stop 1 + end if + + print *, "All Boozer chartmap coordinate tests passed!" + +contains + + subroutine nc_check(status) + integer, intent(in) :: status + if (status /= NF90_NOERR) then + print *, trim(nf90_strerror(status)) + error stop 2 + end if + end subroutine nc_check + + subroutine run_roundtrip_checks(path, nerrors) + character(len=*), intent(in) :: path + integer, intent(inout) :: nerrors + + class(coordinate_system_t), allocatable :: cs + real(dp) :: u(3), u_back(3), u_cyl(3), x(3), xcyl(3) + real(dp) :: du_theta, du_zeta, max_dphi + real(dp) :: zeta_period, dphi_geom + integer :: ir, it, iz, ierr_local + real(dp), parameter :: rho_vals(3) = [0.12_dp, 0.45_dp, 0.82_dp] + real(dp), parameter :: theta_vals(4) = [0.1_dp, 1.4_dp, 3.0_dp, 4.9_dp] + real(dp), parameter :: zeta_fracs(4) = [0.0_dp, 0.17_dp, 0.41_dp, 0.73_dp] + real(dp), parameter :: tol_u = 5.0e-6_dp + + max_dphi = 0.0_dp + call make_chartmap_coordinate_system(cs, path) + + select type (ccs => cs) + type is (chartmap_coordinate_system_t) + if (ccs%zeta_convention /= BOOZER) then + print *, " FAIL: expected BOOZER zeta convention, got ", & + ccs%zeta_convention + nerrors = nerrors + 1 + return + end if + + zeta_period = TWOPI/real(ccs%num_field_periods, dp) + + 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)*zeta_period] + 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) + + dphi_geom = angular_diff(xcyl(2), u(3), TWOPI) + max_dphi = max(max_dphi, dphi_geom) + + call ccs%from_cart(x, u_back, ierr_local) + if (ierr_local /= 0) then + print *, " FAIL: from_cart ierr=", ierr_local, " at u=", u + nerrors = nerrors + 1 + cycle + end if + + du_theta = angular_diff(u_back(2), u(2), TWOPI) + du_zeta = angular_diff(u_back(3), u(3), zeta_period) + if (abs(u_back(1) - u(1)) > tol_u .or. du_theta > tol_u .or. & + du_zeta > tol_u) then + print *, " FAIL: from_cart roundtrip mismatch" + print *, " u =", u + print *, " u_back=", u_back + nerrors = nerrors + 1 + cycle + end if + + call ccs%from_cyl(xcyl, u_cyl, ierr_local) + if (ierr_local /= chartmap_from_cyl_ok) then + print *, " FAIL: from_cyl ierr=", ierr_local, " at u=", u + nerrors = nerrors + 1 + cycle + end if + + du_theta = angular_diff(u_cyl(2), u(2), TWOPI) + du_zeta = angular_diff(u_cyl(3), u(3), zeta_period) + if (abs(u_cyl(1) - u(1)) > tol_u .or. du_theta > tol_u .or. & + du_zeta > tol_u) then + print *, " FAIL: from_cyl roundtrip mismatch" + print *, " u =", u + print *, " u_cyl =", u_cyl + nerrors = nerrors + 1 + end if + end do + end do + end do + + if (max_dphi < 5.0e-2_dp) then + print *, " FAIL: synthetic Boozer chartmap is too close to cylindrical" + nerrors = nerrors + 1 + else + print *, " PASS: Boozer chartmap uses non-cylindrical toroidal parameter" + end if + class default + print *, " FAIL: chartmap file did not load as chartmap type" + nerrors = nerrors + 1 + end select + end subroutine run_roundtrip_checks + + pure real(dp) function angular_diff(a, b, period) + real(dp), intent(in) :: a, b, period + angular_diff = abs(modulo(a - b + 0.5_dp*period, period) - 0.5_dp*period) + end function angular_diff + + subroutine write_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 = 9, ntheta = 11, nzeta = 10, nfp = 3 + real(dp), parameter :: r0 = 180.0_dp + real(dp), parameter :: a = 45.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.08_dp*rho_val*sin(theta_val) + & + 0.02_dp*sin(2.0_dp*zeta_val) + rmaj = r0 + a*rho_val*cos(theta_val) + & + 0.04_dp*a*rho_val*cos(theta_val - 2.0_dp*zeta_val) + zpos = a*rho_val*sin(theta_val) + & + 0.02_dp*a*rho_val*sin(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_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_boozer_file + +end program test_chartmap_boozer_coordinates diff --git a/test/source/test_chartmap_validator.f90 b/test/source/test_chartmap_validator.f90 index 247981bc..cc3cd6b7 100644 --- a/test/source/test_chartmap_validator.f90 +++ b/test/source/test_chartmap_validator.f90 @@ -21,6 +21,16 @@ program test_chartmap_validator print *, " PASS: validator accepted a known-good chartmap file" end if + call write_good_boozer_file("chartmap_good_boozer.nc") + call validate_chartmap_file("chartmap_good_boozer.nc", ierr, message) + if (ierr /= 0) then + print *, " FAIL: validator rejected a Boozer chartmap file" + print *, trim(message) + nerrors = nerrors + 1 + else + print *, " PASS: validator accepted a Boozer chartmap file" + end if + call write_missing_x_file("chartmap_missing_x.nc") call validate_chartmap_file("chartmap_missing_x.nc", ierr, message) if (ierr == 0 .or. index(message, "missing variable x") == 0) then @@ -189,6 +199,83 @@ subroutine write_good_file(filename) call nc_check(nf90_close(ncid)) end subroutine write_good_file + subroutine write_good_boozer_file(filename) + character(len=*), intent(in) :: filename + integer :: ncid + integer :: dim_rho, dim_theta, dim_zeta + integer :: var_rho, var_theta, var_zeta, var_x, var_y, var_z, & + var_num_field_periods + integer, parameter :: nrho = 3, ntheta = 5, nzeta = 6, nfp = 2 + real(dp), parameter :: r0 = 150.0_dp + real(dp), parameter :: a = 35.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.08_dp*rho_val*sin(theta_val) + & + 0.02_dp*sin(2.0_dp*zeta_val) + rmaj = r0 + a*rho_val*cos(theta_val) + & + 0.04_dp*a*rho_val*cos(theta_val - 2.0_dp*zeta_val) + zpos = a*rho_val*sin(theta_val) + & + 0.02_dp*a*rho_val*sin(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(filename), NF90_NETCDF4, ncid)) + call nc_check(nf90_put_att(ncid, NF90_GLOBAL, "zeta_convention", "boozer")) + 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_num_field_periods)) + 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_num_field_periods, nfp)) + call nc_check(nf90_close(ncid)) + end subroutine write_good_boozer_file + subroutine write_missing_units_file(filename) character(len=*), intent(in) :: filename integer :: ncid