Draft:Super droplet method

• Comment: Currently too much of this is a how-to/textbook, and we do not have those, see WP:NOTTEXTBOOK. Look at some other pages such as Broyden's method, QM/MM, Anderson acceleration and others. I thnk the whole #SDM Monte-Carlo Algorithm section should be heavily trimmed, perhaps removed. Ldm1954 (talk) 13:19, 26 September 2025 (UTC)

---

In mathematical modeling of aerosols, clouds and precipitation, Super Droplet Method (SDM) is a Monte-Carlo approach for representing collisions and coalescence of particles in atmospheric fluid dynamics simulations. The method and its name was introduced in a 2007 arXiv e-print by Shin-ichiro Shima et al.. (preceded by a 2006 patent application and followed by 2008 RIMS Kôkyûroku and 2009 QJRMS papers).

File:Super droplet method teaser.gif

Evolution of mass density distribution and the underlying population of 2¹¹ super-particles in a Monte-Carlo SDM simulation employing additive coagulation kernel and exponential initial particle mass spectrum, for which the Safronov-Golovin analytical solution to the Smoluchowski coagulation equations can be used as a reference (animation based on Fig. 2 from Shima et al. 2007, generated using the basic Python implementation of SDM provided in this article).

SDM algorithm is a probabilistic alternative to the deterministic model of the process embodied in the Smoluchowski coagulation equations. Among the key characteristics of SDM is that it is not subject to the "curse of dimensionality" that hampers application of other methods when multiple particle attributes need to be resolved in a simulation. The algorithm is embarrassingly parallel, has linear time complexity, constant state vector size (number conservation of simulated particles during collisions) and zero numerical diffusion.

SDM in the cloud microphysics model taxonomy

Particle-based approaches, including SDM, for simulating the dynamics of sizes of particles in clouds are considered as one of three distinct cloud microphysics modelling paradigms, which can be classified as either Eulerian or Lagriangian in terms of formulation of dynamics in particle attribute space:

Lagrangian microphysics approaches:

• particle-based (in zero spatial dimensions as in cloud parcel models, analogous to moving-sectional/method-of-characteristics/method-of-lines approaches)

Eulerian microphysics approaches:

• bulk (in which moment based description of particle attribute distribution is used, with assumed size/mass distributions),
• bin (in which nonparametric histogram-based description of particle attribute space is used, possibly involving moment-based description within each histogram bin)

The term "super-droplet" approach/method has been used either in reference to the particular Monte-Carlo algorithm (even if not used to model clouds), or more broadly in reference to the particle-based approach for modeling of atmospheric clouds (even if neglecting coalescence processes). Applications of particle-based methods in atmospheric modelling, including for cloud microphysics modelling, and including Monte-Carlo techniques and super-particle nomenclature, predate SDM. Analogous Monte-Carlo particle-based methods (or particle swarm methods) have also been used for modelling accretion in astrophysical context.

The canonical version of SDM follows the "all or nothing" algorithm presented by Shima et al.. The SDM collision algorithm is a Gillespie algorithm (see also Kinetic Monte Carlo and Dynamic Monte Carlo methods) in the sense that the expected number of droplet collisions which is represented by superdroplets over a given time in a certain volume matches the expected number of real droplet collisions, although the variance can be substantially higher.

SDM Monte-Carlo Algorithm for Coagulation of Particles

The algorithm formulation presented herein follows the notation of Shima et al., but uses pseudocode rather than mathematical notation. The code snippets are valid Python (used to generate the animation above) and refer to NumPy and SciPy components.

Super-particle state

SDM models evolution of a particulate system composed of $n_s$ super droplets, with $i$-th super droplet representing a multiplicity ${\xi }_i$ of particles. Each super particle carries a set of attributes, among which there are extensive attributes, such as mass $m_i$, as well as auxiliary attributes which do not change upon aggregation, such as position in space.

Well-mixed control volume

The considered particle-laden volume is split into cells. Cell volumes should be small enough to consider them well mixed, and hence to assume that any particle in the cell can collide with any other with a probability dependent only on particle extensive attributes. In subsequent formulation of the algorithm, a single cell of volume $\Delta v$ is considered.

Attribute sampling

SDM allows for arbitrary initial sampling of the particle attribute distribution. Among the commonly used methods, there is inverse transform sampling which yields uniform initial multiplicities across the population. Such constant-multiplicity sampling implies evaluation of the quantile function at random locations chosen uniformly from the 0...1 range as in the listing below:

import numpy as np
from scipy import stats

def sample(*,
    rng: np.random.Generator,
    dist: stats.rv_continuous,
    norm: float,
    n_s: int,
):
    state = np.empty(
        shape=n_s,
        dtype=[('ξ', np.int64), ('m', np.float64)]
    )
    state['ξ'] =  round(norm / n_s)
    state['m'] =  dist.ppf(rng.uniform(0, 1, n_s))
    return state

