Coagulation project documentation

Continuation of the M.Sc. thesis Investigating the Feasibility of utilizing Stochastic Sampling of Colliding Circumstellar Dust Particle Pairs for Computational Cost Optimization of the Smoluchowski Coagulation Equation’s Numerical Integration (V. Mader, Heidelberg University, 2024; supervisor C. P. Dullemond): a Rust re-implementation of the local dust coagulation-fragmentation model, a suite of time integrators up to an own Radau IIA, a systematic study of Monte Carlo kernel sampling, an interactive web visualization, and a joint paper.

Where to start

1. Vision and research goals: goals, research questions, hypotheses, paper pitch.
2. Thesis review: what the thesis established and the findings (F1 to F10) that change the method.
3. Roadmap: milestones M0 to M11 with tasks and acceptance criteria.

Contents

Conventions

• IDs: goals Gn, research questions RQn, hypotheses Hn, findings Fn, milestones Mn and tasks Mn.m, experiments En, questions Qn, risks Xn. IDs are stable; retired items are struck through, not renumbered.
• Status lines: every planning document starts with its status and version; bump the version when content changes materially.
• Decisions: a choice with alternatives becomes an ADR (Proposed until Vincent accepts it); smaller scope decisions go into the decision log.
• Diagrams: Mermaid in Markdown, so they render on GitHub and in most editors and stay diffable.
• Language: English; symbols as in notation.md; SI units.
• Sources: ../related/ (project root) holds symlinks to the thesis (msc-thesis/), its code (msc-thesis.py/), the thesis working folder with papers (Write Master's Thesis./), and the paper inbox (Papers/). They are read-only: nothing there is edited, and what we need is rebuilt in this repo, mostly from scratch in Rust.

Vision and research goals

Status: v0.2 (2026-11-10), scope settled in the grill session (decision log in 06). Iterate in place; record decisions as ADRs.

1. One-paragraph summary

Solving the Smoluchowski coagulation-fragmentation equation on a grid costs O(N³) per right-hand side once fragmentation is included, and O(N^{3p}) when dust particles carry p properties (mass, filling factor, composition). That cost is why almost all disk-scale dust evolution codes describe grains by mass alone. Dullemond proposed estimating the double sum by Monte Carlo sampling of colliding pairs, while keeping the grid, detailed balance, and implicit time stepping. The M.Sc. thesis (Mader 2024) tested a variant of this and found a qualitatively correct evolution but poor accuracy. Our review shows that the thesis estimator was biased and its proposal distribution was not the one documented (01-thesis-review.md). This project rebuilds the full pipeline in Rust as a clean, tested, and fast research code. It then answers the open question properly: when, and by how much, does unbiased and variance-reduced kernel sampling with implicit time integration beat full summation, for one and for two particle properties? The result is a methods paper with Dullemond, plus an interactive web tool that makes every component of the model visible and switchable.

2. Goals

Non-goals for the paper: global radial transport (1D disk), gas evolution, radiative transfer, GPU kernels. The architecture stays ready for a non-local model (ADR 0015): radial transport and a global disk are the most likely follow-up, so nothing in the local code may assume a single cell.

3. Research questions

4. Working hypotheses

5. What is new compared with prior work

6. Paper pitch (to refine with Dullemond)

Working title options:

1. Unbiased Monte Carlo sampling of the coagulation kernel for implicit dust evolution with multiple particle properties
2. Stochastic kernel summation for the Smoluchowski equation: bias, variance, and implicit time integration

Draft abstract (to be rewritten once results exist):

Grid-based solutions of the coagulation-fragmentation equation scale as N³ per time step and as N^{3p} when grains carry p properties. We study a Monte Carlo estimator of the collision sum that samples colliding pairs, preserves mass conservation for every sample, and remains compatible with implicit time integration. We show that sampling without reweighting, as used in earlier tests, underestimates all collision rates, and we derive the variance-optimal proposal distribution. Combined with a hybrid exact/sampled split and sparse Jacobians in a Radau IIA integrator, the method reproduces full-kernel steady states of growth and fragmentation to X % at Y times lower cost for p = 1 and Z times lower cost for a mass and filling-factor model.

Candidate venues: A&A (where the group publishes, e.g. Beutel and Dullemond 2023), or J. Comput. Phys. if the numerical analysis dominates. Decide with Dullemond.

Planned key figures: (1) model overview diagram; (2) bias vs variance decomposition of the thesis setup; (3) proposal comparison; (4) noise floor vs S and Δt; (5) work-precision diagrams p = 1 and p = 2; (6) break-even surface over (N_m, N_φ, accuracy); (7) steady-state distributions with and without bouncing.

7. Working principles

• Reproduce before improving. Every new method is compared against the faithful baseline, and the baseline against the thesis.
• One switch per physical assumption. Every component (velocity term, outcome model, boundary policy, estimator, integrator) is a config toggle, so its effect can be isolated and visualized.
• Correctness is tested, not assumed. Analytic solutions, parity with coag_py/scipy, conservation laws, and statistical tests of estimators run in CI.
• Decisions are written down. Architecture and method choices go into ADRs with alternatives and consequences.
• Results are reproducible. Every output carries its config hash, git revision, and RNG seed.

8. Settled scope (2026-11-10)

Critical review of the M.Sc. thesis and its code

Purpose: record what the thesis established, what it left open, and what in its method or code must change before the results can go into a paper. Every finding cites the source so it can be checked. Paths are relative to the project root (related/ holds the symlinks).

Sources reviewed in full: thesis TeX (related/msc-thesis/src/chapters/*.tex), thesis Python (related/msc-thesis.py/src), Dullemond’s coag_py (related/msc-thesis.py/lib/coag_py), Dullemond’s notes Dust evolution with binning methods (19 pp.).

1. What the thesis did

flowchart LR
    A[Disk model<br/>irradiated power law,<br/>r0 = 10 au midplane] --> B[Dust kinematics<br/>BR, RD, AZ, TU]
    B --> C[Collision rates<br/>sigma dv]
    C --> D[Outcome probabilities<br/>MB-smoothed v_frag = 1 m/s,<br/>no bouncing]
    D --> E[Kernel K_kij<br/>Kovetz-Olund + cancellation,<br/>MRN fragmentation]
    E --> F[Integration<br/>scipy Radau per macro step]
    E --> G[Sampled kernel<br/>pairs without replacement]
    G --> F
    F --> H[Diagnostics<br/>mass error, accuracy error,<br/>RMS dN/dt, sparsity]

2. What is solid and carries over

1. The kernel construction (Kovetz and Olund plus cancellation handling) matches coag_py.create_coag_kernel and conserves mass per collision to about 10⁻¹². It becomes the reference for M4 parity tests.
2. The physical model is a reasonable, documented “standard model” of local dust evolution at one radius (close to Birnstiel 2010 and DustPy). Good enough as the testbed for a numerical-methods paper.
3. The diagnostics (mass error, RMS dN/dt, sparsity, sampled-pair maps) are the right ones to keep; some need better definitions (F9).
4. The main conclusion (sampling pays off when the particle state has more than one property) is plausible and is exactly the argument Dullemond’s notes build toward. It has not been tested yet. Testing it is the core of the paper.

3. Findings that change the method

Severity: A invalidates a quantitative result, B biases or weakens a result, C presentation or citation.

F1 (A): the sampling estimator is biased low; this explains the “slower evolution”

Thesis ch. 4, “Sampling of the Total Kernel”, Observation 3, attributes the slower approach to equilibrium at low ρ_sample to “effectively lower collision rates”. That is a direct consequence of the estimator, not a property of sampling in general:

• Pairs are drawn without replacement (np.random.choice(..., replace=False), related/msc-thesis.py/src/models/kernel/sampled.py:94-95, config allow_duplicate_sampling = false).
• The sampled kernel is the full kernel restricted to the drawn pairs, with no reweighting (sampled.py:69-70, Kernel(..., ijs=ijs)).
• Hence E[F̂_k] = Σ_ij π_ij K_kij N_i N_j with inclusion probability π_ij ≤ 1. Every term is scaled down by its inclusion probability, so the rate of every process is underestimated. A Horvitz-Thompson weight 1/π_ij, or importance sampling with replacement and weight 1/(S p_ij), removes the bias.
• The with-replacement branch (sampled.py:71-72) divides by P_ij but not by S, and it is applied to a kernel built from the set of unique pairs, so multiplicities are lost. It is not an unbiased estimator either, and the thesis runs did not use it.
• Dullemond’s notes (Eq. 1.70) and his 2021 prototype (coag_py/mc_sampling/coag.py, sample_kernel_sum) already use the unbiased form K/(F_i F_j) with draws with replacement. The thesis diverged from it.

Consequence: the accuracy error in thesis Fig. 106 mixes bias (systematic lag) and variance (noise). Only the variance part is intrinsic to sampling. The headline negative result of the thesis is partly an artifact and should be re-measured with an unbiased estimator. This is the strongest single argument for a follow-up paper.

F2 (A): the documented sampling probability is not what the runs used

The thesis defines P_ij ∝ ΔM^net_ij N_i N_j m_i m_j with ΔM^net_ij = Σ_k m_k |K_kij|. In the code path used for all sampling figures (examples 104, 106, 107 import SolverV2):

• SolverV2.__init__ sets self.W_ij = weights(self.mg, self.kernels[0]) (related/msc-thesis.py/src/models/solver/solver_v2.py:91).
• weights() (solver_v2.py:223-226) calls np.sum([...]) without axis=0, which reduces to a scalar.
• SampledKernel uses this scalar (sampled.py:39-44), so P_ij ∝ N_i N_j m_i m_j. The ΔM^net factor is a constant and drops out after normalization.

So the runs sampled by N_i N_j m_i m_j, i.e. by the product of bin masses ρ_i ρ_j, not by the documented relevance measure. The per-pair throughput (which vanishes for pairs that are cut at the grid edge, and is tiny where Δv is small) never entered. This may explain the “peaks at the upper edge” (Observation 4), since pairs near m_max carry most mass but partially have zero kernel.

F3 (B): the documented proposal is not the variance-optimal one

Even as documented, P_ij ∝ ΔM^net N_i N_j m_i m_j has an extra m_i m_j factor justified only heuristically. For the unbiased IS estimator of the vector F = Σ_ij K_·ij N_i N_j, the proposal that minimizes E‖F̂ − F‖ in a weighted norm ‖·‖w is p_ij ∝ ‖K·ij‖_w N_i N_j (standard result for vector-valued IS; derivation in 04-experiments.md, E7). With the mass-weighted L1 norm this is exactly p_ij ∝ ΔM^net_ij N_i N_j, i.e. the thesis’ relevance measure without m_i m_j. The m_i m_j factor double counts mass and oversamples big-big pairs.

F4 (B): one frozen sample per macro step, with long steps

Each of the 200 macro steps calls scipy Radau once over [t_n, t_{n+1}] with a kernel sampled once at t_n (solver_v2.py:126-155, N_subst = 1). On a log time grid the late macro steps span about 10⁸ s, so a single random draw drives the solution for a long interval and MC noise does not average out within that interval. Noise averaging needs either more, shorter steps, or several independent samples per step (averaged), or resampling per integrator step. This is a design variable to study (E9), not a fixed choice.

F5 (B): integrator tolerances and order

• scipy Radau is the 3-stage Radau IIA method of order 5 (also stated in the react_0d.py docstring). The thesis abstract, ch. 3, and ch. 5 say “4th order”. Fix in any write-up.
• solve_react_0d_equation passes no rtol/atol, so scipy defaults apply (rtol = 10⁻³, atol = 10⁻⁶). N_k is in m⁻³ and spans many decades (rough estimate: about 3·10² m⁻³ for monomers at 10 au, many orders less for the largest bins). An absolute tolerance of 10⁻⁶ m⁻³ means bins below that are not error-controlled. Needs a scaled state or per-component atol (ADR 0003). To verify in M6.
• The professor’s react_0d.py is a wrapper around scipy Radau with an analytic Jacobian plus a Newton-iterated backward Euler. It is not a custom Radau implementation. Our own Radau IIA (ADR 0008) is new work.

F6 (B): molecular viscosity in the radial gas velocity

The thesis (ch. 2, “Gas Viscosity” and “Radial Gas Velocity”) and the code (related/msc-thesis.py/src/models/disk/disk.py:54,72, physics.viscosity(u_th, lambda_mfp)) use the molecular viscosity ν = ½ u_th λ_mfp in the Lynden-Bell and Pringle formula for u_g. Birnstiel et al. (2010) use the turbulent viscosity ν_t = α c_s H_p. ν_mol is smaller than ν_t by the Reynolds number (a rough estimate for this model gives Re ≈ 10⁵ at 10 au), so u_g is underestimated by about five orders of magnitude. Impact on Δv is small (the drag term of radial drift scales with u_g and differs between grains only through St), but it is physically wrong and must be fixed in M2 with a toggle for both choices so the effect can be shown.

F7 (B): negative densities are clipped and the mass is redistributed silently

fix_negative_densities (related/msc-thesis.py/src/functions/solver/__init__.py, called at solver_v2.py:157) zeroes negative bins and rescales all positive bins to restore total mass. This conserves mass by construction, so the reported mass error after such a fix does not measure the integrator’s conservation. Negative values are a symptom (step too large, noisy estimator, tolerance) and must be reported, not hidden (global rule: no fallbacks that hide bugs). New code: count and report negativity events, and fail in tests.

F8 (B): kernel boundary and fragmentation choices

• Pairs with m_i + m_j beyond the grid are skipped entirely for coagulation (kernel.py, near_upper_bound, k_h >= N_m), including their loss term. Particles in the top bins become inert, which creates the “traffic jam” the thesis describes. Make the policy explicit and configurable (ADR 0005).
• MRN fragments extend up to m_i + m_j. Dullemond’s notes argue for independent shattering of each body up to max(i, j), and add cratering for extreme mass ratios. Both are in scope for M4.
• Bouncing is absent (P_bounce = 0). It is part of the target model (Dominik and Dullemond 2024) and changes the equilibrium substantially.

F9 (C): diagnostic definitions

• Δ_stab is defined as a global relative mass error in ch. 3 but as a sum of per-bin relative errors in ch. 4. Use the global one.
• Δ_acc divides by ρ^c_i(t), which is zero in many bins (the thesis footnote mentions gaps). Replace by a mass-weighted norm that is well defined (see 03-verification.md).
• assert N.all() >= 0 (sampled.py:49) compares a boolean with 0 and is always true.

F10 (C): text and citation errors to fix before reuse

4. What the thesis left open (and Dullemond’s notes flag as TODO)

5. A gap nobody has addressed yet: the linear solve

For p particle properties the flattened state has N = N_m N_φ ⋯ entries, the full kernel sum costs O(N³) with fragmentation, and so does a dense LU of the N × N Jacobian in any implicit method. Sampling reduces the RHS cost to O(S · L) (S draws, L recipe entries per draw) but not the linear algebra. An implicit sampled scheme only wins if the Jacobian is also treated cheaply:

• Sparse Jacobian assembled from the same sampled recipes (O(S · L) non-zeros), solved by sparse LU or by Krylov (GMRES/BiCGStab) with a cheap preconditioner.
• Approximate Jacobians, which W-methods (Rosenbrock-W) and simplified-Newton Radau tolerate by design.
• Jacobian-free Newton-Krylov using sampled RHS evaluations.

This turns “does sampling pay off” from an RHS question into a whole-scheme question, and it is a natural, novel section of the paper (E13).

6. Summary for the plan

1. Re-implement the thesis model faithfully (with F6 to F8 fixed behind toggles) and reproduce its figures as a regression baseline.
2. Replace the biased estimator by unbiased IS / HT estimators, use the variance-optimal proposal, and re-measure the cost-accuracy trade-off with work-normalized efficiency.
3. Add variance reduction, the implicit sampled form, and cheap Jacobian strategies.
4. Move to two particle properties, where the payoff is expected, and compare against full summation, tensor-train or other fast deterministic methods where feasible, and the group’s RPMC.

Roadmap: milestones, tasks, acceptance criteria

Status: v0.2 (2026-11-10), updated with the grill decisions. Durations are rough estimates for part-time work and will be revised after M1.

Conventions: milestone Mn, task Mn.m. A milestone is done when all its acceptance criteria (AC) are checked by an automated test or a recorded experiment, just check is green, and the docs reflect the result. Experiments En are specified in 04-experiments.md, tests in 03-verification.md, crates in 05-architecture.md.

1. Overview

flowchart TD
    M0[M0 Plan and literature] --> M1[M1 Workspace, core types, CI]
    M1 --> M2[M2 Disk model]
    M2 --> M3[M3 Kinematics, rates, outcomes]
    M1 --> M4a[M4a Kernel framework + analytic kernels]
    M3 --> M4b[M4b Physical kernels:<br/>coag, frag, bouncing]
    M4a --> M4b
    M4a --> M5[M5 Explicit integrators<br/>+ analytic benchmarks]
    M5 --> M6[M6 Implicit integrators,<br/>own Radau IIA]
    M4b --> M6
    M6 --> R[R Thesis reproduction<br/>baseline]
    R --> M7[M7 Unbiased sampling,<br/>implicit sampled form]
    M7 --> M8[M8 Variance reduction,<br/>noise, Jacobians, cost]
    M8 --> M9[M9 Second property:<br/>filling factor]
    M0 --> V1[M10a Web app: design]
    V1 --> V2[M10b Views ship with<br/>each milestone M2 to M9]
    M2 --> V2
    V2 --> V3[M10c Web app: integration,<br/>sampling, comparisons]
    M8 --> V3
    M4b --> P0[M11a Short note to Dullemond<br/>with website URL]
    M7 --> P1[M11b Full pitch:<br/>bias result, paper proposal]
    M9 --> P2[M11c Paper draft,<br/>reproducibility package]
    P1 --> P2

gantt
    title Indicative timeline (part-time)
    dateFormat YYYY-MM-DD
    axisFormat %b %y
    section Foundations
    M0 Plan and literature          :done, m0, 2026-11-10, 14d
    M1 Workspace and core           :m1, after m0, 14d
    section Physics
    M2 Disk model                   :m2, after m1, 14d
    M3 Kinematics and outcomes      :m3, after m2, 21d
    M4a Kernel framework            :m4a, after m1, 14d
    M4b Physical kernels            :m4b, after m3, 21d
    section Numerics
    M5 Explicit integrators         :m5, after m4a, 14d
    M6 Implicit integrators, Radau  :m6, after m4b, 28d
    R Thesis reproduction           :r, after m6, 14d
    section Research
    M7 Unbiased sampling            :m7, after r, 28d
    M8 Variance, noise, cost        :m8, after m7, 42d
    M9 Second property              :m9, after m8, 42d
    section Visualization
    M10a Design                     :m10a, after m0, 28d
    M10b Views per milestone        :m10b, after m2, 280d
    M10c Integration and sampling   :m10c, after m8, 28d
    section Paper
    M11a Short note to Dullemond    :milestone, m11a, after m4b, 0d
    M11b Full pitch                 :milestone, m11b, after m7, 0d
    M11c Draft and package          :m11c, after m9, 42d

2. Milestones

M0: Plan and literature (this document set)

• M0.1 Read thesis, code, coag_py, Dullemond notes, local papers. Done.
• M0.2 Thesis review with findings (01-thesis-review.md). Done.
• M0.3 Literature map, glossary, notation. Done.
• M0.4 Vision, roadmap, verification, experiments, architecture, ADR set. Done (drafts).
• M0.5 Iterate on the plan with Vincent (grill session, 27 decisions in the decision log); ADRs 0001 to 0016 accepted. Done.
• M0.6 Decide audience and licence, repo hosting. Done: private until submission, open-source with the paper (06-risks-and-open-questions.md, Q1).

AC: Vincent has accepted the vision, roadmap, and ADRs. Done 2026-11-10.

M1: Workspace, core types, CI

• M1.1 Done. git init in repo/ (done with this plan; related/ stays outside and read-only), Cargo workspace per ADR 0002, justfile (check, test, bench, run), lints, deny.toml, project AGENTS.md.
• M1.2 Done (constants moved to M2.1, where the disk model needs them and the thesis-vs-CODATA parity question arises). coag-core: Real trait (f64, double-double; f128 behind a nightly feature, ADR 0016), physical constants (CODATA, IAU), Axis (linear, logarithmic; edges, centers, widths, index lookup), ParticleGrid with flattened multi-index (ADR 0004), State (N per cell and bin; one cell for now, ADR 0015) with mass and number moments.
• M1.3 Done. Config schema with serde + TOML, every toggle explicit, config hash (ADR 0010).
• M1.4 Done (coag-sim: .npy writer and reader, manifest.json with SHA-256 config hash). Output writer and run manifest (ADR 0011).
• M1.5 Done (ChaCha8 keyed by the seed, injective stream ids, golden values pinned). Deterministic RNG streams (ADR 0009).
• M1.6 Repository and CI done (vincentmader/coagulation, .github/workflows/check.yml); docs site built and tested locally (site/docs), deploy pending (DNS points to mk1, Caddy site file not yet in vps-infra). Private GitHub repository coagulation; CI runs just check; a deploy job publishes a docs page rendered from doc/ to https://coag.mader.xyz (static site on mk11, Caddy site file in vps-infra, per hetzner-deploy; DNS is set up). Local just check before each commit stays the primary gate.

AC: axis round-trip index(value(k)) = k for 10⁶ random k; geometric centers and q_m match the thesis grid to 10⁻¹⁵; manifest written and parsed back; same seed gives bit-identical draws on macOS and wasm32; f64 and double-double builds of coag-core pass the same tests; https://coag.mader.xyz serves the rendered docs.

M2: Disk model

• M2.1 Physical constants (CODATA, IAU; thesis values behind a toggle for parity). Star, gas surface density power law with normalization to M_gas (thesis) and alternative normalization (Σ at 1 au), radial grid.
• M2.2 Temperature (irradiated, optional floor), c_s, Ω_K, H_p, ρ_g(r, z), n_g, λ_mfp, u_th, ν_mol, ν_t, Re, pressure gradient, u_n, u_g with both viscosity choices (finding F6).
• M2.3 LocalConditions evaluated at (r₀, z₀) and passed down as plain data (no lower crate calls the disk model, ADR 0015); the cell position is a config key; analytic derivatives where available instead of finite differences.
• M2.4 Parity with thesis Disk/DiskRegion and with Dullemond’s diskmodel.py (for σ = 1).
• M2.5 Export radial profiles for the web app and for figures (thesis Figs. temperature, sound speed, gas density 2×2).
• M2.6 Disk view live on coag.mader.xyz (first app view, after the M10a design is approved).

AC: all profiles match thesis Python to 10⁻¹² relative (same inputs, molecular viscosity); u_g with ν_t matches an independent closed form for a power-law disk; profile figures regenerated.

M3: Dust kinematics, collision rates, outcomes

• M3.1 Particle properties per bin: a(m, φ), St (Epstein; Stokes regime as toggle), stopping time.
• M3.2 Relative velocities as separate, toggleable terms: Brownian, radial drift, azimuthal drift, turbulence (Ormel and Cuzzi 2007, all three regimes, cross checked with coag_py/dv_rel.py), differential settling (fixed), combination rule (quadrature; option for a velocity PDF).
• M3.3 Cross sections (compact; porous form for M9) and C_ij.
• M3.4 Outcome models: hard threshold, Maxwell-Boltzmann fragmentation probability (Stammler and Birnstiel 2022), bouncing with the Dominik and Dullemond (2024) local recipe implemented verbatim (ambiguities logged for the first contact) and validated against one of their local figures, v_frag scaling options.
• M3.5 Parity with thesis relative_velocity, collision_rate, collision_outcome_probabilities.
• M3.6 DustPy (uv, project-local) as an external oracle at one radius with transport off; matched settings documented; the DustPy-like robustness setup defined.
• M3.7 Δv and collision views live on the website.

AC: each Δv term matches thesis Python to 10⁻¹² (thesis settings); turbulent Δv reproduces Ormel and Cuzzi limiting cases; P^coag + P^bounce + P^frag = 1 to round-off for all pairs; Δv and C_ij maps regenerated (thesis Figs. relative velocity, cross section, rate coefficient).

M4a: Kernel framework and analytic kernels

• M4a.1 Sparse collision recipes (ADR 0005): per pair a list of (k, weight); assembly into RHS, Jacobian, dense K_kij export.
• M4a.2 Kovetz and Olund mapping with near-cancellation handling; upper boundary policies (drop pair, as thesis; pile up; outflow with accounting).
• M4a.3 Analytic test kernels: constant, sum, product.
• M4a.4 Diagnostics: per-pair mass defect, ΔX^rel, sparsity, ΔM^net_ij.
• M4a.5 f64 vs double-double comparison of per-pair mass defects and cancellation handling (E3).

AC: dense export equals coag_py.create_coag_kernel to 10⁻¹⁴ relative on its test grids (symmetric and asymmetric forms); per-pair mass defect < 10 ε_mach (m_i + m_j) for all pairs with cancellation handling; with handling disabled, the 10²⁰-span test from the coag_py docstring shows the breakdown.

M4b: Physical kernels

• M4b.1 Coagulation kernel R^coag X^coag.
• M4b.2 Fragmentation: pulverization; power-law shattering with configurable slope (MRN -11/6, Kun and Herrmann -3.27), upper cut-off m_i + m_j (thesis) or independent shattering up to max(i, j) (notes Eqs. 1.37 to 1.39); cratering with reverse Kovetz and Olund and its cancellation handling; bouncing with mass loss.
• M4b.3 Total kernel = Σ over outcomes; component kernels kept separately for visualization.
• M4b.4 Parity with thesis Kernel for coag, frag (MRN), total; thesis kernel figures (101, 102, 34) regenerated.
• M4b.5 Kernel and recipe-inspector views live; this triggers the short note to Dullemond (M11a).

AC: kernel parity to 10⁻¹² with thesis settings; ΔX^rel ≤ 10⁻¹² for coag, frag, total (thesis claim reproduced); sparsity 3 % / 30 % at N_m = 100 reproduced.

M5: Explicit integrators and analytic benchmarks

• M5.1 Integrator trait (ADR 0008); explicit Euler and Dormand-Prince 5(4) with step control (Heun and RK4 dropped: no experiment needs them).
• M5.2 Explicit stability limit estimator (CFL-like, Laibe and Lombart 2021).
• M5.3 Benchmarks: constant, sum, product kernels against analytic solutions (Ohtsuki 1990; coag_py examples), including coarse-grid artificial acceleration; an aggregation-fragmentation benchmark with known steady state (Sekiyama et al. 2017) if feasible.

AC: observed convergence orders 1 and 5 within ±0.2 on the constant kernel; Euler on the thesis linear-grid example reproduces thesis Fig. 112; analytic solutions matched to the discretization error documented by coag_py tests.

M6: Implicit integrators, own Radau IIA

• M6.1 Backward Euler with Newton (matches coag_py backwardeuler, niter), linearized one-step form (Brauer 2008 App. C).
• M6.2 One Rosenbrock-W method (tolerates approximate Jacobians; needed for E13).
• M6.3 Radau IIA order 5 after Hairer and Wanner: transformed stage system, simplified Newton, step-size control and error estimate, Jacobian reuse, dense output; analytic Jacobian from recipes.
• M6.4 Scaled variables and per-component tolerances (ADR 0003); positivity reporting instead of silent clipping (F7).
• M6.5 Parity with scipy Radau (via the Python bridge, ADR 0013) on the full thesis kernel.
• M6.6 Stiffness study: eigenvalue spectrum of J along the trajectory; step sizes chosen by each method.
• M6.7 Integration view live (trajectories, mass error, step sizes).

AC: Radau matches scipy Radau to 10⁻⁸ relative at rtol = 10⁻¹⁰; convergence order 5 on a smooth test; steady state of the thesis full model reached with mass error ≤ 10⁻¹² and no negativity events at default tolerances.

R: Thesis reproduction baseline

• R.1 Re-run thesis chapter 3 integrations (pure coag lin/log, pure frag, full model) and chapter 4 sampling runs with the thesis estimator (without replacement, scalar weights, as found in F1/F2) and with the documented proposal.
• R.2 Write results/baseline/ with all figures and a short report: what matches, what differs and why.

AC: thesis Figs. 103, 104, 106, 107, 109 reproduced qualitatively (same features); quantitative deviations explained.

M7: Unbiased sampling and implicit sampled form

• M7.1 Estimator trait (ADR 0007): full sum; thesis estimator (baseline); IS with replacement; HT without replacement (e.g. Poisson or conditional Poisson sampling); 1D (index) sampling for hit-and-stick.
• M7.2 Proposals: uniform, ρ_i ρ_j, F_i F_j factorized (Dullemond, W_i = 1, m_i, m_i²), ΔM^net N_i N_j, ‖K_·ij‖ N_i N_j; alias-method sampler.
• M7.3 Detailed balance per sample: each draw adds its full recipe (ADR 0006).
• M7.4 Implicit sampled kernel (Eq. 1.71) as a sparse quadratic kernel, usable by all implicit integrators.
• M7.5 Statistical test harness: replicas, bias with confidence intervals, variance, CRN.
• M7.6 Experiments E5, E6, E7, E8.

AC: for all unbiased estimators, |bias| < 3 standard errors over 10³ replicas on fixed N; mass conservation per step at round-off for every sampled step; E6 answers RQ1 with a figure.

Checkpoint M11b: after M7, the full pitch to Dullemond: bias finding, corrected estimator, first plots, paper proposal (German draft per email-writing, Vincent approves and sends). Sampling views go live on the website at the same time.

M8: Variance reduction, noise, Jacobians, cost

• M8.1 Hybrid split estimator (exact top-q pairs + sampled remainder), stratified sampling by k-blocks or by outcome type, control variates (previous-step sum, coarse-grid sum), common random numbers across replicas and stages.
• M8.2 Sample handling in integrators: frozen per step, averaged batches, resampling per step with step-size control decoupled from noise.
• M8.3 Noise floor and stability at steady state (E11, E12); balance-preserving sampling for the growth-fragmentation equilibrium.
• M8.4 Jacobian strategies: dense exact, sparse from sampled recipes, approximate (dominant pairs), JFNK with GMRES; preconditioners (E13).
• M8.5 Cost model: counted recipe evaluations, flops, wall-clock; work-precision diagrams (E14).
• M8.6 Performance pass: Rayon parallelism over pairs, SIMD-friendly layouts, benchmarks (criterion).

AC: η (work-normalized efficiency) reported for every estimator and integrator combination at N_m ∈ {50, 100, 200, 400}; the best combination beats full summation in η at some N_m, or a documented negative result with reasons.

Go/no-go gate for M9 (criteria fixed 2026-11-10, before any results): go if the best p = 1 scheme reaches η ≥ 2 × full summation at N_m = 200 at the 1 % target, or if the measured cost model projects ≥ 5 × for p = 2 at N_m × N_φ = 100 × 20. Otherwise publish p = 1 as a methods paper (bias analysis, optimal proposal, Jacobian study) and move p = 2 to follow-up work. The decision and its numbers are written into this file. Compute for the p = 2 sweeps is estimated here; paid server time needs Vincent’s OK.

M9: Second particle property (filling factor)

• M9.1 Filling factor grid, φ-dependent radius, St, cross section.
• M9.2 Minimal porosity physics behind a swappable PorosityModel trait, so richer physics (static compression after Kataoka et al. 2013, Ormel et al. 2007 regimes) can be added later: the Okuzumi et al. (2012) hit-and-stick volume rule, one published collisional compaction rule, compact fragments; φ enters radius, St, and cross section.
• M9.3 2D recipes (bilinear Kovetz and Olund in (m, V) or (m, φ)) with conservation of mass and the chosen second moment.
• M9.4 Full 2D sum as reference at modest resolution; sampled schemes at production resolution (E15).
• M9.5 Compare with φ(m, t) as dependent variable (Schöll et al. 2026) where the comparison is meaningful.

AC: 2D full and sampled results agree within noise at N_m × N_φ = 50 × 20; break-even surface over (N_m, N_φ, accuracy) computed; H5 decided.

M10: Interactive web visualization

See 07-visualization.md. Public at https://coag.mader.xyz from M2 on. Design before code (global rule).

• M10a Requirements, information architecture, and a design in Claude Design (or similar), in parallel with M1 and M2; Vincent approves.
• M10b Views ship with their milestones (disk M2, Δv and collisions M3, kernels M4b, integration M6, sampling M7), each with a screenshot review per ui-design.
• M10c Comparison views and run browser for experiment outputs (M8).
• M10d Shared Impressum and privacy pages for the mader.xyz family (static-website), linked from the app before the first public view.

AC: every config toggle from M2 to M8 reachable in the UI; live recomputation for N_m ≤ 100 under 1 s for kernels; screenshot review passed.

M11: Paper and outreach

See 08-paper-and-outreach.md.

• M11a Short note to Dullemond once the kernel views are live (after M4b): the thesis work continues, here is the URL, a paper is planned. German draft, Vincent approves and sends.
• M11b Full pitch after M7 (see the M7 checkpoint).
• M11c Paper draft (LaTeX or Typst per journal), figures from experiment outputs only, reproducibility package (tagged release, configs, scripts, data DOI e.g. Zenodo).
• M11d Internal review loop with co-authors; submission.

AC: every figure regenerates from a tagged commit with one command; co-authors sign off.

3. Task tracking

Until the repo has an issue tracker, track tasks in this file: mark finished tasks with “Done (commit abc1234)”. When GitHub is set up (Q1), convert each milestone into a GitHub milestone and each task into an issue with the same ID.

Verification and validation

How we know each layer is right. Tests live next to the code (unit), in crates/*/tests (integration), and in experiments/ (studies whose outputs are checked in CI at small size). Policy details follow the global testing skill.

1. Test pyramid

flowchart BT
    U[Unit tests<br/>axes, constants, formulas,<br/>single recipes] --> P[Property tests<br/>conservation per pair,<br/>symmetry, positivity]
    P --> O[Oracle and parity tests<br/>thesis Python, coag_py,<br/>scipy Radau, analytic solutions]
    O --> S[Statistical tests<br/>estimator bias and variance<br/>over replicas]
    S --> E[Experiment regression<br/>small-size versions of E1 to E16<br/>with stored reference outputs]

2. Oracles

Fixtures are generated data from deterministic configs, not real observational data.

3. Invariants checked automatically

Statistical tests use fixed seeds, so they are deterministic; their thresholds are chosen so a correct estimator fails with probability below 10⁻⁶ per test, documented next to each test.

4. Error metrics (replacing the thesis definitions)

Let ρ_k = m_k N_k and ρ = Σ ρ_k.

• Mass error: δ_M(t) = (ρ(t) - ρ(0)) / ρ(0). (Thesis chapter 3 definition.)
• Accuracy error vs reference (full kernel, same integrator, tight tolerances): mass-weighted L1 distance of the normalized distributions, e_acc(t) = ½ Σ_k | ρ_k(t)/ρ(t) − ρ^ref_k(t)/ρ^ref(t) |, which lies in [0, 1], is defined when bins are empty, and equals the fraction of mass “in the wrong bins”. Report also the log-space L2 error on bins above a mass floor for the shape at the small-mass tail.
• Time-lag diagnostic: mass-weighted mean log-mass ⟨ln m⟩(t) vs reference; separates a systematic lag (bias, F1) from noise.
• Steady-state noise floor: time-RMS of ‖dN/dt‖ (mass weighted) over the last decade in time, and the spread of ρ_k across replicas.
• Estimator diagnostics on fixed N: bias B_k, variance V_k, and the mass-weighted mean squared error MSE = Σ_k m_k² (B_k² + V_k) / ρ².
• Cost: recipe entries evaluated, kernel entries touched, Jacobian non-zeros, LU/Krylov flops, wall-clock (single thread, release build, fixed machine noted in the manifest).
• Work-normalized efficiency: η = 1 / (MSE × cost); the comparison metric for estimators. For full runs, work-precision diagrams (e_acc at t_end vs cost).
• Accuracy targets: e_acc = 1 % (headline), 10 % (quick look), 0.1 % where affordable; each at transient times and at steady state.

5. Integrator verification

6. Physics validation

7. Continuous integration budget

just check must finish in under 5 minutes on the M1 Pro: unit, property, parity, and statistical tests at small sizes (N_m ≤ 50, R ≤ 10³). Full studies run through heavy and are not part of check.

Experiments and the theory behind them

Each experiment has an ID, the question it answers, a config under experiments/En_<slug>/config.toml, a script that writes outputs with a run manifest, and one or more figures. Sizes given are the study sizes; CI runs a reduced version (see 03-verification.md).

1. Estimator theory (basis for E5 to E13)

Fix the state N (one evaluation of the RHS). The target is the vector

F_k = Σ_(i,j) K_kij N_i N_j = Σ_(i,j) a_ij v^(ij)_k, with a_ij = N_i N_j R_ij and v^(ij)_k = X_kij,

where v^(ij) is the per-collision recipe of pair (i, j) (its gain/loss vector, ADR 0005). Write T_ij = a_ij v^(ij) for the term vector.

Thesis estimator (without replacement, no weights). Drawing a set 𝒮 of distinct pairs with inclusion probabilities π_ij and summing T_ij over 𝒮 gives E[F̂] = Σ π_ij T_ij. With π_ij < 1 every process is slowed down by its own factor π_ij, so the bias is B = −Σ (1 − π_ij) T_ij. This is finding F1.

Horvitz-Thompson (HT). F̂ = Σ_{(i,j)∈𝒮} T_ij / π_ij is unbiased for any design with π_ij > 0 wherever T_ij ≠ 0.

Importance sampling with replacement (IS). Draw S pairs i.i.d. from p, F̂ = (1/S) Σ_s T_{(ij)s} / p{(ij)_s}. Unbiased; covariance (1/S)(Σ T T^⊤/p − F F^⊤). This is Dullemond’s Eq. 1.70 with p_ij = F_i F_j.

Variance-optimal proposal. For a weighted norm ‖x‖²_w = Σ_k w_k x_k², E‖F̂ − F‖²_w = (1/S)(Σ_ij ‖T_ij‖²_w / p_ij − ‖F‖²_w). Minimizing over p with Σp = 1 (Lagrange, Cauchy-Schwarz) gives

p*_ij ∝ ‖T_ij‖_w = N_i N_j R_ij ‖v^(ij)‖_w.

With the mass weighting w_k = m_k² (error in mass per bin) and replacing L2 by L1 (a common robust choice) this is p_ij ∝ N_i N_j Σ_k m_k |K_kij| = N_i N_j ΔM^net_ij: the relevance measure the thesis wrote down, without the extra m_i m_j factor (F3). Cost per draw differs between pairs (recipe length L_ij is 3 to 4 for coagulation, up to N_m for shattering), so the cost-aware optimum is p_ij ∝ ‖T_ij‖ / sqrt(L_ij). E7 tests this.

Hybrid estimator. Sort pairs by p*_ij; sum the top set 𝒟 holding a fraction q of the proposal mass exactly, sample the rest. Unbiased, variance from the tail only. For peaked distributions (the typical case, notes section 1.3.1) a small 𝒟 removes most variance. E12.

Implicit sampled form. Replacing a_ij = N_i N_j R_ij by a_ij evaluated at the unknown end-of-step state, but with weights from the start-of-step proposal (Dullemond Eq. 1.71), gives a sparse quadratic kernel K̃ with one entry block per drawn pair. Unbiased at the start of the step; first-order bias in Δt that vanishes as the step shrinks. E10 measures it.

Detailed balance per sample. Because each draw contributes a full recipe vector v^(ij) whose mass moment is zero, every sampled RHS conserves mass exactly, independently of which pairs were drawn. Mass conservation does not depend on estimator quality. ADR 0006.

2. Experiment catalogue

Physics and baseline

Estimators on a fixed state (no time integration)

Estimators inside time integration

Scaling and comparison

3. Experiment hygiene

• Headline experiments (E9 to E15) run in two physics setups (corrected thesis model, DustPy-like), each with and without bouncing.

• Every run writes a manifest (config hash, git revision, seed, machine, timings). No figure in the paper without a manifest.

• Experiments that use replicas use independent RNG streams per replica (ADR 0009) so results are reproducible and parallel.

• Heavy runs go through heavy (global rules) and are logged in experiments/LOG.md with date, command, duration, and outcome.

Software architecture

Status: v0.2 (2026-11-10). The binding decisions are the ADRs in ../adr; this document shows how they fit together. Crate versions are chosen and pinned in M1 after checking the current releases (global rule against outdated APIs); none are fixed here.

1. Layering

The model is a pipeline of pure stages. Each stage is a crate with a narrow interface, so it can be tested alone, swapped by config, and visualized on its own.

flowchart LR
    subgraph core [coag-core]
        C1[constants]
        C2[Axis, ParticleGrid]
        C3[State N_k]
        C4[Config types]
    end
    subgraph phys [physics]
        D[coag-disk<br/>star, gas, DiskRegion]
        K[coag-dust<br/>a, St, dv terms,<br/>C_ij, outcomes]
    end
    subgraph num [numerics]
        R[coag-kernel<br/>collision recipes,<br/>K assembly, Jacobian]
        S[coag-sampling<br/>proposals, estimators,<br/>sampled kernels]
        I[coag-integrate<br/>Integrator trait,<br/>Euler ... Radau IIA]
    end
    subgraph app [applications]
        X[coag-sim<br/>runs, diagnostics,<br/>outputs, manifest]
        CLI[coag-cli]
        W[coag-web<br/>Leptos + wasm]
        PY[coag-py<br/>PyO3 bridge]
    end
    core --> D --> K --> R
    core --> R
    R --> S
    R --> I
    S --> I
    I --> X
    S --> X
    X --> CLI
    X --> W
    X --> PY

Dependency rule: arrows point from lower to higher layers only. coag-kernel knows nothing about disks (it takes rate matrices); coag-integrate knows nothing about dust (it takes an Rhs); coag-sampling depends only on recipes. This is what makes analytic test kernels, the thesis kernel, and 2D kernels interchangeable.

2. Workspace layout (planned for M1)

repo/
  AGENTS.md  CLAUDE.md  README.md  LICENSE  justfile  Cargo.toml  deny.toml
  crates/
    coag-core/        constants, axes, grids, state, config schema, errors
    coag-disk/        disk model and DiskRegion
    coag-dust/        particle properties, relative velocities, rates, outcomes
    coag-kernel/      recipes, Kovetz-Olund, fragmentation models, analytic kernels
    coag-sampling/    proposals, alias sampler, estimators, sampled kernels
    coag-integrate/   integrators, step control, linear solvers, Jacobian strategies
    coag-sim/         run orchestration, diagnostics, outputs, manifest
    coag-cli/         binary: run, sweep, inspect
    coag-py/          PyO3 bindings (oracles, paper figures)
  web/                Leptos SPA (after design approval, M10)
  experiments/        E1 ... E16: config.toml, script, README, LOG.md
  tests/fixtures/     thesis and coag_py reference arrays (+ generator scripts)
  paper/              manuscript and figure scripts (M11)
  doc/                this documentation

../related/ (outside the repo) holds read-only reference material (thesis, its code, papers, coag_py). Nothing there is edited, and nothing is copied by default: what we need is rebuilt here, mostly from scratch in Rust.

3. Core abstractions

classDiagram
    class Axis {
        +edges() &[f64]
        +center(k) f64
        +width(k) f64
        +index_below(x) Option~usize~
    }
    class ParticleGrid {
        +axes: Vec~Axis~
        +len() usize
        +flat(multi) usize
        +multi(flat) MultiIndex
        +mass(flat) f64
    }
    class State {
        +n: Vec~f64~
        +mass_density(grid) f64
        +number_density() f64
    }
    class RateMatrices {
        +collision: Matrix
        +by_outcome: OutcomeMatrices
    }
    class Recipe {
        +pair: PairIndex
        +rate: f64
        +entries: Vec~RecipeEntry~
        +mass_defect(grid) f64
    }
    class RecipeSet {
        +recipes: Vec~Recipe~
        +rhs(state, out)
        +jacobian(state) SparseMatrix
        +dense() Tensor3
    }
    class Proposal {
        <<trait>>
        +weights(state, recipes) Vec~f64~
    }
    class Estimator {
        <<trait>>
        +draw(state, rng) SampledKernel
    }
    class SampledKernel {
        +recipes: Vec~WeightedRecipe~
        +rhs(state, out)
        +jacobian(state) SparseMatrix
    }
    class Rhs {
        <<trait>>
        +eval(t, y, out)
        +jacobian(t, y) JacobianRepr
    }
    class Integrator {
        <<trait>>
        +step(rhs, t, y, dt) StepResult
    }
    ParticleGrid o-- Axis
    RecipeSet o-- Recipe
    RecipeSet ..|> Rhs
    SampledKernel ..|> Rhs
    Estimator --> Proposal
    Estimator --> SampledKernel
    Integrator --> Rhs

Key choices, each backed by an ADR:

• State is N_k (particles per volume per bin) on a flattened multi-index, so p = 1 and p = 2 share all numerics (ADR 0004).
• Recipes are the unit of physics: one per pair (i, j) and outcome, a short list of (k, weight) with a guaranteed zero mass moment. RHS, Jacobian, dense tensor, sampled kernels, and visualizations are all assembled from recipes (ADR 0005, ADR 0006).
• Estimators return a SampledKernel, which implements the same Rhs trait as the full RecipeSet. Integrators never know whether they integrate a full or a sampled kernel (ADR 0007).
• Integrators are own implementations behind one trait; scipy is an oracle, not a dependency (ADR 0008).
• Cells: the state is cells × particle grid (one cell for now); physics below coag-disk receives LocalConditions and never calls the disk model, so radial transport can be added later by operator splitting or a coupled Jacobian (ADR 0015).
• Precision: numerics are generic over Real (f64 default, double-double, nightly f128) with an f64-specialized fast path (ADR 0016).
• Jacobian representation is an enum (dense, sparse, approximate, matrix-free) so E13 can swap strategies (ADR 0014).

4. One time step with sampling

sequenceDiagram
    participant Sim as coag-sim Runner
    participant Est as Estimator
    participant RNG as RngStream
    participant Int as Integrator (Radau IIA)
    participant SK as SampledKernel
    participant Diag as Diagnostics
    Sim->>Est: draw(state_n, recipes)
    Est->>RNG: stream(run, replica, step)
    Est->>Est: proposal weights p_ij, alias table
    Est-->>Sim: SampledKernel (weighted recipes, implicit form)
    Sim->>Int: step(SK, t_n, N_n, dt)
    loop simplified Newton
        Int->>SK: eval(t, y)
        Int->>SK: jacobian(t, y) (reused)
    end
    Int-->>Sim: N_n+1, error estimate, dt_next
    Sim->>Diag: mass error, negativity, cost counters
    Diag-->>Sim: record (manifest + arrays)

The sample is frozen within the step by default (deterministic ODE step, clean step-size control). E9 compares alternatives; the Runner owns that policy, not the integrator.

5. Configuration and toggles

One TOML file describes a run. Every physical and numerical choice is explicit; defaults reproduce the thesis (ADR 0010). Sketch:

[grid.mass]
min = 1e-15          # kg
max = 1e-8
bins = 100
scale = "log"

[disk]
star_mass = 1.0      # solar masses
radius = 10.0        # au
surface_density = { profile = "power_law", exponent = 0.8, normalization = "disk_mass", disk_to_star = 0.01 }
gas_velocity_viscosity = "turbulent"   # "molecular" reproduces the thesis (finding F6)

[dust.relative_velocity]
brownian = true
radial_drift = true
azimuthal_drift = true
turbulence = "ormel_cuzzi_2007"
settling = false
combine = "quadrature"

[collisions]
outcome = "maxwell_boltzmann"          # "threshold", "dominik_dullemond_2024"
v_frag = 1.0                           # m/s
bouncing = false
fragments = { model = "power_law", slope = -1.8333, upper = "sum" }   # "max_body", "pulverization"
cratering = false
upper_boundary = "drop_pair"           # "pile_up", "outflow"
cancellation_handling = true

[sampling]
estimator = "importance"               # "full", "thesis", "horvitz_thompson", "hybrid"
proposal = "optimal_mass_l1"
samples = 2000
hold = "per_step"                      # "per_macro_step", "batch_average"
implicit_form = true

[numerics]
precision = "f64"                     # "double_double", "f128" (nightly only)

[integrator]
method = "radau_iia5"
rtol = 1e-6
atol_scaled = 1e-12
jacobian = "sparse_sampled"

[time]
start = 1.0          # s
end = 3e10
outputs = 200
scale = "log"

[run]
seed = 42
replicas = 1

6. Performance plan

7. Web application (M10, after design approval)

flowchart LR
    subgraph browser [Browser]
        UI[Leptos SPA<br/>panels, toggles, plots]
        WASM[coag core in wasm<br/>disk, dv, kernels,<br/>small integrations]
        WK[Web Worker<br/>long computations]
    end
    subgraph server [mk11: coag.mader.xyz]
        ST[Static files<br/>Caddy]
        RES[Experiment outputs<br/>npy + JSON manifests]
    end
    UI <--> WASM
    UI <--> WK
    WK --> WASM
    UI -->|fetch| RES
    ST --> UI

Small problems (N_m ≤ 100, kernels and short integrations) run live in wasm; large studies are computed offline with coag-cli and loaded from files. No backend server is needed. Details in 07-visualization.md.

8. Error handling

• Library crates return Result with typed errors (thiserror); the CLI reports them (anyhow at the binary edge only).
• Invalid physics (negative densities beyond tolerance, mass error above bound, non-finite values) is an error with context, never silently corrected (finding F7).
• Debug builds assert conservation per recipe; release builds check it per step at negligible cost.

Risks and open questions

Open questions for Vincent

Risks

Decisions log

Accepted decisions move from here into ADRs. Rejected options stay in the ADR’s “Alternatives” section.

Interactive visualization: requirements

Status: requirements draft. No frontend code before a design is approved (global ui-design rule); M10a produces that design from this document, in parallel with M1 and M2; views then ship with their milestones.

1. Purpose

A browser tool in which every layer of the model can be inspected, every component toggled, and every integrator and estimator compared side by side. Two uses: (1) our own research instrument while developing M2 to M9; (2) a demonstration for Dullemond and the group, and possibly figures for talks. It is public at https://coag.mader.xyz from the first view on, so it can be checked live and shared early.

2. Views

3. Interaction principles

• Every toggle maps 1:1 to a config key; the current UI state can be exported as a config.toml and imported again, so anything seen in the browser can be rerun with coag-cli.
• Linked views: clicking a pair (i, j) in any matrix highlights the same pair in Δv, C_ij, P, kernel, and sampling views.
• Long computations run in a Web Worker; the UI never blocks.
• Every plot can be exported as SVG/PNG plus the data array.

4. Technical constraints

• Stack: Leptos SPA, core crates compiled to wasm (ADR 0012); plotting library chosen per the data-visualization skill in M10a (heat maps and animated lines with 10⁴ points must stay smooth).
• Deterministic results identical to the CLI for the same config and seed (same RNG, ADR 0009).
• Hosting: static files on mk11 behind the shared Caddy (site file in vps-infra/mk11/caddy/sites/), deployed by GitHub Actions (hetzner-deploy); no backend. DNS for coag.mader.xyz is set up.
• Public and English only (documented exception to the German/English default: the audience is physicists). No analytics, no cookies, all assets bundled, so the privacy policy stays minimal; shared Impressum and privacy pages for the mader.xyz family are linked from the app.

5. Acceptance

• All toggles from M2 to M8 reachable.
• Kernel recomputation at N_m = 100 below 1 s in the browser on the M1 Pro.
• Screenshot review per ui-design passed (symmetric spacing, consistent typography, readable color scales in light and dark mode).

Paper and outreach plan

1. Collaboration

• Primary contact: Prof. Dr. Cornelis P. Dullemond (ITA, Heidelberg), thesis supervisor and author of the method’s original formulation (notes section 1.3, coag_py MC prototype). Co-authorship is natural: the method idea is his, the systematic study and implementation are ours.
• Possible further co-authors / contacts in or near the group: M. Beutel (RPMC, comparison), authors of Schöll et al. 2026 (porosity, fractal dimension). Only via Dullemond.
• Timing, in two steps: (1) a short note once the kernel views are live on https://coag.mader.xyz (after M4b): the thesis work continues, here is the URL, a paper is planned. (2) The full pitch after M7: the bias finding (F1) with one figure, the corrected estimator matching the full kernel, the plan for p = 2, and the proposal to write a paper together. Ask for his view on scope (porosity recipe, RPMC comparison, journal), on open points of the bouncing recipe (M3.4), and on whether coag_py may be cited as code or redistributed.
• Form: German email (email-writing skill). Drafted by Claude, approved and sent by Vincent. Register (“Sie” or “du”) per Vincent’s past correspondence with him.

2. Paper structure (draft)

1. Introduction: multi-property dust evolution and its cost; grid vs particle methods; contribution.
2. The coagulation-fragmentation equation on a grid: recipes, detailed balance, kernels used (coag, shattering, cratering, bouncing).
3. Monte Carlo estimation of the collision sum: estimators (thesis, HT, IS, hybrid), bias of unweighted sampling, optimal proposal, mass conservation per sample.
4. Time integration with sampled kernels: implicit form, Radau IIA and W-methods, Jacobian strategies, noise at steady state.
5. Results for p = 1: bias/variance, proposals, noise floor, work-precision, break-even.
6. Results for p = 2 (mass, filling factor): cost surface, accuracy, physical results.
7. Comparison with other approaches (RPMC, low-rank/TT, DG).
8. Discussion and outlook (radial transport, GPUs, more properties). Appendices: kernel construction details, Radau implementation notes, test suite.

3. Reproducibility package

• Tagged release of the repository with the exact commit used for the paper.
• experiments/ configs and scripts for every figure; one command per figure (just figure <id>).
• Output data archived with a DOI (Zenodo), including manifests.
• Code licence decided at release (Q1).

4. Writing tooling

LaTeX with the journal class (A&A provides aa.cls), since co-authors and the journal expect it; figures from Python (matplotlib) via the PyO3 bridge (ADR 0013). Typst only for internal notes.

5. Side outputs

• A short technical note or blog post on the web tool (optional, after the paper).
• A talk at the group meeting (ask Dullemond).

Architecture Decision Records

Format and rules: ADR 0001. Template: 0000-template.md.

ADR 0001: Record architecture decisions

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: doc/README.md

Context

The project will run for months or years, with gaps, several AI agents, and later co-authors who need to understand why the code looks the way it does. Method choices (estimators, integrators, units) are research decisions whose reasons must survive.

Decision

Record every significant architecture or method decision as an ADR in doc/adr/, numbered sequentially, using 0000-template.md (a light MADR variant). Statuses: Proposed, Accepted, Superseded, Rejected. ADRs are never edited in substance after acceptance; a new ADR supersedes an old one and both link to each other. The index is README.md.

Alternatives considered

• Decisions only in commit messages: not discoverable.
• A single design document: loses the history of alternatives.

Consequences

• Small overhead per decision.
• Agents and co-authors can read the reasons before changing code (project AGENTS.md will point here).

ADR 0002: Rust Cargo workspace with layered crates

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: 05-architecture.md; ADR 0012, 0013

Context

The thesis code is Python with dense N³ arrays and Python loops, which limits resolution and makes the p = 2 study impossible. The global default stack is Rust. The same core must run natively (studies), in wasm (web app), and from Python (oracles, figures).

Decision

A Cargo workspace in repo/ (the git root; a private GitHub repository coagulation until submission) with the crates coag-core, coag-disk, coag-dust, coag-kernel, coag-sampling, coag-integrate, coag-sim, coag-cli, coag-py, and later web/. Dependencies point only from lower to higher layers as drawn in the architecture document. Library crates are #![forbid(unsafe_code)] unless a measured hotspot justifies otherwise (then an ADR). publish = false until the licence question (Q1) is decided. Workspace lints and edition per the rust skill, versions pinned in M1 after checking current releases.

Alternatives considered

• One crate: simpler at first, but couples physics to numerics and slows compile times; splitting later is a bigger refactor.
• Python + Numba (as in early thesis versions): fast to write, but no wasm path, weaker typing, and the cost study needs control over memory layout.
• C++/Fortran (common in astro codes): no advantage over Rust for this project and against the default stack.

Consequences

• Clear test boundaries per layer; analytic kernels and 2D kernels plug in without touching integrators.
• Some boilerplate for crate interfaces.
• Co-authors who only use Python need the coag-py bridge (ADR 0013).

ADR 0003: SI units at interfaces, scaled variables inside integrators

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: Finding F5; notation.md

Context

The thesis uses SI, Dullemond’s code uses cgs. N_k spans many decades (about 3·10² m⁻³ for monomers at 10 au, far less for large bins) and time spans 1 s to 10¹¹ s. scipy defaults (rtol 10⁻³, atol 10⁻⁶) were used without scaling, so small bins were not error controlled.

Decision

All public APIs, configs, and outputs use SI units (kg, m, s, K), with documented convenience units only in config input (au, solar mass, year), converted at parse time. Integrators work on scaled variables: y_k = N_k m_k / ρ₀ (mass fractions per bin) by default, time in units of a reference collision time; tolerances are given in scaled form (rtol, atol_scaled). Conversions live in coag-core and are tested for round trips. Comparisons with coag_py convert cgs at the fixture boundary.

Alternatives considered

• cgs throughout (astro convention): conflicts with the thesis and its figures, invites mixed-unit bugs at the thesis parity boundary.
• A typed units crate (uom): strong guarantees but heavy generics in hot loops and awkward in wasm/PyO3 interfaces; reconsider if unit bugs appear.
• Unscaled N_k in the integrator: reproduces the thesis tolerance problem.

Consequences

• Error control becomes meaningful for every bin (mass-weighted).
• Every interface must say its unit; doc comments state units for each quantity.
• Parity tests with the thesis must use the thesis’ unscaled tolerances when reproducing its behavior.

ADR 0004: State is N_k on a flattened multi-property index

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: Dullemond notes section 1.3.1; M9

Context

The Smoluchowski equation is solved for N_k (particles per volume per bin) in both the thesis and coag_py. The p = 2 study (mass and filling factor) needs the same numerics on a 2D grid; Dullemond’s notes propose flattening (i, ε) into one index.

Decision

The solver state is a contiguous Vec<f64> of N over a ParticleGrid, which owns one Axis per property and maps between a flat index and the multi-index (row-major, mass fastest-varying is decided in M1 by benchmark). Kernels, estimators, integrators, and diagnostics only see flat indices; physics crates use multi-indices and per-property values. Moments (number, mass, second property) are computed by the grid.

Alternatives considered

• Separate code paths for 1D and 2D: duplicates numerics and makes the comparison unfair.
• State n(m) (density per unit mass): needs Δm factors everywhere; N_k is what the discrete kernel acts on.
• Mass densities ρ_k as state: natural for plots, but the kernel is bilinear in N; we scale to mass fractions only inside the integrator (ADR 0003).

Consequences

• p = 1 and p = 2 share all numerics and tests.
• Index arithmetic is centralized and tested.
• Recipes for p = 2 need a 2D Kovetz and Olund mapping (bilinear), designed in M9.

ADR 0005: Kernels are built from sparse per-pair collision recipes

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: Findings F8; 04-experiments.md section 1

Context

The thesis stores K_kij as dense N_m³ arrays (7.5·10⁶ doubles at N_m = 200; impossible for p = 2: (N_m N_φ)³). Hit-and-stick touches 3 to 4 bins per pair, fragmentation up to N_m. Sampling needs per-pair units with zero mass moment, and the Jacobian is assembled from the same pairs.

Decision

The kernel is a RecipeSet: for every unordered pair (i ≤ j) and every outcome with non-zero rate, a Recipe { pair, rate, entries: [(k, w)] } where w are the per-collision changes of N_k (loss -1 at i and j, Kovetz and Olund gains, fragment distribution, cratering remnant). Each recipe is constructed so Σ_k m_k w_k = 0 to round-off (checked in tests). RHS, Jacobian, the dense tensor (export for plots and parity), and sampled kernels are assembled from recipes. Boundary behavior when m_i + m_j exceeds the grid is an explicit policy: drop_pair (thesis), pile_up, outflow (mass leaves and is accounted).

Alternatives considered

• Dense tensor: simple, but memory and cost scale as N³ even for sparse physics, and a sampled kernel would still need per-pair grouping.
• Generic sparse tensor (COO/CSF): loses the per-pair grouping that detailed balance and sampling need.
• Computing recipes on the fly per evaluation: saves memory, but costs repeated work; possible later as an option for very large p = 2 grids.

Consequences

• Memory O(N_pairs · L) instead of O(N³).
• Per-pair conservation becomes a unit-testable property.
• The visualization can show a recipe for a clicked pair directly.
• Changing rates (e.g. time-dependent disk) requires rebuilding rates, not entries; recipes store rates separately from entries for that reason.

ADR 0006: Detailed balance per sample

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: Dullemond notes section 1.3.4; ADR 0005

Context

Mass conservation to round-off is the property that keeps long integrations stable. A sampled RHS must not trade conservation for speed. Dullemond’s notes require that each sampled pair contributes to all k it affects.

Decision

A drawn pair always contributes its complete recipe (all k entries, loss and gain) with one common weight. Sampling over k inside a fragmentation recipe (to reduce cost for long recipes) is allowed only if the sampled sub-recipe is itself mass conserving (e.g. a sampled fragment set renormalized to the pair’s mass). Every sampled step is checked for mass conservation like a full step.

Alternatives considered

• Independent samples per k (one option in the notes): breaks conservation per step; kept only as a comparison variant in E8, never as default.
• Post-hoc rescaling to restore mass (as fix_negative_densities did): hides errors (finding F7).

Consequences

• Mass conservation is independent of estimator quality, so errors in experiments are purely estimator bias and variance.
• Recipe length L enters the cost per draw; the cost-aware proposal accounts for it.

ADR 0007: Unbiased estimators as the default, thesis estimator as baseline

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: Findings F1 to F4; RQ1 to RQ3; E5 to E8

Context

The thesis estimator samples without replacement and does not reweight, so it underestimates every rate by its inclusion probability (F1). The runs also used a proposal different from the documented one (F2). Dullemond’s notes and prototype use an unbiased importance-sampling form.

Decision

An Estimator trait returns a SampledKernel that implements Rhs. Implementations: Full (exact), Thesis (faithful baseline, without replacement, no weights, configurable weights including the scalar bug), ImportanceSampling (with replacement, weights 1/(S p)), HorvitzThompson (without replacement, weights 1/π), Hybrid (exact dominant set plus sampled tail). Proposals are a separate trait (uniform, ρ_iρ_j, factorized F_iF_j with weight W, ΔM^net N_iN_j, optimal ‖T‖ with cost awareness). Sampling uses the alias method. Default for new studies: Hybrid or ImportanceSampling with the optimal proposal, decided after E7 and E12.

Alternatives considered

• Keep only the thesis estimator and tune ρ_sample: cannot remove the bias.
• Only factorized index sampling (Dullemond’s sample2d): cheaper proposals, but cannot express pair-level relevance such as Δv-dependent rates; kept as a proposal option.

Consequences

• Statistical tests (bias, variance) become part of CI.
• The thesis results are reproducible as a baseline and the bias can be quantified (E6).

ADR 0008: Own integrators behind one trait, scipy as oracle

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: Finding F5; M5, M6; E4, E9, E10, E13

Context

The thesis used scipy Radau through coag_py.react_0d, one call per macro step, with default tolerances. Studying sampled kernels needs control over stages, Newton iterations, Jacobian reuse, and step-size control (noise must not drive step rejection). Rust ODE crates exist, but none is known to offer Radau IIA with the hooks we need; this is re-checked in M6 before writing code.

Decision

Integrators are implemented in coag-integrate behind trait Integrator { fn step(&mut self, rhs: &dyn Rhs, t, y, dt) -> StepResult } with a separate step controller. Methods: explicit Euler, Dormand-Prince 5(4), backward Euler with Newton, linearized implicit Euler (Brauer 2008), one Rosenbrock-W method, Radau IIA order 5 following Hairer and Wanner (transformed stage system, simplified Newton, embedded error estimate). The Jacobian is requested from the Rhs in a representation chosen by config (ADR 0014). scipy Radau is used only as a test oracle (ADR 0013).

Alternatives considered

• scipy via Python bridge in production: no wasm, little control, Python overhead in the loop.
• Wrap Fortran RADAU5: fast and proven, but foreign build chain, no wasm, still limited control.
• An existing Rust crate: adopt if one passes the M6 parity tests and exposes Jacobian and stage hooks; otherwise use it as an extra oracle.

Consequences

• Significant implementation and testing effort in M6 (order tests, scipy parity, step statistics).
• Full control for the research questions RQ4 and RQ6.
• Note for write-ups: Radau IIA (3 stages) has order 5, not 4.

ADR 0009: Reproducible random number streams

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: E6 to E16; ADR 0012

Context

Estimator studies need many independent replicas, parallel execution, and identical results in the CLI and the browser. The thesis used the global np.random state without seeds.

Decision

Use a seedable, platform-independent generator from the rand ecosystem (a counter-based or PCG/ChaCha generator, chosen in M1 by speed and wasm support). Each random stream is derived deterministically from (run seed, replica, time step, purpose) so results do not depend on thread scheduling or on the number of draws in other streams. The seed and derivation scheme are recorded in the run manifest.

Alternatives considered

• One global RNG: results depend on execution order and thread count.
• OS entropy: not reproducible.

Consequences

• Replicas run in parallel with identical results to serial runs.
• Common random numbers across estimator variants become possible (variance reduction in comparisons).

ADR 0010: One TOML config per run, every assumption a toggle

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: 05-architecture.md section 5; 07-visualization.md

Context

The thesis config (config.toml) already exposed many switches but not all (e.g. the viscosity choice, boundary policy, estimator details were code edits). The research needs to isolate each assumption, and the web UI needs a 1:1 mapping to config keys.

Decision

Each run is described by one TOML file deserialized with serde into typed structs with deny_unknown_fields. Every physical and numerical choice is an explicit key with a documented default; defaults reproduce the thesis setup, and fixes of thesis behavior are separate keys. Convenience units in input are converted at parse time (ADR 0003). The canonical serialized config is hashed and stored in the run manifest. The web UI reads and writes the same format.

Alternatives considered

• Hard-coded settings per experiment script: not reproducible from outputs.
• Command line flags only: unwieldy for dozens of options.

Consequences

• Any figure can be regenerated from its manifest.
• Adding a physics option means adding a config key, a default, a test, and a UI control.

ADR 0011: Outputs as npy arrays plus JSON manifest

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: 03-verification.md; ADR 0010, 0012, 0013

Context

Outputs must be readable by numpy/matplotlib (figures), by the browser (web app), and by our own tests, and must carry provenance. Sizes are moderate (N_t × N up to about 10⁷ doubles per run).

Decision

Each run writes a directory with one .npy file per array (little-endian f64, readable by numpy without extra dependencies and parsable in a few lines in Rust/wasm) and a manifest.json with config, config hash, git revision, seed, machine, timings, and diagnostics summaries. Large studies may later add Arrow/Parquet if needed (new ADR).

Alternatives considered

• HDF5: standard in astro, but a heavy C dependency and poor wasm support.
• Parquet/Arrow: good for tables, less natural for N-d arrays; overkill now.
• Plain text (as in the thesis np.savetxt): slow, large, loses precision.

Consequences

• Simple tooling everywhere; easy to archive on Zenodo.
• Many small files per run; acceptable.

ADR 0012: Web frontend: Leptos SPA with the core compiled to wasm

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: 07-visualization.md; global ui-design rule

Context

Goal G5 asks for an interactive tool in which every model component can be toggled and visualized. The default web stack is a Leptos SPA. The physics and kernel code should be identical to the CLI to avoid divergent results.

Decision

The web app in web/ is a Leptos SPA that links the core crates compiled to wasm32 and runs heavier work in a Web Worker. Precomputed experiment outputs are loaded as static files (ADR 0011). Hosting: a public static site at https://coag.mader.xyz on mk11 (shared Caddy), deployed by GitHub Actions; English only. No frontend code is written before a design is approved (M10a); the design sprint runs in parallel with M1 and M2, and views then ship with their milestones. The plotting library is chosen in M10a per the data-visualization skill.

Alternatives considered

• Python dashboard (Dash, Panel, Streamlit): fast to build, but needs a server, duplicates computation in Python, and diverges from the Rust core.
• Plain TypeScript app with a wasm module: possible, but against the default stack and splits the codebase.
• Desktop app (Tauri): no advantage for sharing with the group.

Consequences

• Same numbers in browser and CLI for the same config and seed.
• Core crates must stay wasm-compatible (no threads required, no file system in core, feature-gated Rayon).

ADR 0013: Python only for figures and test oracles

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: Global stack rule (Python only where its ecosystem is the reason); ADR 0008

Context

Paper figures are expected in matplotlib quality, and scipy plus the thesis and coag_py code are the oracles for parity tests. Neither justifies Python in the simulation loop.

Decision

A thin PyO3 crate coag-py exposes run, kernel, and estimator entry points to Python. Python (managed with uv, linted per the python skill) is used for: fixture generation from the thesis and coag_py code, scipy Radau parity checks, DustPy runs as an external oracle, and figure scripts for the paper. The simulation itself never depends on Python.

Alternatives considered

• Matplotlib-free figures in Rust (plotters): possible, but slower to iterate to publication quality.
• Calling Python from Rust tests: fragile in CI; fixtures are generated once and committed instead.

Consequences

• Two toolchains in the repo (cargo, uv), both driven from the justfile.
• Fixtures pin the oracle behavior at a known commit of the thesis code.

ADR 0014: Linear algebra with a pluggable Jacobian representation

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: RQ6, E13; ADR 0008

Context

Implicit steps need solves with I − h γ J. For p = 1 dense LU on N ≤ 800 is cheap; for p = 2 (N up to about 4000 to 16000) dense LU costs as much as the full kernel sum (thesis review section 5). Sampled kernels yield sparse Jacobians.

Decision

Jacobians are represented by an enum: Dense, Sparse (CSR, assembled from recipes), Approximate (from a dominant-pair subset), MatrixFree (Jacobian-vector products by recipe sums). Solvers: dense LU, sparse LU, GMRES/BiCGStab with simple preconditioners (diagonal, ILU(0) on the dominant-pair Jacobian). Use the faer crate for dense and sparse factorizations if it fits after the M6 check of its current API and wasm support; otherwise nalgebra for dense plus own Krylov code.

Alternatives considered

• Dense only: makes the p = 2 cost study meaningless.
• Binding LAPACK/SuiteSparse: fast, but C dependencies and no wasm.

Consequences

• E13 can compare strategies under one interface.
• Some own Krylov and preconditioner code to test (against dense solves on small problems).

ADR 0015: Ready for non-local models without building them

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: vision section 2 (non-goals); ADR 0004; grill decision Q19

Context

The paper studies a local model at one radius (r₀ = 10 au). Radial transport and a global disk model are the most likely follow-up projects. The most expensive refactor in codes of this kind is physics that reads global disk state directly, so that adding a spatial dimension touches every layer.

Decision

• The solver state is organized as cells × particle grid. Today the number of cells is one.
• coag-disk evaluates the disk model at a cell position and returns a plain LocalConditions value (T, ρ_g, c_s, H_p, Ω_K, α, ν, pressure gradient, u_g, u_n, …). Everything below it (kinematics, rates, recipes, estimators, integrators) receives LocalConditions and never calls the disk model.
• Rhs implementations and the run loop operate on a list of cells, so a transport operator can be added later, either by operator splitting or as additional terms in a coupled Jacobian.
• The cell position (radius, height) is a config key from the start.
• No transport, advection, or diffusion code is written for the paper.

Alternatives considered

• Build a 1D radial model now: multiplies cost by N_r and adds a second research question to a paper that already has several.
• Ignore non-locality: fastest now, but locks physics to a single global point and forces a broad refactor later.

Consequences

• Small overhead (one extra loop level, one plain-data struct).
• Recipes and rate matrices are built per cell; for one cell this is identical to the local model.
• A later transport milestone needs a new ADR on splitting vs coupled integration.

ADR 0016: Switchable floating-point precision

• Status: Accepted (2026-11-10)
• Date: 2026-11-10
• Deciders: Vincent Mader
• Related: findings F5, F9; experiments E3, E4; ADR 0003, 0014; grill decision Q21

Context

The thesis reports mass errors around 10⁻¹² and relies on analytic cancellation handling because f64 cannot represent m_i + m_j when m_i ≪ m_j. Comparing f64 against an extended precision answers whether that error is a round-off floor or a method error. As of November 2026 Rust’s native f128 is nightly-only, with stabilization targeted for late 2026 at the earliest. Double-double arithmetic (about 106 bits of mantissa) is available on stable Rust, in pure Rust, and works in wasm.

Decision

• Numerics are generic over a Real trait defined in coag-core: axes, recipes, kernel assembly, RHS, explicit integrators, and the implicit integrators with a small own dense LU.
• Supported types: f64 (default), double-double (the stable extended type; the crate is chosen in M1 after checking current releases), and native f128 behind a nightly-only cargo feature until it stabilizes.
• The production path stays specialized for f64: faer-based linear algebra, alias-table sampling, and probabilities are f64 only.
• Precision is a config key ([numerics] precision), recorded in the run manifest.
• E3 and E4 compare f64 against double-double for per-pair mass defects, cancellation handling, and long-term mass error.

Alternatives considered

• f64 only: simplest, but leaves the round-off question unanswered.
• Arbitrary precision (MPFR via rug): C dependency, no wasm, much slower than needed.
• Native f128 only: blocked on nightly.

Consequences

• Generic code in the core crates (trait bounds, constants via Real::from_f64).
• Extended-precision runs are limited to small grids (own dense LU is O(N³), double-double arithmetic is several times slower than f64, to be measured in M1).
• Tests run for both f64 and double-double in CI on small sizes.
• Found in M1: twofloat documents its exp, ln, and powf as not accurate to full double-double precision, so the Real trait offers only arithmetic, sqrt, and abs. Grids are defined in f64 and promoted exactly, so f64 and extended runs use bit-identical bin masses; physics that needs transcendental functions computes in f64.
• Found in M1: on nightly 1.101 for aarch64-apple-darwin, f128::sqrt returns 0. Real::sqrt for f128 uses two Newton steps from the f64 root instead (tested to f128 accuracy).

ADR NNNN:

• Status: Accepted (2026-11-10) | Accepted | Superseded by NNNN | Rejected
• Date: YYYY-MM-DD
• Deciders: Vincent Mader
• Related: findings, experiments, other ADRs

Context

What problem forces a decision, with the facts that matter.

Decision

What we do, stated so it can be checked.

Alternatives considered

Each alternative with the reason it lost.

Consequences

What becomes easier, what becomes harder, what we must now do.

Literature and source map

Every source this project builds on, with its role, where it lives locally (if anywhere), and its reading status. Paths are relative to this file; related/ is the folder of symlinks at the project root.

Bibliographic details were checked against publisher or ADS records on 2026-11-10.

Status: read (studied for this plan), skimmed (abstract and relevant sections), to read (needed before the milestone named in brackets), cited (only known from the thesis bibliography).

1. Primary sources of this project

2. Coagulation physics and dust evolution in disks

3. Porosity and additional particle properties (the regime where sampling should pay off)

4. Numerical methods for the coagulation equation

5. Stochastic and reduced-cost methods (positioning of the paper)

6. Context references from the thesis introduction (not needed for the paper)

andrews_2007, kitamura_2002, burton_2013, natta_2006, kant_1755, laplace_1799, newton_1687, courant_friedrichs_lewy_1928, schumann_1940, ALMA/ESO/Hubble image credits. Full entries: references.bib.

7. Local material that turned out not to be relevant

Glossary

Terms as used in this project. Where the thesis, Dullemond’s notes, and the code use different words for the same thing, all are listed and the project term is in bold. Symbols are defined in notation.md.

Astrophysical setting

Collision physics

Smoluchowski equation and discretization

Time integration

Monte Carlo kernel sampling

Other Monte Carlo families (for positioning)

Software terms (this project)

Notation

One symbol table for the project, with the equivalents in the thesis, in Dullemond’s notes and coag_py, and the identifier planned for the Rust code. Units are SI (ADR 0003).

Indices

Grids

Disk and gas

Dust

Rates, kernels, estimators

Diagnostics