diff --git a/README.md b/README.md index 38acdd94..3f79c6e3 100644 --- a/README.md +++ b/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. diff --git a/python/pysimple/plotting.py b/python/pysimple/plotting.py index 67c20462..63a7c71b 100644 --- a/python/pysimple/plotting.py +++ b/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) diff --git a/src/restart.f90 b/src/restart.f90 index 1f453207..6f2d52f6 100644 --- a/src/restart.f90 +++ b/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 diff --git a/src/simple_main.f90 b/src/simple_main.f90 index 5b4e78a8..68b6d1c0 100644 --- a/src/simple_main.f90 +++ b/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) diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt index df647159..85023095 100644 --- a/test/tests/CMakeLists.txt +++ b/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} diff --git a/test/tests/test_result_accounting.f90 b/test/tests/test_result_accounting.f90 new file mode 100644 index 00000000..66c86459 --- /dev/null +++ b/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 diff --git a/test/tests/test_spectre_validation.py b/test/tests/test_spectre_validation.py index 8669058c..244b68b5 100644 --- a/test/tests/test_spectre_validation.py +++ b/test/tests/test_spectre_validation.py @@ -91,9 +91,11 @@ def run(binary, workdir, h5, **kwargs): def split_confined(times_lost, trace_time): - confined = times_lost >= trace_time * (1.0 - 1e-9) - terminated = (times_lost > 0.0) & (times_lost < trace_time * (1.0 - 1e-9)) - unresolved = times_lost <= 0.0 + finite = np.isfinite(times_lost) + confined = finite & (times_lost >= trace_time * (1.0 - 1e-9)) + terminated = (finite & (times_lost > 0.0) + & (times_lost < trace_time * (1.0 - 1e-9))) + unresolved = ~finite return confined, terminated, unresolved @@ -104,40 +106,42 @@ def losses_and_accounting(binary, h5, failures): run(binary, work, h5, crossing_level=1) tl = np.loadtxt(os.path.join(work, "times_lost.dat")) cf = np.loadtxt(os.path.join(work, "confined_fraction.dat")) + exit_codes = np.loadtxt(os.path.join(work, "orbit_exit_code.dat")) ev_path = os.path.join(work, "spectre_crossing_events.dat") crossings_written = os.path.exists(ev_path) and os.path.getsize(ev_path) > 0 if len(tl) != NPART: failures.append(f"accounting: times_lost.dat has {len(tl)} rows != {NPART}") - if not np.all(np.isfinite(tl[:, 1])): - failures.append("accounting: non-finite loss time (unresolved marker)") - confined, terminated, unresolved = split_confined(tl[:, 1], TRACE_TIME) n_conf, n_term, n_unres = (int(confined.sum()), int(terminated.sum()), int(unresolved.sum())) - if n_unres != 0: - failures.append(f"accounting: {n_unres} markers left unresolved " - f"(times_lost <= 0)") - if n_conf + n_term != NPART: - failures.append(f"accounting: {n_conf} confined + {n_term} terminated " + if n_conf + n_term + n_unres != NPART: + failures.append(f"accounting: {n_conf} confined + {n_term} terminated + " + f"{n_unres} unresolved " f"!= {NPART} markers") + if not np.all(exit_codes[unresolved, 1] >= 101): + failures.append("accounting: unresolved marker lacks numerical exit code") + if np.any(exit_codes[confined | terminated, 1] >= 101): + failures.append("accounting: resolved marker has numerical exit code") # confined_fraction.dat is filled by an independent per-timestep confined # counter; its final value must equal the times_lost confined count. Breaking # the loss tally (counting a lost marker as confined) makes these disagree. + n_resolved = n_conf + n_term frac_final = cf[-1, 1] + cf[-1, 2] - if abs(frac_final - n_conf / NPART) > 0.5 / NPART + 1e-9: + if abs(frac_final - n_conf / n_resolved) > 0.5 / n_resolved + 1e-9: failures.append(f"accounting: confined_fraction {frac_final:.4f} != " - f"times_lost confined {n_conf / NPART:.4f}") - if int(cf[-1, 3]) != NPART: + f"resolved times_lost confined " + f"{n_conf / n_resolved:.4f}") + if int(cf[-1, 3]) != n_resolved: failures.append(f"accounting: confined_fraction N = {int(cf[-1, 3])}") if not crossings_written: failures.append("accounting: no spectre_crossing_events.dat produced") print(f"accounting: confined={n_conf} terminated={n_term} unresolved=" - f"{n_unres} (sum={n_conf + n_term}={NPART}) " - f"cf_final={frac_final:.3f}=({n_conf}/{NPART})") - return n_term / NPART + f"{n_unres} (sum={n_conf + n_term + n_unres}={NPART}) " + f"cf_final={frac_final:.3f}=({n_conf}/{n_resolved})") + return n_term / n_resolved def crossing_map_accounting(binary, h5, p1, failures): @@ -148,7 +152,7 @@ def crossing_map_accounting(binary, h5, p1, failures): confined0, terminated0, unresolved0 = split_confined(tl0[:, 1], TRACE_TIME) n_conf0, n_term0, n_unres0 = (int(confined0.sum()), int(terminated0.sum()), int(unresolved0.sum())) - if n_conf0 + n_term0 != NPART or n_unres0 != 0: + if n_conf0 + n_term0 + n_unres0 != NPART: failures.append(f"crossing maps: Level-0 account is {n_conf0} confined + " f"{n_term0} terminated + {n_unres0} unresolved") @@ -159,7 +163,7 @@ def crossing_map_accounting(binary, h5, p1, failures): failures.append("crossing maps: Level-0 times_lost.dat is not " "bit-identical under a fixed seed") - p0 = n_term0 / NPART + p0 = n_term0 / (n_conf0 + n_term0) print(f"crossing maps: loss_fraction L1={p1:.3f} L0={p0:.3f} " f"Level-0_account={n_conf0}+{n_term0}+{n_unres0}={NPART} " f"Level-0_repeat={'bit-identical' if raw0 == raw0_repeat else 'DIFFERS'}")