diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 24c28c40..418723af 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -14,6 +14,7 @@ field/field_coils.f90 field/field_vmec.f90 field/field_splined.f90 + field/field_boozer_chartmap.f90 field/vmec_field_eval.f90 field/field_newton.F90 field.F90 diff --git a/src/boozer_converter.F90 b/src/boozer_converter.F90 index 1e0b1090..e316b9b3 100644 --- a/src/boozer_converter.F90 +++ b/src/boozer_converter.F90 @@ -19,6 +19,7 @@ module boozer_sub public :: vmec_to_boozer, boozer_to_vmec public :: delthe_delphi_BV public :: reset_boozer_batch_splines + public :: load_boozer_from_chartmap ! Constants real(dp), parameter :: TWOPI = 2.0_dp*3.14159265358979_dp @@ -877,6 +878,160 @@ subroutine reset_boozer_batch_splines if (allocated(delt_delp_B_grid)) deallocate (delt_delp_B_grid) end subroutine reset_boozer_batch_splines + subroutine load_boozer_from_chartmap(filename) + !> Populate module-level Boozer batch splines from an extended chartmap + !> NetCDF file, bypassing the VMEC-based compute_boozer_data path. + use vector_potentail_mod, only: torflux, ns, hs + use new_vmec_stuff_mod, only: nper, ns_A, ns_s, ns_tp + use boozer_coordinates_mod, only: ns_s_B, ns_tp_B, ns_B, n_theta_B, & + n_phi_B, hs_B, h_theta_B, h_phi_B, & + use_B_r, use_del_tp_B + use netcdf + + character(len=*), intent(in) :: filename + + integer :: ncid, status, dimid, varid + integer :: n_rho, n_theta, n_zeta, nfp_file + real(dp) :: torflux_val + real(dp), allocatable :: rho(:), theta(:), zeta(:) + real(dp), allocatable :: A_phi_arr(:), B_theta_arr(:), B_phi_arr(:) + real(dp), allocatable :: Bmod_arr(:, :, :) + real(dp), allocatable :: y_aphi(:, :), y_bcovar(:, :), y_bmod(:, :, :, :) + real(dp) :: s_min, s_max, rho_min, rho_max + integer :: i + integer :: spline_order + integer :: order_3d(3) + logical :: periodic_3d(3) + real(dp) :: x_min_3d(3), x_max_3d(3) + + call reset_boozer_batch_splines + + ! Open file + status = nf90_open(trim(filename), nf90_nowrite, ncid) + if (status /= nf90_noerr) then + print *, "load_boozer_from_chartmap: cannot open ", trim(filename) + error stop + end if + + ! Read dimensions + call nc_check(nf90_inq_dimid(ncid, "rho", dimid), "dim rho") + call nc_check(nf90_inquire_dimension(ncid, dimid, len=n_rho), "len rho") + call nc_check(nf90_inq_dimid(ncid, "theta", dimid), "dim theta") + call nc_check(nf90_inquire_dimension(ncid, dimid, len=n_theta), "len theta") + call nc_check(nf90_inq_dimid(ncid, "zeta", dimid), "dim zeta") + call nc_check(nf90_inquire_dimension(ncid, dimid, len=n_zeta), "len zeta") + + ! Read coordinate grids + allocate (rho(n_rho), theta(n_theta), zeta(n_zeta)) + call nc_check(nf90_inq_varid(ncid, "rho", varid), "var rho") + call nc_check(nf90_get_var(ncid, varid, rho), "get rho") + call nc_check(nf90_inq_varid(ncid, "theta", varid), "var theta") + call nc_check(nf90_get_var(ncid, varid, theta), "get theta") + call nc_check(nf90_inq_varid(ncid, "zeta", varid), "var zeta") + call nc_check(nf90_get_var(ncid, varid, zeta), "get zeta") + + ! Read scalar attributes and variables + call nc_check(nf90_get_att(ncid, nf90_global, "torflux", torflux_val), & + "att torflux") + call nc_check(nf90_inq_varid(ncid, "num_field_periods", varid), & + "var num_field_periods") + call nc_check(nf90_get_var(ncid, varid, nfp_file), & + "get num_field_periods") + + ! Read 1D profiles + allocate (A_phi_arr(n_rho), B_theta_arr(n_rho), B_phi_arr(n_rho)) + call nc_check(nf90_inq_varid(ncid, "A_phi", varid), "var A_phi") + call nc_check(nf90_get_var(ncid, varid, A_phi_arr), "get A_phi") + call nc_check(nf90_inq_varid(ncid, "B_theta", varid), "var B_theta") + call nc_check(nf90_get_var(ncid, varid, B_theta_arr), "get B_theta") + call nc_check(nf90_inq_varid(ncid, "B_phi", varid), "var B_phi") + call nc_check(nf90_get_var(ncid, varid, B_phi_arr), "get B_phi") + + ! Read 3D Bmod: NetCDF dims (zeta, theta, rho) map to Fortran + ! (rho, theta, zeta) via NF90 dimension reversal + allocate (Bmod_arr(n_rho, n_theta, n_zeta)) + call nc_check(nf90_inq_varid(ncid, "Bmod", varid), "var Bmod") + call nc_check(nf90_get_var(ncid, varid, Bmod_arr), "get Bmod") + + call nc_check(nf90_close(ncid), "close") + + ! Set global parameters used by splint_boozer_coord + torflux = torflux_val + nper = nfp_file + + ! Set boozer_coordinates_mod parameters + ns_s_B = 5 + ns_tp_B = 5 + ns_B = n_rho + n_theta_B = n_theta + n_phi_B = n_zeta + hs_B = rho(2) - rho(1) + h_theta_B = theta(2) - theta(1) + h_phi_B = zeta(2) - zeta(1) + use_B_r = .false. + use_del_tp_B = .false. + + ! Set vector_potentail_mod parameters for A_phi spline + rho_min = rho(1) + rho_max = rho(n_rho) + s_min = rho_min**2 + s_max = rho_max**2 + ns = n_rho + hs = (s_max - s_min) / real(n_rho - 1, dp) + ns_A = 5 + + ! Build A_phi batch spline over s + spline_order = ns_A + allocate (y_aphi(n_rho, 1)) + y_aphi(:, 1) = A_phi_arr + call construct_batch_splines_1d(s_min, s_max, y_aphi, spline_order, & + .false., aphi_batch_spline) + aphi_batch_spline_ready = .true. + deallocate (y_aphi) + + ! Build B_theta, B_phi batch spline over rho_tor + spline_order = ns_s_B + allocate (y_bcovar(n_rho, 2)) + y_bcovar(:, 1) = B_theta_arr + y_bcovar(:, 2) = B_phi_arr + call construct_batch_splines_1d(rho(1), rho(n_rho), y_bcovar, spline_order, & + .false., bcovar_tp_batch_spline) + bcovar_tp_batch_spline_ready = .true. + deallocate (y_bcovar) + + ! Build Bmod 3D batch spline over (rho_tor, theta_B, phi_B) + order_3d = [ns_s_B, ns_tp_B, ns_tp_B] + periodic_3d = [.false., .true., .true.] + x_min_3d = [rho(1), theta(1), zeta(1)] + x_max_3d = [rho(n_rho), theta(n_theta), zeta(n_zeta)] + + allocate (y_bmod(n_rho, n_theta, n_zeta, 1)) + y_bmod(:, :, :, 1) = Bmod_arr + call construct_batch_splines_3d(x_min_3d, x_max_3d, y_bmod, order_3d, & + periodic_3d, bmod_br_batch_spline) + bmod_br_batch_spline_ready = .true. + bmod_br_num_quantities = 1 + deallocate (y_bmod) + + print *, 'Loaded Boozer splines from chartmap: ', trim(filename) + print *, ' nfp=', nfp_file, ' ns=', n_rho, ' ntheta=', n_theta, & + ' nphi=', n_zeta + print *, ' torflux=', torflux_val + + contains + + subroutine nc_check(stat, loc) + integer, intent(in) :: stat + character(len=*), intent(in) :: loc + if (stat /= nf90_noerr) then + print *, "load_boozer_from_chartmap: NetCDF error at ", trim(loc), & + ": ", trim(nf90_strerror(stat)) + error stop + end if + end subroutine nc_check + + end subroutine load_boozer_from_chartmap + subroutine build_boozer_aphi_batch_spline use vector_potentail_mod, only: ns, hs, sA_phi use new_vmec_stuff_mod, only: ns_A diff --git a/src/field.F90 b/src/field.F90 index f1dabdbb..60f7698d 100644 --- a/src/field.F90 +++ b/src/field.F90 @@ -8,6 +8,9 @@ module field use field_vmec, only: vmec_field_t, create_vmec_field use field_coils, only: coils_field_t, create_coils_field use field_splined, only: splined_field_t, create_splined_field + use field_boozer_chartmap, only: boozer_chartmap_field_t, & + create_boozer_chartmap_field, & + is_boozer_chartmap implicit none @@ -42,6 +45,14 @@ subroutine field_clone(source, dest) class default error stop 'field_clone: Allocation failure (splined)' end select + type is (boozer_chartmap_field_t) + allocate (boozer_chartmap_field_t :: dest) + select type (dest) + type is (boozer_chartmap_field_t) + dest = source + class default + error stop 'field_clone: Allocation failure (boozer_chartmap)' + end select class default error stop 'field_clone: Unsupported field type' end select @@ -59,6 +70,7 @@ subroutine field_from_file(filename, field) type(coils_field_t) :: raw_coils type(splined_field_t), allocatable :: splined_coils type(vmec_field_t) :: vmec_field + type(boozer_chartmap_field_t), allocatable :: bc_temp integer :: file_type, ierr character(len=2048) :: message @@ -78,14 +90,19 @@ subroutine field_from_file(filename, field) call create_vmec_field(vmec_field) call field_clone(vmec_field, field) case (refcoords_file_chartmap) - print *, & - 'field_from_file: chartmap NetCDF is a coordinate system file,', & - ' not a field:' - print *, ' filename = ', trim(filename) - print *, & - 'Set coord_input to this chartmap file and set field_input to a', & - ' VMEC wout.' - error stop + if (is_boozer_chartmap(filename)) then + call create_boozer_chartmap_field(filename, bc_temp) + call move_alloc(bc_temp, field) + else + print *, & + 'field_from_file: chartmap NetCDF is a coordinate system ', & + 'file, not a field:' + print *, ' filename = ', trim(filename) + print *, & + 'Set coord_input to this chartmap file and set field_input', & + ' to a VMEC wout.' + error stop + end if case (refcoords_file_unknown) print *, 'field_from_file: Unknown NetCDF file type: ', trim(filename) error stop diff --git a/src/field/field_boozer_chartmap.f90 b/src/field/field_boozer_chartmap.f90 new file mode 100644 index 00000000..17cb27af --- /dev/null +++ b/src/field/field_boozer_chartmap.f90 @@ -0,0 +1,300 @@ +module field_boozer_chartmap + !> Boozer field from extended chartmap NetCDF file. + !> + !> Reads an extended chartmap file containing Boozer coordinate geometry + !> (X, Y, Z on a Boozer angle grid) plus magnetic field data: + !> - A_phi(rho) : toroidal vector potential (1D) + !> - B_theta(rho) : covariant poloidal B component (1D, surface function) + !> - B_phi(rho) : covariant toroidal B component (1D, surface function) + !> - Bmod(rho,theta,zeta) : field magnitude (3D) + !> - torflux : total toroidal flux (scalar attribute) + !> + !> This type serves as a magnetic_field_t for the non-canonical path and + !> provides data to populate the Boozer batch splines for the symplectic path. + + use, intrinsic :: iso_fortran_env, only: dp => real64 + use field_base, only: magnetic_field_t + use interpolate, only: BatchSplineData1D, BatchSplineData3D, & + construct_batch_splines_1d, construct_batch_splines_3d, & + evaluate_batch_splines_1d_der2, & + evaluate_batch_splines_3d, & + destroy_batch_splines_1d, destroy_batch_splines_3d + use netcdf + + implicit none + + private + public :: boozer_chartmap_field_t, create_boozer_chartmap_field, & + is_boozer_chartmap + + type, extends(magnetic_field_t) :: boozer_chartmap_field_t + type(BatchSplineData1D) :: aphi_spline + type(BatchSplineData1D) :: bcovar_spline + type(BatchSplineData3D) :: bmod_spline + real(dp) :: torflux = 0.0_dp + integer :: nfp = 1 + integer :: ns = 0 + integer :: ntheta = 0 + integer :: nphi = 0 + real(dp) :: hs = 0.0_dp + real(dp) :: h_theta = 0.0_dp + real(dp) :: h_phi = 0.0_dp + character(len=512) :: filename = '' + logical :: initialized = .false. + contains + procedure :: evaluate => boozer_chartmap_evaluate + final :: boozer_chartmap_cleanup + end type boozer_chartmap_field_t + +contains + + function is_boozer_chartmap(filename) result(is_bc) + character(len=*), intent(in) :: filename + logical :: is_bc + integer :: ncid, status, ival + + is_bc = .false. + status = nf90_open(trim(filename), nf90_nowrite, ncid) + if (status /= nf90_noerr) return + + status = nf90_get_att(ncid, nf90_global, "boozer_field", ival) + if (status == nf90_noerr .and. ival == 1) is_bc = .true. + + status = nf90_close(ncid) + end function is_boozer_chartmap + + subroutine create_boozer_chartmap_field(filename, field) + use libneo_coordinates, only: coordinate_system_t, & + make_chartmap_coordinate_system + + character(len=*), intent(in) :: filename + type(boozer_chartmap_field_t), allocatable, intent(out) :: field + + integer :: ncid, status, dimid, varid + integer :: n_rho, n_theta, n_zeta, nfp_file + real(dp) :: torflux_val + real(dp), allocatable :: rho(:), theta(:), zeta(:) + real(dp), allocatable :: A_phi_arr(:), B_theta_arr(:), B_phi_arr(:) + real(dp), allocatable :: Bmod_arr(:, :, :) + real(dp), allocatable :: s_grid(:) + real(dp), allocatable :: y_aphi(:, :), y_bcovar(:, :), y_bmod(:, :, :, :) + real(dp) :: s_min, s_max, rho_min, rho_max + real(dp) :: h_s, h_theta_val, h_phi_val + integer :: i + integer, parameter :: spline_order_1d = 5 + integer, parameter :: spline_order_3d(3) = [5, 5, 5] + logical, parameter :: periodic_3d(3) = [.false., .true., .true.] + real(dp) :: x_min_3d(3), x_max_3d(3) + class(coordinate_system_t), allocatable :: cs + + allocate (field) + + ! Open file + status = nf90_open(trim(filename), nf90_nowrite, ncid) + call check_nc(status, "open " // trim(filename)) + + ! Read dimensions + call check_nc(nf90_inq_dimid(ncid, "rho", dimid), "inq_dim rho") + call check_nc(nf90_inquire_dimension(ncid, dimid, len=n_rho), "get_dim rho") + call check_nc(nf90_inq_dimid(ncid, "theta", dimid), "inq_dim theta") + call check_nc(nf90_inquire_dimension(ncid, dimid, len=n_theta), "get_dim theta") + call check_nc(nf90_inq_dimid(ncid, "zeta", dimid), "inq_dim zeta") + call check_nc(nf90_inquire_dimension(ncid, dimid, len=n_zeta), "get_dim zeta") + + ! Read coordinate arrays + allocate (rho(n_rho), theta(n_theta), zeta(n_zeta)) + call check_nc(nf90_inq_varid(ncid, "rho", varid), "inq_var rho") + call check_nc(nf90_get_var(ncid, varid, rho), "get_var rho") + call check_nc(nf90_inq_varid(ncid, "theta", varid), "inq_var theta") + call check_nc(nf90_get_var(ncid, varid, theta), "get_var theta") + call check_nc(nf90_inq_varid(ncid, "zeta", varid), "inq_var zeta") + call check_nc(nf90_get_var(ncid, varid, zeta), "get_var zeta") + + ! Read scalar attributes and variables + call check_nc(nf90_get_att(ncid, nf90_global, "torflux", torflux_val), & + "get_att torflux") + call check_nc(nf90_inq_varid(ncid, "num_field_periods", varid), & + "inq_var num_field_periods") + call check_nc(nf90_get_var(ncid, varid, nfp_file), & + "get_var num_field_periods") + + ! Read 1D profiles + allocate (A_phi_arr(n_rho), B_theta_arr(n_rho), B_phi_arr(n_rho)) + + call check_nc(nf90_inq_varid(ncid, "A_phi", varid), "inq_var A_phi") + call check_nc(nf90_get_var(ncid, varid, A_phi_arr), "get_var A_phi") + + call check_nc(nf90_inq_varid(ncid, "B_theta", varid), "inq_var B_theta") + call check_nc(nf90_get_var(ncid, varid, B_theta_arr), "get_var B_theta") + + call check_nc(nf90_inq_varid(ncid, "B_phi", varid), "inq_var B_phi") + call check_nc(nf90_get_var(ncid, varid, B_phi_arr), "get_var B_phi") + + ! Read 3D Bmod (stored as zeta, theta, rho in NetCDF → read into rho, theta, zeta) + allocate (Bmod_arr(n_rho, n_theta, n_zeta)) + call check_nc(nf90_inq_varid(ncid, "Bmod", varid), "inq_var Bmod") + call read_3d_reordered(ncid, varid, n_rho, n_theta, n_zeta, Bmod_arr) + + call check_nc(nf90_close(ncid), "close") + + ! Grid parameters + rho_min = rho(1) + rho_max = rho(n_rho) + h_theta_val = theta(2) - theta(1) + h_phi_val = zeta(2) - zeta(1) + + ! Build A_phi spline over s = rho^2 (matching VMEC/boozer_converter convention) + allocate (s_grid(n_rho)) + do i = 1, n_rho + s_grid(i) = rho(i)**2 + end do + s_min = s_grid(1) + s_max = s_grid(n_rho) + h_s = (s_max - s_min) / real(n_rho - 1, dp) + + allocate (y_aphi(n_rho, 1)) + y_aphi(:, 1) = A_phi_arr + call construct_batch_splines_1d(s_min, s_max, y_aphi, spline_order_1d, & + .false., field%aphi_spline) + + ! Build B_theta, B_phi spline over rho_tor + allocate (y_bcovar(n_rho, 2)) + y_bcovar(:, 1) = B_theta_arr + y_bcovar(:, 2) = B_phi_arr + call construct_batch_splines_1d(rho_min, rho_max, y_bcovar, spline_order_1d, & + .false., field%bcovar_spline) + + ! Build Bmod 3D spline over (rho_tor, theta_B, phi_B) + x_min_3d = [rho_min, theta(1), zeta(1)] + x_max_3d = [rho_max, theta(n_theta), zeta(n_zeta)] + + allocate (y_bmod(n_rho, n_theta, n_zeta, 1)) + y_bmod(:, :, :, 1) = Bmod_arr + call construct_batch_splines_3d(x_min_3d, x_max_3d, y_bmod, & + spline_order_3d, periodic_3d, & + field%bmod_spline) + + ! Store metadata + field%torflux = torflux_val + field%nfp = nfp_file + field%ns = n_rho + field%ntheta = n_theta + field%nphi = n_zeta + field%hs = h_s + field%h_theta = h_theta_val + field%h_phi = h_phi_val + field%filename = filename + + ! Set up coordinate system from the same chartmap file + call make_chartmap_coordinate_system(cs, filename) + allocate (field%coords, source=cs) + + field%initialized = .true. + + print *, 'Loaded Boozer chartmap field from ', trim(filename) + print *, ' nfp=', nfp_file, ' ns=', n_rho, ' ntheta=', n_theta, & + ' nphi=', n_zeta + print *, ' torflux=', torflux_val + end subroutine create_boozer_chartmap_field + + subroutine boozer_chartmap_evaluate(self, x, Acov, hcov, Bmod, sqgBctr) + !> Evaluate field at x = (rho, theta_B, phi_B) in chartmap coordinates. + !> Input x(1) = rho = sqrt(s), matching the chartmap coordinate system. + !> Returns covariant components in (rho, theta_B, phi_B) coordinates. + class(boozer_chartmap_field_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: Acov(3), hcov(3), Bmod + real(dp), intent(out), optional :: sqgBctr(3) + + real(dp) :: rho_tor, s, theta_B, phi_B + real(dp) :: A_phi_val, A_theta_val + real(dp) :: y1d_aphi(1), dy1d_aphi(1), d2y1d_aphi(1) + real(dp) :: y1d_bc(2), dy1d_bc(2), d2y1d_bc(2) + real(dp) :: x_3d(3), y_bmod(1) + real(dp) :: B_theta_val, B_phi_val, Bmod_val + real(dp), parameter :: twopi = 8.0_dp * atan(1.0_dp) + + rho_tor = x(1) + theta_B = x(2) + phi_B = x(3) + s = rho_tor**2 + + ! A_theta = torflux * s (linear, by definition) + ! In rho coordinates: A_theta(rho) = torflux * rho^2 + ! Covariant component transforms: A_rho = A_s * ds/drho = A_s * 2*rho + ! But A_s = 0 (gauge), so A_rho = 0 regardless. + A_theta_val = self%torflux * s + + ! A_phi from 1D spline at s + call evaluate_batch_splines_1d_der2(self%aphi_spline, s, y1d_aphi, & + dy1d_aphi, d2y1d_aphi) + A_phi_val = y1d_aphi(1) + + ! B_theta, B_phi from 1D spline at rho_tor + call evaluate_batch_splines_1d_der2(self%bcovar_spline, rho_tor, y1d_bc, & + dy1d_bc, d2y1d_bc) + B_theta_val = y1d_bc(1) + B_phi_val = y1d_bc(2) + + ! Bmod from 3D spline at (rho_tor, theta_B, phi_B) + x_3d(1) = rho_tor + x_3d(2) = modulo(theta_B, twopi) + x_3d(3) = modulo(phi_B, twopi / real(self%nfp, dp)) + + call evaluate_batch_splines_3d(self%bmod_spline, x_3d, y_bmod) + Bmod_val = y_bmod(1) + + Bmod = Bmod_val + + ! Covariant vector potential in (rho, theta, phi) coordinates. + ! A_theta and A_phi are the same in (s, theta, phi) and (rho, theta, phi) + ! since the angular coordinates are unchanged. A_rho = 0. + Acov(1) = 0.0_dp + Acov(2) = A_theta_val + Acov(3) = A_phi_val + + ! Covariant unit field direction: h = B/|B| + ! B_theta and B_phi are covariant and independent of the radial + ! coordinate choice (they refer to the angular basis vectors). + hcov(1) = 0.0_dp + hcov(2) = B_theta_val / Bmod_val + hcov(3) = B_phi_val / Bmod_val + + if (present(sqgBctr)) then + error stop "sqgBctr not implemented for boozer_chartmap_field_t" + end if + end subroutine boozer_chartmap_evaluate + + subroutine boozer_chartmap_cleanup(self) + type(boozer_chartmap_field_t), intent(inout) :: self + + if (self%initialized) then + call destroy_batch_splines_1d(self%aphi_spline) + call destroy_batch_splines_1d(self%bcovar_spline) + call destroy_batch_splines_3d(self%bmod_spline) + self%initialized = .false. + end if + end subroutine boozer_chartmap_cleanup + + subroutine read_3d_reordered(ncid, varid, n1, n2, n3, arr) + !> Read a NetCDF variable with dims (zeta, theta, rho) into Fortran + !> array (n1=rho, n2=theta, n3=zeta). NF90 reverses dimension order + !> so the data lands directly in (rho, theta, zeta) layout. + integer, intent(in) :: ncid, varid, n1, n2, n3 + real(dp), intent(out) :: arr(n1, n2, n3) + + call check_nc(nf90_get_var(ncid, varid, arr), "get_var 3D") + end subroutine read_3d_reordered + + subroutine check_nc(status, location) + integer, intent(in) :: status + character(len=*), intent(in) :: location + + if (status /= nf90_noerr) then + print *, "NetCDF error at ", trim(location), ": ", & + trim(nf90_strerror(status)) + error stop "NetCDF operation failed" + end if + end subroutine check_nc + +end module field_boozer_chartmap diff --git a/src/field/field_can_boozer.f90 b/src/field/field_can_boozer.f90 index bc2ea7c8..93b948ea 100644 --- a/src/field/field_can_boozer.f90 +++ b/src/field/field_can_boozer.f90 @@ -4,8 +4,32 @@ module field_can_boozer implicit none + public :: integ_to_ref_boozer_chartmap, ref_to_integ_boozer_chartmap + contains + subroutine integ_to_ref_boozer_chartmap(xinteg, xref) + !> Convert integrator coords (s, theta_B, phi_B) to + !> chartmap reference coords (rho, theta_B, phi_B). + real(dp), intent(in) :: xinteg(3) + real(dp), intent(out) :: xref(3) + + xref(1) = sqrt(max(xinteg(1), 0.0_dp)) + xref(2) = xinteg(2) + xref(3) = xinteg(3) + end subroutine integ_to_ref_boozer_chartmap + + subroutine ref_to_integ_boozer_chartmap(xref, xinteg) + !> Convert chartmap reference coords (rho, theta_B, phi_B) to + !> integrator coords (s, theta_B, phi_B). + real(dp), intent(in) :: xref(3) + real(dp), intent(out) :: xinteg(3) + + xinteg(1) = xref(1)**2 + xinteg(2) = xref(2) + xinteg(3) = xref(3) + end subroutine ref_to_integ_boozer_chartmap + subroutine evaluate_boozer(f, r, th_c, ph_c, mode_secders) type(field_can_t), intent(inout) :: f real(dp), intent(in) :: r, th_c, ph_c diff --git a/src/field/field_splined.f90 b/src/field/field_splined.f90 index 4c5655ee..8909a25f 100644 --- a/src/field/field_splined.f90 +++ b/src/field/field_splined.f90 @@ -30,6 +30,7 @@ module field_splined destroy_batch_splines_3d use field_base, only: magnetic_field_t use field_vmec, only: vmec_field_t + use field_boozer_chartmap, only: boozer_chartmap_field_t use libneo_coordinates, only: coordinate_system_t, vmec_coordinate_system_t, & chartmap_coordinate_system_t, RHO_TOR, RHO_POL, & PSI_TOR_NORM, PSI_POL_NORM @@ -461,6 +462,14 @@ subroutine evaluate_at_ref_coords(source, ref_coords, x_scaled, x_ref, x_cart, & end select needs_scaling = .true. return + type is (boozer_chartmap_field_t) + ! Boozer chartmap evaluates directly in (rho, theta_B, phi_B), + ! which matches the chartmap reference coordinate system. + call source%evaluate(x_ref, Acov, hcov, Bmod) + x_cart = 0d0 + e_cov = 0d0 + needs_scaling = .true. + return class default continue end select diff --git a/src/field_can.f90 b/src/field_can.f90 index cd7efbd6..d58607fd 100644 --- a/src/field_can.f90 +++ b/src/field_can.f90 @@ -7,11 +7,13 @@ module field_can_mod identity_transform, field_can_t use field_can_test, only : evaluate_test use field_can_flux, only : evaluate_flux, integ_to_ref_flux, ref_to_integ_flux -use field_can_boozer, only : evaluate_boozer, integ_to_ref_boozer, ref_to_integ_boozer +use field_can_boozer, only : evaluate_boozer, integ_to_ref_boozer, ref_to_integ_boozer, & + integ_to_ref_boozer_chartmap, ref_to_integ_boozer_chartmap use field_can_meiss, only : init_meiss, evaluate_meiss, & integ_to_ref_meiss, ref_to_integ_meiss use field_can_albert, only : evaluate_albert, init_albert, integ_to_ref_albert, & ref_to_integ_albert +use field_boozer_chartmap, only : boozer_chartmap_field_t implicit none @@ -124,7 +126,8 @@ end function id_from_name subroutine init_field_can(field_id, field_noncan) use get_can_sub, only : get_canonical_coordinates, get_canonical_coordinates_with_field - use boozer_sub, only : get_boozer_coordinates, get_boozer_coordinates_with_field + use boozer_sub, only : get_boozer_coordinates, get_boozer_coordinates_with_field, & + load_boozer_from_chartmap use field_can_meiss, only : get_meiss_coordinates use field_can_albert, only : get_albert_coordinates @@ -140,7 +143,15 @@ subroutine init_field_can(field_id, field_noncan) case (CANFLUX) call get_canonical_coordinates_with_field(field_noncan) case (BOOZER) - call get_boozer_coordinates_with_field(field_noncan) + select type (f => field_noncan) + type is (boozer_chartmap_field_t) + call load_boozer_from_chartmap(f%filename) + ! Chartmap ref coords use rho=sqrt(s), integrator uses s + integ_to_ref => integ_to_ref_boozer_chartmap + ref_to_integ => ref_to_integ_boozer_chartmap + class default + call get_boozer_coordinates_with_field(field_noncan) + end select case (MEISS) call get_meiss_coordinates case (ALBERT) diff --git a/src/simple_main.f90 b/src/simple_main.f90 index cd7d8842..9f23ac34 100644 --- a/src/simple_main.f90 +++ b/src/simple_main.f90 @@ -41,17 +41,19 @@ module simple_main subroutine main use params, only: read_config, netcdffile, ns_s, ns_tp, multharm, & integmode, params_init, swcoll, generate_start_only, & - isw_field_type, & + isw_field_type, field_input, startmode, & ntestpart, ntimstep, coord_input use timing, only: init_timer, print_phase_time use magfie_sub, only: TEST, VMEC, init_magfie use samplers, only: init_starting_surf use version, only: simple_version + use field_boozer_chartmap, only: is_boozer_chartmap implicit none character(256) :: config_file type(tracer_t) :: norb + logical :: chartmap_mode ! Print version on startup print '(A,A)', 'SIMPLE version ', simple_version @@ -85,12 +87,32 @@ subroutine main call print_phase_time('Collision initialization completed') end if + chartmap_mode = .false. + if (len_trim(field_input) > 0) then + chartmap_mode = is_boozer_chartmap(field_input) + end if + if (isw_field_type == TEST) then ! TEST field uses analytic tokamak - no VMEC needed for sampling call init_magfie(TEST) call print_phase_time('TEST field initialization completed') call sample_particles_test_field call print_phase_time('TEST field particle sampling completed') + else if (chartmap_mode) then + ! Boozer chartmap: no VMEC, use Boozer magfie for field line tracing. + if (startmode /= 2) then + error stop 'Boozer chartmap requires startmode=2 '// & + '(load particles from start.dat)' + end if + + call init_magfie(isw_field_type) + call print_phase_time('Boozer magfie initialization completed') + + call init_starting_surf + call print_phase_time('Starting surface initialization completed') + + call sample_particles + call print_phase_time('Particle sampling completed') else call init_magfie(VMEC) call print_phase_time('VMEC magnetic field initialization completed') @@ -126,6 +148,7 @@ end subroutine main subroutine init_field(self, vmec_file, ans_s, ans_tp, amultharm, aintegmode) use field_base, only: magnetic_field_t use field, only: field_from_file + use field_boozer_chartmap, only: boozer_chartmap_field_t, is_boozer_chartmap use timing, only: print_phase_time use magfie_sub, only: TEST, CANFLUX, VMEC, BOOZER, MEISS, ALBERT, & REFCOORDS, set_magfie_refcoords_field @@ -140,13 +163,43 @@ subroutine init_field(self, vmec_file, ans_s, ans_tp, amultharm, aintegmode) integer, intent(in) :: ans_s, ans_tp, amultharm, aintegmode class(magnetic_field_t), allocatable :: field_temp character(:), allocatable :: vmec_equilibrium_file + logical :: use_boozer_chartmap self%integmode = aintegmode + ! Check if field_input is a Boozer chartmap (no VMEC needed) + use_boozer_chartmap = .false. + if (len_trim(field_input) > 0) then + use_boozer_chartmap = is_boozer_chartmap(field_input) + end if + ! TEST field is analytic - no VMEC or field files needed if (isw_field_type == TEST) then self%fper = twopi ! Full torus for analytic tokamak call print_phase_time('TEST field mode - no input files required') + else if (use_boozer_chartmap) then + ! Boozer chartmap: file-based, no VMEC initialization needed + call init_reference_coordinates(coord_input) + call print_phase_time('Reference coordinate system '// & + 'initialization completed') + + block + use libneo_coordinates, only: chartmap_coordinate_system_t + select type (cs => ref_coords) + type is (chartmap_coordinate_system_t) + self%fper = twopi / real(cs%num_field_periods, dp) + class default + self%fper = twopi + end select + end block + + call init_stl_wall_if_enabled(coord_input) + call print_phase_time('STL wall initialization completed') + + if (self%integmode >= 0) then + call field_from_file(field_input, field_temp) + call print_phase_time('Boozer chartmap field loading completed') + end if else vmec_equilibrium_file = select_vmec_equilibrium_file(vmec_file, & field_input, & diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt index d249240c..46af006e 100644 --- a/test/tests/CMakeLists.txt +++ b/test/tests/CMakeLists.txt @@ -331,6 +331,26 @@ add_test(NAME test_array_utils COMMAND test_array_utils.x) add_test(NAME test_chartmap_metadata COMMAND test_chartmap_metadata.x) set_tests_properties(test_chartmap_metadata PROPERTIES LABELS "unit") + # Boozer chartmap test: generate test file with a compiled NetCDF tool, then run Fortran test + add_executable(generate_test_boozer_chartmap.x generate_test_boozer_chartmap.c) + target_link_libraries(generate_test_boozer_chartmap.x netcdf m) + add_custom_command( + OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/test_boozer_chartmap.nc + COMMAND $ ${CMAKE_CURRENT_BINARY_DIR}/test_boozer_chartmap.nc + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + DEPENDS generate_test_boozer_chartmap.x + COMMENT "Generating test Boozer chartmap NetCDF" + ) + add_custom_target(generate_test_boozer_chartmap + DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/test_boozer_chartmap.nc) + add_executable(test_boozer_chartmap.x test_boozer_chartmap.f90) + target_link_libraries(test_boozer_chartmap.x simple) + add_dependencies(test_boozer_chartmap.x generate_test_boozer_chartmap) + add_test(NAME test_boozer_chartmap COMMAND test_boozer_chartmap.x) + set_tests_properties(test_boozer_chartmap PROPERTIES + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + LABELS "unit") + if(SIMPLE_ENABLE_CGAL) add_executable(test_stl_wall_intersection.x test_stl_wall_intersection.f90) target_link_libraries(test_stl_wall_intersection.x simple) diff --git a/test/tests/generate_test_boozer_chartmap.c b/test/tests/generate_test_boozer_chartmap.c new file mode 100644 index 00000000..34f8e9ef --- /dev/null +++ b/test/tests/generate_test_boozer_chartmap.c @@ -0,0 +1,164 @@ +#include +#include +#include +#include +#include + +#define CHECK_NC(call, context) \ + do { \ + int status_ = (call); \ + if (status_ != NC_NOERR) { \ + fprintf(stderr, "NetCDF error in %s: %s\n", context, \ + nc_strerror(status_)); \ + return 1; \ + } \ + } while (0) + +int main(int argc, char **argv) +{ + enum + { + nfp = 2, + n_rho = 17, + n_theta = 33, + n_zeta = 33 + }; + + const double b0 = 2.0e4; + const double r0 = 150.0; + const double a = 50.0; + const double iota = 0.5; + const double twopi = 2.0 * acos(-1.0); + const double epsilon = a / r0; + const double torflux = 0.5 * b0 * a * a; + const char *filename = argc > 1 ? argv[1] : "test_boozer_chartmap.nc"; + + double rho[n_rho]; + double theta[n_theta]; + double zeta[n_zeta]; + double a_phi[n_rho]; + double b_theta[n_rho]; + double b_phi[n_rho]; + double *x = NULL; + double *y = NULL; + double *z = NULL; + double *bmod = NULL; + int ncid; + int dim_rho, dim_theta, dim_zeta; + int dims_3d[3]; + int var_rho, var_theta, var_zeta, var_x, var_y, var_z, var_aphi, var_btheta; + int var_bphi, var_bmod, var_nfp; + + x = malloc((size_t)n_rho * n_theta * n_zeta * sizeof(*x)); + y = malloc((size_t)n_rho * n_theta * n_zeta * sizeof(*y)); + z = malloc((size_t)n_rho * n_theta * n_zeta * sizeof(*z)); + bmod = malloc((size_t)n_rho * n_theta * n_zeta * sizeof(*bmod)); + if (x == NULL || y == NULL || z == NULL || bmod == NULL) { + fprintf(stderr, "Allocation failure\n"); + free(x); + free(y); + free(z); + free(bmod); + return 1; + } + + for (int ir = 0; ir < n_rho; ++ir) { + rho[ir] = (double)ir / (double)(n_rho - 1); + a_phi[ir] = -torflux * iota * rho[ir] * rho[ir]; + b_theta[ir] = b0 * epsilon * iota; + b_phi[ir] = b0 * r0; + } + rho[0] = 1.0e-6; + a_phi[0] = -torflux * iota * rho[0] * rho[0]; + + for (int it = 0; it < n_theta; ++it) { + theta[it] = twopi * (double)it / (double)n_theta; + } + + for (int iz = 0; iz < n_zeta; ++iz) { + zeta[iz] = (twopi / (double)nfp) * (double)iz / (double)n_zeta; + } + + for (int iz = 0; iz < n_zeta; ++iz) { + const double zeta_val = zeta[iz]; + for (int it = 0; it < n_theta; ++it) { + const double theta_val = theta[it]; + const double cos_theta = cos(theta_val); + const double sin_theta = sin(theta_val); + for (int ir = 0; ir < n_rho; ++ir) { + const double rho_val = rho[ir]; + const double r_minor = a * rho_val; + const double radius = r0 + r_minor * cos_theta; + const size_t idx = ((size_t)iz * n_theta + (size_t)it) * n_rho + + (size_t)ir; + x[idx] = radius * cos(zeta_val); + y[idx] = radius * sin(zeta_val); + z[idx] = r_minor * sin_theta; + bmod[idx] = b0 / (1.0 + rho_val * epsilon * cos_theta); + } + } + } + + CHECK_NC(nc_create(filename, NC_NETCDF4, &ncid), "create file"); + CHECK_NC(nc_def_dim(ncid, "rho", n_rho, &dim_rho), "def dim rho"); + CHECK_NC(nc_def_dim(ncid, "theta", n_theta, &dim_theta), "def dim theta"); + CHECK_NC(nc_def_dim(ncid, "zeta", n_zeta, &dim_zeta), "def dim zeta"); + + dims_3d[0] = dim_zeta; + dims_3d[1] = dim_theta; + dims_3d[2] = dim_rho; + + CHECK_NC(nc_def_var(ncid, "rho", NC_DOUBLE, 1, &dim_rho, &var_rho), "def rho"); + CHECK_NC(nc_def_var(ncid, "theta", NC_DOUBLE, 1, &dim_theta, &var_theta), "def theta"); + CHECK_NC(nc_def_var(ncid, "zeta", NC_DOUBLE, 1, &dim_zeta, &var_zeta), "def zeta"); + CHECK_NC(nc_def_var(ncid, "x", NC_DOUBLE, 3, dims_3d, &var_x), "def x"); + CHECK_NC(nc_def_var(ncid, "y", NC_DOUBLE, 3, dims_3d, &var_y), "def y"); + CHECK_NC(nc_def_var(ncid, "z", NC_DOUBLE, 3, dims_3d, &var_z), "def z"); + CHECK_NC(nc_def_var(ncid, "A_phi", NC_DOUBLE, 1, &dim_rho, &var_aphi), "def A_phi"); + CHECK_NC(nc_def_var(ncid, "B_theta", NC_DOUBLE, 1, &dim_rho, &var_btheta), "def B_theta"); + CHECK_NC(nc_def_var(ncid, "B_phi", NC_DOUBLE, 1, &dim_rho, &var_bphi), "def B_phi"); + CHECK_NC(nc_def_var(ncid, "Bmod", NC_DOUBLE, 3, dims_3d, &var_bmod), "def Bmod"); + CHECK_NC(nc_def_var(ncid, "num_field_periods", NC_INT, 0, NULL, &var_nfp), + "def num_field_periods"); + + CHECK_NC(nc_put_att_text(ncid, var_x, "units", 2, "cm"), "put x units"); + CHECK_NC(nc_put_att_text(ncid, var_y, "units", 2, "cm"), "put y units"); + CHECK_NC(nc_put_att_text(ncid, var_z, "units", 2, "cm"), "put z units"); + CHECK_NC(nc_put_att_text(ncid, NC_GLOBAL, "rho_convention", 7, "rho_tor"), + "put rho_convention"); + CHECK_NC(nc_put_att_text(ncid, NC_GLOBAL, "zeta_convention", 6, "boozer"), + "put zeta_convention"); + CHECK_NC(nc_put_att_double(ncid, NC_GLOBAL, "rho_lcfs", NC_DOUBLE, 1, + (const double[]){1.0}), + "put rho_lcfs"); + CHECK_NC(nc_put_att_int(ncid, NC_GLOBAL, "boozer_field", NC_INT, 1, + (const int[]){1}), + "put boozer_field"); + CHECK_NC(nc_put_att_double(ncid, NC_GLOBAL, "torflux", NC_DOUBLE, 1, &torflux), + "put torflux"); + + CHECK_NC(nc_enddef(ncid), "enddef"); + + CHECK_NC(nc_put_var_double(ncid, var_rho, rho), "put rho"); + CHECK_NC(nc_put_var_double(ncid, var_theta, theta), "put theta"); + CHECK_NC(nc_put_var_double(ncid, var_zeta, zeta), "put zeta"); + CHECK_NC(nc_put_var_double(ncid, var_x, x), "put x"); + CHECK_NC(nc_put_var_double(ncid, var_y, y), "put y"); + CHECK_NC(nc_put_var_double(ncid, var_z, z), "put z"); + CHECK_NC(nc_put_var_double(ncid, var_aphi, a_phi), "put A_phi"); + CHECK_NC(nc_put_var_double(ncid, var_btheta, b_theta), "put B_theta"); + CHECK_NC(nc_put_var_double(ncid, var_bphi, b_phi), "put B_phi"); + CHECK_NC(nc_put_var_double(ncid, var_bmod, bmod), "put Bmod"); + CHECK_NC(nc_put_var_int(ncid, var_nfp, (const int[]){nfp}), "put num_field_periods"); + CHECK_NC(nc_close(ncid), "close"); + + printf("Generated %s\n", filename); + printf(" torflux=%14.6e, B0=%8.1f, R0=%8.1f\n", torflux, b0, r0); + printf(" nfp=%d, nrho=%d, ntheta=%d, nzeta=%d\n", nfp, n_rho, n_theta, n_zeta); + + free(x); + free(y); + free(z); + free(bmod); + return 0; +} diff --git a/test/tests/test_boozer_chartmap.f90 b/test/tests/test_boozer_chartmap.f90 new file mode 100644 index 00000000..7459d174 --- /dev/null +++ b/test/tests/test_boozer_chartmap.f90 @@ -0,0 +1,133 @@ +program test_boozer_chartmap + use, intrinsic :: iso_fortran_env, only: dp => real64 + use field_boozer_chartmap, only: boozer_chartmap_field_t, & + create_boozer_chartmap_field, & + is_boozer_chartmap + + implicit none + + type(boozer_chartmap_field_t), allocatable :: field + real(dp) :: x(3), Acov(3), hcov(3), Bmod + character(len=256) :: filename + logical :: is_bc + integer :: nfail + + nfail = 0 + filename = 'test_boozer_chartmap.nc' + + ! Test detection + is_bc = is_boozer_chartmap(filename) + if (.not. is_bc) then + print *, 'FAIL: is_boozer_chartmap returned false' + nfail = nfail + 1 + else + print *, 'PASS: is_boozer_chartmap detected file correctly' + end if + + ! Test loading + call create_boozer_chartmap_field(filename, field) + + if (.not. field%initialized) then + print *, 'FAIL: field not initialized' + error stop + end if + + if (field%nfp /= 2) then + print *, 'FAIL: nfp =', field%nfp, ' expected 2' + nfail = nfail + 1 + else + print *, 'PASS: nfp = 2' + end if + + if (field%torflux <= 0.0_dp) then + print *, 'FAIL: torflux =', field%torflux, ' expected > 0' + nfail = nfail + 1 + else + print *, 'PASS: torflux =', field%torflux + end if + + ! Test field evaluation at mid-radius, theta=0, phi=0 + ! x(1) = rho = sqrt(s), so rho=0.5 means s=0.25 + x = [0.5_dp, 0.0_dp, 0.0_dp] ! rho=0.5, theta=0, phi=0 + call field%evaluate(x, Acov, hcov, Bmod) + + ! Bmod should be positive and reasonable (analytic tokamak ~ B0) + if (Bmod <= 0.0_dp) then + print *, 'FAIL: Bmod =', Bmod, ' expected > 0' + nfail = nfail + 1 + else + print *, 'PASS: Bmod =', Bmod + end if + + ! A_theta = torflux * s = torflux * rho^2 = torflux * 0.25 + if (abs(Acov(2) - field%torflux * 0.25_dp) > 1.0e-8_dp * abs(Acov(2))) then + print *, 'FAIL: Acov(2) =', Acov(2), ' expected', field%torflux * 0.25_dp + nfail = nfail + 1 + else + print *, 'PASS: Acov(2) = torflux * rho^2' + end if + + ! hcov(1) should be 0 (no radial component in Boozer) + if (abs(hcov(1)) > 1.0e-15_dp) then + print *, 'FAIL: hcov(1) =', hcov(1), ' expected 0' + nfail = nfail + 1 + else + print *, 'PASS: hcov(1) = 0' + end if + + ! hcov(2) and hcov(3) should be B_theta/Bmod and B_phi/Bmod + if (abs(hcov(2)) < 1.0e-15_dp .or. abs(hcov(3)) < 1.0e-15_dp) then + print *, 'FAIL: hcov(2) or hcov(3) is zero' + nfail = nfail + 1 + else + print *, 'PASS: hcov(2) =', hcov(2), ' hcov(3) =', hcov(3) + end if + + ! Test at another point: rho=0.7, theta=1, phi=0.5 + x = [0.7_dp, 1.0_dp, 0.5_dp] + call field%evaluate(x, Acov, hcov, Bmod) + + if (Bmod <= 0.0_dp) then + print *, 'FAIL: Bmod at second point =', Bmod + nfail = nfail + 1 + else + print *, 'PASS: Bmod at second point =', Bmod + end if + + ! Test integ_to_ref / ref_to_integ coordinate transforms + block + use field_can_boozer, only: integ_to_ref_boozer_chartmap, & + ref_to_integ_boozer_chartmap + real(dp) :: x_ref(3), x_integ(3), x_roundtrip(3), err + + x_ref = [0.7_dp, 2.5_dp, 1.3_dp] ! rho, theta, phi + call ref_to_integ_boozer_chartmap(x_ref, x_integ) + + ! x_integ(1) should be rho^2 = 0.49 + if (abs(x_integ(1) - 0.49_dp) > 1.0e-14_dp) then + print *, 'FAIL: ref_to_integ s =', x_integ(1), ' expected 0.49' + nfail = nfail + 1 + else + print *, 'PASS: ref_to_integ rho->s correct' + end if + + call integ_to_ref_boozer_chartmap(x_integ, x_roundtrip) + err = maxval(abs(x_roundtrip - x_ref)) + if (err > 1.0e-14_dp) then + print *, 'FAIL: coord transform roundtrip error =', err + nfail = nfail + 1 + else + print *, 'PASS: coord transform roundtrip < 1e-14' + end if + end block + + ! Summary + print *, '' + if (nfail == 0) then + print *, 'All tests passed.' + else + print *, nfail, ' tests failed.' + error stop 'Test failures' + end if + +end program test_boozer_chartmap