diff --git a/Makefile b/Makefile index 8555a5ee..d6a6defa 100644 --- a/Makefile +++ b/Makefile @@ -1,7 +1,7 @@ CONFIG ?= Release FLAGS ?= BUILD_DIR := build -GOLDEN_LIBNEO_REF := 18b91aabffcc3c8833648f01a057e278809c15f2 +GOLDEN_LIBNEO_REF := e451e0a42f5228d85b0479de40aa35db8b48239d # Prevent ambient shell env from silently changing which libneo is fetched. # Pass the ref explicitly via: make ... LIBNEO_REF= diff --git a/src/field/field_can_meiss.f90 b/src/field/field_can_meiss.f90 index adca27c9..1254a2c1 100644 --- a/src/field/field_can_meiss.f90 +++ b/src/field/field_can_meiss.f90 @@ -72,7 +72,11 @@ module field_can_meiss type(BatchSplineData3D) :: spl_transform_batch logical :: transform_splines_initialized = .false. -integer, parameter :: order(3) = [5, 5, 5] ! Restored to original 5th order for accuracy +integer, parameter :: default_spline_order = 5 +integer, parameter :: default_ode_step_limit = 1000000 +integer :: order(3) = [default_spline_order, default_spline_order, & + default_spline_order] +integer :: ode_step_limit = default_ode_step_limit logical, parameter :: periodic(3) = [.False., .True., .True.] ! Coordinate scaling for s <-> r transform (default: r = sqrt(s)) @@ -141,7 +145,7 @@ subroutine rh_can_wrapper(r_c, z, dz, context) end subroutine rh_can_wrapper subroutine init_meiss(field_noncan_, n_r_, n_th_, n_phi_, rmin, rmax, thmin, thmax, & - scaling, transformation_relerr_) + scaling, transformation_relerr_, spline_order_, ode_step_limit_) use, intrinsic :: ieee_arithmetic, only: ieee_is_finite use, intrinsic :: iso_fortran_env, only: int64 use new_vmec_stuff_mod, only : nper @@ -151,6 +155,8 @@ subroutine init_meiss(field_noncan_, n_r_, n_th_, n_phi_, rmin, rmax, thmin, thm real(dp), intent(in), optional :: rmin, rmax, thmin, thmax class(coordinate_scaling_t), intent(in), optional :: scaling real(dp), intent(in), optional :: transformation_relerr_ + integer, intent(in), optional :: spline_order_ + integer, intent(in), optional :: ode_step_limit_ call cleanup_meiss call field_clone(field_noncan_, field_noncan) @@ -172,11 +178,19 @@ subroutine init_meiss(field_noncan_, n_r_, n_th_, n_phi_, rmin, rmax, thmin, thm n_th = default_n_th n_phi = default_n_phi transformation_relerr = default_transformation_relerr + order = default_spline_order + ode_step_limit = default_ode_step_limit xmin = default_xmin xmax = [1d0, twopi, twopi/nper] if (present(n_r_)) n_r = n_r_ if (present(n_th_)) n_th = n_th_ if (present(n_phi_)) n_phi = n_phi_ + if (present(spline_order_)) order = spline_order_ + if (present(ode_step_limit_)) ode_step_limit = ode_step_limit_ + if (any(order < 3) .or. any(order > 5)) then + error stop 'canonical map spline order must be between 3 and 5' + end if + if (ode_step_limit < 1) error stop 'canonical map ODE step limit must be positive' if (any([n_r, n_th, n_phi] < order + 1)) then error stop 'canonical map grid dimensions must exceed the spline order' end if @@ -324,11 +338,22 @@ subroutine ref_to_integ_meiss(xref, xinteg) end subroutine ref_to_integ_meiss -subroutine get_meiss_coordinates +subroutine get_meiss_coordinates(ierr) + integer, intent(out), optional :: ierr + integer :: status + #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT print *, 'field_can_meiss.init_transformation' #endif - call init_transformation + call init_transformation(status) + if (status /= 0) then + call cleanup_meiss + if (present(ierr)) then + ierr = status + return + end if + error stop 'Meiss transformation ODE failed' + end if #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT print *, 'field_can_meiss.spline_transformation' @@ -339,12 +364,21 @@ subroutine get_meiss_coordinates print *, 'field_can_meiss.init_canonical_field_components' #endif call init_canonical_field_components + if (present(ierr)) ierr = 0 end subroutine get_meiss_coordinates -subroutine init_transformation +subroutine init_transformation(ierr) + integer, intent(out), optional :: ierr + integer :: status + call init_transformation_arrays() - call compute_transformation() + call compute_transformation(status) + if (present(ierr)) then + ierr = status + else if (status /= 0) then + error stop 'Meiss transformation ODE failed' + end if end subroutine init_transformation subroutine init_transformation_arrays() @@ -358,30 +392,39 @@ subroutine init_transformation_arrays() chi_gauge = 0.0_dp end subroutine init_transformation_arrays -subroutine compute_transformation() +subroutine compute_transformation(ierr) !> Compute transformation data via integration (expensive operation) + integer, intent(out) :: ierr integer :: i_ctr + integer, allocatable :: slice_status(:) logical :: parallel_slices i_ctr = 0 + allocate(slice_status(n_phi), source=0) parallel_slices = .true. select type (field_noncan) type is (spectre_field_t) parallel_slices = .false. end select - !$omp parallel private(i_ctr) if(parallel_slices) - !$omp do - do i_ctr = 1, n_phi - call compute_phi_slice(i_ctr) - enddo - !$omp end do - !$omp end parallel + if (parallel_slices) then + !$omp parallel do private(i_ctr) + do i_ctr = 1, n_phi + call compute_phi_slice(i_ctr, slice_status(i_ctr)) + end do + !$omp end parallel do + else + do i_ctr = 1, n_phi + call compute_phi_slice(i_ctr, slice_status(i_ctr)) + end do + end if + ierr = maxval(slice_status) end subroutine compute_transformation -subroutine compute_phi_slice(i_phi) +subroutine compute_phi_slice(i_phi, ierr) !> Compute transformation for single phi slice integer, intent(in) :: i_phi + integer, intent(out) :: ierr real(dp) :: y(2) integer :: i_th @@ -389,6 +432,7 @@ subroutine compute_phi_slice(i_phi) call print_progress(i_phi) #endif + ierr = 0 do i_th = 1, n_th lam_phi(1, i_th, i_phi) = 0d0 chi_gauge(1, i_th, i_phi) = 0d0 @@ -397,7 +441,8 @@ subroutine compute_phi_slice(i_phi) if (is_hr_zero_on_slice(i_th, i_phi)) then call set_identity_slice(i_th, i_phi) else - call integrate_radial_slice(i_th, i_phi, y) + call integrate_radial_slice(i_th, i_phi, y, ierr) + if (ierr /= 0) return endif enddo end subroutine compute_phi_slice @@ -413,14 +458,21 @@ subroutine set_identity_slice(i_th, i_phi) enddo end subroutine set_identity_slice -subroutine integrate_radial_slice(i_th, i_phi, y) +subroutine integrate_radial_slice(i_th, i_phi, y, ierr) !> Integrate along radial direction for given (i_th, i_phi) integer, intent(in) :: i_th, i_phi real(dp), intent(inout) :: y(2) + integer, intent(out) :: ierr integer :: i_r - + + ierr = 0 do i_r = 2, n_r - call integrate(i_r, i_th, i_phi, y) + call integrate(i_r, i_th, i_phi, y, ierr) + if (ierr /= 0) then + write (*, '(A,I0,A,I0,A,I0)') 'Meiss ODE failed at r index ', i_r, & + ', theta index ', i_th, ', phi index ', i_phi + return + end if enddo end subroutine integrate_radial_slice @@ -442,11 +494,12 @@ subroutine print_progress(i_phi) !$omp end critical end subroutine print_progress -subroutine integrate(i_r, i_th, i_phi, y) +subroutine integrate(i_r, i_th, i_phi, y, ierr) use odeint_allroutines_sub, only: odeint_allroutines integer, intent(in) :: i_r, i_th, i_phi real(dp), dimension(2), intent(inout) :: y + integer, intent(out) :: ierr integer, parameter :: ndim = 2 real(dp) :: r1, r2 @@ -473,7 +526,9 @@ subroutine integrate(i_r, i_th, i_phi, y) #endif context = grid_indices_t(i_th, i_phi) - call odeint_allroutines(y, ndim, context, r1, r2, relaxed_relerr, rh_can_wrapper) + call odeint_allroutines(y, ndim, context, r1, r2, relaxed_relerr, & + rh_can_wrapper, step_limit=ode_step_limit, ierr=ierr) + if (ierr /= 0) return lam_phi(i_r, i_th, i_phi) = y(1) chi_gauge(i_r, i_th, i_phi) = y(2) @@ -599,14 +654,13 @@ end function is_hr_zero_on_slice subroutine init_canonical_field_components real(dp) :: xcan(3) - real(dp), dimension(:,:,:), allocatable :: Ath, Aphi, hth, hphi, Bmod real(dp), dimension(:,:,:,:), allocatable :: y_batch - real(dp) :: xref(3), Acov(3), hcov(3), lam, dlam(3), chi, dchi(3) - integer :: i_r, i_th, i_phi, dims(3) + real(dp) :: xref(3), Acov(3), hcov(3), Bmod, lam, dlam(3), chi, dchi(3) + integer :: i_r, i_th, i_phi - allocate(Ath(n_r,n_th,n_phi), Aphi(n_r,n_th,n_phi)) - allocate(hth(n_r,n_th,n_phi), hphi(n_r,n_th,n_phi)) - allocate(Bmod(n_r,n_th,n_phi)) + if (allocated(lam_phi)) deallocate(lam_phi) + if (allocated(chi_gauge)) deallocate(chi_gauge) + allocate(y_batch(n_r, n_th, n_phi, 5)) do i_phi=1,n_phi do i_th=1,n_th @@ -634,39 +688,32 @@ subroutine init_canonical_field_components real(dp) :: x_integ_local(3), x_ref_local(3), jac_local(3) x_integ_local = xref call coord_scaling%inverse(x_integ_local, x_ref_local, jac_local) - call fld%evaluate(x_ref_local, Acov, hcov, Bmod(i_r, i_th, i_phi)) + call fld%evaluate(x_ref_local, Acov, hcov, Bmod) end block type is (splined_field_t) ! Splined field already in scaled coordinates - call fld%evaluate(xref, Acov, hcov, Bmod(i_r, i_th, i_phi)) + call fld%evaluate(xref, Acov, hcov, Bmod) type is (spectre_field_t) ! SPECTRE chart is its own radial coordinate (identity scaling) call bias_spectre_upper_endpoint(fld, xref) - call fld%evaluate(xref, Acov, hcov, Bmod(i_r, i_th, i_phi)) + call fld%evaluate(xref, Acov, hcov, Bmod) class default error stop "init_canonical_field_components: unsupported field type" end select ! Note: We only use theta/phi components of Acov/hcov below - Ath(i_r, i_th, i_phi) = Acov(2) + Acov(3)*dlam(2) - dchi(2) - Aphi(i_r, i_th, i_phi) = Acov(3)*(1.0d0 + dlam(3)) - dchi(3) - hth(i_r, i_th, i_phi) = hcov(2) + hcov(3)*dlam(2) - hphi(i_r, i_th, i_phi) = hcov(3)*(1.0d0 + dlam(3)) + y_batch(i_r, i_th, i_phi, 1) = & + Acov(2) + Acov(3)*dlam(2) - dchi(2) + y_batch(i_r, i_th, i_phi, 2) = & + Acov(3)*(1.0d0 + dlam(3)) - dchi(3) + y_batch(i_r, i_th, i_phi, 3) = hcov(2) + hcov(3)*dlam(2) + y_batch(i_r, i_th, i_phi, 4) = hcov(3)*(1.0d0 + dlam(3)) + y_batch(i_r, i_th, i_phi, 5) = Bmod end do end do end do - ! Construct batch spline for all 5 field components: [Ath, Aph, hth, hph, Bmod] - dims = shape(Ath) - allocate(y_batch(dims(1), dims(2), dims(3), 5)) - - y_batch(:,:,:,1) = Ath - y_batch(:,:,:,2) = Aphi - y_batch(:,:,:,3) = hth - y_batch(:,:,:,4) = hphi - y_batch(:,:,:,5) = Bmod - call construct_batch_splines_3d(xmin, xmax, y_batch, order, periodic, spl_field_batch) batch_splines_initialized = .true. end subroutine init_canonical_field_components diff --git a/src/field/field_can_spectre.f90 b/src/field/field_can_spectre.f90 index afa7933f..cf127f54 100644 --- a/src/field/field_can_spectre.f90 +++ b/src/field/field_can_spectre.f90 @@ -35,14 +35,15 @@ module field_can_spectre spectre_volumes, spectre_mvol, set_spectre_volume_lock, & set_spectre_construction_grid - !> Per-volume Meiss construction grid. r, th default to the historical hardcoded - !> 48; phi = -1 selects the automatic policy below. set_spectre_construction_grid + !> Per-volume Meiss construction grid. phi = -1 selects the automatic policy + !> below. set_spectre_construction_grid !> threads the namelist knob (params: spectre_ncon_r/th/phi) in before the build, !> kept as a setter because field_can_* cannot use params (dependency cycle). - integer :: ncon_r = 48, ncon_th = 48, ncon_phi = -1 + integer :: ncon_r = 48, ncon_th = 48, ncon_phi = -1, ncon_order = 5 + integer :: ncon_ode_step_limit = 1000000 + real(dp) :: ncon_ode_relerr = 1.0e-2_dp - !> Automatic phi resolution. NPHI_FULL is the historical hardcoded value, used - !> for fields carrying toroidal harmonics (bit-identical). Axisymmetric SPECTRE + !> Automatic phi resolution. Axisymmetric SPECTRE !> fields (all n = 0) are phi-invariant, so the construction only needs a minimal !> periodic-quintic phi grid; AXISYM_NPHI is the dominant memory saver for tokamak !> cases (arrays and splines scale linearly in n_phi). It stays above the quintic @@ -103,23 +104,34 @@ integer function active_volume(r) result(lvol) end if end function active_volume - subroutine set_spectre_construction_grid(n_r, n_th, n_phi) + subroutine set_spectre_construction_grid(n_r, n_th, n_phi, spline_order, & + ode_step_limit, ode_relerr) !> Set the per-volume Meiss construction grid before the build (n_phi <= 0 !> selects the automatic phi policy). Called from init_spectre_field. integer, intent(in) :: n_r, n_th, n_phi + integer, intent(in), optional :: spline_order + integer, intent(in), optional :: ode_step_limit + real(dp), intent(in), optional :: ode_relerr ncon_r = n_r ncon_th = n_th ncon_phi = n_phi + ncon_order = 5 + if (present(spline_order)) ncon_order = spline_order + ncon_ode_step_limit = 1000000 + if (present(ode_step_limit)) ncon_ode_step_limit = ode_step_limit + ncon_ode_relerr = 1.0e-2_dp + if (present(ode_relerr)) ncon_ode_relerr = ode_relerr end subroutine set_spectre_construction_grid - subroutine construct_spectre_coordinates(field_noncan) + subroutine construct_spectre_coordinates(field_noncan, ierr) !> Build one Meiss chart per SPECTRE volume on r in [lvol-1, lvol]. The axis !> volume starts at AXIS_OFFSET so the construction never touches the r = 0 !> coordinate singularity. class(magnetic_field_t), intent(in) :: field_noncan + integer, intent(out), optional :: ierr - integer :: lvol, nphi_eff + integer :: lvol, nphi_eff, meiss_ierr real(dp) :: rmin, rmax logical :: axisym @@ -141,10 +153,11 @@ subroutine construct_spectre_coordinates(field_noncan) nphi_eff = NPHI_FULL end if - print '(A,I0,A,I0,A,L1,A,I0,A,L1)', & + print '(A,I0,A,I0,A,L1,A,I0,A,L1,A,I0,A,I0,A,ES9.2)', & ' spectre_construction_grid: n_r=', ncon_r, ' n_th=', ncon_th, & ' axisym=', axisym, ' n_phi=', nphi_eff, & - ' phi_auto=', ncon_phi <= 0 + ' phi_auto=', ncon_phi <= 0, ' order=', ncon_order, & + ' ode_steps=', ncon_ode_step_limit, ' ode_relerr=', ncon_ode_relerr allocate (spectre_volumes(spectre_mvol)) @@ -154,12 +167,24 @@ subroutine construct_spectre_coordinates(field_noncan) rmax = real(lvol, dp) call init_meiss(field_noncan, ncon_r, ncon_th, nphi_eff, rmin, rmax, & - 0.0_dp, twopi, identity_scaling_t()) - call get_meiss_coordinates + 0.0_dp, twopi, identity_scaling_t(), & + transformation_relerr_=ncon_ode_relerr, & + spline_order_=ncon_order, & + ode_step_limit_=ncon_ode_step_limit) + call get_meiss_coordinates(meiss_ierr) + if (meiss_ierr /= 0) then + call cleanup_spectre + if (present(ierr)) then + ierr = lvol + return + end if + error stop 'SPECTRE canonical construction ODE failed' + end if call harvest_meiss_volume(spectre_volumes(lvol)) end do call cleanup_meiss + if (present(ierr)) ierr = 0 end subroutine construct_spectre_coordinates diff --git a/src/params.f90 b/src/params.f90 index 1ca156be..3ec9b21f 100644 --- a/src/params.f90 +++ b/src/params.f90 @@ -146,17 +146,18 @@ module params logical :: spectre_sbeg_is_toroidal_flux = .false. !> SPECTRE per-volume Meiss construction grid (n_r, n_th, n_phi). Each volume - !> allocates rank-3 arrays plus quintic batch splines of this size, so the - !> total peak memory scales as Mvol*n_r*n_th*n_phi. spectre_ncon_r/th default - !> to the historical hardcoded 48 (bit-identical). spectre_ncon_phi = -1 means - !> automatic: the historical 32 for fields with toroidal harmonics - !> (bit-identical), and a minimal phi grid for axisymmetric fields (all n = 0), + !> allocates rank-3 arrays plus batch splines of this size, so the + !> total peak memory scales as Mvol*n_r*n_th*n_phi. spectre_ncon_phi = -1 means + !> automatic: 32 for fields with toroidal harmonics and a minimal phi grid + !> for axisymmetric fields (all n = 0), !> the dominant memory saver for tokamak cases. A positive value forces that - !> phi count verbatim, bypassing the axisymmetric clamp (raise it back to 32 - !> when high-order symplectic convergence needs the full phi resolution). + !> phi count verbatim, bypassing the axisymmetric clamp. integer :: spectre_ncon_r = 48 integer :: spectre_ncon_th = 48 integer :: spectre_ncon_phi = -1 + integer :: spectre_ncon_order = 5 + integer :: spectre_ncon_ode_max_steps = 1000000 + real(dp) :: spectre_ncon_ode_relerr = 1.0e-2_dp namelist /config/ notrace_passing, nper, npoiper, ntimstep, ntestpart, & trace_time, num_surf, sbeg, phibeg, thetabeg, contr_pp, & @@ -176,7 +177,8 @@ module params batch_size, ran_seed, reuse_batch, field_input, coord_input, & wall_input, wall_units, integ_coords, crossing_level, & spectre_sbeg_is_toroidal_flux, spectre_ncon_r, spectre_ncon_th, & - spectre_ncon_phi, & + spectre_ncon_phi, spectre_ncon_order, spectre_ncon_ode_max_steps, & + spectre_ncon_ode_relerr, & output_results_netcdf, & output_error, output_orbits_macrostep, & ! callback macrostep_time_grid, checkpoint_interval, restart diff --git a/src/samplers.f90 b/src/samplers.f90 index d66b69bb..befb8e4b 100644 --- a/src/samplers.f90 +++ b/src/samplers.f90 @@ -204,6 +204,48 @@ subroutine sample_spectre_surface(zstart) call save_starting_points(zstart) end subroutine sample_spectre_surface + subroutine init_spectre_start_bounds(zstart) + use params, only: bmin, bmax, bmod00 + use magfie_sub, only: magfie + + real(dp), dimension(:, :), intent(in) :: zstart + + integer, parameter :: nrho = 8, ntheta = 64, nzeta = 64 + integer :: ir, it, iz, nr + real(dp) :: rho_lo, rho_hi, rho, theta, zeta, bmod, bsum + real(dp) :: sqrtg, bder(3), hcovar(3), hctrvr(3), hcurl(3) + + rho_lo = minval(zstart(1, :)) + rho_hi = maxval(zstart(1, :)) + nr = merge(1, nrho, abs(rho_hi - rho_lo) <= & + epsilon(1.0d0)*max(1.0d0, abs(rho_lo), abs(rho_hi))) + bmin = huge(1.0d0) + bmax = -huge(1.0d0) + bsum = 0.0d0 + do ir = 1, nr + if (nr == 1) then + rho = rho_lo + else + rho = rho_lo + (rho_hi - rho_lo)*real(ir - 1, dp)/real(nr - 1, dp) + end if + do it = 1, ntheta + theta = twopi*(real(it, dp) - 0.5d0)/real(ntheta, dp) + do iz = 1, nzeta + zeta = twopi*(real(iz, dp) - 0.5d0)/real(nzeta, dp) + call magfie([rho, theta, zeta], bmod, sqrtg, bder, hcovar, & + hctrvr, hcurl) + bmin = min(bmin, bmod) + bmax = max(bmax, bmod) + bsum = bsum + bmod + end do + end do + end do + bmod00 = bsum/real(nr*ntheta*nzeta, dp) + + print *, 'init_spectre_start_bounds: rho = ', rho_lo, rho_hi, & + ' bmod00 = ', bmod00, ' bmin = ', bmin, ' bmax = ', bmax + end subroutine init_spectre_start_bounds + subroutine sample_surface_fieldline(zstart) real(dp), dimension(:, :), intent(inout) :: zstart diff --git a/src/simple_main.f90 b/src/simple_main.f90 index 5f5fca4e..79f318ab 100644 --- a/src/simple_main.f90 +++ b/src/simple_main.f90 @@ -59,7 +59,8 @@ subroutine main ntestpart, ntimstep, coord_input, restart use timing, only: init_timer, print_phase_time use magfie_sub, only: TEST, VMEC, SPECTRE, init_magfie - use samplers, only: init_starting_surf, sample_spectre_surface + use samplers, only: init_starting_surf, sample_spectre_surface, & + init_spectre_start_bounds use version, only: simple_version use field_boozer_chartmap, only: is_boozer_chartmap use field, only: is_spectre_file @@ -129,6 +130,7 @@ subroutine main call sample_spectre_surface(zstart) else call sample_particles(.true.) + if (startmode == 2) call init_spectre_start_bounds(zstart) end if call print_phase_time('SPECTRE particle sampling completed') @@ -377,7 +379,9 @@ subroutine init_spectre_field(self) use new_vmec_stuff_mod, only: nper, rmajor use util, only: twopi use params, only: field_input, spectre_ncon_r, spectre_ncon_th, & - spectre_ncon_phi + spectre_ncon_phi, spectre_ncon_order, & + spectre_ncon_ode_max_steps + use params, only: spectre_ncon_ode_relerr use timing, only: print_phase_time use orbit_symplectic_base, only: sympl_rmax @@ -423,7 +427,9 @@ subroutine init_spectre_field(self) ! extended linearly, so iterates out there stay finite. sympl_rmax = real(sf%data%Mvol + 1, dp) call set_spectre_construction_grid(spectre_ncon_r, spectre_ncon_th, & - spectre_ncon_phi) + spectre_ncon_phi, spectre_ncon_order, & + spectre_ncon_ode_max_steps, & + spectre_ncon_ode_relerr) call init_field_can(SPECTRE, sf) call print_phase_time('SPECTRE per-volume canonical construction completed') end if diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt index 0c8a2612..df647159 100644 --- a/test/tests/CMakeLists.txt +++ b/test/tests/CMakeLists.txt @@ -886,6 +886,7 @@ add_test(NAME test_spectre_sympl_volume WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) set_tests_properties(test_spectre_sympl_volume PROPERTIES LABELS "integration;python" + RUN_SERIAL TRUE TIMEOUT 300) # SPECTRE Level-0 interface crossing map with event location (#443). tok2vol is @@ -926,6 +927,7 @@ add_test(NAME test_spectre_sympl_crossing WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) set_tests_properties(test_spectre_sympl_crossing PROPERTIES LABELS "integration;python" + RUN_SERIAL TRUE TIMEOUT 300) # SPECTRE construction-grid memory knob and axisymmetric phi auto-clamp (#442): @@ -952,6 +954,7 @@ add_test(NAME test_spectre_validation WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) set_tests_properties(test_spectre_validation PROPERTIES LABELS "integration;python" + RUN_SERIAL TRUE TIMEOUT 300) add_subdirectory(field_can) diff --git a/test/tests/field_can/test_field_can_meiss.f90 b/test/tests/field_can/test_field_can_meiss.f90 index a7f113c5..02821f8c 100644 --- a/test/tests/field_can/test_field_can_meiss.f90 +++ b/test/tests/field_can/test_field_can_meiss.f90 @@ -12,8 +12,9 @@ program test_field_can_meiss use field_can_meiss, only: init_meiss, get_meiss_coordinates, & spline_transformation, init_canonical_field_components, & xmin, xmax, h_r, h_phi, h_th, ah_cov_on_slice, n_r, n_phi, n_th, & - lam_phi, chi_gauge, transformation_relerr, ref_to_integ_meiss, & + spl_transform_batch, transformation_relerr, ref_to_integ_meiss, & integ_to_ref_meiss + use interpolate, only: evaluate_batch_splines_3d use params, only: canonical_grid_nr, canonical_grid_ntheta, & canonical_grid_nphi, canonical_ode_relerr, field_input, coord_input, & read_config @@ -153,6 +154,7 @@ subroutine write_transformation(filename) integer :: funit integer :: i_r, i_th, i_phi real(dp) :: r, th, phi + real(dp) :: transformation(2) open (newunit=funit, file=filename, status='unknown') write (funit, *) '#', ' r', ' phi', ' th', ' lam_phi', ' chi_gauge' @@ -163,8 +165,9 @@ subroutine write_transformation(filename) phi = xmin(3) + h_phi*(i_phi - 1) do i_r = 1, n_r r = xmin(1) + h_r*(i_r - 1) - write (funit, *) r, phi, th, lam_phi(i_r, i_th, i_phi), & - chi_gauge(i_r, i_th, i_phi) + call evaluate_batch_splines_3d(spl_transform_batch, & + [r, th, phi], transformation) + write (funit, *) r, phi, th, transformation end do end do end do diff --git a/test/tests/test_spectre_construction_grid.py b/test/tests/test_spectre_construction_grid.py index e2ad7f9f..261b96a6 100644 --- a/test/tests/test_spectre_construction_grid.py +++ b/test/tests/test_spectre_construction_grid.py @@ -30,6 +30,7 @@ # Construction identity: confined-marker p^2 drift. The coarse-grid Meiss chart # stays at ~2.5e-10 (measured); a broken construction drifts to O(0.1). IDENTITY_TOL = 1.0e-6 +ORDER3_IDENTITY_TOL = 1.0e-4 # Clamp must cut phi well below the historical 32 to save the phi memory factor. CLAMP_PHI_MAX = 16 # The clamp must cut peak RSS to a clear fraction of the unclamped phi=32 build. @@ -98,12 +99,14 @@ def run_config(binary, h5, integmode=3, ntestpart=2, trace_time=1.0e-6, def parse_grid(out): m = re.search(r"spectre_construction_grid: n_r=(\d+) n_th=(\d+) " - r"axisym=([TF]) n_phi=(\d+) phi_auto=([TF])", out) + r"axisym=([TF]) n_phi=(\d+) phi_auto=([TF]) order=(\d+) " + r"ode_steps=(\d+) ode_relerr=\s*([0-9.E+-]+)", out) if not m: raise RuntimeError("spectre_construction_grid line missing from stdout") return dict(n_r=int(m.group(1)), n_th=int(m.group(2)), axisym=m.group(3) == "T", n_phi=int(m.group(4)), - phi_auto=m.group(5) == "T") + phi_auto=m.group(5) == "T", order=int(m.group(6)), + ode_steps=int(m.group(7)), ode_relerr=float(m.group(8))) def main(): @@ -113,6 +116,10 @@ def main(): # Default: axisymmetric tok2vol auto-clamps the phi grid. rss_clamped, out_default, _ = run_config(binary, tok2vol) grid = parse_grid(out_default) + if (grid["n_r"], grid["n_th"], grid["order"]) != (48, 48, 5): + failures.append(f"defaults: unexpected construction settings {grid}") + if not np.isclose(grid["ode_relerr"], 1.0e-2): + failures.append(f"defaults: ODE tolerance {grid['ode_relerr']:.3e}") if not grid["axisym"]: failures.append("clamp: tok2vol not detected as axisymmetric") if not grid["phi_auto"]: @@ -126,7 +133,8 @@ def main(): # Red proof: forcing spectre_ncon_phi = 32 bypasses the clamp and restores # the historical full-resolution build; peak RSS jumps by ~the phi ratio. rss_full, out_full, _ = run_config(binary, tok2vol, - extra=("spectre_ncon_phi = 32",)) + extra=("spectre_ncon_phi = 32", + "spectre_ncon_order = 5")) grid_full = parse_grid(out_full) if grid_full["phi_auto"] or grid_full["n_phi"] != 32: failures.append(f"red proof: explicit phi=32 not honored " @@ -148,6 +156,17 @@ def main(): print(f"memory: coarse(16x16, phi={parse_grid(out_coarse)['n_phi']})=" f"{rss_coarse:.0f} MB (default {rss_clamped:.0f} MB)") + rss_cubic, out_cubic, _ = run_config( + binary, tok2vol, + extra=("spectre_ncon_phi = 32", "spectre_ncon_order = 3")) + if parse_grid(out_cubic)["order"] != 3: + failures.append("order: cubic construction setting not honored") + if not rss_cubic < rss_full: + failures.append(f"order: cubic RSS {rss_cubic:.0f} MB not below " + f"quintic {rss_full:.0f} MB") + print(f"memory: cubic({grid['n_r']}x{grid['n_th']}, phi=32)=" + f"{rss_cubic:.0f} MB (quintic {rss_full:.0f} MB)") + # Construction identity: a coarse-grid chart still conserves the # guiding-center energy p^2 for confined low-energy markers. _, out_id, data = run_config( @@ -165,6 +184,23 @@ def main(): print(f"identity: coarse-grid confined={int(confined.sum())} " f"max|p^2-1|={drift:.3e}") + _, _, data_cubic = run_config( + binary, tok2vol, ntestpart=16, trace_time=1.0e-5, ntimstep=50, + npoiper2=1024, relerr="1d-12", face_al=50.0, + extra=("spectre_ncon_r = 16", "spectre_ncon_th = 16", + "spectre_ncon_order = 3")) + confined_cubic = data_cubic[:, 1] >= 1.0e-5 * (1.0 - 1e-9) + if confined_cubic.sum() < 8: + failures.append(f"order identity: only {int(confined_cubic.sum())} " + "confined markers") + else: + drift_cubic = float(np.max(np.abs(data_cubic[confined_cubic, 8] ** 2 - 1.0))) + if not drift_cubic < ORDER3_IDENTITY_TOL: + failures.append(f"order identity: cubic p^2 drift {drift_cubic:.3e} " + f">= {ORDER3_IDENTITY_TOL:.0e}") + print(f"identity: cubic confined={int(confined_cubic.sum())} " + f"max|p^2-1|={drift_cubic:.3e}") + if failures: print("\nFAILURES:") for f in failures: diff --git a/test/tests/test_spectre_freeboundary.f90 b/test/tests/test_spectre_freeboundary.f90 index a63e7c9f..af31f192 100644 --- a/test/tests/test_spectre_freeboundary.f90 +++ b/test/tests/test_spectre_freeboundary.f90 @@ -3,7 +3,9 @@ program test_spectre_freeboundary use, intrinsic :: ieee_arithmetic, only: ieee_is_finite use field_spectre, only: spectre_field_t, create_spectre_field use field_can_spectre, only: construct_spectre_coordinates, & - set_spectre_construction_grid, spectre_mvol, cleanup_spectre + set_spectre_construction_grid, spectre_mvol, cleanup_spectre, & + evaluate_spectre_can + use field_can_base, only: field_can_t use interface_crossing, only: apply_crossing, crossing_info_t, & CROSSING_LEVEL0, CROSS_CROSSING, CROSS_LOSS use magfie_sub, only: set_magfie_spectre_field, init_magfie, SPECTRE, M_TO_CM @@ -26,10 +28,14 @@ program test_spectre_freeboundary call set_magfie_spectre_field(field) call init_magfie(SPECTRE) - call set_spectre_construction_grid(8, 8, 8) - call construct_spectre_coordinates(field) + call set_spectre_construction_grid(8, 8, 8, 5, 1, 1.0e-4_dp) + call construct_spectre_coordinates(field, ierr) + if (ierr == 0) error stop 'ODE step-limit failure was not catchable' + call set_spectre_construction_grid(8, 8, 8, 5, 1000000, 1.0e-4_dp) + call construct_spectre_coordinates(field, ierr) + if (ierr /= 0) error stop 'canonical construction failed' if (spectre_mvol /= 3) error stop 'canonical construction omitted vacuum volume' - call cleanup_spectre + call check_repeated_construction(field) call check_field_and_gc(field) call check_vacuum_fo(field) @@ -39,6 +45,34 @@ program test_spectre_freeboundary contains + subroutine check_repeated_construction(field) + type(spectre_field_t), intent(in) :: field + type(field_can_t) :: first, current + integer :: i, failures + + call evaluate_spectre_can(first, 1.5_dp, 0.7_dp, 0.3_dp, 0) + call cleanup_spectre + call construct_spectre_coordinates(field, ierr) + if (ierr /= 0) error stop 'repeated canonical construction failed' + call evaluate_spectre_can(current, 1.5_dp, 0.7_dp, 0.3_dp, 0) + if (first%Bmod /= current%Bmod .or. first%Ath /= current%Ath .or. & + first%Aph /= current%Aph .or. any(first%dAth /= current%dAth) .or. & + any(first%dAph /= current%dAph)) then + error stop 'repeated canonical construction changed the chart' + end if + + failures = 0 + !$omp parallel do private(current) reduction(+:failures) + do i = 1, 64 + call evaluate_spectre_can(current, 1.5_dp, 0.7_dp, 0.3_dp, 0) + if (.not. ieee_is_finite(current%Bmod) .or. & + current%Bmod /= first%Bmod) failures = failures + 1 + end do + !$omp end parallel do + if (failures /= 0) error stop 'threaded canonical evaluation changed the field' + call cleanup_spectre + end subroutine check_repeated_construction + subroutine check_field_and_gc(field) type(spectre_field_t), intent(in) :: field type(crossing_info_t) :: info diff --git a/test/tests/test_spectre_startmode2.py b/test/tests/test_spectre_startmode2.py index 5768f86a..0def2ad7 100644 --- a/test/tests/test_spectre_startmode2.py +++ b/test/tests/test_spectre_startmode2.py @@ -16,14 +16,18 @@ ]) -def config(field): +def config(field, generate_start_only): return f"""&config field_input = '{field}' integ_coords = 6 integmode = 0 ntestpart = {len(STARTS)} startmode = 2 -generate_start_only = .True. +generate_start_only = {'.True.' if generate_start_only else '.False.'} +trace_time = 1d-7 +ntimstep = 5 +npoiper2 = 64 +relerr = 1d-8 spectre_ncon_r = 8 spectre_ncon_th = 8 spectre_ncon_phi = 8 @@ -31,12 +35,12 @@ def config(field): """ -def main(): - binary, field = map(Path, sys.argv[1:]) +def run(binary, field, generate_start_only): with tempfile.TemporaryDirectory() as tmp: work = Path(tmp) np.savetxt(work / "start.dat", STARTS, fmt="%.17g") - (work / "simple.in").write_text(config(field.resolve())) + (work / "simple.in").write_text( + config(field.resolve(), generate_start_only)) proc = subprocess.run([binary.resolve(), "simple.in"], cwd=work, capture_output=True, text=True, timeout=120) if proc.returncode != 0: @@ -44,7 +48,18 @@ def main(): actual = np.loadtxt(work / "start.dat") if not np.array_equal(actual, STARTS): raise AssertionError(f"explicit starts were replaced:\n{actual}") - print("SPECTRE startmode=2 explicit starts PASS") + if not generate_start_only: + times_lost = np.atleast_2d(np.loadtxt(work / "times_lost.dat")) + if not np.all(np.isfinite(times_lost[:, 2])): + raise AssertionError( + f"non-finite trapped-passing diagnostics:\n{times_lost}") + + +def main(): + binary, field = map(Path, sys.argv[1:]) + run(binary, field, True) + run(binary, field, False) + print("SPECTRE startmode=2 explicit starts and diagnostics PASS") if __name__ == "__main__":