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
| ID | Goal | Measure of done |
|---|---|---|
| G1 | Faithful, modular re-implementation of the thesis model (disk, kinematics, outcomes, kernels, integrators) in Rust. | Parity tests against thesis Python and coag_py pass to stated tolerances; thesis figures reproduced. |
| G2 | Integrator suite from explicit Euler to an own Radau IIA, behind one trait. | Convergence orders verified on analytic kernels; Radau matches scipy to 10⁻⁸ relative at equal tolerances. |
| G3 | Correct and efficient MC sampling: unbiased estimators, optimal proposals, variance reduction, implicit sampled form, detailed balance per sample. | Bias consistent with zero at 3σ over 10³ replicas; mass error at round-off; documented efficiency gains. |
| G4 | Answer the scaling question for p = 1 and p = 2 (mass, filling factor). | Work-precision diagrams and a break-even N where sampling wins (or a proof that it does not). |
| G5 | Interactive web visualization of every model layer and integrator, with toggles for each physical component. | Public at https://coag.mader.xyz (mk11), growing with every milestone from M2 on; design approved first. |
| G6 | Paper with C. P. Dullemond (and possibly group members), with a reproducibility package. | Draft accepted by co-authors; submitted to A&A or a computational journal. |
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
| ID | Question | Experiments |
|---|---|---|
| RQ1 | How much of the thesis accuracy error is bias and how much is variance? | E5, E6 |
| RQ2 | Which proposal p_ij minimizes variance per unit cost, and how close do cheap approximations (ρ_i ρ_j, ΔM^net N_i N_j, factorized F_i F_j) come to it? | E7 |
| RQ3 | Does sampling the same pairs for all k (one recipe per draw, in detailed balance) beat independent samples per k? (Dullemond’s open “COMMENT ON WHICH WORKS BEST”.) | E8 |
| RQ4 | How should samples interact with the integrator: frozen per step, averaged over several draws, resampled per stage, or used in the implicit form of Eq. 1.71? What does each do to stability, step-size control, and the noise floor at steady state? | E9, E10, E11 |
| RQ5 | Does MC noise destabilize the coagulation-fragmentation balance at the small-grain end, and can a hybrid estimator (exact sum over dominant pairs, sampled tail) or a balance-preserving sampling remove the problem? | E11, E12 |
| RQ6 | What does the implicit linear solve cost with sampled kernels, and can sparse or approximate Jacobians (W-methods, simplified Newton, JFNK) keep the whole scheme sub-cubic? | E13 |
| RQ7 | Where is the break-even point between full summation and sampling, as a function of N_m, N_φ, required accuracy, and physics (with or without fragmentation and bouncing)? | E14, E15 |
| RQ8 | How does the method compare with representative particle MC (RPMC), low-rank / tensor-train methods, and DG schemes at equal accuracy? | E16 (scope to be set with Dullemond) |
4. Working hypotheses
| ID | Hypothesis | Falsified if |
|---|---|---|
| H1 | The slower evolution in the thesis is caused by the missing 1/π weights; an unbiased estimator removes the lag. | Unbiased runs still lag the full-kernel solution by more than their noise. |
| H2 | p_ij ∝ ‖K_·ij‖ N_i N_j (mass-weighted) is near-optimal and clearly better than ρ_i ρ_j. | Measured variance ratio < 1.5 in favor of the optimal proposal. |
| H3 | A hybrid estimator (exact over the pairs holding a fraction q of the proposal mass, sampled rest) gives most of the variance reduction at a small cost overhead. | Efficiency η does not improve over plain IS for any q. |
| H4 | With implicit stepping, MC noise acts as a bounded forcing; the steady-state noise floor scales as S^{-1/2} and is independent of Δt once Δt is large. | Noise floor grows with Δt or the implicit solve diverges. |
| H5 | For p = 1 sampling wins only at large N_m or low accuracy targets; for p = 2 it wins at practical resolutions (N_m ≈ 100, N_φ ≈ 20 to 50) for 1 % mass-weighted accuracy, provided the Jacobian is treated sparsely or approximately. | Break-even N lies beyond what any study would use. |
5. What is new compared with prior work
| Prior work | What it does | What this project adds |
|---|---|---|
| Thesis (Mader 2024) | Biased pair sampling, scipy Radau, mass only | Unbiased estimators, optimal proposals, variance reduction, cost model including the linear solve, two properties |
Dullemond notes and coag_py prototype (2021) | Formulates IS (Eq. 1.70) and implicit form (Eq. 1.71); explicit/implicit demos on constant/linear kernels | Systematic study on physical kernels with fragmentation and bouncing, answers the listed TODOs, scalable implementation |
| DustPy, Birnstiel 2010, Brauer 2008 | Full summation, implicit, mass only | Sub-cubic RHS and Jacobian, multi-property path |
| Okuzumi 2009, Schöll et al. 2026 | Second property as averaged or dependent variable φ(m, t) | Keeps φ as an independent dimension at affordable cost |
| Ormel 2007, Zsom and Dullemond 2008, Beutel and Dullemond 2023, Gurrutxaga et al. 2026 | Particle MC, explicit in time | Grid-based MC that keeps implicit time stepping and deterministic mass conservation |
| Matveev et al. 2015, 2016 | Low-rank / tensor-train deterministic speed-up | Different trade-off (no rank assumption); compare honestly, possibly combine |
| Lombart and Laibe 2021, Lombart et al. 2024 | High-order DG, few bins | Orthogonal (fewer bins vs cheaper sum); discuss combination |
6. Paper pitch (to refine with Dullemond)
Working title options:
- Unbiased Monte Carlo sampling of the coagulation kernel for implicit dust evolution with multiple particle properties
- 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)
| Topic | Decision |
|---|---|
| Spatial scope | Local model at one radius for the paper; ready for non-local models (ADR 0015). |
| Reference physics | Corrected thesis model (findings F6 to F8 fixed) as primary setup; a DustPy-like setup as robustness check; the exact thesis model only in the reproduction (R). |
| Bouncing | Dominik and Dullemond (2024) local recipe, verbatim; headline experiments run with and without bouncing. |
| Second property | Filling factor φ on a log grid, minimal porosity physics behind a swappable trait, richer physics later. |
| Precision | f64 default, double-double and nightly f128 switchable (ADR 0016). |
| Metrics | Work-normalized efficiency η on snapshots; work-precision diagrams at e_acc = 1 % (headline), 10 %, and 0.1 % where affordable, transient and steady state. |
| Comparisons | Literature numbers and reproduced test problems; RPMC via Dullemond if possible; no own competitor implementations unless a referee asks. |
| Go/no-go for p = 2 | Criteria fixed in the roadmap (M8) before any results. |
| Web | Public at coag.mader.xyz from the first view, English only, views ship with their milestones. |
| Contact | Short note to Dullemond once kernel views are live (after M4b); full pitch after M7. |