POTATO: link libneo instead of symlinking a sibling checkout

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)
@@ -21,6 +21,33 @@ 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()
+# 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

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,10 +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

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

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,17 +26,19 @@ 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
 !
   integer :: i,iter,nsearch,ierr,ndummy,kx,kxc,k
   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,16 +170,19 @@ 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
         if(farr(i).eq.0.d0) cycle
         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