Require valid states for SPECTRE physical losses

src/interface_crossing.f90

@@ -14,7 +14,7 @@ module interface_crossing
     !> the tangential sheet-drift kick and the drift-order v_par term that Level 0
     !> omits. Both levels are energy-exact in the crossing and reflection branch.
 
-    use, intrinsic :: iso_fortran_env, only: dp => real64
+    use, intrinsic :: iso_fortran_env, only: dp => real64, int64
     use magfie_sub, only: magfie
     use parmot_mod, only: ro0
 
@@ -25,7 +25,7 @@ module interface_crossing
     public :: crossing_log_reset, crossing_log_record, crossing_log_write, &
         crossing_log_count, crossing_log_count_type
     public :: CROSSING_LEVEL0, CROSSING_LEVEL1, CROSS_CROSSING, CROSS_REFLECTION, &
-        CROSS_LOSS, CROSS_STOP, CROSS_SHEET, CROSS_RECOVERY
+        CROSS_LOSS, CROSS_STOP, CROSS_SHEET, CROSS_RECOVERY, CROSS_INVALID
 
     integer, parameter :: CROSSING_LEVEL0 = 0
     integer, parameter :: CROSSING_LEVEL1 = 1
@@ -41,6 +41,9 @@ module interface_crossing
     !> iface and both volume fields are zero so post-processing cannot confuse
     !> the recovery map with a discontinuity crossing.
     integer, parameter :: CROSS_RECOVERY = 6
+    !> A candidate event whose physical state or field values are unusable.
+    !> Callers must recover from the pre-step state; this is never a loss.
+    integer, parameter :: CROSS_INVALID = 7
 
     !> Inner cutoff for the innermost volume: rho_g = 0 is the coordinate axis
     !> where sqrt(g) = 0, so the marker reflects trivially before reaching it
@@ -120,15 +123,22 @@ subroutine apply_crossing(y_iface, iface, direction, mvol, level, y_out, info)
         type(crossing_info_t), intent(out) :: info
 
         real(dp) :: rho_face, rho_home, rho_target
-        real(dp) :: bmod_home, bmod_target, perp_inv, mu, vpar, radicand
+        real(dp) :: bmod_home, bmod_target, perp_inv, mu, vpar, radicand, pitch
+
+        info%event_type = CROSS_INVALID
+        y_out = y_iface
 
         if (level /= CROSSING_LEVEL0 .and. level /= CROSSING_LEVEL1) then
             error stop 'apply_crossing: unknown crossing_level'
         end if
+        if ((direction /= -1 .and. direction /= 1) .or. iface < 0 .or. &
+                iface > mvol .or. .not. all_finite(y_iface) .or. &
+                y_iface(4) <= 0.0_dp .or. abs(y_iface(5)) > 1.0_dp) return
 
         rho_face = real(iface, dp)
-        y_out = y_iface
-        vpar = y_iface(4)*y_iface(5)
+        pitch = max(-1.0_dp, min(1.0_dp, y_iface(5)))
+        y_out(5) = pitch
+        vpar = y_iface(4)*pitch
 
         info%iface = iface
         info%theta = y_iface(2)
@@ -149,14 +159,16 @@ subroutine apply_crossing(y_iface, iface, direction, mvol, level, y_out, info)
         end if
 
         bmod_home = bmod_at([rho_home, y_iface(2), y_iface(3)])
+        if (.not. finite_value(bmod_home) .or. bmod_home <= 0.0_dp) return
         info%bmod_home = bmod_home
         info%bmod_target = bmod_home
 
         ! perp_inv = v_perp^2/|B| = z(4)^2 (1 - z(5)^2)/|B| is the perpendicular
         ! invariant the RK45 path stores; mu = perp_inv/2 is the magnetic moment
         ! held fixed across the interface, so the whole map runs on this invariant.
-        perp_inv = y_iface(4)**2*(1.0_dp - y_iface(5)**2)/bmod_home
+        perp_inv = y_iface(4)**2*max(1.0_dp - pitch**2, 0.0_dp)/bmod_home
         mu = 0.5_dp*perp_inv
+        if (.not. finite_value(mu) .or. mu < 0.0_dp) return
         info%mu = mu
 
         if (direction == 1 .and. iface == mvol) then
