diff --git a/POTATO/CMakeLists.txt b/POTATO/CMakeLists.txt index 094dd649..617d7413 100644 --- a/POTATO/CMakeLists.txt +++ b/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 ) diff --git a/POTATO/SRC/CMakeLists.txt b/POTATO/SRC/CMakeLists.txt index 40381191..82ae5d78 100644 --- a/POTATO/SRC/CMakeLists.txt +++ b/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 diff --git a/POTATO/SRC/resonance_mode_bounds_mod.f90 b/POTATO/SRC/resonance_mode_bounds_mod.f90 new file mode 100644 index 00000000..635d0840 --- /dev/null +++ b/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 diff --git a/POTATO/SRC/resonant_int.f90 b/POTATO/SRC/resonant_int.f90 index 6052e307..f4ab8cd7 100644 --- a/POTATO/SRC/resonant_int.f90 +++ b/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)) diff --git a/POTATO/SRC/test_signed_toroidal_bound.f90 b/POTATO/SRC/test_signed_toroidal_bound.f90 new file mode 100644 index 00000000..24f13393 --- /dev/null +++ b/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