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

Concern	Approach	When
Kernel assembly O(N_m² L)	Recipes built once per rate change; Rayon over pairs	M4
Full RHS	Loop over recipes, contiguous arrays, no allocation per call	M4
Sampling	Alias method O(1) per draw after O(P) setup; reuse tables while the proposal is stale-tolerant	M7
Jacobian and solve	Dense LU for small N; sparse LU or GMRES with preconditioner for large N	M6, M8
Parallel replicas	Independent RNG streams, Rayon over replicas	M7
Benchmarks	criterion benchmarks per crate; results recorded in experiments, not in CI	M8

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.

Keyboard shortcuts

Coagulation: Monte Carlo kernel sampling for the Smoluchowski equation