From d96bb0dc62015666e68e05ae538e1d011e0ed8cb Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Mon, 22 Jun 2026 15:50:01 +0200 Subject: [PATCH 01/10] POTATO: trace edge bananas across the separatrix Add the orbit-side support for guiding-center bananas whose tips cross the last closed flux surface, complementing the libneo field_eq_mod::allow_sol flag. - potato.in gains edge_extension (default .false.); tt.f90 sets allow_sol from it, so orbits may enter the scrape-off layer in the extended equilibrium field. - New itest_type=4 single-orbit trace: starts one guiding-center orbit from (orbit_Rstart, orbit_Zstart) with pitch orbit_lambda and integrates it over one bounce via find_bounce with write_orb, writing the trajectory to fort.100. This is the minimal check that an edge banana now closes (ierr=0). - CMakeLists: drop the global -std=f2008/-Wfatal-errors, which leaked into the fetched fortplot (needs F2018 assumed-rank) and vode (legacy DO termination), breaking the build. The deps now compile with the default standard; POTATO sources are unaffected. --- POTATO/CMakeLists.txt | 2 +- POTATO/SRC/potato_input_mod.f90 | 19 +++++++++++++- POTATO/SRC/tt.f90 | 45 ++++++++++++++++++++++++++++++--- 3 files changed, 61 insertions(+), 5 deletions(-) diff --git a/POTATO/CMakeLists.txt b/POTATO/CMakeLists.txt index 2d1be317..9cb4e80d 100644 --- a/POTATO/CMakeLists.txt +++ b/POTATO/CMakeLists.txt @@ -5,7 +5,7 @@ enable_language(Fortran) set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/../cmake") set(CMAKE_Fortran_MODULE_DIRECTORY ${PROJECT_BINARY_DIR}) -add_compile_options(-cpp -g -fbacktrace -Wfatal-errors -std=f2008) +add_compile_options(-cpp -g -fbacktrace) add_compile_options(-O3 -march=native -mtune=native) include(Util) diff --git a/POTATO/SRC/potato_input_mod.f90 b/POTATO/SRC/potato_input_mod.f90 index a72c8e5c..75dae622 100644 --- a/POTATO/SRC/potato_input_mod.f90 +++ b/POTATO/SRC/potato_input_mod.f90 @@ -48,12 +48,24 @@ module potato_input_mod double precision :: toten_plot = 1d0 double precision :: perpinv_plot = 4.5d-5 + ! Single-orbit trace (itest_type=4): start point on the poloidal plane and + ! pitch cosine of one guiding-center orbit, traced over one bounce period. + double precision :: orbit_Rstart = 210d0 ! start major radius [cm] + double precision :: orbit_Zstart = 0d0 ! start height [cm] + double precision :: orbit_lambda = 0.4d0 ! pitch cosine v_par/v at start + ! Canonical frequency plot double precision :: enkin_over_temp = 1d0 ! Profile input file character(len=256) :: profile_file = 'profile_poly.in' + ! Edge extension: when .true. the orbit integrator accepts guiding-center + ! points outside the separatrix (scrape-off layer), so bananas whose tips + ! cross the last closed flux surface can be traced in the extended + ! equilibrium field. Sets field_eq_mod::allow_sol. Default .false. + logical :: edge_extension = .false. + namelist /potato_nml/ & itest_type, E_alpha, A_alpha, Z_alpha, & rho_pol, rho_pol_max, scalfac_energy, scalfac_efield, & @@ -63,7 +75,8 @@ module potato_input_mod adaptive_jperp, npoi_init, nlagr_sampling, eps_sampling, & itermax_sampling, & toten_plot, perpinv_plot, enkin_over_temp, & - profile_file + profile_file, edge_extension, & + orbit_Rstart, orbit_Zstart, orbit_lambda contains @@ -126,6 +139,10 @@ subroutine print_potato_input(iunit) write(iunit, '(A,ES12.5)') ' perpinv_plot = ', perpinv_plot write(iunit, '(A,ES12.5)') ' enkin_over_temp = ', enkin_over_temp write(iunit, '(A,A)') ' profile_file = ', trim(profile_file) + write(iunit, '(A,L1)') ' edge_extension = ', edge_extension + write(iunit, '(A,ES12.5)') ' orbit_Rstart = ', orbit_Rstart + write(iunit, '(A,ES12.5)') ' orbit_Zstart = ', orbit_Zstart + write(iunit, '(A,ES12.5)') ' orbit_lambda = ', orbit_lambda write(iunit, '(A)') '================================' end subroutine print_potato_input diff --git a/POTATO/SRC/tt.f90 b/POTATO/SRC/tt.f90 index 842cd956..279d420d 100644 --- a/POTATO/SRC/tt.f90 +++ b/POTATO/SRC/tt.f90 @@ -2,7 +2,7 @@ program classify_orbits ! use parmot_mod, only : rmu,ro0 - use orbit_dim_mod, only : write_orb,iunit1,numbasef + use orbit_dim_mod, only : write_orb,iunit1,numbasef,neqm use get_matrix_mod, only : iclass use global_invariants, only : dtau,toten,perpinv,cE_ref,Phi_eff use bounds_fixpoints_mod, only : nregions @@ -18,7 +18,10 @@ program classify_orbits rho_pol, rho_pol_max, scalfac_energy, & scalfac_efield, Rmax_orbit, ntimstep, & npoicut, m_min, m_max, n_tor, & - toten_plot, perpinv_plot, profile_file + toten_plot, perpinv_plot, profile_file, & + edge_extension, & + orbit_Rstart, orbit_Zstart, orbit_lambda + use field_eq_mod, only : allow_sol ! implicit none ! @@ -30,10 +33,14 @@ program classify_orbits double precision,parameter :: ev=1.6022d-12 ! logical :: classes_talk,plot_orbits,compute_equibrium_profiles,compute_resonant_torque,plot_poicut + logical :: trace_single_orbit ! integer :: ifdir_type,ierr,m,iunit,i,log_u double precision :: v0,bmod_ref double precision, dimension(:,:), allocatable :: resint + double precision, dimension(neqm) :: z_single + double precision :: taub_single,delphi_single,extraset_single(1) + external :: find_bounce, velo ! ! @@ -43,9 +50,13 @@ program classify_orbits call read_potato_input('potato.in') call print_potato_input(6) if (is_logging_enabled()) call print_potato_input(log_u) +! +! Allow orbits to cross the separatrix into the scrape-off layer when requested: + allow_sol = edge_extension ! iunit=71 ! + trace_single_orbit=.false. select case(itest_type) case(1) call tee_message('Plot orbits') @@ -62,6 +73,12 @@ program classify_orbits plot_orbits=.false. compute_equibrium_profiles=.false. compute_resonant_torque=.true. + case(4) + call tee_message('Trace single orbit') + plot_orbits=.false. + compute_equibrium_profiles=.false. + compute_resonant_torque=.false. + trace_single_orbit=.true. case default call tee_message('unknown test') call close_logging() @@ -150,7 +167,29 @@ program classify_orbits ! ! call test_eqmagprofs !WARNING: contains "stop" inside, comment this test out to continue ! - +!....................................... +! +! Trace a single guiding-center orbit over one bounce period and write it. +! With edge_extension=.true. the orbit may cross the separatrix into the SOL. +! + if(trace_single_orbit) then + call tee_message('Tracing single orbit over one bounce') + z_single(1)=orbit_Rstart ! R [cm] + z_single(2)=0.d0 ! phi + z_single(3)=orbit_Zstart ! Z [cm] + z_single(4)=1.d0 ! p, monoenergetic at E_alpha + z_single(5)=orbit_lambda ! pitch cosine v_par/v + write_orb=.true. + iunit1=100 + extraset_single=0.d0 + call find_bounce(0,velo,dtau,z_single,taub_single,delphi_single, & + extraset_single,ierr) + close(iunit1) + write_orb=.false. + write(*,'(A,ES14.6,A,ES14.6,A,I0)') ' single orbit: taub=',taub_single, & + ' delphi=',delphi_single,' ierr=',ierr + endif +! From a45630cb7df25c49d637351ee9efc4fdef076ed3 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Mon, 22 Jun 2026 16:06:47 +0200 Subject: [PATCH 02/10] POTATO: frequency radial scan (itest_type=5) Scan freq_n midplane start radii from freq_Rmin to freq_Rmax at fixed pitch orbit_lambda; at each, trace one bounce via find_bounce and record the bounce frequency omega_b = 2 pi v0 / taub and toroidal precession omega_phi = delphi v0 / taub against the flux label rho_pol (from psi at the start point). Writes freq_scan.dat. With edge_extension the scan extends across the separatrix into the near-SOL where banana tips still turn. This is the POTATO side of the rung-2 frequency benchmark against thin-orbit NEO-RT. --- POTATO/SRC/potato_input_mod.f90 | 13 ++++++- POTATO/SRC/tt.f90 | 62 ++++++++++++++++++++++++++++++--- 2 files changed, 70 insertions(+), 5 deletions(-) diff --git a/POTATO/SRC/potato_input_mod.f90 b/POTATO/SRC/potato_input_mod.f90 index 75dae622..16cbadd5 100644 --- a/POTATO/SRC/potato_input_mod.f90 +++ b/POTATO/SRC/potato_input_mod.f90 @@ -54,6 +54,13 @@ module potato_input_mod double precision :: orbit_Zstart = 0d0 ! start height [cm] double precision :: orbit_lambda = 0.4d0 ! pitch cosine v_par/v at start + ! Frequency radial scan (itest_type=5): trace one bounce at each of freq_n + ! midplane start radii from freq_Rmin to freq_Rmax (pitch orbit_lambda), + ! and write rho_pol, omega_b, omega_phi to freq_scan.dat. + double precision :: freq_Rmin = 178d0 ! inner start radius [cm] + double precision :: freq_Rmax = 221d0 ! outer start radius [cm], into SOL + integer :: freq_n = 40 ! number of surfaces + ! Canonical frequency plot double precision :: enkin_over_temp = 1d0 @@ -76,7 +83,8 @@ module potato_input_mod itermax_sampling, & toten_plot, perpinv_plot, enkin_over_temp, & profile_file, edge_extension, & - orbit_Rstart, orbit_Zstart, orbit_lambda + orbit_Rstart, orbit_Zstart, orbit_lambda, & + freq_Rmin, freq_Rmax, freq_n contains @@ -143,6 +151,9 @@ subroutine print_potato_input(iunit) write(iunit, '(A,ES12.5)') ' orbit_Rstart = ', orbit_Rstart write(iunit, '(A,ES12.5)') ' orbit_Zstart = ', orbit_Zstart write(iunit, '(A,ES12.5)') ' orbit_lambda = ', orbit_lambda + write(iunit, '(A,ES12.5)') ' freq_Rmin = ', freq_Rmin + write(iunit, '(A,ES12.5)') ' freq_Rmax = ', freq_Rmax + write(iunit, '(A,I0)') ' freq_n = ', freq_n write(iunit, '(A)') '================================' end subroutine print_potato_input diff --git a/POTATO/SRC/tt.f90 b/POTATO/SRC/tt.f90 index 279d420d..4ca11d4e 100644 --- a/POTATO/SRC/tt.f90 +++ b/POTATO/SRC/tt.f90 @@ -20,8 +20,10 @@ program classify_orbits npoicut, m_min, m_max, n_tor, & toten_plot, perpinv_plot, profile_file, & edge_extension, & - orbit_Rstart, orbit_Zstart, orbit_lambda - use field_eq_mod, only : allow_sol + orbit_Rstart, orbit_Zstart, orbit_lambda, & + freq_Rmin, freq_Rmax, freq_n + use field_eq_mod, only : allow_sol, psi_axis, psi_sep + use field_sub, only : psif ! implicit none ! @@ -33,14 +35,18 @@ program classify_orbits double precision,parameter :: ev=1.6022d-12 ! logical :: classes_talk,plot_orbits,compute_equibrium_profiles,compute_resonant_torque,plot_poicut - logical :: trace_single_orbit + logical :: trace_single_orbit,freq_scan ! integer :: ifdir_type,ierr,m,iunit,i,log_u double precision :: v0,bmod_ref double precision, dimension(:,:), allocatable :: resint double precision, dimension(neqm) :: z_single double precision :: taub_single,delphi_single,extraset_single(1) - external :: find_bounce, velo + external :: find_bounce, velo, magfie +! frequency radial scan (itest_type=5): + integer :: ifreq + double precision :: Rstep,omega_b,omega_phi,s_norm,rhopol + double precision :: bmod_d,sqrtg_d,bder_d(3),hcovar_d(3),hctrvr_d(3),hcurl_d(3),x_d(3) ! ! @@ -57,6 +63,7 @@ program classify_orbits iunit=71 ! trace_single_orbit=.false. + freq_scan=.false. select case(itest_type) case(1) call tee_message('Plot orbits') @@ -79,6 +86,12 @@ program classify_orbits compute_equibrium_profiles=.false. compute_resonant_torque=.false. trace_single_orbit=.true. + case(5) + call tee_message('Frequency radial scan') + plot_orbits=.false. + compute_equibrium_profiles=.false. + compute_resonant_torque=.false. + freq_scan=.true. case default call tee_message('unknown test') call close_logging() @@ -190,6 +203,47 @@ program classify_orbits ' delphi=',delphi_single,' ierr=',ierr endif ! +!....................................... +! +! Frequency radial scan: at each midplane start radius trace one bounce and +! record bounce frequency omega_b = 2 pi v0 / taub and toroidal precession +! frequency omega_phi = delphi v0 / taub against the flux label rho_pol. +! + if(freq_scan) then + call tee_message('Frequency radial scan') + write_orb=.false. + open(iunit,file='freq_scan.dat',status='replace',action='write') + write(iunit,'(A)') '# R_start[cm] rho_pol omega_b[1/s] omega_phi[1/s] '// & + 'taub delphi ierr' + Rstep=(freq_Rmax-freq_Rmin)/dble(max(freq_n-1,1)) + do ifreq=1,freq_n + z_single(1)=freq_Rmin+Rstep*dble(ifreq-1) + z_single(2)=0.d0 + z_single(3)=orbit_Zstart + z_single(4)=1.d0 + z_single(5)=orbit_lambda +! flux label of the start point: evaluate the field there, then normalize psi. + x_d=z_single(1:3) + call magfie(x_d,bmod_d,sqrtg_d,bder_d,hcovar_d,hctrvr_d,hcurl_d) + s_norm=(psif-psi_axis)/(psi_sep-psi_axis) + rhopol=sqrt(max(s_norm,0.d0)) + extraset_single=0.d0 + call find_bounce(0,velo,dtau,z_single,taub_single,delphi_single, & + extraset_single,ierr) + if(ierr.eq.0 .and. taub_single.gt.0.d0) then + omega_b=2.d0*pi*v0/taub_single + omega_phi=delphi_single*v0/taub_single + write(iunit,'(6ES16.7,I3)') z_single(1),rhopol,omega_b,omega_phi, & + taub_single,delphi_single,ierr + else + write(iunit,'(6ES16.7,I3)') z_single(1),rhopol,0.d0,0.d0, & + taub_single,delphi_single,ierr + endif + enddo + close(iunit) + call tee_message('Frequency radial scan done -> freq_scan.dat') + endif +! From 00a2e00279fec9be4a0eb4848b46a01b2f8590ea Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Mon, 22 Jun 2026 16:20:01 +0200 Subject: [PATCH 03/10] test: thin-orbit frequency radial scan for the rung-2 benchmark test_freq_scan loops the flux surface s, calling neort init at each, and writes the bounce frequency Om_th and toroidal precession Om_ph at fixed pitch (eta = (1-lambda^2)/Bmin, i.e. v_par/v=0.4 at the midplane) and fixed energy (10 keV D) to freq_scan_neort.dat, with E x B set to zero. This is the thin-orbit reference curve compared against POTATO itest_type=5. --- test/test_freq_scan.f90 | 49 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) create mode 100644 test/test_freq_scan.f90 diff --git a/test/test_freq_scan.f90 b/test/test_freq_scan.f90 new file mode 100644 index 00000000..205e5647 --- /dev/null +++ b/test/test_freq_scan.f90 @@ -0,0 +1,49 @@ +program test_freq_scan + ! Thin-orbit canonical frequencies vs flux surface for the rung-2 benchmark + ! against POTATO. At each s the bounce frequency Om_th and the toroidal + ! precession Om_ph are evaluated at a fixed pitch (eta = 0.84/Bmin, i.e. + ! v_par/v = 0.4 at the local low-field-side midplane, matching POTATO's + ! orbit_lambda) and a fixed energy. No E x B (Om_tE = 0) to match POTATO's + ! zero-potential profile, so Om_ph is the magnetic precession alone. + use iso_fortran_env, only: dp => real64 + use do_magfie_mod, only: s, inp_swi, do_magfie_init + use neort, only: init + use neort_freq, only: Om_th, Om_ph + use driftorbit, only: mph, mth, vth, etadt, Om_tE, sign_vpar, magdrift + + implicit none + + integer, parameter :: ns = 30 + real(dp), parameter :: smin = 0.05_dp, smax = 0.95_dp + real(dp), parameter :: v_10keV = 9.787e7_dp ! 10 keV deuteron, cm/s + real(dp), parameter :: lambda0 = 0.4_dp ! v_par/v at the midplane + + integer :: i, fid + real(dp) :: v, eta, Omth, dOmthdv, dOmthdeta, Omph, dOmphdv, dOmphdeta + + inp_swi = 9 ! ASDEX Upgrade Boozer format + magdrift = .true. + s = smin + call do_magfie_init("in_file") + + mth = 1 + mph = 2 + vth = v_10keV + v = v_10keV + Om_tE = 0.0_dp ! no E x B, magnetic precession only + sign_vpar = 1.0_dp + + open(newunit=fid, file='freq_scan_neort.dat') + write(fid, '(A)') '# s rho_tor omega_b[1/s] omega_phi[1/s]' + do i = 1, ns + s = smin + (smax - smin)*real(i - 1, dp)/real(ns - 1, dp) + call init + eta = (1.0_dp - lambda0**2)*etadt ! = 0.84/Bmin, lambda=0.4 at midplane + call Om_th(v, eta, Omth, dOmthdv, dOmthdeta) + call Om_ph(v, eta, Omph, dOmphdv, dOmphdeta) + write(fid, '(4ES16.7)') s, sqrt(s), Omth, Omph + end do + close(fid) + print *, 'freq_scan_neort.dat written, ns =', ns + +end program test_freq_scan From 248b6f926f12322359209a7f48e0d95830caf0af Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Mon, 22 Jun 2026 20:36:10 +0200 Subject: [PATCH 04/10] test: freq scan outputs iota and uses 60 surfaces Emit iota per surface so the rho_pol map can be built from NEO-RT's own q (q = 1/iota) instead of a gfile q-profile, and densify the scan to 60 flux surfaces for a smoother edge curve. --- test/test_freq_scan.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_freq_scan.f90 b/test/test_freq_scan.f90 index 205e5647..be9d7ee3 100644 --- a/test/test_freq_scan.f90 +++ b/test/test_freq_scan.f90 @@ -6,14 +6,14 @@ program test_freq_scan ! orbit_lambda) and a fixed energy. No E x B (Om_tE = 0) to match POTATO's ! zero-potential profile, so Om_ph is the magnetic precession alone. use iso_fortran_env, only: dp => real64 - use do_magfie_mod, only: s, inp_swi, do_magfie_init + use do_magfie_mod, only: s, inp_swi, do_magfie_init, iota use neort, only: init use neort_freq, only: Om_th, Om_ph use driftorbit, only: mph, mth, vth, etadt, Om_tE, sign_vpar, magdrift implicit none - integer, parameter :: ns = 30 + integer, parameter :: ns = 60 real(dp), parameter :: smin = 0.05_dp, smax = 0.95_dp real(dp), parameter :: v_10keV = 9.787e7_dp ! 10 keV deuteron, cm/s real(dp), parameter :: lambda0 = 0.4_dp ! v_par/v at the midplane @@ -34,14 +34,14 @@ program test_freq_scan sign_vpar = 1.0_dp open(newunit=fid, file='freq_scan_neort.dat') - write(fid, '(A)') '# s rho_tor omega_b[1/s] omega_phi[1/s]' + write(fid, '(A)') '# s rho_tor omega_b[1/s] omega_phi[1/s] iota' do i = 1, ns s = smin + (smax - smin)*real(i - 1, dp)/real(ns - 1, dp) call init eta = (1.0_dp - lambda0**2)*etadt ! = 0.84/Bmin, lambda=0.4 at midplane call Om_th(v, eta, Omth, dOmthdv, dOmthdeta) call Om_ph(v, eta, Omph, dOmphdv, dOmphdeta) - write(fid, '(4ES16.7)') s, sqrt(s), Omth, Omph + write(fid, '(5ES16.7)') s, sqrt(s), Omth, Omph, iota end do close(fid) print *, 'freq_scan_neort.dat written, ns =', ns From 117c253da49d236f562e1f4c07f93834c0baf305 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 23 Jun 2026 17:03:23 +0200 Subject: [PATCH 05/10] Tighten POTATO low-energy resonance handling Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> --- POTATO/CMakeLists.txt | 12 + POTATO/SRC/CMakeLists.txt | 3 +- POTATO/SRC/box_counting.f90 | 122 ++++---- POTATO/SRC/find_all_roots.f90 | 32 ++- POTATO/SRC/potato_input_mod.f90 | 4 +- POTATO/SRC/resonant_int.f90 | 91 ++++-- POTATO/SRC/sub_potato.f90 | 473 +++++++++++++++++++++----------- POTATO/SRC/tt.f90 | 5 +- 8 files changed, 466 insertions(+), 276 deletions(-) diff --git a/POTATO/CMakeLists.txt b/POTATO/CMakeLists.txt index 9cb4e80d..3ef7a2e0 100644 --- a/POTATO/CMakeLists.txt +++ b/POTATO/CMakeLists.txt @@ -21,6 +21,18 @@ FetchContent_Declare( ) FetchContent_MakeAvailable(fortplot) +# fortnum provides the variable-order Adams ODE integrator with event/root +# detection used by find_bounce to locate the bounce return (replacing the +# fixed-step Poincare-cut crossing search). Same pin as the NEO-RT migration. +if(NOT TARGET fortnum) + FetchContent_Declare( + fortnum + GIT_REPOSITORY https://github.com/lazy-fortran/fortnum.git + GIT_TAG 92de6e949a772cfffc73bb5295fe5e2b056b9c18 + ) + FetchContent_MakeAvailable(fortnum) +endif() + add_subdirectory(SRC) add_executable (potato.x diff --git a/POTATO/SRC/CMakeLists.txt b/POTATO/SRC/CMakeLists.txt index e3cfcd68..a552da6f 100644 --- a/POTATO/SRC/CMakeLists.txt +++ b/POTATO/SRC/CMakeLists.txt @@ -37,8 +37,9 @@ add_library(potato_base neo2_reader.f90 ) find_or_fetch(vode) -target_link_libraries(potato_base PUBLIC +target_link_libraries(potato_base PUBLIC vode +fortnum BLAS::BLAS LAPACK::LAPACK HDF5::HDF5 diff --git a/POTATO/SRC/box_counting.f90 b/POTATO/SRC/box_counting.f90 index b32bfdde..4d2385fb 100644 --- a/POTATO/SRC/box_counting.f90 +++ b/POTATO/SRC/box_counting.f90 @@ -24,7 +24,7 @@ end subroutine timestep_vode subroutine time_in_box(z, cnt, sbox, taub, tau) ! Returns time spent in boxes - use dvode_f90_m, only: vode_opts, set_normal_opts, dvode_f90, get_stats + use dvode_f90_m, only: vode_opts, set_opts, dvode_f90, get_stats use field_sub, only: psif use field_eq_mod, only: psi_axis, psi_sep use orbit_dim_mod, only: neqm @@ -37,23 +37,18 @@ subroutine time_in_box(z, cnt, sbox, taub, tau) real(8), intent(in) :: taub ! Bounce time real(8), intent(out) :: tau(cnt) ! Time in each box - real(8) :: delphi - real(8) :: sprev, snext ! Previous and next flux radius box + real(8) :: smid real(8) :: y(neqm) - integer(4) :: nmax ! Maximum loop iterations - - integer(4) :: k + integer(4) :: k, nsample real(8) :: ti - real(8) :: atol(neqm), rtol, tout, rstats(22) - integer(4) :: itask, istate, istats(31), numevents + real(8) :: atol(neqm), rtol, tout + integer(4) :: itask, istate, method_flag type (vode_opts) :: options real(8) :: bmod, phi_elec, s - real(8) :: sold, told, yold(neqm) - integer(4) :: sind, sind0 ! s index - - integer(4) :: jroots(2) + real(8) :: sold, told + integer(4) :: sind ! s index external timestep_vode @@ -63,83 +58,58 @@ subroutine time_in_box(z, cnt, sbox, taub, tau) atol = 1d-13 itask = 1 istate = 1 - numevents = 2 - options = set_normal_opts(abserr_vector=atol, relerr=rtol, nevents=numevents) + method_flag = 10 + options = set_opts(method_flag=method_flag, abserr_vector=atol, & + relerr=rtol, mxstep=2000) - nmax = 3*size(sbox) + nsample = min(256, max(64, size(sbox))) tau = 0d0 ti = 0d0 call get_bmod_and_Phi(z(1:3), bmod, phi_elec) s = abs((psif-psi_axis)/(psi_sep-psi_axis)) - sind = size(sbox)+1 - sprev = -1e5 - snext = 1e5 - do k = 1,size(sbox) - if(sbox(k)>s) then - sind = k - snext = sbox(k) - if(k>1) sprev = sbox(k-1) - exit - end if - enddo - sind0 = sind - + sold = s told = 0d0 - do k = 1,nmax - yold = y - sold = s - tout = taub + do k = 1,nsample + tout = taub*dble(k)/dble(nsample) call dvode_f90( & - timestep_vode, neqm, y, ti, tout, itask, istate, options, & - g_fcn = sroots ) - if (istate == 2) exit - if (istate == 3) then - tau(sind) = tau(sind) + ti-told - told = ti - call get_stats(rstats, istats, numevents, jroots) - if (jroots(2).ne.0) then - sind = sind + 1 ! moving outwards - sprev = snext - if (sind == size(sbox)+1) then - snext = 1e5 - else - snext = sbox(sind) - end if - end if - if (jroots(1).ne.0) then - sind = sind - 1 ! moving inwards - snext = sprev - if (sind == 1) then - sprev = -1e5 - else - sprev = sbox(sind-1) - end if - end if - elseif (istate == -1) then - print *, 'Error in VODE: k =', k - istate = 1 - exit + timestep_vode, neqm, y, ti, tout, itask, istate, options) + if (istate == -1) then + print *, 'time_in_box: VODE exceeded mxstep at k =', k, & + ' ti =', ti, ' taub =', taub + tau = 0d0 + return + endif + if (istate /= 2) then + print *, 'time_in_box: VODE failed with istate =', istate, & + ' k =', k, ' ti =', ti, ' taub =', taub + tau = 0d0 + return end if - end do - - tau(sind) = tau(sind) + taub-told -contains - - subroutine sroots(neqext, t, yext, ng, gout) - ! For finding roots between boxes + call get_bmod_and_Phi(y(1:3), bmod, phi_elec) - integer, intent(in) :: neqext, ng - real(8), intent(in) :: t, yext(neqext) - real(8), intent(out) :: gout(ng) + s = abs((psif-psi_axis)/(psi_sep-psi_axis)) + smid = 0.5d0*(sold+s) + sind = radial_box_index(smid) + tau(sind) = tau(sind) + ti-told + sold = s + told = ti + end do - call get_bmod_and_Phi(yext(1:3), bmod, phi_elec) +contains - s = abs((psif-psi_axis)/(psi_sep-psi_axis)) + integer function radial_box_index(sval) + real(8), intent(in) :: sval + integer :: ibox - gout(1) = s - sprev - gout(2) = s - snext - end subroutine sroots + radial_box_index = size(sbox) + do ibox = 1,size(sbox) + if (sval <= sbox(ibox)) then + radial_box_index = ibox + exit + endif + enddo + end function radial_box_index end subroutine time_in_box diff --git a/POTATO/SRC/find_all_roots.f90 b/POTATO/SRC/find_all_roots.f90 index ed6dfa08..09dec90a 100644 --- a/POTATO/SRC/find_all_roots.f90 +++ b/POTATO/SRC/find_all_roots.f90 @@ -5,7 +5,8 @@ ! module find_all_roots_mod logical :: customgrid=.false. - integer :: nroots, nsearch_min=100, ncustom, niter=100 + logical :: fail_fast=.false. + integer :: nroots, nsearch_min=100, ncustom, niter=100, max_roots_abort=0 double precision :: relerr_allroots=1.d-12 double precision, dimension(:), allocatable :: xcustom,roots end module find_all_roots_mod @@ -25,7 +26,7 @@ subroutine find_all_roots(fun,x1in,x2in,ierr) ! use find_all_roots_mod, only : nsearch_min,niter,relerr_allroots, & customgrid,ncustom,xcustom, & - nroots,roots + nroots,roots,fail_fast,max_roots_abort ! implicit none ! @@ -33,9 +34,11 @@ subroutine find_all_roots(fun,x1in,x2in,ierr) double precision :: x,x1in,x2in,f,df,hx,dx,errdist,xb,xe,dfb,dfe,xxtr,fxtr double precision :: x1,x2 double precision, dimension(:), allocatable :: xarr,farr,dfarr,dummy1d + logical :: abort_search external :: fun ! ierr=0 + abort_search=.false. x1=x1in x2=x2in ! @@ -137,6 +140,7 @@ subroutine find_all_roots(fun,x1in,x2in,ierr) x=(farr(i)*xarr(i-1)-farr(i-1)*xarr(i))/(farr(i)-farr(i-1)) ! call addroot + if(abort_search) exit ! endif else @@ -166,6 +170,7 @@ subroutine find_all_roots(fun,x1in,x2in,ierr) x=(fxtr*xarr(i-1)-farr(i-1)*xxtr)/(fxtr-farr(i-1)) ! call addroot + if(abort_search) exit ! endif if(fxtr*farr(i).le.0.d0) then @@ -173,9 +178,11 @@ subroutine find_all_roots(fun,x1in,x2in,ierr) x=(farr(i)*xxtr-fxtr*xarr(i))/(farr(i)-fxtr) ! call addroot + if(abort_search) exit ! endif endif + if(abort_search) exit enddo ! !------------ @@ -189,6 +196,12 @@ subroutine addroot ! Root adjustment by Newton method ! double precision :: xx +! + if(max_roots_abort.gt.0 .and. nroots.ge.max_roots_abort) then + ierr=2 + abort_search=.true. + return + endif ! ! extend storage arragy for a new root: nroots=nroots+1 @@ -222,9 +235,22 @@ subroutine addroot enddo ! if(iter.gt.niter) then + ierr=1 + if(fail_fast) then + abort_search=.true. + nroots=nroots-1 + if(nroots.eq.0) then + deallocate(roots) + else + dummy1d=roots(1:nroots) + deallocate(roots) + allocate(roots(nroots)) + roots=dummy1d + endif + return + endif print *,'No convergence of Newton in find_all_roots: error = ',sngl(abs(dx)), & ' tolerance = ',sngl(errdist) - ierr=1 endif ! roots(nroots)=x diff --git a/POTATO/SRC/potato_input_mod.f90 b/POTATO/SRC/potato_input_mod.f90 index 16cbadd5..19718e10 100644 --- a/POTATO/SRC/potato_input_mod.f90 +++ b/POTATO/SRC/potato_input_mod.f90 @@ -43,6 +43,7 @@ module potato_input_mod integer :: nlagr_sampling = 3 double precision :: eps_sampling = 1d-2 integer :: itermax_sampling = 5 + logical :: clip_resonance_classes = .true. ! Orbit plotting double precision :: toten_plot = 1d0 @@ -80,7 +81,7 @@ module potato_input_mod m_min, m_max, n_tor, & nenerg, thermen_max, nbox, & adaptive_jperp, npoi_init, nlagr_sampling, eps_sampling, & - itermax_sampling, & + itermax_sampling, clip_resonance_classes, & toten_plot, perpinv_plot, enkin_over_temp, & profile_file, edge_extension, & orbit_Rstart, orbit_Zstart, orbit_lambda, & @@ -143,6 +144,7 @@ subroutine print_potato_input(iunit) write(iunit, '(A,I0)') ' nlagr_sampling = ', nlagr_sampling write(iunit, '(A,ES12.5)') ' eps_sampling = ', eps_sampling write(iunit, '(A,I0)') ' itermax_sampling = ', itermax_sampling + write(iunit, '(A,L1)') ' clip_resonance_classes = ', clip_resonance_classes write(iunit, '(A,ES12.5)') ' toten_plot = ', toten_plot write(iunit, '(A,ES12.5)') ' perpinv_plot = ', perpinv_plot write(iunit, '(A,ES12.5)') ' enkin_over_temp = ', enkin_over_temp diff --git a/POTATO/SRC/resonant_int.f90 b/POTATO/SRC/resonant_int.f90 index a9ba3dc1..93462702 100644 --- a/POTATO/SRC/resonant_int.f90 +++ b/POTATO/SRC/resonant_int.f90 @@ -134,44 +134,57 @@ subroutine integrate_class_resonances ! ! Computes sum over resonances $x=x^{res}_{(\bm,k)}$ in Eq.(104) for a given class $k$ ! - use find_all_roots_mod, only : customgrid,ncustom,xcustom,nroots,roots - use get_matrix_mod, only : relmargin,iclass + use find_all_roots_mod, only : customgrid,ncustom,xcustom,nroots,roots, & + relerr_allroots,fail_fast,max_roots_abort + use get_matrix_mod, only : iclass use form_classes_doublecount_mod, only : ifuntype,R_class_beg,R_class_end,sigma_class use resint_mod, only : nmodes,marr,narr,twopim2,rm3,delint_mode,respoints_jp - use orbit_dim_mod, only : neqm - use global_invariants, only : toten,perpinv,cE_ref,Phi_eff - use sample_matrix_mod, only : npoi,xarr + use orbit_dim_mod, only : neqm,fb_max_abs_delphi,fb_max_tau + use global_invariants, only : toten,perpinv,cE_ref,Phi_eff,dtau +! xbeg,xend are the search bounds set (and clipped to the |Delta_phi_b| cap) by +! the preceding sample_class_doublecount, so the root search runs over exactly +! the sampled interval. + use sample_matrix_mod, only : npoi,xarr,xbeg,xend,eps use logging_mod, only : tee_message ! implicit none ! double precision, parameter :: pi=3.14159265358979d0, twopi=2.d0*pi, & pi32_over4m=-0.25d0*pi**1.5d0 + integer, parameter :: max_roots_per_class_mode=32 ! - integer :: mode,iroot,ierr - double precision :: relmargin_loc,widthclass,xbeg,xend + integer :: mode,iroot,ierr,max_roots_abort_save double precision :: rescond,dresconddx,psiast,dpsiastdx,taub,delphi double precision :: one_res,sigma,delta_R,Rst,xi,dxi_dx,dpsiast_dRst,absHn2 double precision :: toten_loc,perpinv_loc,fmaxw,A1ast,A2ast double precision :: dens,temp,ddens,dtemp,phi_elec,dPhi_dpsi + double precision :: relerr_allroots_save,fb_max_abs_delphi_save,fb_max_tau_save double precision, dimension(neqm) :: z + character(len=256) :: msg + logical :: fail_fast_save ! toten_loc=toten perpinv_loc=perpinv ! sigma=sigma_class(iclass) delta_R=R_class_end(iclass)-R_class_beg(iclass) -!old=> relmargin_loc=1.d-8 -!old=> widthclass=1.d0 - relmargin_loc=relmargin !<=new - widthclass=abs(R_class_end(iclass)/R_class_beg(iclass)-1.d0) !<=new -! - call classbounds(ifuntype(iclass),relmargin_loc,widthclass,xbeg,xend) ! - customgrid=.true. - ncustom=npoi - allocate(xcustom(ncustom)) - xcustom=xarr + ncustom=count(xarr.ge.xbeg .and. xarr.le.xend) + customgrid=ncustom.ge.2 + if(customgrid) then + allocate(xcustom(ncustom)) + xcustom=pack(xarr,xarr.ge.xbeg .and. xarr.le.xend) + endif + relerr_allroots_save=relerr_allroots + relerr_allroots=max(1.d-8,1.d-2*eps) + fail_fast_save=fail_fast + fail_fast=.true. + max_roots_abort_save=max_roots_abort + max_roots_abort=max_roots_per_class_mode+1 + fb_max_abs_delphi_save=fb_max_abs_delphi + fb_max_tau_save=fb_max_tau + fb_max_abs_delphi=1.5d0*twopi*dble(maxval(abs(marr)))/dble(abs(narr(1))) + fb_max_tau=200.d0*dtau ! ! do mode=1,nmodes @@ -180,13 +193,25 @@ subroutine integrate_class_resonances delint_mode(mode)=0.d0 ! call find_all_roots(get_rescond,xbeg,xend,ierr) +! + if(ierr.eq.2 .or. nroots.gt.max_roots_per_class_mode) then + write(msg,'(A,ES12.4,A,ES12.4,A,I0,A,I0,A,I0,A,I0)') & + 'skipped oscillatory roots: toten=',toten,' perpinv=',perpinv, & + ' iclass=',iclass,' mode=',mode,' m=',marr(mode), & + ' nroots=',nroots + call tee_message(trim(msg)) + respoints_jp(mode,iclass)%nrespoi=0 + respoints_jp(mode,iclass)%toten_res=toten + respoints_jp(mode,iclass)%perpinv_res=perpinv + cycle + endif ! if(ierr.ne.0) then - call tee_message( & - 'integrate_class_resonances: error in find_all_roots') - customgrid=.false. - deallocate(xcustom) - return + call tee_message('integrate_class_resonances: skipped unpolished roots') + respoints_jp(mode,iclass)%nrespoi=0 + respoints_jp(mode,iclass)%toten_res=toten + respoints_jp(mode,iclass)%perpinv_res=perpinv + cycle endif ! respoints_jp(mode,iclass)%nrespoi=nroots @@ -254,7 +279,12 @@ subroutine integrate_class_resonances enddo ! customgrid=.false. - deallocate(xcustom) + relerr_allroots=relerr_allroots_save + fail_fast=fail_fast_save + max_roots_abort=max_roots_abort_save + fb_max_abs_delphi=fb_max_abs_delphi_save + fb_max_tau=fb_max_tau_save + if(allocated(xcustom)) deallocate(xcustom) ! !------------ contains @@ -421,6 +451,7 @@ subroutine resonant_torque double precision :: omdens,trapez_fac,perpinv_max double precision :: torque_int,torque_int_loc double precision :: xjperp,xenerg,totxint,step_energ + double precision :: w_energ,dw_energ double precision :: time_beg,time_end double precision :: dens, temp, ddens, dtemp character(len=256) :: msg @@ -483,8 +514,10 @@ subroutine resonant_torque ! torquebox=0.d0 ! - step_energ=toten_range/dble(nenerg) !integration step over total energy -!step_energ=0.22521463755624047d0 +! Energy grid uniform in velocity w = sqrt(toten - toten_min) ~ v, so the scan +! starts at low energy while still covering the same maximum v/vth range as +! NEO-RT. The energy integral keeps its value through d(toten) = 2 w dw. + dw_energ=sqrt(toten_range)/dble(nenerg) ! if(adaptive_jperp) then open(1901,file='subint_ofH0int_104_vsJperp_fromresp.dat') @@ -493,12 +526,10 @@ subroutine resonant_torque open(1902,file='subint_ofH0int_104_vsJperp_equi.dat') endif ! -! do ienerg=1,nenerg do ienerg=2,nenerg -! do ienerg= 20,20 !10,10 !20,20 !<=fix energy for debugging - xenerg=(dble(ienerg)-0.5d0)/dble(nenerg) - toten=toten_min+toten_range*xenerg -!toten = -3.1211921097605737d0 + w_energ=(dble(ienerg)-0.5d0)*dw_energ + toten=toten_min+w_energ**2 + step_energ=2.d0*w_energ*dw_energ write(msg, '(A,ES22.14)') 'toten = ', toten call tee_message(trim(msg)) ! diff --git a/POTATO/SRC/sub_potato.f90 b/POTATO/SRC/sub_potato.f90 index b4d77f56..1aeeb58c 100644 --- a/POTATO/SRC/sub_potato.f90 +++ b/POTATO/SRC/sub_potato.f90 @@ -81,6 +81,18 @@ module orbit_dim_mod logical :: write_orb=.false. integer :: iunit1=100,next,numbasef double precision :: Rorb_max +! Orbit-integrator tolerances used by find_bounce. 6 significant digits is far +! tighter than the ~10% code-to-code benchmark and the eps~1d-3 class grid need, +! while keeping orbit integration (the dominant cost) affordable. The class +! sampler relaxes these further for its grid build (see sample_class_doublecount). + double precision :: fb_relerr=1d-6, fb_abserr=1d-8 + double precision :: fb_max_abs_delphi=huge(1d0) + double precision :: fb_max_tau=huge(1d0) +! Keep the outer X-point -> rho_pol-boundary classes (clip their search interval +! to the |Delta_phi_b| resonance cap instead of dropping them) so resonances are +! found out to the edge, matching NEO-RT. Default off: keeping them integrates +! many long-period orbits and needs the find_bounce low-energy speed-up first. + logical :: clip_resonance_classes=.false. end module orbit_dim_mod ! !------------------------------------------------------ @@ -111,194 +123,173 @@ subroutine find_bounce(next,velo_ext,dtau_in,z_eqm,taub,delphi,extraset,ierr) ! extraset(next) - extra integrals along the orbit (inout) ! ierr - error flag, 0 = success, 1 = orbit left domain (output) ! - use orbit_dim_mod, only : neqm,write_orb,iunit1,Rorb_max + use orbit_dim_mod, only : neqm,write_orb,iunit1,Rorb_max,fb_relerr, & + fb_abserr,fb_max_abs_delphi,fb_max_tau use field_eq_mod, only : ierrfield use logging_mod, only : log_message + use fortnum_ode, only : ODE_EVENT_RISING + use fortnum_ode_vode, only : vode_state_t, vode_init, vode_integrate_to + use fortnum_status, only : fortnum_status_t, FORTNUM_OK ! implicit none ! -! nousecut - way to close the orbit during primary search: -! .true. - without using Poincare cut, can be used for general Phi distributions -! .false. - with using the cut, valid for Phi=Phi(psi) only -! logical, parameter :: nousecut=.true. - logical, parameter :: nousecut=.false. +! relative/absolute error of orbit integrator (module-settable; the sampler +! relaxes them for speed, see orbit_dim_mod): + double precision :: relerr + double precision :: abserr +! Step budget per orbit. A confined orbit closes in O(1e2..1e3) variable-order +! Adams steps; even the longest kept orbits (|Delta_phi_b| up to ~1.5x the +! resonance cap, tens of toroidal transits) close well under this. An orbit on +! the trapped-passing separatrix has a diverging bounce time and never closes, so +! cap the budget and skip it (ierr=1) instead of grinding to 50000 steps -- those +! non-closing grazers are the bulk of the sampler cost and are discarded by the +! class clip anyway. + integer, parameter :: max_orbit_steps=10000 +! Finite non-writing chunks let resonance sampling reject non-closing tails once +! their toroidal displacement is beyond the largest requested resonance. + double precision, parameter :: nonwrite_chunk_factor=20.d0 ! -! maximum number of Newton iterations for closing the orbit: - integer, parameter :: niter=20 -! -! relative error of orbit integrator: - double precision, parameter :: relerr=1d-10 !8 - integer, intent(in) :: next double precision, intent(in) :: dtau_in double precision, dimension(neqm), intent(in) :: z_eqm double precision, dimension(next), intent(inout) :: extraset double precision, intent(out) :: taub, delphi integer, intent(out) :: ierr - - logical :: firstpass - integer :: ndim, iter - double precision :: dtau - double precision :: dL2_pol,dL2_pol_start,dtau_newt,r_prev,z_prev - double precision :: tau0,RNorm,ZNorm,vnorm,dnorm,vel_pol,dL2_pol_min - double precision :: dZ_dR,sign_delZ,Z_tmp - double precision, dimension(neqm+next) :: z,z_start,vz - character(len=256) :: msg +! + integer :: ndim, nleg + double precision :: t_root, dt_leg, tau_max, vel_pol, RNorm, ZNorm + double precision, dimension(neqm+next) :: z_start, vz + double precision, dimension(neqm+next) :: atol + double precision, dimension(:), allocatable :: yout + type(vode_state_t) :: vstate + type(fortnum_status_t) :: status + logical :: root_found ! external velo_ext ! ndim = neqm+next ierr = 0 -! - z(1:neqm)=z_eqm - if(next.gt.0) then - z(neqm+1:ndim)=extraset - endif - Rorb_max=z(1) -! -! dtau=dtau_in/max(z(4),1d-3) - dtau=dtau_in -! -! Primary search: -! - z_start=z -! - call velo_ext(dtau,z,vz) -! -! unit 2D vector in the direction of the guiding center velocity in RZ-plane: -! - vel_pol=sqrt(vz(1)**2+vz(3)**2) - RNorm=vz(1)/vel_pol - ZNorm=vz(3)/vel_pol -! - tau0=0.d0 -! - if(write_orb) write (iunit1,*) z(1:neqm),vz(5) -! -! first step: -! - call odeint_allroutines(z,ndim,tau0,dtau,relerr,velo_ext) - if(ierrfield.ne.0) then - write(msg,'(A,2ES14.6,A,2ES14.6)') & - 'find_bounce: orbit left domain at R,Z=', & - z(1),z(3),' started from R,Z=',z_start(1),z_start(3) - write(*,'(A)') trim(msg) - call log_message(trim(msg)) - ierr = 1 + relerr = fb_relerr + abserr = fb_abserr + atol = abserr +! + z_start(1:neqm)=z_eqm + if(next.gt.0) z_start(neqm+1:ndim)=extraset + Rorb_max=z_start(1) +! +! Bounce detection by RETURN TO THE START POINT, not by a Poincare-cut crossing. +! The cut crossing fails for deeply trapped orbits (a small banana never reaches +! the cut) -- exactly what stalled the old fixed-step search at low v/vth. The +! orbit closes when its displacement projected on the start velocity direction +! returns through zero moving in the start direction (a rising zero). This is +! the original start-point closing condition and is well defined for trapped, +! deeply trapped and passing orbits alike. fortnum's variable-order Adams +! integrator locates the root on its own interpolant, so taub is exact and no +! separate Newton closing is needed. + call velo_ext(0.d0, z_start, vz) + vel_pol = sqrt(vz(1)**2 + vz(3)**2) + RNorm = vz(1)/vel_pol + ZNorm = vz(3)/vel_pol +! + call vode_init(vstate, ndim, 0.d0, z_start) + vstate%max_steps = max_orbit_steps + dt_leg = merge(dtau_in, nonwrite_chunk_factor*dtau_in, write_orb) + tau_max = min(2000.d0*dtau_in,fb_max_tau) +! + if(write_orb) write (iunit1,*) z_start(1:neqm), 0.d0 +! +! Step off the start plane once without event detection (the start sits exactly +! on it: a rising zero at t=0 that must not be mistaken for the return). + call vode_integrate_to(rhs_w, vstate, dtau_in, relerr, atol, yout, status) + if(status%code /= FORTNUM_OK .or. ierrfield /= 0) then + call report_left(z_start, yout) return endif -! - taub=dtau - if(nousecut) then -! initialize sqrt(2) of the poloidal length of the step - dL2_pol=2.d0*(z(1)-z_start(1))**2+(z(3)-z_start(3))**2 - else -! initialize Poincare cut crossing check -! - call get_poicut(z(1),Z_tmp,dZ_dR) -! - sign_delZ=sign(1.d0,z(3)-Z_tmp) - firstpass=.true. - dL2_pol=2.d0*(z(1)-z_start(1))**2+(z(3)-z_start(3))**2 - endif -! - if(write_orb) then -! - call velo_ext(dtau,z,vz) -! - write (iunit1,*) z(1:neqm),vz(5) - endif -! - do - r_prev=z(1) - z_prev=z(3) -! - call odeint_allroutines(z,ndim,tau0,dtau,relerr,velo_ext) - if(ierrfield.ne.0) then - write(msg,'(A,2ES14.6,A,2ES14.6)') & - 'find_bounce: orbit left domain at R,Z=', & - z(1),z(3),' started from R,Z=',z_start(1),z_start(3) - write(*,'(A)') trim(msg) - call log_message(trim(msg)) - ierr = 1 + Rorb_max=max(Rorb_max, yout(1)) + if(write_orb) call write_orbit_point(yout) +! +! Integrate to the return. write_orb dumps the orbit, advancing in dtau_in-sized +! legs and writing each; otherwise one far leg suffices since the event stops it +! at the bounce. + root_found=.false. + do nleg=1,2*max_orbit_steps + call vode_integrate_to(rhs_w, vstate, vstate%tn+dt_leg, relerr, atol, & + yout, status, event=return_event, & + event_dir=ODE_EVENT_RISING, t_root=t_root, & + root_found=root_found) + if(status%code /= FORTNUM_OK .or. ierrfield /= 0) then + call report_left(z_start, yout) return endif -! - taub=taub+dtau - if(nousecut) then -! check if poloidal distance to the starting point is smaller than -! sqrt(2) of the poloidal length of the step: - dL2_pol=2.d0*((z(1)-r_prev)**2+(z(3)-z_prev)**2) - dL2_pol_start=(z(1)-z_start(1))**2+(z(3)-z_start(3))**2 - if(dL2_pol_start.lt.dL2_pol) exit - else -! check if Poincare cut has been crossed -! - call get_poicut(z(1),Z_tmp,dZ_dR) -! - if(sign_delZ*(z(3)-Z_tmp).lt.0.d0) then - if(firstpass) then -! first crossing (continue integration) - firstpass=.false. - sign_delZ=-sign_delZ - else -! second crossing (stop integration) - exit - endif - endif + Rorb_max=max(Rorb_max, yout(1)) + if(write_orb) call write_orbit_point(yout) + if(root_found) exit + if(abs(yout(2)-z_start(2)).gt.fb_max_abs_delphi) then + call report_left(z_start, yout) + return endif -! - if(write_orb) then -! - call velo_ext(dtau,z,vz) -! - write (iunit1,*) z(1:neqm),vz(5) + if(.not.write_orb) then + if(vstate%tn.ge.tau_max) then + call report_left(z_start, yout) + return + endif endif - Rorb_max=max(Rorb_max,z(1)) enddo + if(.not.root_found) then + call report_left(z_start, yout) + return + endif ! -! End primary search -! -! Newton adjustment -! - do iter=1,niter + taub = t_root + if(next.gt.0) extraset=yout(neqm+1:ndim) + delphi = yout(2)-z_start(2) ! - call velo_ext(dtau,z,vz) + if(write_orb) write (iunit1,*) 'NaN NaN NaN NaN NaN NaN' ! - vnorm=vz(1)*RNorm+vz(3)*ZNorm - dnorm=(z_start(1)-z(1))*RNorm+(z_start(3)-z(3))*ZNorm - if(dnorm**2.lt.dL2_pol*relerr) exit - dtau_newt=dnorm/vnorm + contains ! - call odeint_allroutines(z,ndim,tau0,dtau_newt,relerr,velo_ext) - if(ierrfield.ne.0) then + subroutine rhs_w(t_, y_, dydt_, ctx_) + double precision, intent(in) :: t_ + double precision, intent(in) :: y_(:) + double precision, intent(out) :: dydt_(:) + class(*), intent(in), optional :: ctx_ + if(present(ctx_)) continue + call velo_ext(t_, y_, dydt_) + end subroutine rhs_w +! + function return_event(t_, y_, ctx_) result(g) + double precision, intent(in) :: t_ + double precision, intent(in) :: y_(:) + class(*), intent(in), optional :: ctx_ + double precision :: g + if(present(ctx_)) continue + if(t_ < 0.d0) continue + g = (y_(1)-z_start(1))*RNorm + (y_(3)-z_start(3))*ZNorm + end function return_event +! + subroutine write_orbit_point(y_) + double precision, intent(in) :: y_(:) + call velo_ext(0.d0, y_, vz) + write (iunit1,*) y_(1:neqm), vz(5) + end subroutine write_orbit_point +! + subroutine report_left(zs, zc) + double precision, intent(in) :: zs(:), zc(:) + character(len=256) :: msg +! Orbits grazing the X-point / trapped-passing separatrix do not close. During +! phase-space sampling this is a normal, expected outcome (the class is clipped +! to the closing region), and there are many such orbits, so the per-orbit +! string format + stdout write + log_message dominated the runtime (top self-time +! in perf). Emit the diagnostic only in explicit orbit-dump mode. + if(write_orb) then write(msg,'(A,2ES14.6,A,2ES14.6)') & - 'find_bounce: orbit left domain at R,Z=', & - z(1),z(3),' started from R,Z=',z_start(1),z_start(3) + 'find_bounce: orbit not closed from R,Z=', zs(1), zs(3), & + ' last R,Z=', zc(1), zc(3) write(*,'(A)') trim(msg) call log_message(trim(msg)) - ierr = 1 - return endif -! - taub=taub+dtau_newt - enddo -! - if(next.gt.0) then - extraset=z(neqm+1:ndim) - endif -! - if(write_orb) then -! - call velo_ext(dtau,z,vz) -! - write (iunit1,*) z(1:neqm),vz(5) - write (iunit1,*) 'NaN NaN NaN NaN NaN NaN' - endif -! -! End Newton adjustment -! - delphi=z(2)-z_start(2) + ierr = 1 + end subroutine report_left ! end subroutine find_bounce ! @@ -2395,7 +2386,7 @@ subroutine get_matrix_doublecount ! of class parameter x for adaptive refinement of interpolation grid over x ! use sample_matrix_mod, only : n1,n2,x,amat - use orbit_dim_mod, only : neqm,next + use orbit_dim_mod, only : neqm,next,clip_resonance_classes use get_matrix_mod, only : iclass use global_invariants, only : dtau,toten,perpinv use form_classes_doublecount_mod, only : ifuntype,R_class_beg,R_class_end,sigma_class @@ -2421,8 +2412,17 @@ subroutine get_matrix_doublecount call starter_doublecount(toten,perpinv,sigma,Rst, & psiast,dpsiast_dRst,z,ierr) ! +! No valid orbit here (start off the cut or negative kinetic energy). With the +! class clip on, mark an out-of-cap sentinel so clip_bounds_to_dphi discards the +! node; otherwise leave the original drop-on-fail behaviour. No per-node print -- +! it floods the sampler. if(ierr.ne.0) then - print *,'get_matrix: error in starter' + if(clip_resonance_classes) then + amat(1,1)=0.d0 + amat(2,1)=1.d30 + amat(3,1)=1.d30 + if(next.gt.0) amat(4:3+next,1)=0.d0 + endif return endif ! @@ -2433,7 +2433,15 @@ subroutine get_matrix_doublecount if(fullbounce) then ! call find_bounce(next,velo,dtau,z,taub,delphi,extraset,ierr) - if(ierr.ne.0) return + if(ierr.ne.0) then + if(clip_resonance_classes) then + amat(1,1)=psiast + amat(2,1)=1.d30 + amat(3,1)=1.d30 + if(next.gt.0) amat(4:3+next,1)=0.d0 + endif + return + endif ! else ! @@ -2447,7 +2455,15 @@ subroutine get_matrix_doublecount extraset=0.d0 ! call find_bounce(next,velo_pphint,dtau,z,taub,delphi,extraset,ierr) - if(ierr.ne.0) return + if(ierr.ne.0) then + if(clip_resonance_classes) then + amat(1,1)=psiast + amat(2,1)=1.d30 + amat(3,1)=1.d30 + if(next.gt.0) amat(4:3+next,1)=0.d0 + endif + return + endif ! endif ! @@ -2471,16 +2487,22 @@ subroutine sample_class_doublecount(iunit,ierr) ! use sample_matrix_mod, only : nlagr,n1,n2,itermax,eps,xbeg,xend, & npoi,xarr,amat_arr - use orbit_dim_mod, only : next,numbasef + use orbit_dim_mod, only : next,numbasef,fb_relerr,fb_abserr, & + fb_max_abs_delphi,fb_max_tau, & + clip_resonance_classes use get_matrix_mod, only : relerror,relmargin,iclass use global_invariants, only : dtau,toten,perpinv,sigma use form_classes_doublecount_mod, only : ifuntype,R_class_beg,R_class_end + use potato_input_mod, only : m_min,m_max,n_tor use cc_mod, only : dowrite ! implicit none ! + double precision, parameter :: twopi=6.28318530717958648d0 integer :: iunit,ierr,i - double precision :: psiastbeg,psiastend,widthclass + double precision :: psiastbeg,psiastend,widthclass,dphi_cap + double precision :: fb_relerr_save,fb_abserr_save + double precision :: fb_max_abs_delphi_save,fb_max_tau_save ! external :: get_matrix_doublecount ! @@ -2506,17 +2528,140 @@ subroutine sample_class_doublecount(iunit,ierr) call classbounds(ifuntype(iclass),relmargin,widthclass,xbeg,xend) ! eps=relerror +! +! The grid interpolates to eps~1d-3, so the orbit integration needs no more than +! that; relax the integrator tolerance for the (cost-dominant) sampling, then +! restore the tight default for the accurate weight integration in pertham. + fb_relerr_save=fb_relerr + fb_abserr_save=fb_abserr + fb_max_abs_delphi_save=fb_max_abs_delphi + fb_max_tau_save=fb_max_tau + fb_relerr=1.d-4 + fb_abserr=1.d-6 + if(clip_resonance_classes) then + fb_max_abs_delphi=1.5d0*twopi*dble(max(abs(m_min),abs(m_max))) & + /dble(abs(n_tor)) + fb_max_tau=200.d0*dtau + dphi_cap=fb_max_abs_delphi + call preclip_class_to_dphi(xbeg,xend,dphi_cap) + endif ! call sample_matrix(get_matrix_doublecount,ierr) ! - if(dowrite) then +! Classes touching the X-point / trapped-passing separatrix carry the outer +! (high-rho) resonances, but there tau_b and Delta_phi_b diverge and the +! interpolant cannot converge within itermax, so the whole class was dropped +! (ierr=2) and its resonances lost. No resonance with |m|<=m_max exists past +! |Delta_phi_b| = 2 pi m_max / n_tor, so re-clip the search interval to the grid +! span where |Delta_phi_b| stays within 1.5x that cap (a buffer that keeps the +! m_max root while cutting the divergent tail) and re-sample the bounded +! sub-interval. If the clipped class still does not converge, its interpolant is +! not reliable enough for per-mode root finding, so leave ierr set and skip the +! class upstream. + if(ierr.ne.0 .and. clip_resonance_classes) then + call clip_bounds_to_dphi(xbeg,xend,dphi_cap) + if(xend.gt.xbeg) then + call sample_matrix(get_matrix_doublecount,ierr) + endif + endif +! + fb_relerr=fb_relerr_save + fb_abserr=fb_abserr_save + fb_max_abs_delphi=fb_max_abs_delphi_save + fb_max_tau=fb_max_tau_save +! + if(dowrite) then print *,'npoi = ',npoi do i=1,npoi write(iunit,*) xarr(i),amat_arr(:,1,i) enddo endif ! - end subroutine sample_class_doublecount + contains +! + subroutine preclip_class_to_dphi(xb,xe,dphi_cap) + use sample_matrix_mod, only : n1,n2,x,amat + implicit none + integer, parameter :: nprobe=9 + integer :: i,ilo,ihi + double precision, intent(inout) :: xb,xe + double precision, intent(in) :: dphi_cap + double precision :: dx + double precision :: xprobe(nprobe) + double precision :: dphi_probe(nprobe) + external :: get_matrix_doublecount +! + if(xe.le.xb) return + if(.not.allocated(amat)) allocate(amat(n1,n2)) +! + dx=(xe-xb)/dble(nprobe-1) + do i=1,nprobe + xprobe(i)=xb+dx*dble(i-1) + x=xprobe(i) + call get_matrix_doublecount + dphi_probe(i)=abs(amat(3,1)) + enddo +! + ilo=0 + ihi=0 + do i=1,nprobe + if(dphi_probe(i).le.dphi_cap) then + if(ilo.eq.0) ilo=i + ihi=i + endif + enddo +! + if(ilo.eq.0) then + xe=xb + return + endif +! + if(ilo.gt.1) ilo=ilo-1 + if(ihi.lt.nprobe) ihi=ihi+1 + xb=xprobe(ilo) + xe=xprobe(ihi) + end subroutine preclip_class_to_dphi +! + end subroutine sample_class_doublecount +! +!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +! + subroutine clip_bounds_to_dphi(xb,xe,dphi_cap) +! +! Shrinks the class search interval [xb,xe] to the span of the (partially +! sampled) grid where |Delta_phi_b| <= dphi_cap, discarding the separatrix tail +! where Delta_phi_b diverges. Uses the grid (xarr) and matrix (amat_arr, row 3 +! = Delta_phi_b) left by the preceding sample_matrix call, so it needs no extra +! orbit integration. The physical (resonance-bearing) part is contiguous, so +! the first and last in-cap grid nodes bracket it. +! + use sample_matrix_mod, only : npoi,xarr,amat_arr +! + implicit none +! + double precision, intent(inout) :: xb,xe + double precision, intent(in) :: dphi_cap + integer :: i,ilo,ihi +! + ilo=0 + ihi=0 + do i=1,npoi + if(abs(amat_arr(3,1,i)).le.dphi_cap) then + if(ilo.eq.0) ilo=i + ihi=i + endif + enddo +! + if(ilo.eq.0) then +! no node within the cap: nothing to resolve, collapse the interval to skip it + xe=xb + return + endif +! + xb=xarr(ilo) + xe=xarr(ihi) +! + end subroutine clip_bounds_to_dphi ! !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! diff --git a/POTATO/SRC/tt.f90 b/POTATO/SRC/tt.f90 index 4ca11d4e..050aaada 100644 --- a/POTATO/SRC/tt.f90 +++ b/POTATO/SRC/tt.f90 @@ -2,7 +2,8 @@ program classify_orbits ! use parmot_mod, only : rmu,ro0 - use orbit_dim_mod, only : write_orb,iunit1,numbasef,neqm + use orbit_dim_mod, only : write_orb,iunit1,numbasef,neqm, & + orbit_clip_resonance_classes => clip_resonance_classes use get_matrix_mod, only : iclass use global_invariants, only : dtau,toten,perpinv,cE_ref,Phi_eff use bounds_fixpoints_mod, only : nregions @@ -20,6 +21,7 @@ program classify_orbits npoicut, m_min, m_max, n_tor, & toten_plot, perpinv_plot, profile_file, & edge_extension, & + input_clip_resonance_classes => clip_resonance_classes, & orbit_Rstart, orbit_Zstart, orbit_lambda, & freq_Rmin, freq_Rmax, freq_n use field_eq_mod, only : allow_sol, psi_axis, psi_sep @@ -59,6 +61,7 @@ program classify_orbits ! ! Allow orbits to cross the separatrix into the scrape-off layer when requested: allow_sol = edge_extension + orbit_clip_resonance_classes = input_clip_resonance_classes ! iunit=71 ! From ebf9f622f4ef7f3ef4b1e90769c5228201a69f52 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 23 Jun 2026 18:22:32 +0200 Subject: [PATCH 06/10] Build fortnum permissively (GNU implicit-typing ext) under f2018 project std --- POTATO/CMakeLists.txt | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/POTATO/CMakeLists.txt b/POTATO/CMakeLists.txt index 377ba173..9a4cdb46 100644 --- a/POTATO/CMakeLists.txt +++ b/POTATO/CMakeLists.txt @@ -32,6 +32,11 @@ if(NOT TARGET fortnum) ) FetchContent_MakeAvailable(fortnum) endif() +# fortnum relies on a GNU implicit-typing extension rejected by f2018; build it +# permissively rather than under the project-wide standard. +if(TARGET fortnum) + target_compile_options(fortnum PRIVATE -std=gnu) +endif() add_subdirectory(SRC) From 2a140a86ee7f6eb96f9bae70c97bca944c8d8a6a Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 23 Jun 2026 18:23:10 +0200 Subject: [PATCH 07/10] Fix sample_matrix DEALLOCATE of unallocated xarr/amat_arr when amat preallocated --- POTATO/SRC/sample_matrix.f90 | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/POTATO/SRC/sample_matrix.f90 b/POTATO/SRC/sample_matrix.f90 index 40c1a33d..fa5e6a71 100644 --- a/POTATO/SRC/sample_matrix.f90 +++ b/POTATO/SRC/sample_matrix.f90 @@ -38,9 +38,12 @@ SUBROUTINE sample_matrix(get_matrix,ierr) h=(xend-xbeg)/npoilag/(1.d0+symm_break) hh=symm_break*h/npoilag ! - if(allocated(amat)) then - DEALLOCATE(amat,xarr,amat_arr) - endif +! amat is the shared sample_matrix_mod array; callers (e.g. the resonance class +! preclip) may have allocated it alone, so deallocate each array on its own +! status rather than gating all three on amat being allocated. + if(allocated(amat)) DEALLOCATE(amat) + if(allocated(xarr)) DEALLOCATE(xarr) + if(allocated(amat_arr)) DEALLOCATE(amat_arr) ! ALLOCATE(amat(n1,n2),xarr(npoi),amat_arr(n1,n2,npoi)) ! From c897294e32d51cdbfe2159f845369a935cd4f7c1 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 23 Jun 2026 20:30:27 +0200 Subject: [PATCH 08/10] Resonance search in canonical frequencies (Omega_b, Omega_phi) Root m*Omega_b + n*Omega_phi = 0 instead of the equivalent Delta_phi_b + 2pi m/n = 0. Same roots (loss-free), but Omega_b = 2pi/tau_b and Omega_phi = Delta_phi_b /tau_b stay finite at the trapped-passing separatrix where tau_b, Delta_phi_b diverge, so the interpolation grid converges instead of exploding. Foundation for the log-factored spline that removes the residual infinite-slope cost. --- POTATO/SRC/resonant_int.f90 | 28 ++++++++++++++++++---------- POTATO/SRC/sub_potato.f90 | 11 +++++++++-- 2 files changed, 27 insertions(+), 12 deletions(-) diff --git a/POTATO/SRC/resonant_int.f90 b/POTATO/SRC/resonant_int.f90 index 93462702..83bc63db 100644 --- a/POTATO/SRC/resonant_int.f90 +++ b/POTATO/SRC/resonant_int.f90 @@ -292,25 +292,33 @@ subroutine integrate_class_resonances ! subroutine get_rescond(x,rescond,dresconddx) ! -! Computes resonance condition $F(x)=\Delta\varphi_b+2\pi m_2/m_3$ -! and its derivative $F^\prime(x)$ for $F(x)=0$ root finding. -! Computes as by-products normalized toroidal momentum $\psi^\ast$, -! bounce time $\tau_b$ and toroidal displacement $\Delta\varphi_b$. +! Resonance condition in the canonical FREQUENCIES, $F(x)=m\,\Omega_b+n\,\Omega_\varphi$, +! and its derivative, for $F(x)=0$ root finding. This is the bounce-drift +! resonance $m\,\Omega_b+n\,\Omega_\varphi=0$ -- the same roots as +! $\Delta\varphi_b+2\pi m/n=0$ (multiply by $\tau_b/2\pi$), but $\Omega_b=2\pi/\tau_b$ +! and $\Omega_\varphi=\Delta\varphi_b/\tau_b$ stay finite at the separatrix where +! $\tau_b\to\infty$, so the interpolated grid (built by get_matrix_doublecount, +! which now stores the frequencies in rows 2,3) converges instead of diverging. +! By-products $\psi^\ast$, $\tau_b$, $\Delta\varphi_b$ are recovered for the weight. ! use sample_matrix_mod, only : n1 ! implicit none ! - double precision :: x,rescond,dresconddx + double precision, parameter :: twopi=6.28318530717958648d0 + double precision :: x,rescond,dresconddx,omega_b,omega_phi,rm2 double precision, dimension(n1) :: vec,dvec ! call interpolate_class_doublecount(x,vec,dvec) ! - psiast=vec(1) ! $\psi^\ast$ - taub=vec(2) ! $\tau_b$ - delphi=vec(3) ! $\Delta\varphi_b$ - rescond=delphi+twopim2/rm3 - dresconddx=dvec(3) + psiast=vec(1) ! $\psi^\ast$ + omega_b=vec(2) ! $\Omega_b=2\pi/\tau_b$ + omega_phi=vec(3) ! $\Omega_\varphi=\Delta\varphi_b/\tau_b$ + taub=twopi/omega_b ! recover $\tau_b$ (finite at the root) + delphi=omega_phi*taub ! recover $\Delta\varphi_b$ + rm2=twopim2/twopi ! poloidal harmonic $m=$ twopim2$/2\pi$ + rescond=rm2*omega_b+rm3*omega_phi ! $m\,\Omega_b+n\,\Omega_\varphi$ + dresconddx=rm2*dvec(2)+rm3*dvec(3) ! end subroutine get_rescond ! diff --git a/POTATO/SRC/sub_potato.f90 b/POTATO/SRC/sub_potato.f90 index 1aeeb58c..35ecbda2 100644 --- a/POTATO/SRC/sub_potato.f90 +++ b/POTATO/SRC/sub_potato.f90 @@ -2395,6 +2395,7 @@ subroutine get_matrix_doublecount ! logical :: fullbounce integer :: ierr,k + double precision, parameter :: twopi_gmd=6.28318530717958648d0 double precision :: psiast,dpsiast_dRst,taub,delphi,xi,dxi_dx,Rst,sigma,delta_R, & tau_fr,dphi_fr double precision, dimension(neqm) :: z @@ -2467,9 +2468,15 @@ subroutine get_matrix_doublecount ! endif ! +! Interpolate the canonical FREQUENCIES, not the bounce time / toroidal shift: +! Omega_b = 2 pi / tau_b and Omega_phi = Delta_phi_b / tau_b stay finite and +! smooth at the separatrix (Omega_b -> 0) where tau_b, Delta_phi_b diverge. The +! resonance m*Omega_b + n*Omega_phi = 0 has the same roots as Delta_phi_b + 2 pi +! m/n = 0 (loss-free), but the polynomial grid converges instead of diverging. +! get_rescond recovers tau_b = 2 pi / Omega_b. amat(1,1)=psiast - amat(2,1)=taub - amat(3,1)=delphi + amat(2,1)=twopi_gmd/taub + amat(3,1)=delphi/taub ! if(next.gt.0) then amat(4:3+next,1)=extraset From a8bbe9bb80c363cb0de58eb3fa768dc37b219ea5 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 23 Jun 2026 20:49:03 +0200 Subject: [PATCH 09/10] Clip class to closing region at the separatrix (frequency form) Mark non-closing orbits Omega_b=0 and clip the class endpoint to the closing boundary x_sep (bisection, only X-point classes whose endpoint does not close). Removes the closing->non-closing discontinuity that made the adaptive grid refine without end, keeping all closing orbits (low-v and high-m alike, no tau cap). Correct (0 fake roots) but the boundary search integrates the slow near-X-point orbits to the cap -- a follow-up will detect the grazing geometrically from the banana-tip flux instead of by integration. --- POTATO/SRC/sub_potato.f90 | 110 +++++++++++++++++++++++--------------- 1 file changed, 67 insertions(+), 43 deletions(-) diff --git a/POTATO/SRC/sub_potato.f90 b/POTATO/SRC/sub_potato.f90 index 35ecbda2..6aa93157 100644 --- a/POTATO/SRC/sub_potato.f90 +++ b/POTATO/SRC/sub_potato.f90 @@ -2418,12 +2418,10 @@ subroutine get_matrix_doublecount ! node; otherwise leave the original drop-on-fail behaviour. No per-node print -- ! it floods the sampler. if(ierr.ne.0) then - if(clip_resonance_classes) then - amat(1,1)=0.d0 - amat(2,1)=1.d30 - amat(3,1)=1.d30 - if(next.gt.0) amat(4:3+next,1)=0.d0 - endif + amat(1,1)=0.d0 + amat(2,1)=0.d0 ! Omega_b=0 : orbit does not close (tau_b -> infinity) + amat(3,1)=0.d0 + if(next.gt.0) amat(4:3+next,1)=0.d0 return endif ! @@ -2435,12 +2433,10 @@ subroutine get_matrix_doublecount ! call find_bounce(next,velo,dtau,z,taub,delphi,extraset,ierr) if(ierr.ne.0) then - if(clip_resonance_classes) then - amat(1,1)=psiast - amat(2,1)=1.d30 - amat(3,1)=1.d30 - if(next.gt.0) amat(4:3+next,1)=0.d0 - endif + amat(1,1)=psiast + amat(2,1)=0.d0 + amat(3,1)=0.d0 + if(next.gt.0) amat(4:3+next,1)=0.d0 return endif ! @@ -2457,12 +2453,10 @@ subroutine get_matrix_doublecount ! call find_bounce(next,velo_pphint,dtau,z,taub,delphi,extraset,ierr) if(ierr.ne.0) then - if(clip_resonance_classes) then - amat(1,1)=psiast - amat(2,1)=1.d30 - amat(3,1)=1.d30 - if(next.gt.0) amat(4:3+next,1)=0.d0 - endif + amat(1,1)=psiast + amat(2,1)=0.d0 + amat(3,1)=0.d0 + if(next.gt.0) amat(4:3+next,1)=0.d0 return endif ! @@ -2545,31 +2539,21 @@ subroutine sample_class_doublecount(iunit,ierr) fb_max_tau_save=fb_max_tau fb_relerr=1.d-4 fb_abserr=1.d-6 - if(clip_resonance_classes) then - fb_max_abs_delphi=1.5d0*twopi*dble(max(abs(m_min),abs(m_max))) & - /dble(abs(n_tor)) - fb_max_tau=200.d0*dtau - dphi_cap=fb_max_abs_delphi - call preclip_class_to_dphi(xbeg,xend,dphi_cap) - endif -! - call sample_matrix(get_matrix_doublecount,ierr) -! -! Classes touching the X-point / trapped-passing separatrix carry the outer -! (high-rho) resonances, but there tau_b and Delta_phi_b diverge and the -! interpolant cannot converge within itermax, so the whole class was dropped -! (ierr=2) and its resonances lost. No resonance with |m|<=m_max exists past -! |Delta_phi_b| = 2 pi m_max / n_tor, so re-clip the search interval to the grid -! span where |Delta_phi_b| stays within 1.5x that cap (a buffer that keeps the -! m_max root while cutting the divergent tail) and re-sample the bounded -! sub-interval. If the clipped class still does not converge, its interpolant is -! not reliable enough for per-mode root finding, so leave ierr set and skip the -! class upstream. - if(ierr.ne.0 .and. clip_resonance_classes) then - call clip_bounds_to_dphi(xbeg,xend,dphi_cap) - if(xend.gt.xbeg) then - call sample_matrix(get_matrix_doublecount,ierr) - endif + dphi_cap=0.d0 ! legacy, unused +! +! The class coordinate runs past the trapped-passing/X-point separatrix into a +! region where the orbit no longer closes (find_bounce fails, Omega_b=0). That +! closing->non-closing step, and the tau_b divergence just before it, are what +! made the adaptive grid refine without end. Clip the endpoint to the closing +! boundary so sample_matrix sees a smooth, all-closing interval. Done with the +! clip flag on or off; only the X-point classes (whose endpoint does not close) +! pay for the boundary search. + call clip_class_to_closing(xbeg,xend) +! + if(xend.gt.xbeg) then + call sample_matrix(get_matrix_doublecount,ierr) + else + ierr=1 endif ! fb_relerr=fb_relerr_save @@ -2585,6 +2569,46 @@ subroutine sample_class_doublecount(iunit,ierr) endif ! contains +! + logical function class_closes(xx) +! True if the orbit at class coordinate xx closes (Omega_b = amat(2,1) > 0). + use sample_matrix_mod, only : n1,n2,x,amat + implicit none + double precision, intent(in) :: xx + external :: get_matrix_doublecount + if(.not.allocated(amat)) allocate(amat(n1,n2)) + x=xx + call get_matrix_doublecount + class_closes=amat(2,1).gt.0.d0 + end function class_closes +! + subroutine clip_class_to_closing(xb,xe) +! Clip the separatrix-side endpoint to the closing boundary. If the endpoint +! orbit already closes the class is smooth -- nothing to do (cheap, the common +! case). Otherwise bisect [closing, non-closing] for the boundary x_sep and set +! the endpoint just inside it, so the sampled interval is all-closing and smooth. + implicit none + double precision, intent(inout) :: xb,xe + double precision :: a,b,xm + integer :: it + if(xe.le.xb) return + if(class_closes(xe)) return ! endpoint closes: smooth class, no clip + if(.not.class_closes(xb)) then ! nothing closes: empty class + xe=xb + return + endif + a=xb ! closes + b=xe ! does not close + do it=1,20 + xm=0.5d0*(a+b) + if(class_closes(xm)) then + a=xm + else + b=xm + endif + enddo + xe=a ! last coordinate that still closes + end subroutine clip_class_to_closing ! subroutine preclip_class_to_dphi(xb,xe,dphi_cap) use sample_matrix_mod, only : n1,n2,x,amat From ef2a1603a89db2df83573c15d68a84b2ae95cea0 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 23 Jun 2026 21:57:00 +0200 Subject: [PATCH 10/10] POTATO: link libneo instead of symlinking a sibling checkout POTATO/SRC carried relative symlinks (../../../libneo/src/...) into a sibling libneo working copy for its field, magfie, spline and kinds/constants sources, tying the build to a machine-specific path and an unpinned, ad-hoc libneo state. Replace them with find_or_fetch(libneo); link LIBNEO::magfie (pulls in the neo core) plus mc_efit for chamb_divB0. binsrc and plag_coeff stay local: POTATO calls them as bare external subroutines, whereas libneo wraps them in *_sub modules. Scope -std=f2018, -Wfatal-errors and -cpp to POTATO's own targets so the fetched libneo, fortplot and fortnum build under their own standards; -cpp in particular makes CMake mis-parse libneo's free-form .f90. Keep -O3 -march=native global so libneo's field hot path stays optimized. All four executables build and link; potato.x classifies orbits and computes bounce frequencies, matching the prior output to ~1e-5 (libneo's -ffast-math vs POTATO's plain -O3). --- POTATO/CMakeLists.txt | 16 +++++++++++++--- POTATO/SRC/CMakeLists.txt | 26 +++++++------------------- POTATO/SRC/amn_mod.f90 | 1 - POTATO/SRC/bdivfree.f90 | 1 - POTATO/SRC/bdivfree_mod.f90 | 1 - POTATO/SRC/chamb_divB0.f90 | 1 - POTATO/SRC/extract_fluxcoord_mod.f90 | 1 - POTATO/SRC/field_c_mod.f90 | 1 - POTATO/SRC/field_divB0.f90 | 1 - POTATO/SRC/field_eq_mod.f90 | 1 - POTATO/SRC/field_mod.f90 | 1 - POTATO/SRC/input_files.f90 | 1 - POTATO/SRC/inthecore_mod.f90 | 1 - POTATO/SRC/libneo_kinds.f90 | 1 - POTATO/SRC/magfie_cyl.f90 | 1 - POTATO/SRC/math_constants.f90 | 1 - POTATO/SRC/period_mod.f90 | 1 - POTATO/SRC/rng.f90 | 1 - POTATO/SRC/spl_three_to_five.f90 | 1 - POTATO/SRC/spline5_RZ.f90 | 1 - POTATO/SRC/theta_rz_mod.f90 | 1 - 21 files changed, 20 insertions(+), 41 deletions(-) delete mode 120000 POTATO/SRC/amn_mod.f90 delete mode 120000 POTATO/SRC/bdivfree.f90 delete mode 120000 POTATO/SRC/bdivfree_mod.f90 delete mode 120000 POTATO/SRC/chamb_divB0.f90 delete mode 120000 POTATO/SRC/extract_fluxcoord_mod.f90 delete mode 120000 POTATO/SRC/field_c_mod.f90 delete mode 120000 POTATO/SRC/field_divB0.f90 delete mode 120000 POTATO/SRC/field_eq_mod.f90 delete mode 120000 POTATO/SRC/field_mod.f90 delete mode 120000 POTATO/SRC/input_files.f90 delete mode 120000 POTATO/SRC/inthecore_mod.f90 delete mode 120000 POTATO/SRC/libneo_kinds.f90 delete mode 120000 POTATO/SRC/magfie_cyl.f90 delete mode 120000 POTATO/SRC/math_constants.f90 delete mode 120000 POTATO/SRC/period_mod.f90 delete mode 120000 POTATO/SRC/rng.f90 delete mode 120000 POTATO/SRC/spl_three_to_five.f90 delete mode 120000 POTATO/SRC/spline5_RZ.f90 delete mode 120000 POTATO/SRC/theta_rz_mod.f90 diff --git a/POTATO/CMakeLists.txt b/POTATO/CMakeLists.txt index 9a4cdb46..a6677b30 100644 --- a/POTATO/CMakeLists.txt +++ b/POTATO/CMakeLists.txt @@ -5,7 +5,7 @@ enable_language(Fortran) set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/../cmake") set(CMAKE_Fortran_MODULE_DIRECTORY ${PROJECT_BINARY_DIR}) -add_compile_options(-cpp -g -fbacktrace -Wfatal-errors -std=f2018) +add_compile_options(-g -fbacktrace) add_compile_options(-O3 -march=native -mtune=native) include(Util) @@ -32,12 +32,22 @@ if(NOT TARGET fortnum) ) FetchContent_MakeAvailable(fortnum) endif() -# fortnum relies on a GNU implicit-typing extension rejected by f2018; build it -# permissively rather than under the project-wide standard. +# fortnum relies on a GNU implicit-typing extension; pin it to the GNU standard. if(TARGET fortnum) target_compile_options(fortnum PRIVATE -std=gnu) endif() +# libneo supplies the field, magfie and spline code POTATO used to carry as its +# own copies. Link the upstream library so fixes (e.g. the stretch_coords +# convex-wall load) propagate instead of drifting in a fork. +find_or_fetch(libneo) + +# Hold POTATO's own sources to f2018, stop at the first diagnostic, and run the C +# preprocessor. Applied after the dependency fetches above so the upstream +# libraries (libneo, fortplot, fortnum) build under their own standards instead +# of inheriting these; -cpp in particular makes CMake mis-parse libneo's .f90. +add_compile_options(-cpp -Wfatal-errors -std=f2018) + add_subdirectory(SRC) add_executable (potato.x diff --git a/POTATO/SRC/CMakeLists.txt b/POTATO/SRC/CMakeLists.txt index c0e875d1..cab99073 100644 --- a/POTATO/SRC/CMakeLists.txt +++ b/POTATO/SRC/CMakeLists.txt @@ -1,25 +1,11 @@ +# The field, magfie, spline and kinds/constants sources live in libneo and are +# linked below (LIBNEO::magfie pulls in libneo's core). Only POTATO-specific +# sources are compiled here. binsrc/plag_coeff stay local: POTATO calls them as +# bare external subroutines, whereas libneo wraps them in *_sub modules. add_library(potato_base - libneo_kinds.f90 logging_mod.f90 potato_input_mod.f90 - math_constants.f90 alpha_lifetime_mod.f90 - chamb_divB0.f90 - magfie_cyl.f90 - period_mod.f90 - input_files.f90 - field_c_mod.f90 - field_mod.f90 - inthecore_mod.f90 - field_eq_mod.f90 - field_divB0.f90 - bdivfree_mod.f90 - amn_mod.f90 - theta_rz_mod.f90 - extract_fluxcoord_mod.f90 - bdivfree.f90 - spl_three_to_five.f90 - spline5_RZ.f90 velo.f90 odeint_allroutines.f plag_coeff.f90 @@ -38,11 +24,13 @@ add_library(potato_base ) find_or_fetch(vode) # vode (dvode_f90_m) uses labeled DO termination, deleted in f2018; build it with -# the permissive legacy standard rather than the project-wide f2018. +# the permissive legacy standard rather than POTATO's f2018. target_compile_options(vode PRIVATE -std=legacy) target_link_libraries(potato_base PUBLIC vode fortnum +LIBNEO::magfie +mc_efit BLAS::BLAS LAPACK::LAPACK HDF5::HDF5 diff --git a/POTATO/SRC/amn_mod.f90 b/POTATO/SRC/amn_mod.f90 deleted file mode 120000 index ebf63216..00000000 --- a/POTATO/SRC/amn_mod.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/magfie/amn_mod.f90 \ No newline at end of file diff --git a/POTATO/SRC/bdivfree.f90 b/POTATO/SRC/bdivfree.f90 deleted file mode 120000 index bdc437df..00000000 --- a/POTATO/SRC/bdivfree.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/magfie/bdivfree.f90 \ No newline at end of file diff --git a/POTATO/SRC/bdivfree_mod.f90 b/POTATO/SRC/bdivfree_mod.f90 deleted file mode 120000 index 207ce50f..00000000 --- a/POTATO/SRC/bdivfree_mod.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/magfie/bdivfree_mod.f90 \ No newline at end of file diff --git a/POTATO/SRC/chamb_divB0.f90 b/POTATO/SRC/chamb_divB0.f90 deleted file mode 120000 index 6bd4433a..00000000 --- a/POTATO/SRC/chamb_divB0.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/MC/chamb_divB0.f90 \ No newline at end of file diff --git a/POTATO/SRC/extract_fluxcoord_mod.f90 b/POTATO/SRC/extract_fluxcoord_mod.f90 deleted file mode 120000 index 460a2739..00000000 --- a/POTATO/SRC/extract_fluxcoord_mod.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/magfie/extract_fluxcoord_mod.f90 \ No newline at end of file diff --git a/POTATO/SRC/field_c_mod.f90 b/POTATO/SRC/field_c_mod.f90 deleted file mode 120000 index 81a304b0..00000000 --- a/POTATO/SRC/field_c_mod.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/magfie/field_c_mod.f90 \ No newline at end of file diff --git a/POTATO/SRC/field_divB0.f90 b/POTATO/SRC/field_divB0.f90 deleted file mode 120000 index ba979e1a..00000000 --- a/POTATO/SRC/field_divB0.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/magfie/field_divB0.f90 \ No newline at end of file diff --git a/POTATO/SRC/field_eq_mod.f90 b/POTATO/SRC/field_eq_mod.f90 deleted file mode 120000 index b4842173..00000000 --- a/POTATO/SRC/field_eq_mod.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/magfie/field_eq_mod.f90 \ No newline at end of file diff --git a/POTATO/SRC/field_mod.f90 b/POTATO/SRC/field_mod.f90 deleted file mode 120000 index 4c8e4f4b..00000000 --- a/POTATO/SRC/field_mod.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/magfie/field_mod.f90 \ No newline at end of file diff --git a/POTATO/SRC/input_files.f90 b/POTATO/SRC/input_files.f90 deleted file mode 120000 index 1fd65563..00000000 --- a/POTATO/SRC/input_files.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/magfie/input_files.f90 \ No newline at end of file diff --git a/POTATO/SRC/inthecore_mod.f90 b/POTATO/SRC/inthecore_mod.f90 deleted file mode 120000 index 6a056a5d..00000000 --- a/POTATO/SRC/inthecore_mod.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/magfie/inthecore_mod.f90 \ No newline at end of file diff --git a/POTATO/SRC/libneo_kinds.f90 b/POTATO/SRC/libneo_kinds.f90 deleted file mode 120000 index b880394d..00000000 --- a/POTATO/SRC/libneo_kinds.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/libneo_kinds.f90 \ No newline at end of file diff --git a/POTATO/SRC/magfie_cyl.f90 b/POTATO/SRC/magfie_cyl.f90 deleted file mode 120000 index a248f9fe..00000000 --- a/POTATO/SRC/magfie_cyl.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/magfie/magfie_cyl.f90 \ No newline at end of file diff --git a/POTATO/SRC/math_constants.f90 b/POTATO/SRC/math_constants.f90 deleted file mode 120000 index b87bd015..00000000 --- a/POTATO/SRC/math_constants.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/math_constants.f90 \ No newline at end of file diff --git a/POTATO/SRC/period_mod.f90 b/POTATO/SRC/period_mod.f90 deleted file mode 120000 index 266c77e8..00000000 --- a/POTATO/SRC/period_mod.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/magfie/period_mod.f90 \ No newline at end of file diff --git a/POTATO/SRC/rng.f90 b/POTATO/SRC/rng.f90 deleted file mode 120000 index 38273697..00000000 --- a/POTATO/SRC/rng.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/MC/rng.f90 \ No newline at end of file diff --git a/POTATO/SRC/spl_three_to_five.f90 b/POTATO/SRC/spl_three_to_five.f90 deleted file mode 120000 index 59516009..00000000 --- a/POTATO/SRC/spl_three_to_five.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/spl_three_to_five.f90 \ No newline at end of file diff --git a/POTATO/SRC/spline5_RZ.f90 b/POTATO/SRC/spline5_RZ.f90 deleted file mode 120000 index 4496bd49..00000000 --- a/POTATO/SRC/spline5_RZ.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/magfie/spline5_RZ.f90 \ No newline at end of file diff --git a/POTATO/SRC/theta_rz_mod.f90 b/POTATO/SRC/theta_rz_mod.f90 deleted file mode 120000 index 91d48b83..00000000 --- a/POTATO/SRC/theta_rz_mod.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../libneo/src/magfie/theta_rz_mod.f90 \ No newline at end of file