Report confinement over resolved orbits

README.md

@@ -239,15 +239,23 @@ Diagnostics for slow convergence of Newton iterations are written in `fort.6601`
 1. Physical time
 2. Confined fraction of passing particles
 3. Confined fraction of trapped particles
-4. Total number of particles
+4. Number of numerically resolved particles used as the denominator
 
-The sum of 2. and 3. yields the overall confined fraction at each time.
+The sum of 2. and 3. yields the overall confined fraction among resolved
+particles at each time. Numerically fatal orbits are excluded from both the
+numerator and denominator. If no orbit is resolved, both fractions are `NaN`.
+`unresolved_fraction.dat` reports the excluded numerical-failure fraction
+against the total initial population; its third column remains the total
+number of particles.
 
 `times_lost.dat` contains the loss time of each particle. Columns are:
 1. Particle index. Corresponds to line number in start.dat .
-2. Time t_loss [s] when the particle is lost. Possible values are: -1, `trace_time`, or any other value between 0 and trace_time.
+2. Time t_loss [s] when the particle is lost. Possible values are: -1, `NaN`,
+   `trace_time`, or any other value between 0 and trace_time.
   * If never lost or classified as regular, maximum tracing time `trace_time` is written.
   * If ignored due to contr_pp, which defines deep passing region as confined (we don consider them anymore), -1 is written.
+  * If tracing ends in a numerical fatal condition, `NaN` is written. Such an
+    orbit is unresolved, not physically lost or confined.
 3. Trapping parameter trap_par that is 1 for deeply trapped, 0 for tp boundary
    and negative for passing. Eq. (3.1) in Accelerated Methods paper.
    Whenever trap_par < contr_pp, particle is not traced and counted as confined.

python/pysimple/plotting.py

@@ -44,7 +44,7 @@
     Column 1: Time [s]
     Column 2: confpart_pass (confined passing fraction)
     Column 3: confpart_trap (confined trapped fraction)
-    Column 4: Total number of particles
+    Column 4: Numerically resolved particles used as the denominator
     Note: Total confined fraction = Column 2 + Column 3
 
 Example
@@ -210,7 +210,7 @@ def load_loss_data(directory: str | Path) -> LossData:
         start_pitch = np.zeros(n_particles)
 
     # Load confined_fraction.dat
-    # Columns: time, confpart_pass, confpart_trap, ntestpart
+    # Columns: time, confpart_pass, confpart_trap, resolved particle count
     conf_file = directory / "confined_fraction.dat"
     if conf_file.exists():
         conf_data = np.loadtxt(conf_file)

src/restart.f90

@@ -2,11 +2,11 @@ module restart_mod
     use, intrinsic :: iso_fortran_env, only: dp => real64
     use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
     use params, only: ntestpart, times_lost, orbit_exit_code, &
-                      boundary_event_radial_residual, boundary_event_time_width, &
-                      trap_par, perp_inv, zend, &
-                      confpart_pass, confpart_trap, ntimstep, kt_macro, &
-                      v0, dtaumin, trace_time, ORBIT_EXIT_COMPLETED, &
-                      ORBIT_EXIT_LCFS
+        boundary_event_radial_residual, boundary_event_time_width, &
+        trap_par, perp_inv, zend, &
+        confpart_pass, confpart_trap, ntimstep, kt_macro, &
+        v0, dtaumin, trace_time, ORBIT_EXIT_COMPLETED, &
+        ORBIT_EXIT_LCFS, ORBIT_EXIT_WALL
     implicit none
     private
 
@@ -64,7 +64,8 @@ subroutine read_restart_data()
                 boundary_event_radial_residual(idx) = radial_residual
                 boundary_event_time_width(idx) = time_width
                 particle_done(idx) = exit_code == ORBIT_EXIT_COMPLETED .or. &
-                    exit_code == ORBIT_EXIT_LCFS
+                    exit_code == ORBIT_EXIT_LCFS .or. &
+                    exit_code == ORBIT_EXIT_WALL
             end do
             close (unit)
         end if

src/simple_main.f90

@@ -1118,8 +1118,8 @@ subroutine trace_orbit_spectre(ipart, z, passing, orbit_traj, orbit_times)
         !> confined for the full trace, reflects at forbidden interfaces, or is
         !> lost at the outermost interface. Every crossing/reflection is logged.
         use spectre_orbit, only: spectre_orbit_state_t, spectre_event_t, &