@@ -177,13 +189,19 @@ subroutine apply_crossing(y_iface, iface, direction, mvol, level, y_out, info)
         if (level == CROSSING_LEVEL1) then
             call level1_map(rho_face, rho_home, rho_target, direction, &
                 bmod_home, mu, vpar, y_iface, y_out, info)
+            if (.not. valid_mapped_crossing(y_out, info)) then
+                info%event_type = CROSS_INVALID
+                y_out = y_iface
+            end if
             return
         end if
 
         bmod_target = bmod_at([rho_target, y_iface(2), y_iface(3)])
+        if (.not. finite_value(bmod_target) .or. bmod_target <= 0.0_dp) return
         info%bmod_target = bmod_target
 
         radicand = vpar**2 - 2.0_dp*mu*(bmod_target - bmod_home)
+        if (.not. finite_value(radicand)) return
         if (radicand >= 0.0_dp) then
             info%event_type = CROSS_CROSSING
             info%vpar_after = sign(sqrt(radicand), vpar)
@@ -202,6 +220,39 @@ subroutine apply_crossing(y_iface, iface, direction, mvol, level, y_out, info)
         end if
     end subroutine apply_crossing
 
+    pure logical function valid_mapped_crossing(y, info)
+        real(dp), intent(in) :: y(5)
+        type(crossing_info_t), intent(in) :: info
+
+        valid_mapped_crossing = all_finite(y) .and. y(4) > 0.0_dp .and. &
+            abs(y(5)) <= 1.0_dp + 100.0_dp*epsilon(1.0_dp) .and. &
+            finite_value(info%mu) .and. info%mu >= 0.0_dp .and. &
+            finite_value(info%bmod_home) .and. info%bmod_home > 0.0_dp .and. &
+            finite_value(info%bmod_target) .and. info%bmod_target > 0.0_dp .and. &
+            (info%event_type == CROSS_CROSSING .or. &
+             info%event_type == CROSS_REFLECTION)
+    end function valid_mapped_crossing
+
+    pure logical function finite_value(x)
+        real(dp), intent(in) :: x
+        integer(int64), parameter :: exponent_mask = int(z'7FF0000000000000', int64)
+        integer(int64) :: bits
+
+        bits = transfer(x, bits)
+        finite_value = iand(bits, exponent_mask) /= exponent_mask
+    end function finite_value
+
+    pure logical function all_finite(x)
+        real(dp), intent(in) :: x(:)
+        integer :: i
+
+        all_finite = .false.
+        do i = 1, size(x)
+            if (.not. finite_value(x(i))) return
+        end do
+        all_finite = .true.
+    end function all_finite
+
     function bmod_at(x) result(bmod)
         real(dp), intent(in) :: x(3)
         real(dp) :: bmod

src/spectre_orbit.f90

@@ -15,7 +15,7 @@ module spectre_orbit
     use odeint_allroutines_sub, only: odeint_allroutines, odeint_has_failed
     use alpha_lifetime_sub, only: velo_can
     use interface_crossing, only: apply_crossing, crossing_info_t, axis_offset, &
-        CROSS_LOSS
+        CROSS_LOSS, CROSS_INVALID
 
     implicit none
     private
@@ -139,6 +139,12 @@ subroutine orbit_timestep_spectre(state, z, dtaumin, relerr, level, ierr, event)
         iface = nint(boundary)
         call apply_crossing(z_hit, iface, direction, state%mvol, level, &
             z, event%info)
+        if (event%info%event_type == CROSS_INVALID) then
+            z = z_start
+            event%occurred = .false.
+            ierr = SPECTRE_FAULT
+            return
+        end if
         event%occurred = .true.
 
         if (event%info%event_type == CROSS_LOSS) then

src/spectre_sympl_orbit.f90

@@ -18,10 +18,10 @@ module spectre_sympl_orbit
         ref_to_integ
     use field_can_spectre, only: set_spectre_volume_lock
     use orbit_symplectic_base, only: symplectic_integrator_t, &
-        symplectic_newton_warning_mode
+        symplectic_newton_warning_mode, SYMPLECTIC_STEP_OUTSIDE_DOMAIN
     use orbit_symplectic, only: orbit_timestep_sympl, orbit_sympl_init
     use interface_crossing, only: apply_crossing, crossing_info_t, &
