Use one NetCDF/HDF5 link stack

DOC/config.md

@@ -15,11 +15,17 @@
   reaches its iteration limit, SIMPLE continues only if the final Newton
   correction is finite and no more than ten times the requested relative
   tolerance. Each occurrence is reported by the corresponding `*_maxit`
-  diagnostic. Set it to `.false.` to stop the affected orbit at its last
-  accepted position. Larger corrections, exterior-domain states, singular
-  linear systems, non-finite values, and unresolved boundary events always
-  stop the orbit. Numerical stops are recorded in `orbit_exit_code` and as
-  `NaN` in `times_lost`; they are unresolved markers, not physical losses.
+  diagnostic. The production RK, symplectic, and full-orbit paths then use the
+  same terminal convention: if bounded recovery cannot resolve a numerical
+  microstep, SIMPLE retains any contiguous accepted prefix and holds only the
+  unresolved remainder, records `warning_step_skip`, advances the clock for
+  the complete interval, and retries from that valid state on the next
+  microstep. A warning hold does not terminate or
+  numerically disqualify the marker: a marker that reaches the requested end
+  time remains a resolved survivor, while a later physical boundary event is
+  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`.
 
 * `canonical_grid_nr`, `canonical_grid_ntheta`, and `canonical_grid_nphi`
   control the Meiss or Albert canonical-map grid. Their defaults are 62, 63,

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
@@ -112,6 +112,9 @@ class LossData:
         Confined trapped fraction at each time.
     trace_time : float
         Maximum tracing time [s].
+    orbit_exit_codes : np.ndarray, optional
+        Per-marker exit codes. Codes 0, 1, and 2 are resolved completion,
+        LCFS loss, and wall loss; codes 101--105 are numerical terminations.
     """
 
     n_particles: int
@@ -127,6 +130,7 @@ class LossData:
     confined_pass: np.ndarray
     confined_trap: np.ndarray
     trace_time: float
+    orbit_exit_codes: Optional[np.ndarray] = None
 
     @property
     def confined_fraction(self) -> np.ndarray:
@@ -135,19 +139,42 @@ def confined_fraction(self) -> np.ndarray:
 
     @property
     def lost_mask(self) -> np.ndarray:
-        """Boolean mask for particles that were lost (not confined)."""
+        """Boolean mask for resolved physical LCFS or wall losses."""
+        if self.orbit_exit_codes is not None:
+            return np.isin(self.orbit_exit_codes, (1, 2))
         return (self.loss_times > 0) & (self.loss_times < self.trace_time)
 
     @property
     def confined_mask(self) -> np.ndarray:
         """Boolean mask for particles that remained confined."""
+        if self.orbit_exit_codes is not None:
+            # Code 3 is SIMPLE's deliberate deep-passing shortcut. Fortran
+            # counts those analytically confined markers in the same resolved
+            # denominator, while skipped_mask keeps the shortcut observable.
+            return np.isin(self.orbit_exit_codes, (0, 3))
         return self.loss_times >= self.trace_time
 
     @property
     def skipped_mask(self) -> np.ndarray:
         """Boolean mask for particles skipped (deep passing, contr_pp)."""
+        if self.orbit_exit_codes is not None:
+            return self.orbit_exit_codes == 3
         return self.loss_times < 0
 
+    @property
+    def resolved_mask(self) -> np.ndarray:
+        """Markers eligible for confinement and loss statistics."""
+        if self.orbit_exit_codes is not None:
+            return np.isin(self.orbit_exit_codes, (0, 1, 2, 3))
+        return self.lost_mask | self.confined_mask
+
+    @property
+    def unresolved_mask(self) -> np.ndarray:
+        """Markers terminated for numerical reasons."""
+        if self.orbit_exit_codes is not None:
+            return (self.orbit_exit_codes >= 101) & (self.orbit_exit_codes <= 105)
+        return ~np.isfinite(self.loss_times)
+
 
 def load_loss_data(directory: str | Path) -> LossData:
     """
@@ -209,28 +236,52 @@ def load_loss_data(directory: str | Path) -> LossData:
         start_phi = np.zeros(n_particles)
         start_pitch = np.zeros(n_particles)
 
-    # Load confined_fraction.dat
-    # Columns: time, confpart_pass, confpart_trap, ntestpart
+    # Load confined_fraction.dat. Its final time is authoritative: numerical
+    # NaNs in times_lost.dat must not poison or extend the requested trace time.
+    # 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)
         # Skip first row (t=0, all confined)
         time_grid = np.abs(conf_data[1:, 0])
         confined_pass = conf_data[1:, 1]
         confined_trap = conf_data[1:, 2]
+        trace_time = float(time_grid[-1]) if len(time_grid) > 0 else 1.0
+    else:
+        finite_times = np.abs(loss_times[np.isfinite(loss_times)])
+        trace_time = float(np.max(finite_times)) if finite_times.size else 1.0
+        time_grid = np.logspace(-5, np.log10(max(trace_time, 1.0e-5)), 100)
+
+    exit_file = directory / "orbit_exit_code.dat"
+    if exit_file.exists():
+        exit_data = np.loadtxt(exit_file)
+        exit_data = np.atleast_2d(exit_data)
+        if exit_data.shape[0] != n_particles:
+            raise ValueError("orbit_exit_code.dat particle count does not match times_lost.dat")
+        orbit_exit_codes = exit_data[:, 1].astype(int)
     else:
-        # Generate from loss times if file not available
-        time_grid = np.logspace(-5, 0, 100)
+        # Backward compatibility for output predating explicit exit codes.
+        orbit_exit_codes = np.full(n_particles, 101, dtype=int)
+        orbit_exit_codes[loss_times < 0.0] = 3
+        orbit_exit_codes[np.isfinite(loss_times) & (loss_times >= trace_time)] = 0
+        orbit_exit_codes[(loss_times > 0.0) & (loss_times < trace_time)] = 1
+
+    if not conf_file.exists():
         confined_pass = np.zeros_like(time_grid)
         confined_trap = np.zeros_like(time_grid)
+        resolved = np.isin(orbit_exit_codes, (0, 1, 2, 3))
+        passing = trap_parameter < 0
+        denominator = int(np.sum(resolved))
         for i, t in enumerate(time_grid):
-            confined = np.abs(loss_times) >= t
-            passing = trap_parameter < 0
-            confined_pass[i] = np.sum(confined & passing) / n_particles
-            confined_trap[i] = np.sum(confined & ~passing) / n_particles
-
-    # Determine trace_time from maximum loss time or last time in grid
-    trace_time = max(np.max(np.abs(loss_times)), time_grid[-1] if len(time_grid) > 0 else 1.0)
+            confined = resolved & (
+                np.isin(orbit_exit_codes, (0, 3)) | (loss_times >= t)
+            )
+            if denominator > 0:
+                confined_pass[i] = np.sum(confined & passing) / denominator
+                confined_trap[i] = np.sum(confined & ~passing) / denominator
+            else:
+                confined_pass[i] = np.nan
+                confined_trap[i] = np.nan
 
     return LossData(
         n_particles=n_particles,
@@ -246,6 +297,7 @@ def load_loss_data(directory: str | Path) -> LossData:
         confined_pass=confined_pass,
         confined_trap=confined_trap,
         trace_time=trace_time,
+        orbit_exit_codes=orbit_exit_codes,
     )
 
 
@@ -281,7 +333,9 @@ def compute_energy_confined_fraction(
     if time_grid is None:
         time_grid = data.time_grid
 
-    e_total = float(data.n_particles)
+    e_total = float(np.sum(data.resolved_mask))
+    if e_total == 0.0:
+        return time_grid, np.full_like(time_grid, np.nan, dtype=float)
 
     if slowing_down_curve is not None:
         times_sd, energy_sd = slowing_down_curve
@@ -295,7 +349,7 @@ def compute_energy_confined_fraction(
 
     energy_fraction = np.zeros_like(time_grid)
     for i, t in enumerate(time_grid):
-        lost_before_t = (data.loss_times > 0) & (data.loss_times < t)
+        lost_before_t = data.lost_mask & (data.loss_times < t)
         e_lost = np.sum(energy_at_loss[lost_before_t])
         energy_fraction[i] = (e_total - e_lost) / e_total
 
@@ -412,7 +466,7 @@ def _plot_kde_density_panel(ax, data: LossData) -> None:
 def _plot_energy_loss_panel(ax, data: LossData) -> None:
     """Plot energy loss distribution vs J_perp panel."""
     jperp_bins, particle_count, energy_lost = compute_energy_loss_distribution(data)
-    total = data.n_particles if data.n_particles > 0 else 1
+    total = max(int(np.sum(data.resolved_mask)), 1)
     energy_frac = energy_lost / total
     particle_frac = particle_count / total
 
@@ -426,7 +480,10 @@ def _plot_energy_loss_panel(ax, data: LossData) -> None:
 
 def _plot_starting_positions_panel(ax, data: LossData) -> None:
     """Plot starting positions colored by loss time panel."""
-    tlost_all = np.maximum(np.abs(data.loss_times), 1e-10)
+    tlost_all = np.where(
+        np.isfinite(data.loss_times), np.abs(data.loss_times), data.trace_time
+    )
+    tlost_all = np.maximum(tlost_all, 1e-10)
     scatter = ax.scatter(data.start_theta, data.start_pitch, c=np.log10(tlost_all), s=1, cmap="viridis", vmin=-5, vmax=0)
     cbar = plt.colorbar(scatter, ax=ax)
     cbar.set_label(r"$\log_{10}(t_{\mathrm{loss}})$")
@@ -650,7 +707,8 @@ def _bin_energy_actual(
 ) -> Tuple[np.ndarray, np.ndarray]:
     """Bin particles and energy by J_perp using actual final_p^2 (vectorized)."""
     k = _compute_bin_indices(data.perp_invariant, hp, nperp)
-    part_distr = np.bincount(k, minlength=nperp).astype(float)
+    part_distr = np.bincount(k, weights=data.resolved_mask.astype(float),
+                            minlength=nperp)
     energy_weights = np.where(data.lost_mask, data.final_p**2, 0.0)
     energ_distr = np.bincount(k, weights=energy_weights, minlength=nperp)
     return part_distr, energ_distr
@@ -665,7 +723,8 @@ def _bin_energy_theoretical(
 ) -> Tuple[np.ndarray, np.ndarray]:
     """Bin particles and energy by J_perp using theoretical slowing-down (vectorized)."""
     k = _compute_bin_indices(data.perp_invariant, hp, nperp)
-    part_distr = np.bincount(k, minlength=nperp).astype(float)
+    part_distr = np.bincount(k, weights=data.resolved_mask.astype(float),
+                            minlength=nperp)
 
     energy_adj = energy_sd - THERMAL_ENERGY_FRACTION
     energy_at_loss = np.interp(data.loss_times, times_sd, energy_adj, left=energy_adj[0], right=energy_adj[-1])
@@ -743,13 +802,14 @@ def plot_energy_loss_vs_jperp(
             _, energ_n_theo = _bin_energy_theoretical(data_nocoll, hp, nperp, times_sd, energy_sd)
 
     part_c_safe = np.where(part_c > 0, part_c, 1)
+    part_n_safe = np.where(part_n > 0, part_n, 1) if part_n is not None else None
     jperp = np.arange(1, nperp + 1) / nperp
 
     curves = [
         (energ_c / part_c_safe, "r-", 2, "COLL (actual)"),
-        (energ_n / part_c_safe if energ_n is not None else None, "b-", 2, "NOCOLL (actual)"),
+        (energ_n / part_n_safe if energ_n is not None else None, "b-", 2, "NOCOLL (actual)"),
         (energ_c_theo / part_c_safe if energ_c_theo is not None else None, "r--", 1.5, "COLL (theoretical)"),
-        (energ_n_theo / part_c_safe if energ_n_theo is not None else None, "b--", 1.5, "NOCOLL (theoretical)"),
+        (energ_n_theo / part_n_safe if energ_n_theo is not None else None, "b--", 1.5, "NOCOLL (theoretical)"),
     ]
 
     max_f = 0.0

src/CMakeLists.txt

@@ -121,7 +121,7 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
 endif()
 
 target_link_libraries(simple PUBLIC
-    netcdf netcdff BLAS::BLAS LAPACK::LAPACK
+    BLAS::BLAS LAPACK::LAPACK
 )
 
 target_link_libraries(simple PUBLIC