Several other sampling methods were discussed and applied in studies employing SDM

Time stepping, candidate pairs and attribute update

Simulation using SDM involves time-stepping loop that advances the system state by $\Delta t$ in each step. SDM is a Monte Carlo algorithm and each step employs random number generation. As exemplified in the basic serial implementation of an SDM step given in the listing below, each involves:

• permuting the super particles to select a random set of n_pair non-overlapping pairs among all n_s super particles in the considered volume $\Delta v$;
• shuffling n_pair random numbers from a uniform distribution (stored into array φ)
• computing the p_ratio of all possible particle unordered pairs to the number of considered candidate pairs n_pair;
• enumerating across all candidate pairs with α numbering the pair, and j,k numbering the super droplets within a pair;
• evaluating the probability p_α of collision within the super-droplet pair α as a product of the coagulation kernel, the timestep-to-cell-volume ratio Δt/Δv, the p_ratio and ξ_j where j is chosen to point to the super droplet of higher multiplicity in the pair;
• evaluating the γ factor, which for p_α<1 is a Boolean flag resulting from comparison of the probability of coagulation with the random number from the 0...1 range; while in the case of p_α≥1, γ is augmented by the number of "certain" coagulations to represent repeated collision events within a single step (note that depending on the value of the ξ_j/ξ_k ratio, the number of repeated events may need be truncated incurring a collision-event deficit);
• updating particle attributes ξ_j, ξ_k, m_j and m_k.

Scenario A: single-collision event (γ=1)

Scenario B: multiple-collision event (γ>1 and ξ[j]≠γξ[k])

Scenario C: equal-ratio collision event (γ≥1 and ξ[j]=γξ[k])

from collections import abc

def sdm(*,
    rng: np.random.Generator,
    ξ: abc.MutableSequence[int],
    m: abc.MutableSequence[float],
    kern: abc.Callable[[float, float], float],
    Δt: float,
    Δv: float,
):
    """ SDM step assuming non-zero multiplicities """
    n_s = len(ξ)
    n_pair = n_s // 2
    
    pairs = rng.permutation(n_s)[: 2 * n_pair]
    φ = rng.uniform(0, 1, n_pair)

    p_ratio = n_s * (n_s - 1) / 2 / n_pair
    for α, (j, k) in enumerate(pairs.reshape(-1, 2)):
        p_jk = kern(m[j], m[k]) * Δt / Δv
        if ξ[j] < ξ[k]:
            j, k = k, j
        p_α = ξ[j] * p_ratio * p_jk
        γ = p_α // 1 + (p_α - p_α // 1) > φ[α]
        if γ != 0:
            γ = min(γ, (ξ[j] / ξ[k]) // 1)
            if ξ[j] - γ * ξ[k] > 0:
                ξ[j] -= γ * ξ[k]
                m[k] += γ * m[j]
            else:
                ξ[j] = ξ[k] // 2
                ξ[k] -= ξ[j]
                m[k] += γ * m[j]
                m[j] = m[k]

Collision Scenarios

The above algorithm can yield the following collision scenarios (see animations above, symbol m denotes an arbitrary mass unit, symbol M denotes the mass attribute of the super-droplet of respective color):

- A: single-collision event (γ=1 )

- B: multiple-collision event (γ > 1 and ξ[j] ≠ γξ[k])

- C: equal-ratio collision event (γ ≥ 1 and ξ[j] = γξ[k]), note: an equal-ratio collision event may also occur within scenario A

Concurrency

The SDM coagulation algorithm does not feature data dependencies across candidate pairs, hence is particularly well suited for SIMT concurrent operation. Implementations of SDM leverage both CPU and GPU multi-threading.

Algorithm development and applications

The SDM algorithm has been applied in diverse modelling studies, covering topics such as:

• mixed-phase clouds
• representation of collisional breakup
• aqueous chemical reactions involving soluble aerosol material
• in-cloud electro-coalescence
• role of clouds in nuclear fallout dynamics
• geo-engineering (cloud brightening) (the reference uses particle-based microphysics representation, but collisions are represented in a different way than SDM)
• simulations of processes occurring in cloud-chamber laboratories
• interaction of cosmic rays with clouds
• dynamics of fogs

Open-source implementations

(links point to Git repositories)

• reusable SDM libraries and packages:
  • libcloudph++ (C++/CUDA)
  • PySDM (Python/Numba/CUDA)
  • CLEO (C++/Kokkos), also integrated into ICON model.
  • Droplets.jl (Julia)
• SDM implementations integrated in other software:
  • SCALE-SDM (Fortran)
  • Pencil Code (Fortran)
  • PALM LES (Fortran)
  • LCM1D (Python)
  • superdroplet (Cython/Numba/C++11/Fortran 2008/Julia)
  • NTLP (FORTRAN)
  • LacmoPy (Python/Numba)
  • McSnow (FORTRAN)
  • Particula (Python)