Restore warning-mode orbit recovery without Newton stalls

DOC/config.md

@@ -26,6 +26,10 @@
   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`.
+  A symplectic failure may use the established adaptive-RK pusher from the last
+  accepted state. That fallback is capped at 10,000 internal steps. If both
+  pushers fail from the same state, warning mode records held intervals for the
+  rest of the trace instead of repeating an identical stalled RK solve.
 
 * `canonical_grid_nr`, `canonical_grid_ntheta`, and `canonical_grid_nphi`
   control the Meiss or Albert canonical-map grid. Their defaults are 62, 63,

src/diag_counters.f90

@@ -16,7 +16,9 @@ module diag_counters
     public :: EVT_NEWTON1_MAXIT, EVT_NEWTON2_MAXIT, EVT_RK_GAUSS_MAXIT, &
         EVT_RK_LOBATTO_MAXIT, EVT_FIXPOINT_MAXIT, EVT_R_NEGATIVE, &
         EVT_FO_LOSS, EVT_FO_FAULT, EVT_MIDPOINT_MAXIT, &
-        EVT_WARNING_STEP_SKIP, N_EVENT
+        EVT_WARNING_STEP_SKIP, EVT_SPECTRE_REF_INVERSE_MAXIT, &
+        EVT_SPECTRE_INVALID_STATE, EVT_SYMPLECTIC_RK_RECOVERY, &
+        EVT_SYMPLECTIC_RESUME, N_EVENT
     public :: diag_counters_init, count_event, diag_counters_total, &
         diag_counters_reset, event_name
 
@@ -33,7 +35,11 @@ module diag_counters
     integer, parameter :: EVT_FO_FAULT = 8
     integer, parameter :: EVT_MIDPOINT_MAXIT = 9
     integer, parameter :: EVT_WARNING_STEP_SKIP = 10
-    integer, parameter :: N_EVENT = 10
+    integer, parameter :: EVT_SPECTRE_REF_INVERSE_MAXIT = 11
+    integer, parameter :: EVT_SPECTRE_INVALID_STATE = 12
+    integer, parameter :: EVT_SYMPLECTIC_RK_RECOVERY = 13
+    integer, parameter :: EVT_SYMPLECTIC_RESUME = 14
+    integer, parameter :: N_EVENT = 14
 
     ! Whole cache lines per thread column keep neighbouring threads from sharing
     ! a line. The event id indexes within a column; STRIDE >= N_EVENT.
@@ -102,6 +108,14 @@ function event_name(id) result(name)
             name = 'midpoint_maxit'
         case (EVT_WARNING_STEP_SKIP)
             name = 'warning_step_skip'
+        case (EVT_SPECTRE_REF_INVERSE_MAXIT)
+            name = 'spectre_ref_inverse_maxit'
+        case (EVT_SPECTRE_INVALID_STATE)
+            name = 'spectre_invalid_state'
+        case (EVT_SYMPLECTIC_RK_RECOVERY)
+            name = 'symplectic_rk_recovery'
+        case (EVT_SYMPLECTIC_RESUME)
+            name = 'symplectic_resume'
         case default
             name = 'unknown'
         end select

src/field/field_can_spectre.f90

@@ -15,7 +15,7 @@ module field_can_spectre
     !> Bmod Tesla->Gauss, covariant h meter->cm, covariant A (flux-like, [B]*L^2)
     !> Tesla*m^2 -> Gauss*cm^2.
 
-    use, intrinsic :: iso_fortran_env, only: dp => real64
+    use, intrinsic :: iso_fortran_env, only: dp => real64, int64
     use field_base, only: magnetic_field_t
     use field_spectre, only: spectre_field_t
     use field_can_base, only: field_can_t, n_field_evaluations, twopi
@@ -26,6 +26,8 @@ module field_can_spectre
                            evaluate_batch_splines_3d_der, &
                            evaluate_batch_splines_3d_der2
     use magfie_sub, only: TESLA_TO_GAUSS, M_TO_CM
+    use diag_counters, only: count_event, EVT_SPECTRE_REF_INVERSE_MAXIT, &
+        EVT_SPECTRE_INVALID_STATE
 
     implicit none
     private
@@ -86,6 +88,17 @@ module field_can_spectre
 
 contains
 
+    elemental logical function ieee_is_finite(value)
+        real(dp), intent(in) :: value
+
+        integer(int64), parameter :: exponent_mask = &
+            int(z'7FF0000000000000', int64)
+        integer(int64) :: bits
+
+        bits = transfer(value, bits)
+        ieee_is_finite = iand(bits, exponent_mask) /= exponent_mask
+    end function ieee_is_finite
+
     subroutine set_spectre_volume_lock(lvol)
         !> lvol > 0 pins field evaluation to that volume; 0 restores dispatch on
         !> int(r).
@@ -345,24 +358,49 @@ subroutine ref_to_integ_spectre(xref, xinteg)
         integer, parameter :: MAX_ITER = 16
 
         integer :: lvol, i
-        real(dp) :: x(3), y(2), dy(3, 2), phi_prev
+        real(dp) :: x(3), y(2), dy(3, 2), phi, phi_next, target
+        real(dp) :: residual, correction, best_phi, best_residual, denom
+
+        if (.not. all(ieee_is_finite(xref))) then
+            xinteg = xref
+            call count_event(EVT_SPECTRE_INVALID_STATE)
+            return
+        end if
 
         lvol = active_volume(xref(1))
         xinteg(1) = xref(1)
         xinteg(2) = modulo(xref(2), twopi)
-        xinteg(3) = modulo(xref(3), twopi)
+        target = modulo(xref(3), twopi)
+        phi = target
+        best_phi = phi
+        best_residual = huge(1.0_dp)
 
         do i = 1, MAX_ITER
             x = [clamp_to_volume(xinteg(1), spectre_volumes(lvol)), xinteg(2), &
-                 xinteg(3)]
+                 modulo(phi, twopi)]
             call evaluate_batch_splines_3d_der(spectre_volumes(lvol)%spl_transform, &
                                                x, y, dy)
-            phi_prev = xinteg(3)
-            xinteg(3) = phi_prev - (phi_prev + y(1) - xref(3))/(1d0 + dy(3, 1))
-            if (abs(xinteg(3) - phi_prev) < TOL) return
+            if (.not. ieee_is_finite(y(1)) .or. &
+                .not. ieee_is_finite(dy(3, 1))) exit
+            residual = modulo(phi + y(1) - target + 0.5_dp*twopi, twopi) - &
+                0.5_dp*twopi
+            if (abs(residual) < best_residual) then
+                best_residual = abs(residual)
+                best_phi = phi
+            end if
+            denom = 1.0_dp + dy(3, 1)
+            if (abs(denom) <= 100.0_dp*epsilon(1.0_dp)) exit
+            correction = max(-0.5_dp*twopi, min(0.5_dp*twopi, residual/denom))
+            phi_next = phi - correction
+            if (.not. ieee_is_finite(phi_next)) exit
+            phi = phi_next
+            if (abs(correction) < TOL) then
+                xinteg(3) = phi + twopi*anint((xref(3) - phi)/twopi)
+                return
+            end if
         end do
-        print *, 'WARNING: ref_to_integ_spectre did not converge after', &
-            MAX_ITER, 'iterations'
+        xinteg(3) = best_phi + twopi*anint((xref(3) - best_phi)/twopi)
+        call count_event(EVT_SPECTRE_REF_INVERSE_MAXIT)
     end subroutine ref_to_integ_spectre