-                                 spectre_state_reset, orbit_timestep_spectre, &
-                                 SPECTRE_OK, SPECTRE_BOUNDARY
+            spectre_state_reset, orbit_timestep_spectre, &
+            SPECTRE_OK, SPECTRE_BOUNDARY
         use magfie_sub, only: spectre_field
         use interface_crossing, only: crossing_log_record
         use params, only: crossing_level
@@ -1133,7 +1133,7 @@ subroutine trace_orbit_spectre(ipart, z, passing, orbit_traj, orbit_times)
 
         type(spectre_orbit_state_t) :: state
         type(spectre_event_t) :: event
-        integer :: it, ktau, ierr_orbit, it_final
+        integer :: it, ktau, ierr_orbit, it_final, it_f
         integer(8) :: kt
         real(dp) :: t_stop
 
@@ -1158,6 +1158,15 @@ subroutine trace_orbit_spectre(ipart, z, passing, orbit_traj, orbit_times)
 
             if (ierr_orbit /= SPECTRE_OK) then
                 it_final = it
+                if (ierr_orbit == SPECTRE_BOUNDARY) then
+                    orbit_exit_code(ipart) = ORBIT_EXIT_LCFS
+                else
+                    orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_EVENT
+                    do it_f = it, ntimstep
+                    !$omp atomic update
+                        unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1
+                    end do
+                end if
                 exit
             end if
 
@@ -1176,8 +1185,11 @@ subroutine trace_orbit_spectre(ipart, z, passing, orbit_traj, orbit_times)
 
         if (ierr_orbit == SPECTRE_BOUNDARY) then
             t_stop = (real(kt, dp) + event%t_frac)*dtaumin/v0
+        else if (ierr_orbit /= SPECTRE_OK) then
+            t_stop = ieee_value(0.0_dp, ieee_quiet_nan)
         else
             t_stop = real(kt, dp)*dtaumin/v0
+            orbit_exit_code(ipart) = ORBIT_EXIT_COMPLETED
         end if
 
 !$omp critical
@@ -1196,8 +1208,9 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, &
         !> interface, matching the RK45 path; CROSS_STOP remains only as the
         !> pathological non-convergence fallback inside the microstepper.
         use spectre_sympl_orbit, only: sympl_spectre_state_t, sympl_spectre_reset, &
-                                       orbit_microstep_sympl_spectre, &
-                                       SYMPL_SPECTRE_OK, SYMPL_SPECTRE_SKIM
+            orbit_microstep_sympl_spectre, &
+            SYMPL_SPECTRE_OK, SYMPL_SPECTRE_LOSS, &
+            SYMPL_SPECTRE_STOP, SYMPL_SPECTRE_SKIM
         use field_can_spectre, only: spectre_mvol, set_spectre_volume_lock
         use params, only: crossing_level
         use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
@@ -1210,7 +1223,7 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, &
         real(dp), intent(out) :: orbit_times(:)
 
         type(sympl_spectre_state_t) :: state
-        integer :: it, ktau, ierr_orbit, it_final
+        integer :: it, ktau, ierr_orbit, it_final, it_f
         integer(8) :: kt
         real(dp) :: t_stop, t_frac
 
@@ -1237,6 +1250,24 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, &
 
             if (ierr_orbit /= SYMPL_SPECTRE_OK) then
                 it_final = it
+                select case (ierr_orbit)
+                case (SYMPL_SPECTRE_LOSS)
+                    orbit_exit_code(ipart) = ORBIT_EXIT_LCFS
+                case (SYMPL_SPECTRE_STOP)
+                    orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_FULL_ORBIT
+                    do it_f = it, ntimstep
+                        !$omp atomic update
+                        unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1
+                    end do
+                case (SYMPL_SPECTRE_SKIM)
+                    orbit_exit_code(ipart) = ORBIT_EXIT_COMPLETED
+                case default
+                    orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_EVENT
+                    do it_f = it, ntimstep
+                        !$omp atomic update
+                        unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1
+                    end do
+                end select
                 exit
             end if
 
@@ -1269,10 +1300,13 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, &
 
         if (ierr_orbit == SYMPL_SPECTRE_OK) then
             t_stop = real(kt, dp)*dtaumin/v0
+            orbit_exit_code(ipart) = ORBIT_EXIT_COMPLETED
         else if (ierr_orbit == SYMPL_SPECTRE_SKIM) then
             ! Mirror-confined at an interior interface: cannot be lost, so record
             ! it as confined (times_lost = trace_time) rather than at its stop.
             t_stop = trace_time
+        else if (ierr_orbit /= SYMPL_SPECTRE_LOSS) then
+            t_stop = ieee_value(0.0_dp, ieee_quiet_nan)
         else
             t_stop = (real(kt, dp) + t_frac)*dtaumin/v0
         end if
