Struct ThermalizationSolver

pub struct ThermalizationSolver {

    pub config: SolverConfig,
    pub cosmo: Cosmology,
    pub grid: FrequencyGrid,
    pub delta_n: Vec<f64>,
    pub electron_temp: ElectronTemperature,
    pub z: f64,
    pub injection: Option<InjectionScenario>,
    pub snapshots: Vec<SolverSnapshot>,
    pub step_count: usize,
    pub disable_dcbr: bool,
    pub coupled_dcbr: bool,
    pub dcbr_scale: f64,
    pub number_conserving: bool,
    pub nc_stride: usize,
    pub accumulated_delta_t: f64,
    pub diag: SolverDiagnostics,
    /* private fields */
}

Full PDE solver for CMB spectral distortions.

Evolves the photon occupation number perturbation Δn(x, z) from z_start down to z_end using a coupled implicit scheme:

• Crank-Nicolson for Kompaneets (frequency redistribution)
• Backward Euler for DC/BR emission (photon-number changing)
• Newton iteration coupling all three simultaneously

Lifecycle

ThermalizationSolver::new(cosmo, grid_config)   // or ::builder(...)
  → set_injection(scenario)                      // optional
  → run_with_snapshots(&[z1, z2, ...])           // evolve the PDE
  → snapshots[i].{mu, y, delta_t, ...}           // read results

The solver can be re-run after calling reset(), which clears Δn and diagnostic counters but keeps the cosmology and grid.

Energy injection

Call set_injection() before running. Without an injection, Δn stays near zero (only adiabatic cooling acts). Multiple runs with different injections require reset() between runs.

Configuration

See SolverConfig for timestep and physics options, and GridConfig for frequency grid options. The builder API (ThermalizationSolver::builder) provides a fluent interface for all configuration.

Field visibility (unstable API)

Many fields are currently pub for use by examples, tests, and diagnostic tooling. These are not part of the stable public API and may become pub(crate) in a future release. Mutating state fields (delta_n, z, electron_temp, snapshots, step_count, accumulated_delta_t) between run_with_snapshots calls can break Newton-solver invariants; prefer the builder and set_injection / reset methods for all changes of consequence.

Fields

config: SolverConfig

Timestepping and physics-toggle configuration.

cosmo: Cosmology

Background cosmology.

grid: FrequencyGrid

Frequency grid (non-uniform in x).

delta_n: Vec<f64>

Current photon occupation-number perturbation Δn(x) on the grid.

electron_temp: ElectronTemperature

Electron temperature state ρ_e = T_e/T_z (perturbative quasi-stationary solve).

z: f64

Current redshift of the integration.

injection: Option<InjectionScenario>

Active injection scenario, if any. Set via Self::set_injection.

snapshots: Vec<SolverSnapshot>

Snapshots collected during the last run. Populated by Self::run_with_snapshots / Self::run.

step_count: usize

Number of timesteps taken in the current run.

disable_dcbr: bool

If true, disable DC/BR processes (Kompaneets only). For diagnostics.

coupled_dcbr: bool

If true (default), couple DC/BR into the Kompaneets Newton iteration instead of operator splitting. Uses IMEX: Crank-Nicolson for Kompaneets

• backward Euler for DC/BR, solved simultaneously. More physically consistent than operator splitting, especially at z > 2×10⁶.

dcbr_scale: f64

Diagnostic: scale factor for DC/BR emission rates. Default 1.0. Set to < 1.0 to test whether rates are too strong.

number_conserving: bool

Subtract the temperature shift component from Δn after each DC/BR step at z > 5×10⁴, enforcing photon number conservation (∫x² Δn dx = 0). This prevents DC/BR-created photons from accumulating as a spurious temperature shift that feeds back into over-thermalization. See Chluba (2013), arXiv:1304.6120. Enabled by default.

nc_stride: usize

Apply NC stripping every nc_stride steps (default 1 = every step). Higher values reduce NC-DC/BR feedback at high z.

accumulated_delta_t: f64

Cumulative δT/T subtracted from Δn by number conservation. Added back to ρ_eq so T_e tracks the true (shifted) reference.

diag: SolverDiagnostics

Diagnostic counters and extrema tracked during evolution.

Implementations

impl ThermalizationSolver

pub fn new(cosmo: Cosmology, grid_config: GridConfig) -> Self

Construct a solver with the given cosmology and frequency grid.

All other state (solver config, injection, flags) is set to defaults; mutate the public fields or call Self::set_injection / Self::set_config before running. For a fluent API that performs validation at build time, use Self::builder instead.

pub fn set_injection( &mut self, scenario: InjectionScenario, ) -> Result<(), String>

Attach an energy-injection scenario, validating it first.

Returns Err if the scenario parameters are unphysical (e.g. negative widths, impossible masses). Collects stimulated-emission warnings into SolverDiagnostics::warnings.

pub fn set_config(&mut self, config: SolverConfig)

Replace the solver configuration and reset the current redshift to the new z_start.

pub fn set_initial_delta_n(&mut self, delta_n: Vec<f64>)

Set an initial photon perturbation Δn(x) for the next PDE run.

This replaces the default Δn = 0 initialization in run_with_snapshots. The perturbation is consumed (taken) on the next call to run_with_snapshots.

Panics if any entry is non-finite — silently passing NaN/Inf into the solver causes a deep panic in the Newton step where the source isn’t obvious. Caught early here so the user sees the real culprit.

pub fn reset(&mut self)

Reset solver state for reuse, keeping grid and recombination cache.

Restores all configuration flags to their defaults: a reset solver behaves identically to a freshly constructed one (aside from the cached grid/recombination tables, which are preserved for speed).

pub fn step(&mut self) -> f64

Advance the solver by a single adaptively-chosen timestep.

Returns the dz taken. Most users should call Self::run or Self::run_with_snapshots instead of stepping manually.

pub fn run_with_snapshots( &mut self, snapshot_redshifts: &[f64], ) -> &[SolverSnapshot]

Integrate from z_start to z_end, recording a snapshot at each requested redshift.

snapshot_redshifts may be given in any order; they are sorted internally. Redshifts above z_start or below z_end are clamped to the initial and final state respectively. The returned slice borrows from self.snapshots; for an owned result, use Self::run_to_result.

pub fn run(&mut self, n_snapshots: usize) -> &[SolverSnapshot]

Run the solver with n_snapshots log-spaced snapshot redshifts between z_start and z_end. Note: with n_snapshots=1 the single snapshot is at z_start (the first log-spaced point), not z_end.

pub fn extract_mu_y_joint(&self) -> (f64, f64)

Extract (μ, y) from the current Δn(x) via the default joint least-squares decomposition (B&F 2022 nonlinear BE; see crate::distortion::decompose_distortion).

pub fn run_to_result(&mut self, z_obs: f64) -> SolverResult

Run the solver and return an owned crate::output::SolverResult with a single snapshot at z_obs.

This is the preferred entry point: the result does not borrow the solver and has built-in JSON/CSV/table serialization. For runs that need multiple intermediate snapshots, call Self::run_with_snapshots directly and inspect Self::snapshots.

pub fn builder(cosmo: Cosmology) -> SolverBuilder

Create a builder for configuring a solver with a fluent API.

Example

use spectroxide::prelude::*;

let mut solver = ThermalizationSolver::builder(Cosmology::planck2018())
    .grid(GridConfig::production())
    .injection(InjectionScenario::SingleBurst {
        z_h: 2e5, delta_rho_over_rho: 1e-5, sigma_z: 100.0,
    })
    .z_range(5e5, 1e3)
    .build()
    .unwrap();