From 78e7b8eb7f5f6a35e455cb4f2897bed586456b16 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 08:31:33 +0200 Subject: [PATCH 1/7] Continue after finite Newton iteration limits --- DOC/config.md | 28 +++++++++--------- src/orbit_symplectic.f90 | 36 ++++++++++-------------- src/orbit_symplectic_base.f90 | 1 - src/params.f90 | 4 +-- src/simple_gpu.f90 | 8 ++---- test/tests/test_newton_solver_status.f90 | 15 +++++----- 6 files changed, 40 insertions(+), 52 deletions(-) diff --git a/DOC/config.md b/DOC/config.md index e2010f72..3fca2184 100644 --- a/DOC/config.md +++ b/DOC/config.md @@ -12,20 +12,20 @@ values to refine event location independently of the nonlinear solve. * `symplectic_newton_warning_mode` defaults to `.true.`. When an implicit solve - reaches its iteration limit, SIMPLE continues only if the final Newton - correction is finite and no more than ten times the requested relative - tolerance. Each occurrence is reported by the corresponding `*_maxit` - diagnostic. The production RK, symplectic, and full-orbit paths then use the - same terminal convention: if bounded recovery cannot resolve a numerical - microstep, SIMPLE retains any contiguous accepted prefix and holds only the - unresolved remainder, records `warning_step_skip`, advances the clock for - the complete interval, and retries from that valid state on the next - microstep. A warning hold does not terminate or - numerically disqualify the marker: a marker that reaches the requested end - time remains a resolved survivor, while a later physical boundary event is - still a loss. Set the option to `.false.` for strict diagnostic runs that end - only the affected marker at the first exhausted recovery and report a - 101--105 `orbit_exit_code` with `NaN` in `times_lost`. + reaches its iteration limit, SIMPLE commits any finite final Newton iterate + and continues the orbit. Each occurrence is reported by the corresponding + `*_maxit` diagnostic. Recursive recovery is reserved for an unusable step, + such as a failed linear solve or non-finite iterate. The production RK, + symplectic, and full-orbit paths use the same terminal convention: if that + recovery cannot resolve a numerical microstep, SIMPLE retains any contiguous + accepted prefix and holds only the unresolved remainder, records + `warning_step_skip`, advances the clock for the complete interval, and + retries from that valid state on the next microstep. A warning hold does not + terminate or numerically disqualify the marker: a marker that reaches the + requested end time remains a resolved survivor, while a later physical + boundary event is still a loss. Set the option to `.false.` for strict + diagnostic runs that end only the affected marker at the first exhausted + recovery and report a 101--105 `orbit_exit_code` with `NaN` in `times_lost`. * `canonical_grid_nr`, `canonical_grid_ntheta`, and `canonical_grid_nphi` control the Meiss or Albert canonical-map grid. Their defaults are 62, 63, diff --git a/src/orbit_symplectic.f90 b/src/orbit_symplectic.f90 index f8d113e1..fd437e96 100644 --- a/src/orbit_symplectic.f90 +++ b/src/orbit_symplectic.f90 @@ -12,8 +12,7 @@ module orbit_symplectic SYMPLECTIC_STEP_MAXITER, SYMPLECTIC_STEP_LINEAR_SOLVE, & SYMPLECTIC_STEP_BOUNDARY, SYMPLECTIC_STEP_EVENT_NOT_CONVERGED, & SYMPLECTIC_STEP_BOUNDARY_LIMITED, boundary_event_fraction_tolerance, & - boundary_event_radial_tolerance, symplectic_newton_warning_mode, & - symplectic_newton_warning_factor + boundary_event_radial_tolerance, symplectic_newton_warning_mode use orbit_symplectic_quasi, only: orbit_timestep_quasi, timestep_expl_impl_euler_quasi, & timestep_impl_expl_euler_quasi, timestep_midpoint_quasi, orbit_timestep_rk45, & timestep_rk_gauss_quasi, timestep_rk_lobatto_quasi @@ -50,20 +49,15 @@ pure logical function finite_newton_iterate(x) finite_newton_iterate = .true. end function finite_newton_iterate -logical function accept_warning_maxiter(x, xlast, tolref, rtol) - ! Continue directly only when the last Newton correction is finite and close - ! to the requested tolerance. Larger finite iterates are not trustworthy - ! states: the warning-mode driver retries them at smaller substeps and, if - ! every recovery is exhausted, advances the clock from the last valid state. - real(dp), intent(in) :: x(:), xlast(:), tolref(:), rtol - - accept_warning_maxiter = .false. - if (.not. symplectic_newton_warning_mode .or. rtol <= 0.0_dp) return - if (.not. finite_newton_iterate(x)) return - if (.not. finite_newton_iterate(xlast)) return - if (.not. finite_newton_iterate(tolref)) return - accept_warning_maxiter = all(abs(x - xlast) <= & - symplectic_newton_warning_factor*rtol*abs(tolref)) +logical function accept_warning_maxiter(x) + ! Warning mode preserves the historical SIMPLE continuation contract: a + ! Newton iteration-limit event is diagnostic, not a rejected step. Commit + ! every finite final iterate and leave recursive recovery for genuinely + ! unusable states such as failed linear solves or non-finite iterates. + real(dp), intent(in) :: x(:) + + accept_warning_maxiter = symplectic_newton_warning_mode .and. & + finite_newton_iterate(x) end function accept_warning_maxiter subroutine advance_symplectic_with_retry(si, f, stepper, status, & @@ -596,7 +590,7 @@ recursive subroutine newton1(si, f, x, maxit, xlast, status) end if else call count_event(EVT_NEWTON1_MAXIT) - if (accept_warning_maxiter(x, xlast, tolref, si%rtol)) & + if (accept_warning_maxiter(x)) & status = SYMPLECTIC_STEP_OK end if end subroutine @@ -691,7 +685,7 @@ recursive subroutine newton2(si, f, x, atol, rtol, maxit, xlast, status) end if else call count_event(EVT_NEWTON2_MAXIT) - if (accept_warning_maxiter(x, xlast, tolref, rtol)) & + if (accept_warning_maxiter(x)) & status = SYMPLECTIC_STEP_OK end if end subroutine @@ -849,7 +843,7 @@ recursive subroutine newton_midpoint(si, f, x, atol, rtol, maxit, xlast, status) end if else call count_event(EVT_MIDPOINT_MAXIT) - if (accept_warning_maxiter(x, xlast, tolref, rtol)) & + if (accept_warning_maxiter(x)) & status = SYMPLECTIC_STEP_OK end if end subroutine @@ -974,7 +968,7 @@ recursive subroutine newton_rk_gauss(si, fs, s, x, atol, rtol, maxit, xlast, sta end if else call count_event(EVT_RK_GAUSS_MAXIT) - if (accept_warning_maxiter(x, xlast, tolref, rtol)) & + if (accept_warning_maxiter(x)) & status = SYMPLECTIC_STEP_OK end if end subroutine newton_rk_gauss @@ -1316,7 +1310,7 @@ recursive subroutine newton_rk_lobatto(si, fs, s, x, atol, rtol, maxit, xlast, s end if else call count_event(EVT_RK_LOBATTO_MAXIT) - if (accept_warning_maxiter(x, xlast, tolref, rtol)) & + if (accept_warning_maxiter(x)) & status = SYMPLECTIC_STEP_OK end if end subroutine newton_rk_lobatto diff --git a/src/orbit_symplectic_base.f90 b/src/orbit_symplectic_base.f90 index 6b4cb2ee..ac3db2cc 100644 --- a/src/orbit_symplectic_base.f90 +++ b/src/orbit_symplectic_base.f90 @@ -27,7 +27,6 @@ module orbit_symplectic_base integer, parameter :: SYMPLECTIC_STEP_EVENT_NOT_CONVERGED = 5 integer, parameter :: SYMPLECTIC_STEP_BOUNDARY_LIMITED = 6 logical :: symplectic_newton_warning_mode = .true. - real(dp), parameter :: symplectic_newton_warning_factor = 10.0_dp real(dp) :: boundary_event_fraction_tolerance = -1d0 real(dp) :: boundary_event_radial_tolerance = -1d0 diff --git a/src/params.f90 b/src/params.f90 index 1f46491a..7aba669b 100644 --- a/src/params.f90 +++ b/src/params.f90 @@ -205,8 +205,8 @@ subroutine read_config(config_file) call reset_seed_if_deterministic if (integmode > 0 .and. symplectic_newton_warning_mode) then - print *, 'WARNING: symplectic integrators accept bounded Newton ', & - 'corrections after maxit; failed steps retry from the last ', & + print *, 'WARNING: symplectic integrators accept finite Newton ', & + 'iterates after maxit; unusable steps retry from the last ', & 'valid state; see maxit/step-skip diagnostics' end if diff --git a/src/simple_gpu.f90 b/src/simple_gpu.f90 index 9cf6b52b..7d6e55a0 100644 --- a/src/simple_gpu.f90 +++ b/src/simple_gpu.f90 @@ -15,8 +15,7 @@ module simple_gpu use orbit_symplectic_base, only: symplectic_integrator_t, & SYMPLECTIC_STEP_OK, SYMPLECTIC_STEP_OUTSIDE_DOMAIN, & SYMPLECTIC_STEP_MAXITER, SYMPLECTIC_STEP_LINEAR_SOLVE, & - SYMPLECTIC_STEP_BOUNDARY_LIMITED, symplectic_newton_warning_mode, & - symplectic_newton_warning_factor + SYMPLECTIC_STEP_BOUNDARY_LIMITED, symplectic_newton_warning_mode use orbit_symplectic_euler1, only: sympl_euler1_newton_iter use orbit_symplectic_euler1, only: sympl_euler1_extrapolate_field, sympl_euler1_advance_angles use boozer_sub, only: boozer_state @@ -106,10 +105,7 @@ subroutine gpu_newton1(si, f, x, xlast, warning_mode, status) else status = SYMPLECTIC_STEP_BOUNDARY_LIMITED end if - else if (warning_mode .and. gpu_finite_iterate(x) .and. & - gpu_finite_iterate(xlast) .and. gpu_finite_iterate(tolref) .and. & - all(abs(x - xlast) <= & - symplectic_newton_warning_factor*si%rtol*abs(tolref))) then + else if (warning_mode .and. gpu_finite_iterate(x)) then status = SYMPLECTIC_STEP_OK end if ! Non-convergence diagnostics (CPU writes fort.6601) are omitted on device. diff --git a/test/tests/test_newton_solver_status.f90 b/test/tests/test_newton_solver_status.f90 index 854ef452..dea01fe6 100644 --- a/test/tests/test_newton_solver_status.f90 +++ b/test/tests/test_newton_solver_status.f90 @@ -263,7 +263,7 @@ subroutine test_newton_warning_mode MIDPOINT, GAUSS1, GAUSS2, GAUSS3, GAUSS4, LOBATTO3] type(symplectic_integrator_t) :: strict_integrator type(field_can_t) :: strict_field - real(dp) :: initial_state(4), accepted(2), previous(2), tolref(2) + real(dp) :: initial_state(4), accepted(2), previous(2) integer :: mode_index, step_status eval_field => evaluate_linear_radial @@ -287,24 +287,23 @@ subroutine test_newton_warning_mode end do previous = [1.0_dp, 2.0_dp] - tolref = [1.0_dp, 2.0_dp] accepted = previous + [5.0e-12_dp, 1.0e-11_dp] symplectic_newton_warning_mode = .true. - if (.not. accept_warning_maxiter(accepted, previous, tolref, 1.0e-12_dp)) then - error stop 'warning mode rejected a bounded Newton correction' + if (.not. accept_warning_maxiter(accepted)) then + error stop 'warning mode rejected a finite Newton iterate' end if accepted(1) = huge(1.0_dp) - if (accept_warning_maxiter(accepted, previous, tolref, 1.0e-12_dp)) then - error stop 'warning mode accepted an unbounded Newton correction' + if (.not. accept_warning_maxiter(accepted)) then + error stop 'warning mode rejected a large finite Newton iterate' end if accepted = previous accepted(1) = ieee_value(0.0_dp, ieee_quiet_nan) - if (accept_warning_maxiter(accepted, previous, tolref, 1.0e-12_dp)) then + if (accept_warning_maxiter(accepted)) then error stop 'warning mode accepted a non-finite Newton correction' end if accepted = previous symplectic_newton_warning_mode = .false. - if (accept_warning_maxiter(accepted, previous, tolref, 1.0e-12_dp)) then + if (accept_warning_maxiter(accepted)) then error stop 'strict mode accepted a Newton max-iteration state' end if end subroutine test_newton_warning_mode From 84d000d55e311b1d07d7259b004ad7bbc00112f8 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 08:58:26 +0200 Subject: [PATCH 2/7] Accept finite boundary-limited Newton iterates --- src/orbit_symplectic.f90 | 105 +++++++++++------------ test/tests/test_newton_solver_status.f90 | 26 +++++- 2 files changed, 73 insertions(+), 58 deletions(-) diff --git a/src/orbit_symplectic.f90 b/src/orbit_symplectic.f90 index fd437e96..c8f3f1ff 100644 --- a/src/orbit_symplectic.f90 +++ b/src/orbit_symplectic.f90 @@ -581,17 +581,16 @@ recursive subroutine newton1(si, f, x, maxit, xlast, status) return end if enddo - if (boundary_limited) then - if (step_boundary_limited .and. & - radial_boundary_reached(x, radial_indices, si%rtol)) then - status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN - else - status = SYMPLECTIC_STEP_BOUNDARY_LIMITED - end if - else - call count_event(EVT_NEWTON1_MAXIT) - if (accept_warning_maxiter(x)) & - status = SYMPLECTIC_STEP_OK + if (boundary_limited .and. step_boundary_limited .and. & + radial_boundary_reached(x, radial_indices, si%rtol)) then + status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN + return + end if + call count_event(EVT_NEWTON1_MAXIT) + if (accept_warning_maxiter(x)) then + status = SYMPLECTIC_STEP_OK + else if (boundary_limited) then + status = SYMPLECTIC_STEP_BOUNDARY_LIMITED end if end subroutine @@ -676,17 +675,16 @@ recursive subroutine newton2(si, f, x, atol, rtol, maxit, xlast, status) return end if enddo - if (boundary_limited) then - if (step_limited .and. & - radial_boundary_reached(x, radial_indices, rtol)) then - status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN - else - status = SYMPLECTIC_STEP_BOUNDARY_LIMITED - end if - else - call count_event(EVT_NEWTON2_MAXIT) - if (accept_warning_maxiter(x)) & - status = SYMPLECTIC_STEP_OK + if (boundary_limited .and. step_limited .and. & + radial_boundary_reached(x, radial_indices, rtol)) then + status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN + return + end if + call count_event(EVT_NEWTON2_MAXIT) + if (accept_warning_maxiter(x)) then + status = SYMPLECTIC_STEP_OK + else if (boundary_limited) then + status = SYMPLECTIC_STEP_BOUNDARY_LIMITED end if end subroutine @@ -834,17 +832,16 @@ recursive subroutine newton_midpoint(si, f, x, atol, rtol, maxit, xlast, status) return end if enddo - if (boundary_limited) then - if (step_limited .and. & - radial_boundary_reached(x, radial_indices, rtol)) then - status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN - else - status = SYMPLECTIC_STEP_BOUNDARY_LIMITED - end if - else - call count_event(EVT_MIDPOINT_MAXIT) - if (accept_warning_maxiter(x)) & - status = SYMPLECTIC_STEP_OK + if (boundary_limited .and. step_limited .and. & + radial_boundary_reached(x, radial_indices, rtol)) then + status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN + return + end if + call count_event(EVT_MIDPOINT_MAXIT) + if (accept_warning_maxiter(x)) then + status = SYMPLECTIC_STEP_OK + else if (boundary_limited) then + status = SYMPLECTIC_STEP_BOUNDARY_LIMITED end if end subroutine @@ -959,17 +956,16 @@ recursive subroutine newton_rk_gauss(si, fs, s, x, atol, rtol, maxit, xlast, sta return end if enddo - if (boundary_limited) then - if (step_limited .and. & - radial_boundary_reached(x, radial_indices, rtol)) then - status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN - else - status = SYMPLECTIC_STEP_BOUNDARY_LIMITED - end if - else - call count_event(EVT_RK_GAUSS_MAXIT) - if (accept_warning_maxiter(x)) & - status = SYMPLECTIC_STEP_OK + if (boundary_limited .and. step_limited .and. & + radial_boundary_reached(x, radial_indices, rtol)) then + status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN + return + end if + call count_event(EVT_RK_GAUSS_MAXIT) + if (accept_warning_maxiter(x)) then + status = SYMPLECTIC_STEP_OK + else if (boundary_limited) then + status = SYMPLECTIC_STEP_BOUNDARY_LIMITED end if end subroutine newton_rk_gauss @@ -1301,17 +1297,16 @@ recursive subroutine newton_rk_lobatto(si, fs, s, x, atol, rtol, maxit, xlast, s return end if enddo - if (boundary_limited) then - if (step_limited .and. & - radial_boundary_reached(x, radial_indices, rtol)) then - status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN - else - status = SYMPLECTIC_STEP_BOUNDARY_LIMITED - end if - else - call count_event(EVT_RK_LOBATTO_MAXIT) - if (accept_warning_maxiter(x)) & - status = SYMPLECTIC_STEP_OK + if (boundary_limited .and. step_limited .and. & + radial_boundary_reached(x, radial_indices, rtol)) then + status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN + return + end if + call count_event(EVT_RK_LOBATTO_MAXIT) + if (accept_warning_maxiter(x)) then + status = SYMPLECTIC_STEP_OK + else if (boundary_limited) then + status = SYMPLECTIC_STEP_BOUNDARY_LIMITED end if end subroutine newton_rk_lobatto diff --git a/test/tests/test_newton_solver_status.f90 b/test/tests/test_newton_solver_status.f90 index dea01fe6..9181f069 100644 --- a/test/tests/test_newton_solver_status.f90 +++ b/test/tests/test_newton_solver_status.f90 @@ -133,7 +133,8 @@ end subroutine retryable_boundary_step end module linear_radial_field_backend program test_newton_solver_status - use, intrinsic :: ieee_arithmetic, only: ieee_quiet_nan, ieee_value + use, intrinsic :: ieee_arithmetic, only: ieee_is_finite, ieee_quiet_nan, & + ieee_value use, intrinsic :: iso_fortran_env, only: dp => real64 use field_can_mod, only: field_can_t, eval_field => evaluate use linear_radial_field_backend, only: basin_limited_step, & @@ -261,8 +262,8 @@ program test_newton_solver_status subroutine test_newton_warning_mode integer, parameter :: modes(8) = [EXPL_IMPL_EULER, IMPL_EXPL_EULER, & MIDPOINT, GAUSS1, GAUSS2, GAUSS3, GAUSS4, LOBATTO3] - type(symplectic_integrator_t) :: strict_integrator - type(field_can_t) :: strict_field + type(symplectic_integrator_t) :: strict_integrator, warning_integrator + type(field_can_t) :: strict_field, warning_field real(dp) :: initial_state(4), accepted(2), previous(2) integer :: mode_index, step_status @@ -284,6 +285,25 @@ subroutine test_newton_warning_mode error stop 'strict max-iteration failure changed the accepted state' end if + warning_field%ro0 = 1.0_dp + warning_field%mu = 1.0_dp + call orbit_sympl_init(warning_integrator, warning_field, initial_state, & + 0.01_dp, 1, 0.0_dp, modes(mode_index)) + warning_integrator%atol = 0.0_dp + symplectic_newton_warning_mode = .true. + call orbit_timestep_sympl(warning_integrator, warning_field, step_status) + if (step_status /= SYMPLECTIC_STEP_OK) then + print *, 'warning mode, status:', modes(mode_index), step_status + error stop 'warning mode did not continue the timestep' + end if + if (.not. all(ieee_is_finite(warning_integrator%z))) then + error stop 'warning mode accepted a non-finite state' + end if + if (modes(mode_index) == EXPL_IMPL_EULER .and. & + all(warning_integrator%z == initial_state)) then + error stop 'warning mode did not commit the timestep' + end if + end do previous = [1.0_dp, 2.0_dp] From 8a08d2d66a17cfa20239bbfc96d86d3799adb5ba Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 09:27:35 +0200 Subject: [PATCH 3/7] Restore historical warning-mode Newton steps --- DOC/config.md | 14 ++-- src/orbit_symplectic.f90 | 127 +++++++++++++++++--------------- src/orbit_symplectic_euler1.f90 | 6 +- src/simple_gpu.f90 | 2 +- 4 files changed, 82 insertions(+), 67 deletions(-) diff --git a/DOC/config.md b/DOC/config.md index 3fca2184..884cec4c 100644 --- a/DOC/config.md +++ b/DOC/config.md @@ -13,12 +13,14 @@ * `symplectic_newton_warning_mode` defaults to `.true.`. When an implicit solve reaches its iteration limit, SIMPLE commits any finite final Newton iterate - and continues the orbit. Each occurrence is reported by the corresponding - `*_maxit` diagnostic. Recursive recovery is reserved for an unusable step, - such as a failed linear solve or non-finite iterate. The production RK, - symplectic, and full-orbit paths use the same terminal convention: if that - recovery cannot resolve a numerical microstep, SIMPLE retains any contiguous - accepted prefix and holds only the unresolved remainder, records + and continues the orbit. This default path uses the historical undamped + Newton correction; the radial solver-basin limiter remains active in strict + mode. Each occurrence is reported by the corresponding `*_maxit` diagnostic. + Recursive recovery is reserved for an unusable step, such as a failed linear + solve or non-finite iterate. The production RK, symplectic, and full-orbit + paths use the same terminal convention: if that recovery cannot resolve a + numerical microstep, SIMPLE retains any contiguous accepted prefix and holds + only the unresolved remainder, records `warning_step_skip`, advances the clock for the complete interval, and retries from that valid state on the next microstep. A warning hold does not terminate or numerically disqualify the marker: a marker that reaches the diff --git a/src/orbit_symplectic.f90 b/src/orbit_symplectic.f90 index c8f3f1ff..54e623fc 100644 --- a/src/orbit_symplectic.f90 +++ b/src/orbit_symplectic.f90 @@ -565,7 +565,8 @@ recursive subroutine newton1(si, f, x, maxit, xlast, status) call eval_field(f, x(1), si%z(2), si%z(3), 2) call get_derivatives2(f, x(2)) call sympl_euler1_newton_iter(si, f, x, tolref, xlast, converged, & - linear_failed, step_boundary_limited) + linear_failed, step_boundary_limited, & + .not. symplectic_newton_warning_mode) boundary_limited = boundary_limited .or. step_boundary_limited if (linear_failed) then @@ -581,16 +582,17 @@ recursive subroutine newton1(si, f, x, maxit, xlast, status) return end if enddo - if (boundary_limited .and. step_boundary_limited .and. & - radial_boundary_reached(x, radial_indices, si%rtol)) then - status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN - return - end if - call count_event(EVT_NEWTON1_MAXIT) - if (accept_warning_maxiter(x)) then - status = SYMPLECTIC_STEP_OK - else if (boundary_limited) then - status = SYMPLECTIC_STEP_BOUNDARY_LIMITED + if (boundary_limited) then + if (step_boundary_limited .and. & + radial_boundary_reached(x, radial_indices, si%rtol)) then + status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN + else + status = SYMPLECTIC_STEP_BOUNDARY_LIMITED + end if + else + call count_event(EVT_NEWTON1_MAXIT) + if (accept_warning_maxiter(x)) & + status = SYMPLECTIC_STEP_OK end if end subroutine @@ -651,10 +653,15 @@ recursive subroutine newton2(si, f, x, atol, rtol, maxit, xlast, status) jinv(3,3) = fjac(1,1)*fjac(2,2) - fjac(1,2)*fjac(2,1) correction = matmul(jinv, fvec)/det - call limit_radial_newton_step(x, correction, radial_indices, step_scale, & - step_limited) - boundary_limited = boundary_limited .or. step_limited - x = x - step_scale*correction + if (symplectic_newton_warning_mode) then + step_limited = .false. + x = x - correction + else + call limit_radial_newton_step(x, correction, radial_indices, step_scale, & + step_limited) + boundary_limited = boundary_limited .or. step_limited + x = x - step_scale*correction + end if !call dgesv(n, 1, fjac, n, pivot, fvec, n, info) ! after solution: fvec = (xold-xnew)_Newton @@ -675,16 +682,17 @@ recursive subroutine newton2(si, f, x, atol, rtol, maxit, xlast, status) return end if enddo - if (boundary_limited .and. step_limited .and. & - radial_boundary_reached(x, radial_indices, rtol)) then - status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN - return - end if - call count_event(EVT_NEWTON2_MAXIT) - if (accept_warning_maxiter(x)) then - status = SYMPLECTIC_STEP_OK - else if (boundary_limited) then - status = SYMPLECTIC_STEP_BOUNDARY_LIMITED + if (boundary_limited) then + if (step_limited .and. & + radial_boundary_reached(x, radial_indices, rtol)) then + status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN + else + status = SYMPLECTIC_STEP_BOUNDARY_LIMITED + end if + else + call count_event(EVT_NEWTON2_MAXIT) + if (accept_warning_maxiter(x)) & + status = SYMPLECTIC_STEP_OK end if end subroutine @@ -805,7 +813,7 @@ recursive subroutine newton_midpoint(si, f, x, atol, rtol, maxit, xlast, status) call solve_newton_system(fjac, fvec, status) if (status /= SYMPLECTIC_STEP_OK) return status = SYMPLECTIC_STEP_MAXITER - if (sympl_rmax > 1d0) then + if (sympl_rmax > 1d0 .or. symplectic_newton_warning_mode) then step_limited = .false. x = x - fvec else @@ -832,16 +840,17 @@ recursive subroutine newton_midpoint(si, f, x, atol, rtol, maxit, xlast, status) return end if enddo - if (boundary_limited .and. step_limited .and. & - radial_boundary_reached(x, radial_indices, rtol)) then - status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN - return - end if - call count_event(EVT_MIDPOINT_MAXIT) - if (accept_warning_maxiter(x)) then - status = SYMPLECTIC_STEP_OK - else if (boundary_limited) then - status = SYMPLECTIC_STEP_BOUNDARY_LIMITED + if (boundary_limited) then + if (step_limited .and. & + radial_boundary_reached(x, radial_indices, rtol)) then + status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN + else + status = SYMPLECTIC_STEP_BOUNDARY_LIMITED + end if + else + call count_event(EVT_MIDPOINT_MAXIT) + if (accept_warning_maxiter(x)) & + status = SYMPLECTIC_STEP_OK end if end subroutine @@ -923,7 +932,7 @@ recursive subroutine newton_rk_gauss(si, fs, s, x, atol, rtol, maxit, xlast, sta if (status /= SYMPLECTIC_STEP_OK) return status = SYMPLECTIC_STEP_MAXITER ! after solution: fvec = (xold-xnew)_Newton - if (sympl_rmax > 1d0) then + if (sympl_rmax > 1d0 .or. symplectic_newton_warning_mode) then step_limited = .false. x = x - fvec else @@ -956,16 +965,17 @@ recursive subroutine newton_rk_gauss(si, fs, s, x, atol, rtol, maxit, xlast, sta return end if enddo - if (boundary_limited .and. step_limited .and. & - radial_boundary_reached(x, radial_indices, rtol)) then - status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN - return - end if - call count_event(EVT_RK_GAUSS_MAXIT) - if (accept_warning_maxiter(x)) then - status = SYMPLECTIC_STEP_OK - else if (boundary_limited) then - status = SYMPLECTIC_STEP_BOUNDARY_LIMITED + if (boundary_limited) then + if (step_limited .and. & + radial_boundary_reached(x, radial_indices, rtol)) then + status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN + else + status = SYMPLECTIC_STEP_BOUNDARY_LIMITED + end if + else + call count_event(EVT_RK_GAUSS_MAXIT) + if (accept_warning_maxiter(x)) & + status = SYMPLECTIC_STEP_OK end if end subroutine newton_rk_gauss @@ -1264,7 +1274,7 @@ recursive subroutine newton_rk_lobatto(si, fs, s, x, atol, rtol, maxit, xlast, s if (status /= SYMPLECTIC_STEP_OK) return status = SYMPLECTIC_STEP_MAXITER ! after solution: fvec = (xold-xnew)_Newton - if (sympl_rmax > 1d0) then + if (sympl_rmax > 1d0 .or. symplectic_newton_warning_mode) then step_limited = .false. x = x - fvec else @@ -1297,16 +1307,17 @@ recursive subroutine newton_rk_lobatto(si, fs, s, x, atol, rtol, maxit, xlast, s return end if enddo - if (boundary_limited .and. step_limited .and. & - radial_boundary_reached(x, radial_indices, rtol)) then - status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN - return - end if - call count_event(EVT_RK_LOBATTO_MAXIT) - if (accept_warning_maxiter(x)) then - status = SYMPLECTIC_STEP_OK - else if (boundary_limited) then - status = SYMPLECTIC_STEP_BOUNDARY_LIMITED + if (boundary_limited) then + if (step_limited .and. & + radial_boundary_reached(x, radial_indices, rtol)) then + status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN + else + status = SYMPLECTIC_STEP_BOUNDARY_LIMITED + end if + else + call count_event(EVT_RK_LOBATTO_MAXIT) + if (accept_warning_maxiter(x)) & + status = SYMPLECTIC_STEP_OK end if end subroutine newton_rk_lobatto diff --git a/src/orbit_symplectic_euler1.f90 b/src/orbit_symplectic_euler1.f90 index 4925d524..098240ae 100644 --- a/src/orbit_symplectic_euler1.f90 +++ b/src/orbit_symplectic_euler1.f90 @@ -48,7 +48,7 @@ subroutine sympl_euler1_jacobian(si, f, x, jac) end subroutine sympl_euler1_jacobian subroutine sympl_euler1_newton_iter(si, f, x, tolref, xlast, converged, & - linear_failed, boundary_limited) + linear_failed, boundary_limited, limit_radial_step) !$acc routine seq type(symplectic_integrator_t), intent(in) :: si type(field_can_t), intent(in) :: f @@ -58,6 +58,7 @@ subroutine sympl_euler1_newton_iter(si, f, x, tolref, xlast, converged, & logical, intent(out) :: converged logical, intent(out) :: linear_failed logical, intent(out) :: boundary_limited + logical, intent(in) :: limit_radial_step real(dp) :: fvec(2), fjac(2, 2), ijac(2, 2) real(dp) :: correction(2), determinant, matrix_scale, step_scale @@ -94,7 +95,8 @@ subroutine sympl_euler1_newton_iter(si, f, x, tolref, xlast, converged, & correction(1) = ijac(1, 1)*fvec(1) + ijac(1, 2)*fvec(2) correction(2) = ijac(2, 1)*fvec(1) + ijac(2, 2)*fvec(2) step_scale = 1d0 - if (sympl_rmax <= 1d0 .and. correction(1) < 0d0) then + if (limit_radial_step .and. sympl_rmax <= 1d0 .and. & + correction(1) < 0d0) then step_scale = min(1d0, & 0.8d0*max(0d0, sympl_rmax - x(1))/(-correction(1))) end if diff --git a/src/simple_gpu.f90 b/src/simple_gpu.f90 index 7d6e55a0..c11cb401 100644 --- a/src/simple_gpu.f90 +++ b/src/simple_gpu.f90 @@ -82,7 +82,7 @@ subroutine gpu_newton1(si, f, x, xlast, warning_mode, status) call eval_field_booz(f, x(1), si%z(2), si%z(3), 2) call get_derivatives2(f, x(2)) call sympl_euler1_newton_iter(si, f, x, tolref, xlast, converged, & - linear_failed, step_boundary_limited) + linear_failed, step_boundary_limited, .not. warning_mode) boundary_limited = boundary_limited .or. step_boundary_limited if (linear_failed) then From cbe3c5fb8d13155112e03db286591def47b240cd Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 13:23:27 +0200 Subject: [PATCH 4/7] Bound Newton warnings only for SPECTRE --- src/orbit_symplectic.f90 | 46 +++++++++++++------ src/orbit_symplectic_base.f90 | 5 ++ src/spectre_sympl_orbit.f90 | 5 +- .../field_can/test_spectre_sympl_crossing.f90 | 3 ++ test/tests/test_newton_solver_status.f90 | 21 +++++++-- 5 files changed, 62 insertions(+), 18 deletions(-) diff --git a/src/orbit_symplectic.f90 b/src/orbit_symplectic.f90 index 54e623fc..7830a640 100644 --- a/src/orbit_symplectic.f90 +++ b/src/orbit_symplectic.f90 @@ -12,7 +12,8 @@ module orbit_symplectic SYMPLECTIC_STEP_MAXITER, SYMPLECTIC_STEP_LINEAR_SOLVE, & SYMPLECTIC_STEP_BOUNDARY, SYMPLECTIC_STEP_EVENT_NOT_CONVERGED, & SYMPLECTIC_STEP_BOUNDARY_LIMITED, boundary_event_fraction_tolerance, & - boundary_event_radial_tolerance, symplectic_newton_warning_mode + boundary_event_radial_tolerance, symplectic_newton_warning_mode, & + symplectic_newton_warning_factor use orbit_symplectic_quasi, only: orbit_timestep_quasi, timestep_expl_impl_euler_quasi, & timestep_impl_expl_euler_quasi, timestep_midpoint_quasi, orbit_timestep_rk45, & timestep_rk_gauss_quasi, timestep_rk_lobatto_quasi @@ -49,15 +50,27 @@ pure logical function finite_newton_iterate(x) finite_newton_iterate = .true. end function finite_newton_iterate -logical function accept_warning_maxiter(x) +logical function accept_warning_maxiter(x, xlast, tolref, rtol, & + accept_unbounded) ! Warning mode preserves the historical SIMPLE continuation contract: a ! Newton iteration-limit event is diagnostic, not a rejected step. Commit ! every finite final iterate and leave recursive recovery for genuinely ! unusable states such as failed linear solves or non-finite iterates. - real(dp), intent(in) :: x(:) - - accept_warning_maxiter = symplectic_newton_warning_mode .and. & - finite_newton_iterate(x) + real(dp), intent(in) :: x(:), xlast(:), tolref(:), rtol + logical, intent(in) :: accept_unbounded + + accept_warning_maxiter = .false. + if (.not. symplectic_newton_warning_mode) return + if (.not. finite_newton_iterate(x)) return + if (accept_unbounded) then + accept_warning_maxiter = .true. + return + end if + if (rtol <= 0.0_dp) return + if (.not. finite_newton_iterate(xlast)) return + if (.not. finite_newton_iterate(tolref)) return + accept_warning_maxiter = all(abs(x - xlast) <= & + symplectic_newton_warning_factor*rtol*abs(tolref)) end function accept_warning_maxiter subroutine advance_symplectic_with_retry(si, f, stepper, status, & @@ -566,7 +579,8 @@ recursive subroutine newton1(si, f, x, maxit, xlast, status) call get_derivatives2(f, x(2)) call sympl_euler1_newton_iter(si, f, x, tolref, xlast, converged, & linear_failed, step_boundary_limited, & - .not. symplectic_newton_warning_mode) + .not. symplectic_newton_warning_mode .or. & + .not. si%accept_unbounded_newton_warning) boundary_limited = boundary_limited .or. step_boundary_limited if (linear_failed) then @@ -591,7 +605,8 @@ recursive subroutine newton1(si, f, x, maxit, xlast, status) end if else call count_event(EVT_NEWTON1_MAXIT) - if (accept_warning_maxiter(x)) & + if (accept_warning_maxiter(x, xlast, tolref, si%rtol, & + si%accept_unbounded_newton_warning)) & status = SYMPLECTIC_STEP_OK end if end subroutine @@ -653,7 +668,8 @@ recursive subroutine newton2(si, f, x, atol, rtol, maxit, xlast, status) jinv(3,3) = fjac(1,1)*fjac(2,2) - fjac(1,2)*fjac(2,1) correction = matmul(jinv, fvec)/det - if (symplectic_newton_warning_mode) then + if (symplectic_newton_warning_mode .and. & + si%accept_unbounded_newton_warning) then step_limited = .false. x = x - correction else @@ -691,7 +707,8 @@ recursive subroutine newton2(si, f, x, atol, rtol, maxit, xlast, status) end if else call count_event(EVT_NEWTON2_MAXIT) - if (accept_warning_maxiter(x)) & + if (accept_warning_maxiter(x, xlast, tolref, rtol, & + si%accept_unbounded_newton_warning)) & status = SYMPLECTIC_STEP_OK end if end subroutine @@ -849,7 +866,8 @@ recursive subroutine newton_midpoint(si, f, x, atol, rtol, maxit, xlast, status) end if else call count_event(EVT_MIDPOINT_MAXIT) - if (accept_warning_maxiter(x)) & + if (accept_warning_maxiter(x, xlast, tolref, rtol, & + si%accept_unbounded_newton_warning)) & status = SYMPLECTIC_STEP_OK end if end subroutine @@ -974,7 +992,8 @@ recursive subroutine newton_rk_gauss(si, fs, s, x, atol, rtol, maxit, xlast, sta end if else call count_event(EVT_RK_GAUSS_MAXIT) - if (accept_warning_maxiter(x)) & + if (accept_warning_maxiter(x, xlast, tolref, rtol, & + si%accept_unbounded_newton_warning)) & status = SYMPLECTIC_STEP_OK end if end subroutine newton_rk_gauss @@ -1316,7 +1335,8 @@ recursive subroutine newton_rk_lobatto(si, fs, s, x, atol, rtol, maxit, xlast, s end if else call count_event(EVT_RK_LOBATTO_MAXIT) - if (accept_warning_maxiter(x)) & + if (accept_warning_maxiter(x, xlast, tolref, rtol, & + si%accept_unbounded_newton_warning)) & status = SYMPLECTIC_STEP_OK end if end subroutine newton_rk_lobatto diff --git a/src/orbit_symplectic_base.f90 b/src/orbit_symplectic_base.f90 index ac3db2cc..ee375110 100644 --- a/src/orbit_symplectic_base.f90 +++ b/src/orbit_symplectic_base.f90 @@ -27,6 +27,7 @@ module orbit_symplectic_base integer, parameter :: SYMPLECTIC_STEP_EVENT_NOT_CONVERGED = 5 integer, parameter :: SYMPLECTIC_STEP_BOUNDARY_LIMITED = 6 logical :: symplectic_newton_warning_mode = .true. + real(dp), parameter :: symplectic_newton_warning_factor = 10.0_dp real(dp) :: boundary_event_fraction_tolerance = -1d0 real(dp) :: boundary_event_radial_tolerance = -1d0 @@ -45,6 +46,10 @@ module orbit_symplectic_base real(dp) :: last_step_fraction = 1d0 real(dp) :: last_event_radial_residual = 0d0 real(dp) :: last_event_fraction_width = 0d0 + ! Generic SIMPLE tracing historically commits every finite Newton + ! iteration-limit state in warning mode. Specialized callers may + ! disable that liberal policy and accept only a near-converged state. + logical :: accept_unbounded_newton_warning = .true. end type symplectic_integrator_t !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc diff --git a/src/spectre_sympl_orbit.f90 b/src/spectre_sympl_orbit.f90 index aed17f93..7eae0df8 100644 --- a/src/spectre_sympl_orbit.f90 +++ b/src/spectre_sympl_orbit.f90 @@ -138,7 +138,7 @@ module spectre_sympl_orbit subroutine sympl_spectre_reset(state, si, mvol, mode, level) type(sympl_spectre_state_t), intent(out) :: state - type(symplectic_integrator_t), intent(in) :: si + type(symplectic_integrator_t), intent(inout) :: si integer, intent(in) :: mvol, mode, level if (mode <= 0) then @@ -156,6 +156,9 @@ subroutine sympl_spectre_reset(state, si, mvol, mode, level) state%mode = mode state%level = level state%dt_std = si%dt + ! SPECTRE discontinuities require a valid interface state. Preserve + ! near-converged warning-mode iterates, but recover gross corrections. + si%accept_unbounded_newton_warning = .false. call set_home(state, si%z(1)) end subroutine sympl_spectre_reset diff --git a/test/tests/field_can/test_spectre_sympl_crossing.f90 b/test/tests/field_can/test_spectre_sympl_crossing.f90 index cfda19a5..67271fb2 100644 --- a/test/tests/field_can/test_spectre_sympl_crossing.f90 +++ b/test/tests/field_can/test_spectre_sympl_crossing.f90 @@ -289,6 +289,9 @@ subroutine trace_marker(im, nrec_got, h_series, h_actual, p_series, t_series, & call init_sympl(si, f, z, dtaumin, dtaumin, relerr, integmode) call sympl_spectre_reset(state, si, spectre_mvol, integmode, & crossing_level) + if (si%accept_unbounded_newton_warning) then + error stop 'SPECTRE did not enable bounded Newton warnings' + end if nrec_got = 0 do k = 1, NSTEP diff --git a/test/tests/test_newton_solver_status.f90 b/test/tests/test_newton_solver_status.f90 index 9181f069..50fa5510 100644 --- a/test/tests/test_newton_solver_status.f90 +++ b/test/tests/test_newton_solver_status.f90 @@ -309,21 +309,34 @@ subroutine test_newton_warning_mode previous = [1.0_dp, 2.0_dp] accepted = previous + [5.0e-12_dp, 1.0e-11_dp] symplectic_newton_warning_mode = .true. - if (.not. accept_warning_maxiter(accepted)) then + if (.not. accept_warning_maxiter(accepted, previous, previous, & + 1.0e-12_dp, .true.)) then error stop 'warning mode rejected a finite Newton iterate' end if accepted(1) = huge(1.0_dp) - if (.not. accept_warning_maxiter(accepted)) then + if (.not. accept_warning_maxiter(accepted, previous, previous, & + 1.0e-12_dp, .true.)) then error stop 'warning mode rejected a large finite Newton iterate' end if + if (accept_warning_maxiter(accepted, previous, previous, & + 1.0e-12_dp, .false.)) then + error stop 'bounded warning policy accepted a gross Newton correction' + end if + accepted = previous + [5.0e-12_dp, 1.0e-11_dp] + if (.not. accept_warning_maxiter(accepted, previous, previous, & + 1.0e-12_dp, .false.)) then + error stop 'bounded warning policy rejected a near-converged iterate' + end if accepted = previous accepted(1) = ieee_value(0.0_dp, ieee_quiet_nan) - if (accept_warning_maxiter(accepted)) then + if (accept_warning_maxiter(accepted, previous, previous, & + 1.0e-12_dp, .true.)) then error stop 'warning mode accepted a non-finite Newton correction' end if accepted = previous symplectic_newton_warning_mode = .false. - if (accept_warning_maxiter(accepted)) then + if (accept_warning_maxiter(accepted, previous, previous, & + 1.0e-12_dp, .true.)) then error stop 'strict mode accepted a Newton max-iteration state' end if end subroutine test_newton_warning_mode From e6d764471b2cd86adf85b547943b01674288af9b Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 14:34:38 +0200 Subject: [PATCH 5/7] Add canonical orbit frequency APIs (#489) ## Risk tier - [ ] T0: docs, comments, small build metadata - [ ] T1: pure refactor, no behavior change intended - [ ] T2: local numerical logic - [x] T3: physics, output behavior, coordinate convention - [ ] T4: parallelism, GPU, dependency, CI, or security-sensitive build logic ## Correctness contract ### Intended behavior change Add a side-effect-free Fortran API and a Python API that compute canonical bounce/transit and mean toroidal angular frequencies for trapped and passing guiding-centre orbits. Support one or multiple completed poloidal cycles, return cycle spreads and explicit status metadata, and expose interpolated tip/toroidal cut tracing through Python. ### Behavior that must not change Existing orbit tracing, classification, NetCDF output, and golden-record results are unchanged. The core frequency routine copies the caller's tracer and does no file I/O. ### Coordinate / unit conventions The source-level Fortran routine accepts SIMPLE integration coordinates. The flat/Python boundary accepts `[s, theta, phi, v/v0, lambda]` in public reference coordinates and performs `ref_to_integ`. Angles and displacements are radians, periods are seconds, and frequencies are angular frequencies in rad/s. `omega_b` is positive; `omega_phi` is signed. ### Numerical invariants Trapped periods use successive negative-to-positive parallel-velocity crossings. Passing periods use successive signed `2*pi` advances of unwrapped poloidal angle. Event time and toroidal displacement are interpolated within the crossing step. Multi-cycle results report means and sample standard deviations. ## Tests added - unit: Fortran invalid-option status and Python argument validation - integration: native trapped/passing, three-cycle frequency computation - system: Python trapped multi-cycle frequency and interpolated tip-cut APIs - golden record: unchanged ## Golden-record impact - [x] unchanged - [ ] changed ## Failure modes considered Invalid options, orbit loss, integrator errors, and maximum-step exhaustion return distinct statuses. Unknown orbit class remains explicit. The cut wrapper rejects unsupported cut types and non-symplectic integration. ## Manual validation Validated trapped and passing VMEC seeds through the native API; validated the Python frequency and cut interfaces against the compiled f90wrap backend. ## Verification - `make test TEST=^test_orbit_frequencies$ CONFIG=Fast VERBOSE=0` - `make test TEST=^test_simple_api$ CONFIG=Fast VERBOSE=0` (24 API cases pass) - Full non-regression suite: 75/76 targets pass after installing the optional plotting dependency. The remaining `test_orbit_macro` failure is pre-existing: it accesses `pysimple._backend`, which current `main` does not export. - `git diff --check` --- DOC/coordinates-and-fields.md | 24 ++- python/CMakeLists.txt | 1 + python/README.md | 45 ++++++ python/pysimple/__init__.py | 214 +++++++++++++++++++++++++- src/CMakeLists.txt | 5 + src/orbit_frequencies.f90 | 200 ++++++++++++++++++++++++ src/simple_main.f90 | 96 ++++++++++++ test/python/CMakeLists.txt | 24 +++ test/python/test_simple_api.py | 89 +++++++++++ test/tests/CMakeLists.txt | 7 + test/tests/test_orbit_frequencies.f90 | 60 ++++++++ 11 files changed, 756 insertions(+), 9 deletions(-) create mode 100644 src/orbit_frequencies.f90 create mode 100644 test/tests/test_orbit_frequencies.f90 diff --git a/DOC/coordinates-and-fields.md b/DOC/coordinates-and-fields.md index 06f58d1a..e4a49b22 100644 --- a/DOC/coordinates-and-fields.md +++ b/DOC/coordinates-and-fields.md @@ -1200,7 +1200,29 @@ because they have no reference-coordinate geometry. | `src/magfie.f90` | Unified field evaluation interface | | `src/magfie_can_boozer.f90` | Boozer/Canflux implementations | -### 10.5 Integration +### 10.5 Canonical Frequency API + +`src/orbit_frequencies.f90` provides the source-level Fortran interface +`compute_canonical_frequencies(tracer, initial_state, options, result)`. +`initial_state` uses SIMPLE integration coordinates in the standard ordering +`[s, theta, phi, v/v0, lambda]`; all angles are radians and remain unwrapped. +The caller must initialize `tracer_t` and set `dtaumin`, `v0`, `relerr`, and +`integmode`. The routine copies the tracer before integration and therefore +does not advance caller-owned orbit state. + +The result contains the positive bounce/transit frequency +`omega_b = 2*pi/period` and the signed mean toroidal frequency +`omega_phi = delta_phi/period`. Periods are seconds, frequencies are rad/s, +and toroidal displacement is radians. Named status and orbit-class constants +are exported by the module. A downstream Fortran target can link the CMake +target `simple` and `use orbit_frequencies`; the target publishes SIMPLE's +build-tree module directory. + +The flat wrapper in `simple_main` accepts the public reference-coordinate +state, performs `ref_to_integ`, and is the boundary used by f90wrap and +`pysimple.compute_canonical_frequencies`. + +### 10.6 Integration | File | Purpose | |------|---------| diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index dc3e29f8..5245a404 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -73,6 +73,7 @@ message(STATUS "Python binary output dir: ${CMAKE_CURRENT_BINARY_DIR}") set(FILES_TO_WRAP simple_main.f90 simple.f90 + orbit_frequencies.f90 samplers.f90 cut_detector.f90 callback.f90 diff --git a/python/README.md b/python/README.md index 183794cb..35d3462d 100644 --- a/python/README.md +++ b/python/README.md @@ -52,6 +52,51 @@ First verify the bindings import cleanly: Complete examples are in ``examples/simple_api.py``, ``examples/classify_fast.py``, and ``examples/classify_fractal.py``. +Canonical orbit frequencies +--------------------------- + +``compute_canonical_frequencies`` traces one particle until it completes the +requested number of poloidal cycles: + +.. code-block:: python + + position = [0.6, 0.0, 0.25 * np.pi, 1.0, 0.3] + frequency = pysimple.compute_canonical_frequencies( + position, integrator="midpoint", n_periods=4 + ) + print(frequency["omega_b"], frequency["omega_phi"]) + +For trapped particles, a cycle is measured between same-direction +``v_parallel = 0`` crossings. For passing particles, it is measured between +successive signed ``2*pi`` advances of the unwrapped poloidal angle. The +toroidal displacement is also unwrapped. Periods are seconds, displacements +are radians, and both frequencies are angular frequencies in rad/s. + +Use ``n_periods=1`` for the inexpensive single-cycle value in an axisymmetric +field. Larger values return the mean and sample spread over several cycles, +which is useful in asymmetric fields. Always inspect ``status`` before using a +result; losses, integrator errors, and the maximum-step limit remain visible +instead of being silently discarded. An orbit that leaves the plasma through +the outer radial boundary reports ``FREQ_ORBIT_LOST``. + +The particle species and energy follow the same defaults as ``simple.in`` +(3.5 MeV alphas). A frequency computed for a different particle needs the +matching overrides in ``init``, for example a 5 keV deuteron: + +.. code-block:: python + + pysimple.init("field.nc", deterministic=True, n_e=1, n_d=2, facE_al=700.0) + +A fast particle in a small equilibrium can be genuinely unconfined; a +``FREQ_ORBIT_LOST`` result for a seed that should be confined usually means +the species/energy overrides are missing. + +The same event machinery is available as ``pysimple.trace_to_cut``. Select +``cut="tip"`` for a same-direction ``v_parallel=0`` section or +``cut="toroidal"`` for a field-period section. The returned state is converted +back to public reference coordinates. This driver currently requires one of +SIMPLE's symplectic integrators. + Legacy script note ------------------ diff --git a/python/pysimple/__init__.py b/python/pysimple/__init__.py index 9987f934..1a9374b1 100644 --- a/python/pysimple/__init__.py +++ b/python/pysimple/__init__.py @@ -63,6 +63,19 @@ "lobatto3": LOBATTO3, } +ORBIT_CLASS_UNKNOWN = 0 +ORBIT_CLASS_TRAPPED = 1 +ORBIT_CLASS_PASSING = 2 + +FREQ_SUCCESS = 0 +FREQ_INVALID_INPUT = 1 +FREQ_ORBIT_LOST = 2 +FREQ_INTEGRATOR_ERROR = 3 +FREQ_MAX_STEPS = 4 + +CUT_TIP = 0 +CUT_TOROIDAL = 1 + # Default spline orders used by Fortran DEFAULT_NS_S = 5 DEFAULT_NS_TP = 5 @@ -71,6 +84,7 @@ # Field type constants (mirroring magfie_sub.f90) _TEST_FIELD_ID = -1 _VMEC_FIELD_ID = 1 +_BOOZER_FIELD_ID = 2 # Direct access to Fortran params module params = _fortran_backend.params @@ -81,6 +95,7 @@ # Module state _initialized = False _current_vmec: str | None = None +_current_chartmap = False _tracer: "_fortran_backend.simple.tracer_t | None" = None _field_key: "tuple | None" = None @@ -89,6 +104,22 @@ def _is_test_field() -> bool: return int(_fortran_backend.velo_mod.isw_field_type) == _TEST_FIELD_ID +def _is_boozer_chartmap(path: str) -> bool: + """Return whether *path* declares the Boozer chart-map file contract.""" + try: + import netCDF4 + except ImportError as exc: + raise ImportError( + "netCDF4 is required to initialize SIMPLE from a NetCDF equilibrium" + ) from exc + + try: + with netCDF4.Dataset(path, "r") as dataset: + return int(dataset.getncattr("boozer_field")) == 1 + except (AttributeError, OSError, TypeError, ValueError): + return False + + def _sampling_rng() -> np.random.Generator: deterministic = bool(getattr(params, "deterministic", False)) seed = 0 if deterministic else None @@ -157,12 +188,23 @@ def init( >>> import pysimple >>> pysimple.init('wout.nc', deterministic=True, ntestpart=1000, trace_time=1e-3) """ - global _initialized, _current_vmec, _tracer, _trace_initialized, _field_key + global _initialized, _current_vmec, _current_chartmap, _tracer + global _trace_initialized, _field_key # Reset trace initialization flag to ensure field pointers are updated _trace_initialized = False vmec_path = str(Path(vmec_file).expanduser().resolve()) + chartmap_mode = _is_boozer_chartmap(vmec_path) + + # A Boozer chart map is already expressed in the canonical coordinates + # required by the symplectic integrators. Select that representation when + # the caller has not explicitly requested a field type. + explicit_field_type = {"isw_field_type", "integ_coords"} & param_overrides.keys() + if chartmap_mode and not explicit_field_type: + _fortran_backend.velo_mod.isw_field_type = _BOOZER_FIELD_ID + if hasattr(params, "integ_coords"): + params.integ_coords = _BOOZER_FIELD_ID # Step 1: Set parameters (replaces read_config without file I/O) params.netcdffile = vmec_path @@ -183,10 +225,14 @@ def init( _fortran_backend.velo_mod.isw_field_type = value_int if hasattr(params, "integ_coords"): params.integ_coords = value_int - elif not hasattr(params, key): + continue + # Fortran namelist names are case-insensitive (e.g. facE_al in + # simple.in), while f90wrap exposes them lowercased. + if not hasattr(params, key): + key = key.lower() + if not hasattr(params, key): raise ValueError(f"Unknown SIMPLE parameter: {key}") - else: - setattr(params, key, value) + setattr(params, key, value) if hasattr(params, "apply_config_aliases"): params.apply_config_aliases() @@ -212,6 +258,18 @@ def init( ) _field_key = field_key + field_type = int(_fortran_backend.velo_mod.isw_field_type) + + # Chart-map inputs do not carry a VMEC wout. Point magfie at the Boozer + # chart-map field before params_init(), whose stevvo setup evaluates B. + if chartmap_mode: + if field_type != _BOOZER_FIELD_ID: + raise ValueError( + "Boozer chart-map input requires isw_field_type=2 " + "(or omit isw_field_type to select it automatically)" + ) + _fortran_backend.magfie_wrapper.wrapper_init_magfie(field_type) + # Step 3: params_init (same as Fortran main()) # This calls reset_seed_if_deterministic() internally! # Also calls reallocate_arrays() which allocates xstart, volstart needed by init_starting_surf @@ -220,12 +278,13 @@ def init( # Step 4: init_magfie - set function pointer for magnetic field evaluation # Use isw_field_type from velo_mod (set via param_overrides above) - field_type = int(_fortran_backend.velo_mod.isw_field_type) - # Step 5: init_starting_surf (required for non-TEST fields) # Match Fortran driver ordering: initialize VMEC magfie for surface setup, # then restore requested field type for tracing. - if field_type != _TEST_FIELD_ID: + if chartmap_mode: + samplers = _fortran_backend.Samplers() + samplers.init_starting_surf() + elif field_type != _TEST_FIELD_ID: _fortran_backend.magfie_wrapper.wrapper_init_magfie(_VMEC_FIELD_ID) samplers = _fortran_backend.Samplers() samplers.init_starting_surf() @@ -235,6 +294,7 @@ def init( _initialized = True _current_vmec = vmec_path + _current_chartmap = chartmap_mode def sample_surface(n_particles: int, s: float) -> np.ndarray: @@ -424,7 +484,8 @@ def _sync_single_surface_state(surface_s: float) -> None: _fortran_backend.params_wrapper.set_sbeg(1, float(surface_s)) field_type = int(_fortran_backend.velo_mod.isw_field_type) - _fortran_backend.magfie_wrapper.wrapper_init_magfie(_VMEC_FIELD_ID) + setup_field_type = field_type if _current_chartmap else _VMEC_FIELD_ID + _fortran_backend.magfie_wrapper.wrapper_init_magfie(setup_field_type) _fortran_backend.Samplers().init_starting_surf() _fortran_backend.magfie_wrapper.wrapper_init_magfie(field_type) @@ -571,6 +632,143 @@ def trace_orbit( _restore_test_field_bounds(test_bounds) +def compute_canonical_frequencies( + position: np.ndarray, + *, + integrator: str | int = MIDPOINT, + n_periods: int = 1, + max_steps: int = 10_000_000, +) -> dict[str, float | int | str]: + """Compute bounce/transit and mean toroidal frequencies for one orbit. + + ``position`` is ``[s, theta, phi, v/v0, lambda]`` in the public + reference coordinates. Set ``n_periods=1`` for a symmetric configuration + or use more periods to quantify cycle variation in an asymmetric field. + Periods are seconds, angles radians, and frequencies rad/s. + """ + if not _initialized: + raise RuntimeError("SIMPLE not initialized. Call pysimple.init() first.") + + if isinstance(integrator, str): + key = integrator.lower() + if key not in _INTEGRATOR_ALIASES: + raise ValueError(f"Unknown integrator: {integrator}") + integrator_code = _INTEGRATOR_ALIASES[key] + else: + integrator_code = int(integrator) + + position = np.ascontiguousarray(position, dtype=np.float64) + if position.shape != (5,): + raise ValueError(f"position must have shape (5,), got {position.shape}") + if n_periods < 1: + raise ValueError("n_periods must be at least one") + if max_steps < 1: + raise ValueError("max_steps must be at least one") + + params.integmode = integrator_code + _ensure_trace_initialized() + _maybe_sync_single_surface_state(position) + + values = np.zeros(6, dtype=np.float64) + metadata = np.zeros(5, dtype=np.int32) + test_bounds = _capture_test_field_bounds() if _is_test_field() else None + + try: + if test_bounds is not None: + _apply_test_field_bounds(float(position[0])) + _simple_main.compute_canonical_frequencies_flat( + _tracer, + position, + int(n_periods), + int(max_steps), + values, + metadata, + ) + finally: + if test_bounds is not None: + _restore_test_field_bounds(test_bounds) + + orbit_names = { + ORBIT_CLASS_UNKNOWN: "unknown", + ORBIT_CLASS_TRAPPED: "trapped", + ORBIT_CLASS_PASSING: "passing", + } + orbit_class = int(metadata[1]) + return { + "period": float(values[0]), + "period_std": float(values[1]), + "delta_phi": float(values[2]), + "delta_phi_std": float(values[3]), + "omega_b": float(values[4]), + "omega_phi": float(values[5]), + "status": int(metadata[0]), + "orbit_class": orbit_names.get(orbit_class, "unknown"), + "orbit_class_code": orbit_class, + "parallel_direction": int(metadata[2]), + "n_periods": int(metadata[3]), + "n_steps": int(metadata[4]), + } + + +def trace_to_cut( + position: np.ndarray, + *, + cut: str | int = "tip", + integrator: str | int = MIDPOINT, + max_events: int = 100, +) -> dict[str, np.ndarray | float | int | str]: + """Trace one orbit to an interpolated tip or toroidal-period cut.""" + if not _initialized: + raise RuntimeError("SIMPLE not initialized. Call pysimple.init() first.") + + cut_aliases = {"any": -1, "tip": CUT_TIP, "toroidal": CUT_TOROIDAL} + if isinstance(cut, str): + key = cut.lower() + if key not in cut_aliases: + raise ValueError(f"Unknown cut: {cut}") + cut_code = cut_aliases[key] + else: + cut_code = int(cut) + if cut_code not in {-1, CUT_TIP, CUT_TOROIDAL}: + raise ValueError(f"Unknown cut: {cut}") + + if isinstance(integrator, str): + key = integrator.lower() + if key not in _INTEGRATOR_ALIASES: + raise ValueError(f"Unknown integrator: {integrator}") + integrator_code = _INTEGRATOR_ALIASES[key] + else: + integrator_code = int(integrator) + if integrator_code <= RK45: + raise ValueError("trace_to_cut requires a symplectic integrator") + if max_events < 1: + raise ValueError("max_events must be at least one") + + position = np.ascontiguousarray(position, dtype=np.float64) + if position.shape != (5,): + raise ValueError(f"position must have shape (5,), got {position.shape}") + + params.integmode = integrator_code + _ensure_trace_initialized() + _maybe_sync_single_surface_state(position) + cut_state = np.zeros(6, dtype=np.float64) + cut_type, status = _simple_main.trace_to_cut_flat( + _tracer, + position, + cut_code, + int(max_events), + cut_state, + ) + cut_names = {CUT_TIP: "tip", CUT_TOROIDAL: "toroidal"} + return { + "state": np.ascontiguousarray(cut_state[:5]), + "parallel_invariant": float(cut_state[5]), + "cut_type": cut_names.get(int(cut_type), "unknown"), + "cut_type_code": int(cut_type), + "status": int(status), + } + + def trace_parallel( positions: np.ndarray, integrator: str | int = MIDPOINT, diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c6ea2aa8..56c09524 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -49,6 +49,7 @@ cut_detector.f90 classification.f90 simple.f90 + orbit_frequencies.f90 parse_ants.f90 profiles.f90 collis_alphas.f90 @@ -92,6 +93,10 @@ list(APPEND SOURCES add_library (simple STATIC ${SOURCES}) + target_include_directories(simple PUBLIC + $ + ) + if(SIMPLE_ENABLE_CGAL) target_link_libraries(simple PUBLIC stl_wall_cgal) endif() diff --git a/src/orbit_frequencies.f90 b/src/orbit_frequencies.f90 new file mode 100644 index 00000000..5d5b93dc --- /dev/null +++ b/src/orbit_frequencies.f90 @@ -0,0 +1,200 @@ +module orbit_frequencies + use, intrinsic :: iso_fortran_env, only: dp => real64 + use simple, only: tracer_t, init_sympl, tstep + use orbit_symplectic_base, only: SYMPLECTIC_STEP_OUTSIDE_DOMAIN, & + SYMPLECTIC_STEP_MAXITER, & + SYMPLECTIC_STEP_BOUNDARY + use util, only: twopi + + implicit none + private + + integer, parameter, public :: ORBIT_CLASS_UNKNOWN = 0 + integer, parameter, public :: ORBIT_CLASS_TRAPPED = 1 + integer, parameter, public :: ORBIT_CLASS_PASSING = 2 + + integer, parameter, public :: FREQ_SUCCESS = 0 + integer, parameter, public :: FREQ_INVALID_INPUT = 1 + integer, parameter, public :: FREQ_ORBIT_LOST = 2 + integer, parameter, public :: FREQ_INTEGRATOR_ERROR = 3 + integer, parameter, public :: FREQ_MAX_STEPS = 4 + + type, public :: frequency_options_t + integer :: n_periods = 1 + integer :: max_steps = 10000000 + end type frequency_options_t + + type, public :: frequency_result_t + integer :: status = FREQ_INVALID_INPUT + integer :: orbit_class = ORBIT_CLASS_UNKNOWN + integer :: parallel_direction = 0 + integer :: n_periods = 0 + integer :: n_steps = 0 + real(dp) :: period = 0.0_dp + real(dp) :: period_std = 0.0_dp + real(dp) :: delta_phi = 0.0_dp + real(dp) :: delta_phi_std = 0.0_dp + real(dp) :: omega_b = 0.0_dp + real(dp) :: omega_phi = 0.0_dp + end type frequency_result_t + + public :: compute_canonical_frequencies + +contains + + subroutine compute_canonical_frequencies(tracer, initial_state, options, result) + type(tracer_t), intent(in) :: tracer + real(dp), intent(in) :: initial_state(5) + type(frequency_options_t), intent(in) :: options + type(frequency_result_t), intent(out) :: result + + type(tracer_t) :: work + real(dp) :: z(5), z_previous(5) + real(dp) :: event_time, event_phi, previous_event_time, previous_event_phi + real(dp) :: target_theta, fraction, step_time + real(dp), allocatable :: periods(:), displacements(:) + integer :: ierr, step, n_events, direction + logical :: have_previous_event, event_found + + result = frequency_result_t() + if (options%n_periods < 1 .or. options%max_steps < 1) return + if (tracer%dtaumin <= 0.0_dp .or. tracer%v0 <= 0.0_dp) return + + allocate (periods(options%n_periods), displacements(options%n_periods)) + periods = 0.0_dp + displacements = 0.0_dp + work = tracer + z = initial_state + step_time = work%dtaumin/work%v0 + + if (work%integmode > 0) then + call init_sympl(work%si, work%f, z, work%dtaumin, work%dtaumin, & + work%relerr, work%integmode) + else + work%dtau = work%dtaumin + end if + + have_previous_event = .false. + n_events = 0 + direction = 0 + target_theta = 0.0_dp + + do step = 1, options%max_steps + z_previous = z + if (work%integmode > 0) then + call tstep(work%si, work%f, z, ierr) + else + call tstep(work, z, ierr) + end if + result%n_steps = step + + if (ierr /= 0) then + ! A step that stops on the outer radial boundary is a physical + ! loss through the plasma edge, exactly like leaving the domain; + ! only genuinely numerical failures report an integrator error. + if (ierr == SYMPLECTIC_STEP_OUTSIDE_DOMAIN & + .or. ierr == SYMPLECTIC_STEP_MAXITER & + .or. ierr == SYMPLECTIC_STEP_BOUNDARY) then + result%status = FREQ_ORBIT_LOST + else + result%status = FREQ_INTEGRATOR_ERROR + end if + return + end if + + event_found = .false. + fraction = 0.0_dp + + if (result%orbit_class == ORBIT_CLASS_UNKNOWN) then + if (z_previous(5) < 0.0_dp .and. z(5) >= 0.0_dp) then + result%orbit_class = ORBIT_CLASS_TRAPPED + result%parallel_direction = 0 + fraction = zero_crossing_fraction(z_previous(5), z(5)) + event_found = .true. + else if (abs(z(2) - initial_state(2)) >= twopi) then + result%orbit_class = ORBIT_CLASS_PASSING + direction = merge(1, -1, z(2) > initial_state(2)) + result%parallel_direction = direction + target_theta = initial_state(2) + real(direction, dp)*twopi + fraction = zero_crossing_fraction(z_previous(2) - target_theta, & + z(2) - target_theta) + event_found = .true. + target_theta = target_theta + real(direction, dp)*twopi + end if + else if (result%orbit_class == ORBIT_CLASS_TRAPPED) then + if (z_previous(5) < 0.0_dp .and. z(5) >= 0.0_dp) then + fraction = zero_crossing_fraction(z_previous(5), z(5)) + event_found = .true. + end if + else + if (direction > 0) then + if (z_previous(2) < target_theta .and. z(2) >= target_theta) then + fraction = zero_crossing_fraction( & + z_previous(2) - target_theta, z(2) - target_theta & + ) + event_found = .true. + end if + else + if (z_previous(2) > target_theta .and. z(2) <= target_theta) then + fraction = zero_crossing_fraction( & + z_previous(2) - target_theta, z(2) - target_theta & + ) + event_found = .true. + end if + end if + if (event_found) target_theta = target_theta + real(direction, dp)*twopi + end if + + if (.not. event_found) cycle + + event_time = (real(step - 1, dp) + fraction)*step_time + event_phi = z_previous(3) + fraction*(z(3) - z_previous(3)) + if (have_previous_event) then + n_events = n_events + 1 + periods(n_events) = event_time - previous_event_time + displacements(n_events) = event_phi - previous_event_phi + if (n_events == options%n_periods) exit + end if + previous_event_time = event_time + previous_event_phi = event_phi + have_previous_event = .true. + end do + + if (n_events /= options%n_periods) then + result%status = FREQ_MAX_STEPS + result%n_periods = n_events + return + end if + + result%n_periods = n_events + result%period = sum(periods)/real(n_events, dp) + result%delta_phi = sum(displacements)/real(n_events, dp) + if (n_events > 1) then + result%period_std = sample_std(periods, result%period) + result%delta_phi_std = sample_std(displacements, result%delta_phi) + end if + result%omega_b = twopi/result%period + result%omega_phi = sum(displacements)/sum(periods) + result%status = FREQ_SUCCESS + end subroutine compute_canonical_frequencies + + pure real(dp) function zero_crossing_fraction(left, right) result(fraction) + real(dp), intent(in) :: left, right + real(dp) :: denominator + + denominator = right - left + if (abs(denominator) <= tiny(denominator)) then + fraction = 0.5_dp + else + fraction = -left/denominator + end if + fraction = max(0.0_dp, min(1.0_dp, fraction)) + end function zero_crossing_fraction + + pure real(dp) function sample_std(values, mean_value) result(std_value) + real(dp), intent(in) :: values(:), mean_value + + std_value = sqrt(sum((values - mean_value)**2)/real(size(values) - 1, dp)) + end function sample_std + +end module orbit_frequencies diff --git a/src/simple_main.f90 b/src/simple_main.f90 index 5b4e78a8..86d0cbc7 100644 --- a/src/simple_main.f90 +++ b/src/simple_main.f90 @@ -909,6 +909,102 @@ subroutine init_counters boundary_event_time_width = -1d0 end subroutine init_counters + subroutine compute_canonical_frequencies_flat(anorb, initial_state, & + n_periods, max_steps, & + values, metadata) + use orbit_frequencies, only: frequency_options_t, frequency_result_t, & + compute_canonical_frequencies + + type(tracer_t), intent(inout) :: anorb + real(dp), intent(in) :: initial_state(5) + integer, intent(in) :: n_periods, max_steps + real(dp), intent(out) :: values(6) + integer, intent(out) :: metadata(5) + + type(frequency_options_t) :: options + type(frequency_result_t) :: result + real(dp) :: z(5) + + call ref_to_integ(initial_state(1:3), z(1:3)) + z(4:5) = initial_state(4:5) + + anorb%dtaumin = dtaumin + anorb%dtau = dtaumin + anorb%v0 = v0 + anorb%relerr = relerr + anorb%integmode = integmode + options%n_periods = n_periods + options%max_steps = max_steps + + call compute_canonical_frequencies(anorb, z, options, result) + + values(1) = result%period + values(2) = result%period_std + values(3) = result%delta_phi + values(4) = result%delta_phi_std + values(5) = result%omega_b + values(6) = result%omega_phi + metadata(1) = result%status + metadata(2) = result%orbit_class + metadata(3) = result%parallel_direction + metadata(4) = result%n_periods + metadata(5) = result%n_steps + end subroutine compute_canonical_frequencies_flat + + subroutine trace_to_cut_flat(anorb, initial_state, requested_cut_type, & + max_events, cut_state, cut_type, status) + use cut_detector, only: cut_detector_t, init_cut_detector => init, & + trace_to_cut + + type(tracer_t), intent(inout) :: anorb + real(dp), intent(in) :: initial_state(5) + integer, intent(in) :: requested_cut_type, max_events + real(dp), intent(out) :: cut_state(6) + integer, intent(out) :: cut_type, status + + type(cut_detector_t) :: detector + real(dp) :: z(5), reference_position(3) + integer :: event_index, ierr + + cut_state = 0.0_dp + cut_type = -1 + status = 1 + if (requested_cut_type < -1 .or. requested_cut_type > 1) return + if (max_events < 1 .or. integmode <= 0) return + + call ref_to_integ(initial_state(1:3), z(1:3)) + z(4:5) = initial_state(4:5) + anorb%dtaumin = dtaumin + anorb%dtau = dtaumin + anorb%v0 = v0 + anorb%relerr = relerr + anorb%integmode = integmode + + call init_sympl(anorb%si, anorb%f, z, dtaumin, dtaumin, relerr, integmode) + call init_cut_detector(detector, anorb%fper, z) + + do event_index = 1, max_events + call trace_to_cut(detector, anorb%si, anorb%f, z, cut_state, & + cut_type, ierr) + if (ierr /= 0) then + status = ierr + return + end if + if (requested_cut_type == -1) exit + if (cut_type == requested_cut_type) exit + end do + + if (requested_cut_type /= -1) then + if (cut_type /= requested_cut_type) then + status = 4 + return + end if + end if + + call integ_to_ref(cut_state(1:3), reference_position) + cut_state(1:3) = reference_position + status = 0 + end subroutine trace_to_cut_flat pure integer function classify_orbit_exit(ierr_orbit, model, mode, & boundary_is_lcfs) integer, intent(in) :: ierr_orbit, model, mode diff --git a/test/python/CMakeLists.txt b/test/python/CMakeLists.txt index ca300c4d..59616cf0 100644 --- a/test/python/CMakeLists.txt +++ b/test/python/CMakeLists.txt @@ -27,11 +27,27 @@ function(add_simple_api_test test_name test_node) ) endfunction() +function(add_chartmap_api_test test_name test_node) + add_test(NAME test_simple_api_${test_name} + COMMAND ${Python_EXECUTABLE} -m pytest + ${CMAKE_CURRENT_SOURCE_DIR}/test_simple_api.py::${test_node} -v + ) + set_tests_properties(test_simple_api_${test_name} PROPERTIES + ENVIRONMENT "PYTHONPATH=${CMAKE_BINARY_DIR}/python:${CMAKE_BINARY_DIR}:${CMAKE_SOURCE_DIR}/python;SIMPLE_TEST_CHARTMAP=${CMAKE_BINARY_DIR}/test/tests/test_boozer_chartmap.nc" + LABELS "python;simple_api;chartmap;slow" + TIMEOUT 120 + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + DEPENDS "generate_test_boozer_chartmap_data" + ) +endfunction() + add_simple_api_test(integration_constants TestConstants::test_integration_constants) add_simple_api_test(default_values TestConstants::test_default_values) add_simple_api_test(init_basic TestInitialization::test_init_basic) add_simple_api_test(init_parameters TestInitialization::test_init_with_parameters) +add_simple_api_test(init_namelist_case TestInitialization::test_init_accepts_namelist_case) add_simple_api_test(init_repeated TestInitialization::test_multiple_reinitializations) +add_chartmap_api_test(init_chartmap TestInitialization::test_chartmap_init_selects_boozer_field) add_simple_api_test(integrator_alias TestIntegratorSelection::test_alias_resolution) add_simple_api_test(sampler_shape TestSampling::test_surface_sampler_shape) add_simple_api_test(sampler_sizes TestSampling::test_sample_different_sizes) @@ -45,6 +61,14 @@ add_simple_api_test(trace_filters TestTraceOrbits::test_confined_and_lost_filter add_simple_api_test(trace_integrators TestTraceOrbits::test_different_integrators) add_simple_api_test(trace_single TestTraceOrbits::test_trace_single_orbit) add_simple_api_test(trace_trajectory TestTraceOrbits::test_trace_with_trajectory) +add_simple_api_test(canonical_frequencies + TestTraceOrbits::test_canonical_frequencies) +add_simple_api_test(canonical_frequency_validation + TestTraceOrbits::test_canonical_frequency_validation) +add_simple_api_test(trace_to_tip + TestTraceOrbits::test_trace_to_tip) +add_simple_api_test(trace_to_cut_validation + TestTraceOrbits::test_trace_to_cut_validation) add_simple_api_test(trace_skipped TestTraceOrbits::test_trace_parallel_preserves_skipped_passing_sentinel) add_simple_api_test(trace_early_loss TestTraceOrbits::test_trace_orbit_uses_backend_loss_time_for_early_loss) diff --git a/test/python/test_simple_api.py b/test/python/test_simple_api.py index 0a718554..5180cc2b 100644 --- a/test/python/test_simple_api.py +++ b/test/python/test_simple_api.py @@ -5,6 +5,7 @@ from __future__ import annotations +import os from pathlib import Path import numpy as np @@ -62,12 +63,49 @@ def test_init_with_parameters(self, vmec_file: str): npoiper2=64 ) + def test_init_accepts_namelist_case(self, vmec_file: str): + """init() should accept namelist-cased names such as facE_al.""" + pysimple.init( + vmec_file, + deterministic=True, + trace_time=1e-4, + ntestpart=1, + facE_al=700.0, + n_d=2, + n_e=1, + ) + assert float(pysimple.params.face_al) == 700.0 + def test_multiple_reinitializations(self, vmec_file: str): """Multiple init() calls should work without crash""" pysimple.init(vmec_file, deterministic=True, trace_time=1e-4) pysimple.init(vmec_file, deterministic=True, trace_time=1e-3) pysimple.init(vmec_file, deterministic=True, trace_time=5e-5) + def test_chartmap_init_selects_boozer_field(self): + """Chart-map initialization must not attempt to load a VMEC wout.""" + chartmap_file = os.environ["SIMPLE_TEST_CHARTMAP"] + pysimple.init( + chartmap_file, + deterministic=True, + ntestpart=1, + npoiper2=64, + trace_time=1e-4, + ) + assert int(pysimple._fortran_backend.velo_mod.isw_field_type) == 2 + result = pysimple.compute_canonical_frequencies( + np.array([0.4, 0.7, 0.1, 1.0, 0.1]), + n_periods=1, + max_steps=10, + ) + # A boundary stop is a physical loss (FREQ_ORBIT_LOST), never an + # unexplained FREQ_INTEGRATOR_ERROR. + assert result["status"] in { + pysimple.FREQ_SUCCESS, + pysimple.FREQ_ORBIT_LOST, + pysimple.FREQ_MAX_STEPS, + } + class TestIntegratorSelection: def test_alias_resolution(self, vmec_file: str): @@ -220,6 +258,57 @@ def test_trace_with_trajectory(self, vmec_file: str): assert result['trajectory'].shape[0] == 5 assert result['times'].shape[0] == result['trajectory'].shape[1] + def test_canonical_frequencies(self, vmec_file: str): + pysimple.init( + vmec_file, + deterministic=True, + ntestpart=1, + npoiper2=1024, + trace_time=1e-3, + ) + trapped = np.array([0.4, 0.7, 0.1, 1.0, 0.1]) + result = pysimple.compute_canonical_frequencies( + trapped, + integrator="midpoint", + n_periods=2, + ) + + assert result["status"] == pysimple.FREQ_SUCCESS + assert result["orbit_class"] == "trapped" + assert result["n_periods"] == 2 + assert result["period"] > 0.0 + assert result["omega_b"] > 0.0 + assert result["omega_b"] == pytest.approx(2.0 * np.pi / result["period"]) + + def test_canonical_frequency_validation(self, vmec_file: str): + pysimple.init(vmec_file, deterministic=True, ntestpart=1) + particle = np.array([0.4, 0.7, 0.1, 1.0, 0.9]) + + with pytest.raises(ValueError, match="n_periods"): + pysimple.compute_canonical_frequencies(particle, n_periods=0) + with pytest.raises(ValueError, match="max_steps"): + pysimple.compute_canonical_frequencies(particle, max_steps=0) + with pytest.raises(ValueError, match="shape"): + pysimple.compute_canonical_frequencies(particle[:4]) + + def test_trace_to_tip(self, vmec_file: str): + pysimple.init(vmec_file, deterministic=True, ntestpart=1, npoiper2=1024) + particle = np.array([0.4, 0.7, 0.1, 1.0, 0.1]) + result = pysimple.trace_to_cut(particle, cut="tip", max_events=10) + + assert result["status"] == 0 + assert result["cut_type"] == "tip" + assert result["state"].shape == (5,) + assert abs(result["state"][4]) < 1e-6 + + def test_trace_to_cut_validation(self, vmec_file: str): + pysimple.init(vmec_file, deterministic=True, ntestpart=1) + particle = np.array([0.4, 0.7, 0.1, 1.0, 0.9]) + with pytest.raises(ValueError, match="Unknown cut"): + pysimple.trace_to_cut(particle, cut="poloidal") + with pytest.raises(ValueError, match="symplectic"): + pysimple.trace_to_cut(particle, integrator="rk45") + def test_trace_parallel_preserves_skipped_passing_sentinel(self, vmec_file: str): """Deep-passing particles skipped by contr_pp must keep loss_time = -1.""" pysimple.init(vmec_file, deterministic=True, trace_time=1e-4, ntestpart=1) diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt index df647159..895625b9 100644 --- a/test/tests/CMakeLists.txt +++ b/test/tests/CMakeLists.txt @@ -350,6 +350,13 @@ if (ENABLE_OPENMP) COMMAND test_sympl_testfield.x WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + add_executable(test_orbit_frequencies.x test_orbit_frequencies.f90) + target_link_libraries(test_orbit_frequencies.x simple) + add_test(NAME test_orbit_frequencies + COMMAND test_orbit_frequencies.x + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + set_tests_properties(test_orbit_frequencies PROPERTIES TIMEOUT 120 LABELS "physics") + # Axis-crossing regression (issue #370): near-axis trapped orbits must # never be lost at the axis and must conserve energy across crossings. add_executable(test_axis_crossing.x test_axis_crossing.f90) diff --git a/test/tests/test_orbit_frequencies.f90 b/test/tests/test_orbit_frequencies.f90 new file mode 100644 index 00000000..21b8fc5d --- /dev/null +++ b/test/tests/test_orbit_frequencies.f90 @@ -0,0 +1,60 @@ +program test_orbit_frequencies + use, intrinsic :: iso_fortran_env, only: dp => real64 + use orbit_frequencies, only: frequency_options_t, frequency_result_t, & + compute_canonical_frequencies, FREQ_SUCCESS, & + FREQ_INVALID_INPUT, ORBIT_CLASS_TRAPPED, & + ORBIT_CLASS_PASSING + use params, only: coord_input, field_input + use simple, only: tracer_t, init_params + use simple_main, only: init_field + + implicit none + + type(tracer_t) :: tracer + type(frequency_options_t) :: options + type(frequency_result_t) :: result + real(dp) :: trapped(5), passing(5) + + field_input = 'wout.nc' + coord_input = 'wout.nc' + call init_field(tracer, 'wout.nc', 5, 5, 3, 3) + call init_params(tracer, 1, 2, 5.0e3_dp, 1024, 1, 1.0e-10_dp) + tracer%integmode = 3 + + options%n_periods = 3 + options%max_steps = 2000000 + + trapped = [0.4_dp, 0.7_dp, 0.1_dp, 1.0_dp, 0.1_dp] + call compute_canonical_frequencies(tracer, trapped, options, result) + call require(result%status == FREQ_SUCCESS, 'trapped frequency status') + call require(result%orbit_class == ORBIT_CLASS_TRAPPED, 'trapped classification') + call require(result%n_periods == options%n_periods, 'trapped period count') + call require(result%period > 0.0_dp, 'trapped positive period') + call require(result%omega_b > 0.0_dp, 'trapped positive omega_b') + + passing = [0.4_dp, 0.7_dp, 0.1_dp, 1.0_dp, 0.9_dp] + call compute_canonical_frequencies(tracer, passing, options, result) + call require(result%status == FREQ_SUCCESS, 'passing frequency status') + call require(result%orbit_class == ORBIT_CLASS_PASSING, 'passing classification') + call require(result%n_periods == options%n_periods, 'passing period count') + call require(result%parallel_direction /= 0, 'passing direction') + call require(result%period > 0.0_dp, 'passing positive period') + call require(result%omega_b > 0.0_dp, 'passing positive omega_b') + + options%n_periods = 0 + call compute_canonical_frequencies(tracer, passing, options, result) + call require(result%status == FREQ_INVALID_INPUT, 'invalid options status') + +contains + + subroutine require(condition, message) + logical, intent(in) :: condition + character(*), intent(in) :: message + + if (.not. condition) then + write (*, '(A)') 'FAILED: '//message + error stop 1 + end if + end subroutine require + +end program test_orbit_frequencies From 48551773c1b682b592112a94f68ce5276e1d602a Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Fri, 17 Jul 2026 07:52:59 +0200 Subject: [PATCH 6/7] Order pyplot module before SIMPLE compilation (#500) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ## Summary - add the explicit target-order dependency required when `diag_albert.f90` consumes `pyplot_module.mod` from the separate `pyplot` target - prevent clean parallel Ninja builds from compiling `diag_albert` before the module file exists ## Evidence The CGAL wall CI job on PR #499 failed in a clean parallel build with: ``` Fatal Error: Cannot open module file ‘pyplot_module.mod’ ``` The existing link edge orders final linking but did not reliably order compilation across targets. This PR adds no source or runtime behavior change. ## Validation - clean `fo build` passed - `fo` passed all stages in 308.3 s --- src/CMakeLists.txt | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 56c09524..dd717fbf 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -103,6 +103,11 @@ list(APPEND SOURCES # Link pyplot library to SIMPLE target_link_libraries(simple PUBLIC pyplot) +# `diag_albert.f90` consumes pyplot_module.mod while `pyplot` is a separate +# target. A link edge alone does not order compilation of the dependent target +# under Ninja, so a clean parallel build can compile diag_albert before the +# module file exists. +add_dependencies(simple pyplot) # Apply SIMPLE-specific compile options target_compile_options(simple PRIVATE ${SIMPLE_COMPILE_OPTIONS}) From 22bd167aed475769c210fbea74f3a0dfe135946e Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Fri, 17 Jul 2026 08:23:31 +0200 Subject: [PATCH 7/7] Accept strict unresolved boundary status --- test/tests/test_newton_solver_status.f90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/test/tests/test_newton_solver_status.f90 b/test/tests/test_newton_solver_status.f90 index 50fa5510..bb33096c 100644 --- a/test/tests/test_newton_solver_status.f90 +++ b/test/tests/test_newton_solver_status.f90 @@ -277,9 +277,13 @@ subroutine test_newton_warning_mode strict_integrator%atol = 0.0_dp symplectic_newton_warning_mode = .false. call orbit_timestep_sympl(strict_integrator, strict_field, step_status) - if (step_status /= SYMPLECTIC_STEP_MAXITER) then + ! The Lobatto path can detect a candidate boundary while exhausting the + ! forced-zero tolerance solve. Both statuses are strict numerical + ! failures, and both must roll back the accepted state. + if (step_status /= SYMPLECTIC_STEP_MAXITER .and. & + step_status /= SYMPLECTIC_STEP_EVENT_NOT_CONVERGED) then print *, 'strict mode, status:', modes(mode_index), step_status - error stop 'strict timestep did not report max iterations' + error stop 'strict timestep did not report a numerical failure' end if if (any(strict_integrator%z /= initial_state)) then error stop 'strict max-iteration failure changed the accepted state'