test: compare signed-mode bounds with tolerance

POTATO/CMakeLists.txt

@@ -88,6 +88,16 @@ add_test(NAME test_find_all_roots_bracketed
    COMMAND test_find_all_roots_bracketed.x
 )
 
+add_executable(test_signed_toroidal_bound.x
+   SRC/test_signed_toroidal_bound.f90
+)
+target_link_libraries(test_signed_toroidal_bound.x
+   potato_base
+)
+add_test(NAME test_signed_toroidal_bound
+   COMMAND test_signed_toroidal_bound.x
+)
+
 add_executable(test_wall_loss.x
    SRC/test_wall_loss.f90
 )

POTATO/SRC/CMakeLists.txt

@@ -34,6 +34,7 @@ add_library(potato_base
   sorting.f90
   binsrc.f90
   bmod_pert.f90
+  resonance_mode_bounds_mod.f90
   resonant_int.f90
   eqmagprofs.f90
   box_counting.f90

POTATO/SRC/resonance_mode_bounds_mod.f90

@@ -0,0 +1,30 @@
+module resonance_mode_bounds_mod
+    implicit none
+
+contains
+
+    pure function resonant_delphi_bound(m_modes, n_modes) result(bound)
+        integer, intent(in) :: m_modes(:), n_modes(:)
+        double precision :: bound
+        double precision, parameter :: pi = 3.14159265358979d0
+
+        ! The signed toroidal mode remains in m*Omega_b+n*Omega_phi=0 and in the
+        ! Fourier phase.  Only this symmetric search extent is a magnitude.  Using
+        ! signed n here made the native n<0 case produce a negative extent and
+        ! silently discard every otherwise valid resonance contribution.
+        bound = 2.d0*pi*(maxval(abs(dble(m_modes))/abs(dble(n_modes))) &
+            + 1.d0/dble(minval(abs(n_modes))))
+    end function resonant_delphi_bound
+
+    pure logical function canonical_flux_outside_lcfs(psi_star, psi_axis, psi_edge)
+        double precision, intent(in) :: psi_star, psi_axis, psi_edge
+
+        ! Outside means beyond the edge in the axis-to-edge flux direction.  The
+        ! earlier psi_star < psi_edge test was valid only when psi decreased from
+        ! axis to edge and rejected every ITER resonance after field_eq selected
+        ! the opposite native GEQDSK gauge.
+        canonical_flux_outside_lcfs = &
+            (psi_star - psi_edge)*(psi_edge - psi_axis) > 0.d0
+    end function canonical_flux_outside_lcfs
+
+end module resonance_mode_bounds_mod

POTATO/SRC/resonant_int.f90

@@ -57,6 +57,7 @@ module resint_mod
 !$omp threadprivate(nperp_max,delint_mode,respoints_jp,respoints_all, &
 !$omp               respoints_all_tmp,respoint,resline_unit,resline_diag_unit, &
 !$omp               resline_unit_is_private,resline_diag_unit_is_private)
+
 end module resint_mod
 !
 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
@@ -170,7 +171,8 @@ subroutine integrate_class_resonances
     use logging_mod,        only : tee_message
     use interp_cache_mod,   only : interp_cache_reset
     use field_sub,          only : psif
-    use field_eq_mod,       only : psi_sep
+    use field_eq_mod,       only : psi_axis,psi_sep
+    use resonance_mode_bounds_mod, only : canonical_flux_outside_lcfs
     use, intrinsic :: ieee_arithmetic, only : ieee_is_finite
     !
     implicit none
@@ -286,7 +288,7 @@ subroutine integrate_class_resonances
             !
             dpsiastdx=dpsiast_dRst*delta_R*dxi_dx !$\difp{\psi^\ast}{x}$
             !
-            if(psiast_res/psi_sep.gt.1.d0) then
+            if(canonical_flux_outside_lcfs(psiast_res,psi_axis,psi_sep)) then
                 ! SOL resonance (rho_pol > 1): the orbit leaves the field domain,
                 ! so pertham/find_bounce cannot close it and would only grind to
                 ! the integrator cap before returning zero.  It is also outside
@@ -495,6 +497,7 @@ subroutine resonant_torque
     use get_matrix_mod,    only : iclass,delphi_max
     use form_classes_doublecount_mod, only : nclasses
     use orbit_dim_mod,     only : numbasef
+    use resonance_mode_bounds_mod, only : resonant_delphi_bound
     use resint_mod,        only : nmodes,marr,narr,delint_mode,respoints_jp,respoints_all,nperp_max, &
         respoints_all_tmp,respoint,resline_unit,resline_diag_unit, &
         resline_unit_is_private,resline_diag_unit_is_private
@@ -545,7 +548,7 @@ subroutine resonant_torque
     ! Bound the class root search to the resonant range: |delphi_b| = 2*pi*|m|/n
     ! at a resonance, so nothing past max|m|/n can resonate.  One n-step margin
     ! keeps the extreme-m root safely inside the trimmed domain.
-    delphi_max=2.d0*pi*(maxval(abs(dble(marr))/dble(narr))+1.d0/dble(minval(narr)))
+    delphi_max=resonant_delphi_bound(marr,narr)
     write(msg, '(A,ES14.6)') &
         'class root search bounded to |delphi_b| <= ', delphi_max
     call tee_message(trim(msg))

POTATO/SRC/test_signed_toroidal_bound.f90

@@ -0,0 +1,27 @@
+program test_signed_toroidal_bound
+    use resonance_mode_bounds_mod, only: resonant_delphi_bound, &
+        canonical_flux_outside_lcfs
+    implicit none
+
+    integer, parameter :: m_modes(7) = [-3, -2, -1, 0, 1, 2, 3]
+    integer, parameter :: n_positive(7) = 3
+    integer, parameter :: n_negative(7) = -3
+    double precision :: positive_bound, negative_bound
+
+    positive_bound = resonant_delphi_bound(m_modes, n_positive)
+    negative_bound = resonant_delphi_bound(m_modes, n_negative)
+
+    if (positive_bound <= 0.d0) error stop "positive n produced nonpositive bound"
+    if (negative_bound <= 0.d0) error stop "negative n produced nonpositive bound"
+    if (abs(positive_bound-negative_bound) > 1d-14) &
+        error stop "search bound depends on toroidal-mode sign"
+
+    if (canonical_flux_outside_lcfs(0.5d0, 0.d0, 1.d0)) &
+        error stop "inside point rejected for increasing flux"
+    if (.not. canonical_flux_outside_lcfs(1.1d0, 0.d0, 1.d0)) &
+        error stop "outside point accepted for increasing flux"
+    if (canonical_flux_outside_lcfs(-0.5d0, 0.d0, -1.d0)) &
+        error stop "inside point rejected for decreasing flux"
+    if (.not. canonical_flux_outside_lcfs(-1.1d0, 0.d0, -1.d0)) &
+        error stop "outside point accepted for decreasing flux"
+end program test_signed_toroidal_bound