simulation.py — RadCluster_2_1 simulation orchestrator.
Responsibilities
----------------
1. Load inputs and initialise physics modules.
2. Dispatch to the appropriate solver mode and physics option.
3. Write timestamped output directory with provenance, CSV, and plots.
Solver modes
------------
full_system → C++ SUNDIALS CVODE BDF, full system, via cpp_bridge.run_cpp_solver
active_window → C++ SUNDIALS CVODE BDF with two independent sliding windows
(SIA + VAC) and OpenMP-parallel RHS. Thread count is auto-set
from N_eq; falls back to single-threaded sliding when OpenMP
is unavailable or only one thread is selected.
Legacy aliases (still accepted): cpp_full → full_system,
sliding_OpenMP → active_window.
Physics options (two orthogonal axes)
-------------------------------------
equations ∈ {'discrete', 'bin_moment'} : per-size ODEs vs. bin-moment grouping
cascade ∈ {'fission', 'fusion'} : He coupling case + cascade spectrum
Combinations map to the canonical physics_option string:
discrete + fission → 'full_CD_fission' (he_mode='case2', Eq. 175)
discrete + fusion → 'full_CD_fusion' (he_mode='case1', Eq. 174)
bin_moment + fission → 'bin_moment_CD_fission' (he_mode='case2')
bin_moment + fusion → 'bin_moment_CD_fusion' (he_mode='case1')
Either form is accepted by RadClusterSimulation. Legacy equations='full_CD'
is silently aliased to 'discrete'.
| py_utils.simulation._remap_bin_moments |
( |
|
mu0_old, |
|
|
|
mu1_old, |
|
|
|
bins_old, |
|
|
|
bins_new |
|
) |
| |
|
protected |
Move (μ₀, μ₁) from one bin partition to another, conserving both sums.
Used when a domain doubling re-derives the geometric bin edges and the
moments have to be expressed on the new partition.
WHY NOT reconstruct-then-reproject. The obvious route — expand the moments
to a per-size c_n with the active shape function, then re-take moments on
the new bins — is NOT conservative. ``distribution_from_moments_hat``
clamps its negative lobe with ``np.maximum(vals, 0.0)``, and that clip adds
mass: for a bin whose mean sits away from the centre the reconstructed sum
exceeds the bin's own μ₀, without bound. Measured on the ⟨100⟩ block of the
I=40000 run (2026-08-15): Σ(reconstructed) / Σμ₀ = 176.9, which is precisely
the ×177 jump in N_100 across that run's doubling, and it took δ_FP from
0.20 to 0.96. The ½⟨111⟩ and vacancy blocks survive the same round trip
only because they are broad and well-centred; the ⟨100⟩ population is
sessile, one-way and sharply peaked, so it sits in the worst case
permanently.
THIS route never builds a per-size array. Each old bin is treated as
uniformly filled and split across the new bins it overlaps, by overlap
fraction w (Σ_K w = 1 because the new bins tile the old bin's span). Both
moments ride the same weights, so Σμ₀ and Σμ₁ are conserved to round-off
by construction — the clip cannot enter because no reconstruction happens.
The one approximation is the within-bin placement of content: splitting μ₁
by width smears a peak across the pieces of a straddled bin. A final clamp
puts each new bin's mean back inside its own range, mirroring the
consistency rule in ``post_process._floor_bin_moments`` (an empty bin
carries no content; μ₁/μ₀ lies in the bin). Both are second-order next to
a factor of 177.
Returns (mu0_new, mu1_new), each length len(bins_new).
Definition at line 57 of file simulation.py.