-        crossing_log_record, CROSS_LOSS, CROSS_STOP, CROSS_RECOVERY, &
+        crossing_log_record, CROSS_LOSS, CROSS_STOP, CROSS_RECOVERY, CROSS_INVALID, &
         CROSS_REFLECTION, CROSS_SHEET
     use spectre_sheet_gc, only: sheet_gc_state_t, sheet_gc_initialize, &
         sheet_gc_advance, sheet_gc_to_y, SHEET_GC_OK
@@ -386,13 +386,26 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, &
                 h_try = budget
                 cycle
             end if
-            used = used + h_land
-            budget = state%dt_std - used
-            call update_landing_stats(resid)
-
             call landed_state(si, f, real(iface, dp), direction, y_iface)
             call apply_crossing(y_iface, iface, direction, state%mvol, &
                 state%level, y_out, info)
+            if (info%event_type == CROSS_INVALID) then
+                si = si0
+                f = f0
+                h_try = 0.5_dp*h_try
+                if (h_try >= h_min_frac*state%dt_std) cycle
+                call recover_local_error(state, si, f, ipart, t_base_sec + &
+                    dt_sec*used/state%dt_std, STOP_LANDING, budget, used, ierr, &
+                    t_frac, resume_gc)
+                if (resume_gc) then
+                    h_try = budget
+                    cycle
+                end if
+                exit
+            end if
+            used = used + h_land
+            budget = state%dt_std - used
+            call update_landing_stats(resid)
             if (info%event_type == CROSS_LOSS) then
                 call crossing_log_record(ipart, &
                     t_base_sec + dt_sec*used/state%dt_std, info)
@@ -559,7 +572,8 @@ subroutine classify_step(state, rho_start, si, ierr_step, boundary, iface, &
             ! The stepper refuses to commit r > sympl_rmax (= Mvol). From the
             ! outermost volume that is the outward-crossing signal, with si%z
             ! still holding the pre-step state.
-            if (nint(state%home_hi) == state%mvol) then
+            if (ierr_step == SYMPLECTIC_STEP_OUTSIDE_DOMAIN .and. &
+                    nint(state%home_hi) == state%mvol) then
                 boundary = .true.
                 iface = state%mvol
                 direction = 1

test/tests/test_interface_crossing_ownership.f90

@@ -18,9 +18,10 @@ end module interface_crossing_ownership_fixture
 
 program test_interface_crossing_ownership
     use, intrinsic :: iso_fortran_env, only: dp => real64
+    use, intrinsic :: ieee_arithmetic, only: ieee_quiet_nan, ieee_value
     use magfie_sub, only: magfie
     use interface_crossing, only: apply_crossing, crossing_info_t, &
-        CROSSING_LEVEL0, CROSS_CROSSING, CROSS_REFLECTION
+        CROSSING_LEVEL0, CROSS_CROSSING, CROSS_REFLECTION, CROSS_INVALID
     use interface_crossing_ownership_fixture, only: stepped_field
 
     implicit none
@@ -45,4 +46,14 @@ program test_interface_crossing_ownership
     if (info%event_type /= CROSS_REFLECTION) error stop 'reflection type'
     if (info%vol_to /= 1) error stop 'reflection ownership'
     if (y_out(1) /= 1.0_dp) error stop 'reflection moved off interface'
+
+    y_in = [1.0_dp, 0.2_dp, 0.3_dp, 1.0_dp, 0.5_dp]
+    y_in(5) = ieee_value(0.0_dp, ieee_quiet_nan)
+    call apply_crossing(y_in, 1, 1, 3, CROSSING_LEVEL0, y_out, info)
+    if (info%event_type /= CROSS_INVALID) &
+        error stop 'non-finite pitch became a physical event'
+    y_in = [3.0_dp, 0.2_dp, 0.3_dp, -1.0_dp, 0.5_dp]
+    call apply_crossing(y_in, 3, 1, 3, CROSSING_LEVEL0, y_out, info)
+    if (info%event_type /= CROSS_INVALID) &
+        error stop 'invalid momentum became an outer loss'
 end program test_interface_crossing_ownership