Draft:Super Droplet Method

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 a 2009 journal paper).

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 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.

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

• TODO: mixed phase
• TODO: breakup
• TODO: turbulent mixing
• TODO: aqueous chemistry
• TODO: electro-coalescence

• TODO: idealised benchmark test cases
• TODO: nuclear fallout
• TODO: geo-engineering (cloud brightening) (the reference uses particle-based microphysics representation, but collisions are represented in a different way than SDM)
• TODO: wildfire simulations
• TODO: cloud-chamber modelling
• TODO: cosmic rays
• TODO: DNS vs. LES
• TODO: contrails
• TODO: fogs
• TODO: development of parameterization for larger-scale models

Open-source implementations

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

References

This page will be placed in the following categories if it is moved to the article namespace.

Categories:

• Cloud and fog physics
• Computational physics
• Monte Carlo methods
• Stochastic simulation