diff --git a/src/interface_crossing.f90 b/src/interface_crossing.f90 index 00ab3ed2..b6eb70bc 100644 --- a/src/interface_crossing.f90 +++ b/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 diff --git a/src/spectre_orbit.f90 b/src/spectre_orbit.f90 index d0187cc8..995f4019 100644 --- a/src/spectre_orbit.f90 +++ b/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 diff --git a/src/spectre_sympl_orbit.f90 b/src/spectre_sympl_orbit.f90 index aed17f93..176d9dc6 100644 --- a/src/spectre_sympl_orbit.f90 +++ b/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 diff --git a/test/tests/test_interface_crossing_ownership.f90 b/test/tests/test_interface_crossing_ownership.f90 index 379f5e7c..38c51fef 100644 --- a/test/tests/test_interface_crossing_ownership.f90 +++ b/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