feat: add SIMPLE Boozer chartmap core path

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

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

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