@@ -1642,12 +1676,14 @@ subroutine write_results
         !> Write the per-particle and confined-fraction result files from the
         !> shared result arrays. progress_monitor calls this every
         !> checkpoint_interval seconds, so a run killed mid-flight keeps its
-        !> last flushed output. Confined fractions are normalised by ntestpart,
-        !> as in the final write, so a partial file is a converging lower bound
-        !> and never exceeds one; particles not yet finished keep their sentinel
-        !> values (times_lost = -1).
-        integer :: i, num_lost, unit
-        real(dp) :: inverse_times_lost_sum, norm
+        !> last flushed output. Confined fractions are conditional on the
+        !> numerically resolved population at each time. A partial file remains
+        !> a lower bound because unfinished particles have not contributed to
+        !> the numerator. Particles not yet finished keep times_lost = -1.
+        use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
+
+        integer :: i, num_lost, resolved_count, unit
+        real(dp) :: inverse_times_lost_sum, norm, pass_fraction, trap_fraction
 
         norm = real(max(ntestpart, 1), dp)
 
@@ -1660,6 +1696,9 @@ subroutine write_results
         num_lost = 0
         inverse_times_lost_sum = 0.0d0
         do i = 1, ntestpart
+            if (orbit_exit_code(i) >= ORBIT_EXIT_NUMERICAL_DOMAIN) then
+                times_lost(i) = ieee_value(0.0_dp, ieee_quiet_nan)
+            end if
             write (unit, *) i, times_lost(i), trap_par(i), zstart(1, i), &
                 perp_inv(i), zend(:, i)
             if (orbit_exit_code(i) == ORBIT_EXIT_LCFS .or. &
@@ -1693,8 +1732,16 @@ subroutine write_results
 
         open (newunit=unit, file='confined_fraction.dat', recl=1024)
         do i = 1, ntimstep
-            write (unit, *) dble(kt_macro(i))*dtaumin/v0, confpart_pass(i)/norm, &
-                confpart_trap(i)/norm, ntestpart
+            resolved_count = ntestpart - unresolved_orbits(i)
+            if (resolved_count > 0) then
+                pass_fraction = confpart_pass(i)/real(resolved_count, dp)
+                trap_fraction = confpart_trap(i)/real(resolved_count, dp)
+            else
+                pass_fraction = ieee_value(0.0_dp, ieee_quiet_nan)
+                trap_fraction = ieee_value(0.0_dp, ieee_quiet_nan)
+            end if
+            write (unit, *) dble(kt_macro(i))*dtaumin/v0, pass_fraction, &
+                trap_fraction, resolved_count
         end do
         close (unit)
 

test/tests/CMakeLists.txt

@@ -821,6 +821,15 @@ set_tests_properties(test_restart_io PROPERTIES
     LABELS "unit"
     TIMEOUT 15)
 
+add_executable(test_result_accounting.x test_result_accounting.f90)
+target_link_libraries(test_result_accounting.x simple)
+add_test(NAME test_result_accounting
+    COMMAND test_result_accounting.x
+    WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_tests_properties(test_result_accounting PROPERTIES
+    LABELS "unit"
+    TIMEOUT 15)
+
 # Restart integration test (TEST field, deterministic)
 add_test(NAME test_restart_integration
     COMMAND ${Python3_EXECUTABLE}

test/tests/test_result_accounting.f90

@@ -0,0 +1,135 @@
+program test_result_accounting
+    use, intrinsic :: ieee_arithmetic, only: ieee_is_nan
+    use, intrinsic :: iso_fortran_env, only: dp => real64, int64
+    use magfie_sub, only: TEST
+    use params, only: ntestpart, ntimstep, kt_macro, dtaumin, v0, &
+        confpart_pass, confpart_trap, unresolved_orbits, &
+        times_lost, orbit_exit_code, trap_par, perp_inv, &
+        zstart, zend, boundary_event_radial_residual, &
+        boundary_event_time_width, isw_field_type, class_plot, &
+        ntcut, ORBIT_EXIT_COMPLETED, ORBIT_EXIT_LCFS, &
+        ORBIT_EXIT_NUMERICAL_MAXITER
+    use simple_main, only: write_results
+    implicit none
+
+    integer :: nerr
+
+    nerr = 0
+    call setup_results()
+    call test_resolved_denominator()
+    call test_zero_resolved()
+    call teardown_results()
+
+    if (nerr > 0) error stop 1
+    print *, 'All result accounting tests passed!'
+
+contains
+
+    subroutine setup_results()
+        ntestpart = 3
+        ntimstep = 1
+        dtaumin = 1.0_dp
+        v0 = 1.0_dp
+        isw_field_type = TEST
+        class_plot = .false.
+        ntcut = 0
+
+        allocate (kt_macro(1), confpart_pass(1), confpart_trap(1))
+        allocate (unresolved_orbits(1), times_lost(3), orbit_exit_code(3))
+        allocate (trap_par(3), perp_inv(3), zstart(5, 3), zend(5, 3))
+        allocate (boundary_event_radial_residual(3))
+        allocate (boundary_event_time_width(3))
+
+        kt_macro = [10_int64]
+        trap_par = 0.0_dp
+        perp_inv = 0.0_dp
+        zstart = 0.0_dp
+        zend = 0.0_dp
+        boundary_event_radial_residual = -1.0_dp
+        boundary_event_time_width = -1.0_dp
+    end subroutine setup_results
+
+    subroutine test_resolved_denominator()
+        integer :: idx, unit
+        real(dp) :: time, pass_fraction, trap_fraction, loss_time
+        real(dp) :: unresolved_fraction
+        integer :: resolved_count, total_count
+
+        confpart_pass = 1.0_dp
+        confpart_trap = 0.0_dp
+        unresolved_orbits = 1
+        times_lost = [10.0_dp, 5.0_dp, 5.0_dp]
+        orbit_exit_code = [ORBIT_EXIT_COMPLETED, ORBIT_EXIT_LCFS, &
+            ORBIT_EXIT_NUMERICAL_MAXITER]
+
+        call write_results()
+
+        open (newunit=unit, file='confined_fraction.dat', status='old')
+        read (unit, *) time, pass_fraction, trap_fraction, resolved_count
+        close (unit)
+        call assert_close(pass_fraction, 0.5_dp, 'resolved confined fraction')
+        call assert_close(trap_fraction, 0.0_dp, 'resolved trapped fraction')
+        call assert_true(resolved_count == 2, 'resolved denominator is two')
+
+        open (newunit=unit, file='times_lost.dat', status='old')
+        read (unit, *) idx, loss_time
+        read (unit, *) idx, loss_time
+        read (unit, *) idx, loss_time
+        close (unit)
+        call assert_true(idx == 3, 'read numerical particle row')
+        call assert_true(ieee_is_nan(loss_time), 'numerical loss time is NaN')
+
+        open (newunit=unit, file='unresolved_fraction.dat', status='old')
+        read (unit, *) time, unresolved_fraction, total_count
+        close (unit)
+        call assert_close(unresolved_fraction, 1.0_dp/3.0_dp, &
+            'unresolved fraction keeps total denominator')
+        call assert_true(total_count == 3, 'unresolved denominator is total')
+    end subroutine test_resolved_denominator
+
+    subroutine test_zero_resolved()
+        integer :: unit, resolved_count
+        real(dp) :: time, pass_fraction, trap_fraction
+
+        confpart_pass = 0.0_dp
+        confpart_trap = 0.0_dp
+        unresolved_orbits = ntestpart
+
+        call write_results()
+
+        open (newunit=unit, file='confined_fraction.dat', status='old')
+        read (unit, *) time, pass_fraction, trap_fraction, resolved_count
+        close (unit)
+        call assert_true(resolved_count == 0, 'zero resolved denominator')
+        call assert_true(ieee_is_nan(pass_fraction), 'zero denominator pass NaN')
+        call assert_true(ieee_is_nan(trap_fraction), 'zero denominator trap NaN')
+    end subroutine test_zero_resolved
+
+    subroutine teardown_results()
+        deallocate (kt_macro, confpart_pass, confpart_trap, unresolved_orbits)
+        deallocate (times_lost, orbit_exit_code, trap_par, perp_inv)
+        deallocate (zstart, zend, boundary_event_radial_residual)
+        deallocate (boundary_event_time_width)
+    end subroutine teardown_results
+
+    subroutine assert_close(actual, expected, message)
+        real(dp), intent(in) :: actual, expected
+        character(*), intent(in) :: message
+
+        if (abs(actual - expected) > 1.0e-12_dp) then
+            print *, 'FAIL:', message, 'expected', expected, 'got', actual
+            nerr = nerr + 1
+        end if
+    end subroutine assert_close
+
+    subroutine assert_true(condition, message)
+        logical, intent(in) :: condition
+        character(*), intent(in) :: message
+
+        if (.not. condition) then
+            print *, 'FAIL:', message
+            nerr = nerr + 1
+        end if
+    end subroutine assert_true
+
+end program test_result_accounting