diff --git a/Makefile b/Makefile index cd71e26a..8555a5ee 100644 --- a/Makefile +++ b/Makefile @@ -1,7 +1,7 @@ CONFIG ?= Release FLAGS ?= BUILD_DIR := build -GOLDEN_LIBNEO_REF := cbb57e125390cd21ea906cc0ec9bff5b28c6fc26 +GOLDEN_LIBNEO_REF := 18b91aabffcc3c8833648f01a057e278809c15f2 # Prevent ambient shell env from silently changing which libneo is fetched. # Pass the ref explicitly via: make ... LIBNEO_REF= diff --git a/app/spectre_poincare.f90 b/app/spectre_poincare.f90 index f44a1283..919a907d 100644 --- a/app/spectre_poincare.f90 +++ b/app/spectre_poincare.f90 @@ -4,15 +4,16 @@ program spectre_poincare !> !> Reads a namelist (default spectre_poincare.in), loads the SPECTRE HDF5 !> equilibrium, seeds field lines per volume at theta = 0.2 with s spread - !> across the volume, and integrates each line with the semi-implicit - !> symplectic Euler map in field_line_spectre. Field periodicity puts every + !> across the volume, and integrates each line with the implicit midpoint + !> map in field_line_spectre. Field periodicity puts every !> section exactly on zeta0 + k*2*pi/Nfp, so no interpolation is needed. !> !> Output per volume: poincare_vol.dat (seed, section, s, theta, R, Z) !> and iota_vol.dat (seed, fitted rotational transform). use, intrinsic :: iso_fortran_env, only: dp => real64 use spectre_reader, only: spectre_data_t, load_spectre, free_spectre - use field_line_spectre, only: field_line_psi, field_line_step, spectre_rz + use field_line_spectre, only: field_line_psi, field_line_step, spectre_rz, & + prepare_field_line_spline, cleanup_field_line_spline implicit none @@ -23,12 +24,14 @@ program spectre_poincare real(dp), parameter :: s_terminate = 1.0_dp - 1.0e-10_dp integer, parameter :: max_volumes = 128 - character(len=1024) :: spectre_file + character(len=1024) :: spectre_file, seed_file real(dp) :: zeta0 integer :: nturns, nsteps_per_period, nseeds_per_volume + integer :: spline_n_s, spline_n_theta, spline_n_zeta, spline_order integer :: volumes(max_volumes) namelist /poincare/ spectre_file, zeta0, nturns, & - nsteps_per_period, nseeds_per_volume, volumes + nsteps_per_period, nseeds_per_volume, volumes, seed_file, spline_n_s, & + spline_n_theta, spline_n_zeta, spline_order type(spectre_data_t) :: data character(len=1024) :: config_file @@ -37,10 +40,15 @@ program spectre_poincare real(dp) :: h spectre_file = 'spectre.h5' + seed_file = '' zeta0 = 0.0_dp nturns = 200 nsteps_per_period = 256 nseeds_per_volume = 4 + spline_n_s = 32 + spline_n_theta = 32 + spline_n_zeta = 32 + spline_order = 3 volumes = -1 config_file = 'spectre_poincare.in' @@ -140,11 +148,14 @@ subroutine trace_volume(data, lvol, h, n_terminated) integer :: nsec, iseed, ksec, istep, uposec, uiota real(dp) :: s, theta, zeta, psi, R, Z, iota - real(dp), allocatable :: zeta_sec(:), theta_sec(:) + real(dp), allocatable :: zeta_sec(:), theta_sec(:), seeds(:, :) character(len=64) :: poi_name, iota_name logical :: terminated nsec = nturns*data%Nfp + call volume_seeds(lvol, seeds) + call prepare_field_line_spline(data, lvol, spline_n_s, spline_n_theta, & + spline_n_zeta, spline_order) allocate (zeta_sec(0:nsec), theta_sec(0:nsec)) write (poi_name, '(A,I0,A)') 'poincare_vol', lvol, '.dat' @@ -154,9 +165,9 @@ subroutine trace_volume(data, lvol, h, n_terminated) write (uposec, '(A)') '# seed section s theta R Z' write (uiota, '(A)') '# seed iota' - do iseed = 1, nseeds_per_volume - s = seed_s(iseed, nseeds_per_volume) - theta = seed_theta + do iseed = 1, size(seeds, 1) + s = seeds(iseed, 1) + theta = seeds(iseed, 2) zeta = zeta0 psi = field_line_psi(data, lvol, s, theta, zeta) terminated = .false. @@ -193,8 +204,45 @@ subroutine trace_volume(data, lvol, h, n_terminated) close (uposec) close (uiota) deallocate (zeta_sec, theta_sec) + call cleanup_field_line_spline end subroutine trace_volume + subroutine volume_seeds(lvol, seeds) + integer, intent(in) :: lvol + real(dp), allocatable, intent(out) :: seeds(:, :) + + integer :: unit, ios, ivol, count, i + real(dp) :: sval, theta + + if (len_trim(seed_file) == 0) then + allocate (seeds(nseeds_per_volume, 2)) + do i = 1, nseeds_per_volume + seeds(i, :) = [seed_s(i, nseeds_per_volume), seed_theta] + end do + return + end if + + open (newunit=unit, file=trim(seed_file), status='old', action='read') + count = 0 + do + read (unit, *, iostat=ios) ivol, sval, theta + if (ios /= 0) exit + if (ivol == lvol - 1) count = count + 1 + end do + rewind (unit) + allocate (seeds(count, 2)) + i = 0 + do + read (unit, *, iostat=ios) ivol, sval, theta + if (ios /= 0) exit + if (ivol == lvol - 1) then + i = i + 1 + seeds(i, :) = [sval, theta] + end if + end do + close (unit) + end subroutine volume_seeds + subroutine write_section(unit, iseed, ksec, s, theta, R, Z) integer, intent(in) :: unit, iseed, ksec real(dp), intent(in) :: s, theta, R, Z diff --git a/src/field_line_spectre.f90 b/src/field_line_spectre.f90 index 59602f44..de70e768 100644 --- a/src/field_line_spectre.f90 +++ b/src/field_line_spectre.f90 @@ -9,20 +9,28 @@ module field_line_spectre !> dtheta/dzeta = -dA_zeta/ds / dA_theta/ds !> dpsi/dzeta = dA_zeta/dtheta - dA_theta/dtheta * dA_zeta/ds/dA_theta/ds !> (the explicit d/dzeta terms cancel exactly). One step advances the pair - !> (theta, psi) by a semi-implicit symplectic Euler map. Both update signs + !> (theta, psi) by an implicit midpoint map. Both update signs !> set the field-line orientation, so a sign flip reverses iota; keep them !> as written. use, intrinsic :: iso_fortran_env, only: dp => real64 use spectre_reader, only: spectre_data_t use spectre_basis, only: spectre_vecpot_t, eval_spectre_vector_potential + use interpolate, only: BatchSplineData3D, construct_batch_splines_3d, & + destroy_batch_splines_3d, & + evaluate_batch_splines_3d_der implicit none private public :: field_line_psi, field_line_step, spectre_rz + public :: prepare_field_line_spline, cleanup_field_line_spline integer, parameter :: max_newton = 20 real(dp), parameter :: newton_tol = 1.0e-13_dp + real(dp), parameter :: fd_scale = epsilon(1.0_dp)**(1.0_dp/3.0_dp) + real(dp), parameter :: twopi = 6.283185307179586476925286766559_dp + type(BatchSplineData3D) :: hotloop_spline + integer :: hotloop_lvol = 0 contains @@ -35,16 +43,62 @@ pure function field_line_psi(data, lvol, s, theta, zeta) result(psi) type(spectre_vecpot_t) :: av - call eval_spectre_vector_potential(data, lvol, s, theta, zeta, av) + call eval_spectre_vector_potential(data, lvol, s, theta, zeta, av, & + second_derivatives=.false.) psi = av%Ath end function field_line_psi - !> Advance (s, theta, zeta, psi) by one symplectic Euler step of size h. - !> The radial coordinate s is solved by Newton from - !> F(s1) = A_theta(s1,theta0,zeta0) - psi0 - h*dpsi/dzeta(s1,theta0,zeta0) - !> with the analytic F' built from the ss and st second derivatives; theta - !> and psi then advance explicitly. s carries the previous radius as the - !> Newton seed and returns the solved radius. + subroutine prepare_field_line_spline(data, lvol, n_s, n_theta, n_zeta, order) + type(spectre_data_t), intent(in) :: data + integer, intent(in) :: lvol, n_s, n_theta, n_zeta, order + + type(spectre_vecpot_t) :: av + real(dp), allocatable :: values(:, :, :, :) + real(dp) :: xmin(3), xmax(3), x(3) + integer :: i, j, k + + call cleanup_field_line_spline + if (min(n_s, n_theta, n_zeta) <= order) then + error stop 'field-line spline grid must exceed its order' + end if + xmin = [-1.0_dp, 0.0_dp, 0.0_dp] + xmax = [1.0_dp, twopi, twopi/real(data%Nfp, dp)] + allocate (values(n_s, n_theta, n_zeta, 2)) + do k = 1, n_zeta + x(3) = grid_coordinate(k, n_zeta, xmin(3), xmax(3)) + do j = 1, n_theta + x(2) = grid_coordinate(j, n_theta, xmin(2), xmax(2)) + do i = 1, n_s + x(1) = grid_coordinate(i, n_s, xmin(1), xmax(1)) + call eval_spectre_vector_potential(data, lvol, x(1), x(2), & + x(3), av, & + second_derivatives=.false.) + values(i, j, k, :) = [av%Ath, av%Azt] + end do + end do + end do + call construct_batch_splines_3d(xmin, xmax, values, [order, order, order], & + [.false., .true., .true.], hotloop_spline) + hotloop_lvol = lvol + end subroutine prepare_field_line_spline + + subroutine cleanup_field_line_spline + if (hotloop_lvol /= 0) call destroy_batch_splines_3d(hotloop_spline) + hotloop_lvol = 0 + end subroutine cleanup_field_line_spline + + pure function grid_coordinate(i, n, xmin, xmax) result(x) + integer, intent(in) :: i, n + real(dp), intent(in) :: xmin, xmax + real(dp) :: x + + x = xmin + (xmax - xmin)*real(i - 1, dp)/real(n - 1, dp) + end function grid_coordinate + + !> Advance the canonical pair (theta, psi) by implicit midpoint. The two + !> midpoint residuals are solved in (s_mid, theta_mid); the endpoint radius + !> is then recovered from A_theta(s1,theta1,zeta1) = psi1. Numerical + !> Jacobians use only values and first derivatives of the vector potential. subroutine field_line_step(data, lvol, h, s, theta, zeta, psi) type(spectre_data_t), intent(in) :: data integer, intent(in) :: lvol @@ -52,46 +106,150 @@ subroutine field_line_step(data, lvol, h, s, theta, zeta, psi) real(dp), intent(inout) :: s, theta, zeta, psi type(spectre_vecpot_t) :: av - real(dp) :: theta0, zeta0, psi0, s1, ratio, flow, fun, dfun, dgds - integer :: iter - logical :: converged + real(dp) :: s0, theta0, zeta0, psi0, smid, theta_delta + logical :: midpoint_ok, radius_ok + s0 = s theta0 = theta zeta0 = zeta psi0 = psi - s1 = s - ratio = 0.0_dp - fun = 0.0_dp + smid = s0 + theta_delta = 0.0_dp + call solve_midpoint(data, lvol, h, s0, theta0, zeta0, psi0, smid, & + theta_delta, av, midpoint_ok) + if (.not. midpoint_ok) then + write (*, '(A,I0,A,2(ES13.5,1X))') & + 'field_line_step: midpoint Newton failed in volume ', lvol, & + ', (s, theta) = ', smid, theta0 + theta_delta + error stop 'field_line_step: midpoint Newton did not converge' + end if + + theta = theta0 + 2.0_dp*theta_delta + psi = 2.0_dp*av%Ath - psi0 + zeta = zeta0 + h + s = 2.0_dp*smid - s0 + call recover_radius(data, lvol, theta, zeta, psi, s, radius_ok) + if (.not. radius_ok) error stop 'field_line_step: endpoint radius failed' + end subroutine field_line_step + + subroutine solve_midpoint(data, lvol, h, s0, theta0, zeta0, psi0, & + s, theta_delta, av, converged) + type(spectre_data_t), intent(in) :: data + integer, intent(in) :: lvol + real(dp), intent(in) :: h, s0, theta0, zeta0, psi0 + real(dp), intent(inout) :: s, theta_delta + type(spectre_vecpot_t), intent(out) :: av + logical, intent(out) :: converged + + real(dp) :: residual(2), jac(2, 2), det, ds, dtheta + integer :: iter + converged = .false. + do iter = 1, max_newton + call midpoint_residual(data, lvol, h, s0, theta0, zeta0, psi0, & + s, theta_delta, residual, av) + if (maxval(abs(residual)) < newton_tol) then + converged = .true. + return + end if + call midpoint_jacobian(data, lvol, h, s0, theta0, zeta0, psi0, & + s, theta_delta, jac) + det = jac(1, 1)*jac(2, 2) - jac(1, 2)*jac(2, 1) + if (abs(det) <= tiny(det)) return + ds = (-residual(1)*jac(2, 2) + residual(2)*jac(1, 2))/det + dtheta = (residual(1)*jac(2, 1) - residual(2)*jac(1, 1))/det + s = s + ds + theta_delta = theta_delta + dtheta + end do + end subroutine solve_midpoint + subroutine midpoint_residual(data, lvol, h, s0, theta0, zeta0, psi0, & + s, theta_delta, residual, av) + type(spectre_data_t), intent(in) :: data + integer, intent(in) :: lvol + real(dp), intent(in) :: h, s0, theta0, zeta0, psi0, s, theta_delta + real(dp), intent(out) :: residual(2) + type(spectre_vecpot_t), intent(out) :: av + + real(dp) :: ratio, flow + + call evaluate_hotloop_av(data, lvol, s, theta0 + theta_delta, & + zeta0 + 0.5_dp*h, av) + ratio = av%dAzt(1)/av%dAth(1) + flow = av%dAzt(2) - av%dAth(2)*ratio + residual(1) = theta_delta + 0.5_dp*h*ratio + residual(2) = av%Ath - psi0 - 0.5_dp*h*flow + end subroutine midpoint_residual + + subroutine midpoint_jacobian(data, lvol, h, s0, theta0, zeta0, psi0, & + s, theta_delta, jac) + type(spectre_data_t), intent(in) :: data + integer, intent(in) :: lvol + real(dp), intent(in) :: h, s0, theta0, zeta0, psi0, s, theta_delta + real(dp), intent(out) :: jac(2, 2) + + type(spectre_vecpot_t) :: av + real(dp) :: plus(2), minus(2), delta + + delta = fd_scale*max(1.0_dp, abs(s)) + call midpoint_residual(data, lvol, h, s0, theta0, zeta0, psi0, & + s + delta, theta_delta, plus, av) + call midpoint_residual(data, lvol, h, s0, theta0, zeta0, psi0, & + s - delta, theta_delta, minus, av) + jac(:, 1) = (plus - minus)/(2.0_dp*delta) + + delta = fd_scale + call midpoint_residual(data, lvol, h, s0, theta0, zeta0, psi0, & + s, theta_delta + delta, plus, av) + call midpoint_residual(data, lvol, h, s0, theta0, zeta0, psi0, & + s, theta_delta - delta, minus, av) + jac(:, 2) = (plus - minus)/(2.0_dp*delta) + end subroutine midpoint_jacobian + + subroutine recover_radius(data, lvol, theta, zeta, psi, s, converged) + type(spectre_data_t), intent(in) :: data + integer, intent(in) :: lvol + real(dp), intent(in) :: theta, zeta, psi + real(dp), intent(inout) :: s + logical, intent(out) :: converged + + type(spectre_vecpot_t) :: av + real(dp) :: residual + integer :: iter + + converged = .false. do iter = 1, max_newton - call eval_spectre_vector_potential(data, lvol, s1, theta0, zeta0, av) - ratio = av%dAzt(1)/av%dAth(1) - flow = av%dAzt(2) - av%dAth(2)*ratio - fun = av%Ath - psi0 - h*flow - if (abs(fun) < newton_tol) then + call evaluate_hotloop_av(data, lvol, s, theta, zeta, av) + residual = av%Ath - psi + if (abs(residual) < newton_tol) then converged = .true. - exit + return end if - dgds = av%d2Azt(2) & - - (av%d2Ath(2)*av%dAzt(1) + av%dAth(2)*av%d2Azt(1))/av%dAth(1) & - + av%dAth(2)*av%dAzt(1)*av%d2Ath(1)/av%dAth(1)**2 - dfun = av%dAth(1) - h*dgds - s1 = s1 - fun/dfun + s = s - residual/av%dAth(1) end do + end subroutine recover_radius - if (.not. converged) then - write (*, '(A,I0,A,ES13.5,A,ES13.5)') & - 'field_line_step: Newton failed to converge in volume ', lvol, & - ', s = ', s1, ', |F| = ', abs(fun) - error stop 'field_line_step: Newton did not converge' - end if + subroutine evaluate_hotloop_av(data, lvol, s, theta, zeta, av) + type(spectre_data_t), intent(in) :: data + integer, intent(in) :: lvol + real(dp), intent(in) :: s, theta, zeta + type(spectre_vecpot_t), intent(out) :: av - theta = theta0 - h*ratio - psi = av%Ath - zeta = zeta0 + h - s = s1 - end subroutine field_line_step + real(dp) :: x(3), values(2), derivatives(3, 2) + + if (hotloop_lvol == lvol) then + x = [s, modulo(theta, twopi), & + modulo(zeta, twopi/real(data%Nfp, dp))] + call evaluate_batch_splines_3d_der(hotloop_spline, x, values, derivatives) + av%Ath = values(1) + av%Azt = values(2) + av%dAth = derivatives(:, 1) + av%dAzt = derivatives(:, 2) + else + call eval_spectre_vector_potential(data, lvol, s, theta, zeta, av, & + second_derivatives=.false.) + end if + end subroutine evaluate_hotloop_av !> Cylindrical R, Z of a point in volume lvol from the SPEC interface Fourier !> blend, mirroring libneo_coordinates_spectre for output only: interior diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt index 1ee57cb2..0c8a2612 100644 --- a/test/tests/CMakeLists.txt +++ b/test/tests/CMakeLists.txt @@ -844,6 +844,16 @@ set_tests_properties(test_spectre_poincare PROPERTIES LABELS "integration;python" TIMEOUT 120) +add_executable(test_field_line_midpoint.x test_field_line_midpoint.f90) +target_link_libraries(test_field_line_midpoint.x simple) +add_test(NAME test_field_line_midpoint + COMMAND test_field_line_midpoint.x + ${libneo_SOURCE_DIR}/test/spectre/data/spectre_test.h5 + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +set_tests_properties(test_field_line_midpoint PROPERTIES + LABELS "unit;spectre" + TIMEOUT 15) + # SPECTRE guiding-center RK45 tracing per volume with boundary-stop (#438). # Uses simple.x with the committed libneo SPECTRE HDF5 fixture; VMEC-free. add_test(NAME test_spectre_gc_rk45 diff --git a/test/tests/test_field_line_midpoint.f90 b/test/tests/test_field_line_midpoint.f90 new file mode 100644 index 00000000..115cc8b4 --- /dev/null +++ b/test/tests/test_field_line_midpoint.f90 @@ -0,0 +1,29 @@ +program test_field_line_midpoint + use, intrinsic :: iso_fortran_env, only: dp => real64 + use field_line_spectre, only: field_line_psi, field_line_step + use spectre_reader, only: spectre_data_t, load_spectre, free_spectre + implicit none + + type(spectre_data_t) :: data + character(len=1024) :: field + real(dp) :: state(4), initial(4), h + integer :: ierr + + call get_command_argument(1, field) + call load_spectre(trim(field), data, ierr) + if (ierr /= 0) error stop 'failed to load SPECTRE fixture' + + state(1:3) = [0.2_dp, 0.3_dp, 0.1_dp] + state(4) = field_line_psi(data, 2, state(1), state(2), state(3)) + initial = state + h = 3.141592653589793238462643383279_dp/(64.0_dp*real(data%Nfp, dp)) + + call field_line_step(data, 2, h, state(1), state(2), state(3), state(4)) + call field_line_step(data, 2, -h, state(1), state(2), state(3), state(4)) + if (maxval(abs(state - initial)) > 2.0e-11_dp) then + error stop 'implicit midpoint step is not reversible' + end if + + call free_spectre(data) + print *, 'SPECTRE field-line midpoint reversibility PASS' +end program test_field_line_midpoint diff --git a/test/tests/test_spectre_poincare.py b/test/tests/test_spectre_poincare.py index 13a65896..a3276991 100644 --- a/test/tests/test_spectre_poincare.py +++ b/test/tests/test_spectre_poincare.py @@ -29,12 +29,15 @@ ENV_REACH = 0.05 -def write_namelist(path, h5, nturns, nsteps, nseeds, volumes=None): +def write_namelist(path, h5, nturns, nsteps, nseeds, volumes=None, + seed_file=None): lines = ["&poincare", f" spectre_file = '{h5}'", " zeta0 = 0.0", f" nturns = {nturns}", f" nsteps_per_period = {nsteps}", f" nseeds_per_volume = {nseeds}"] if volumes is not None: lines.append(" volumes = " + ", ".join(str(v) for v in volumes)) + if seed_file is not None: + lines.append(f" seed_file = '{seed_file}'") lines.append("/") with open(path, "w") as f: f.write("\n".join(lines) + "\n") @@ -156,11 +159,8 @@ def main(): print(f"BDD3 vol{lvol}: iota [{iota.min():.5f},{iota.max():.5f}] " f"ref [{ref_iota.min():.5f},{ref_iota.max():.5f}]") - # BDD 4: step-halving convergence. Symplectic Euler is first order, so - # the final section point shift is O(h): it stays small and halves on - # each doubling (the literal 1e-6 target needs ~1e6 steps/period). A - # loose Newton tolerance would floor the shift and break the halving; - # Newton runs at 1e-13, so the shift keeps shrinking. + # BDD 4: implicit midpoint is second order, so the final-section shift + # falls by about four when the step is halved. pts = {} for nsteps in (256, 512, 1024): run_app(binary, work, h5, nturns=3, nsteps=nsteps, nseeds=4, @@ -175,10 +175,21 @@ def shift(n1, n2): d1, d2 = shift(256, 512), shift(512, 1024) if not d1 < 1.0e-3: failures.append(f"BDD4: coarse shift {d1:.3e} not small") - if not d2 < 0.75 * d1: + if not d2 < 0.40 * d1: failures.append(f"BDD4: shift not converging {d2:.3e} vs {d1:.3e}") print(f"BDD4: shift 256->512 = {d1:.3e}, 512->1024 = {d2:.3e}") + seed_path = os.path.join(work, "seeds.dat") + explicit = np.array([[1, -0.4, 0.1], [1, 0.2, 0.7]]) + np.savetxt(seed_path, explicit, fmt=["%d", "%.17g", "%.17g"]) + run_app(binary, work, h5, nturns=1, nsteps=64, nseeds=1, + volumes=[2], seed_file=seed_path) + initial = load_volume(work, 2) + initial = initial[initial[:, 1] == 0][:, [2, 3]] + if not np.allclose(initial, explicit[:, 1:], rtol=0.0, atol=1.0e-14): + failures.append("BDD5: explicit field-line seeds were replaced") + print("BDD5: explicit field-line seeds preserved") + if failures: print("\nFAILURES:") for f in failures: