PySDM is a package for simulating the dynamics of population of particles immersed in moist air using the particle-based (a.k.a. super-droplet) approach to represent aerosol/cloud/rain microphysics. The package core is a Pythonic high-performance implementation of the Super-Droplet Method (SDM) Monte-Carlo algorithm for representing collisional growth (Shima et al. 2009), hence the name. PySDM has two alternative parallel number-crunching backends available: multi-threaded CPU backend based on Numba and GPU-resident backend built on top of ThrustRTC.
It is worth here to distinguish the dependencies of the PySDM core subpackage
(named simply PySDM
) vs. PySDM_examples
and PySDM_tests
subpackages.
PySDM core subpackage dependencies are all available through PyPI, the key dependencies are Numba and Numpy.
The Numba backend is the default, and features multi-threaded parallelism for multi-core CPUs. It uses the just-in-time compilation technique based on the LLVM infrastructure.
The ThrustRTC backend offers GPU-resident operation of PySDM leveraging the SIMT parallelisation model.
The dependencies of PySDM examples and test subpackages are summarised in
the requirements.txt
file.
Noteworthy, one of the examples (ICMW_2012_case_1
) uses MPyDATA,
a concurently developed sister project to PySDM.
Hints on the installation workflow can be sought in the .travis.yml file
used in the continuous integration workflow of PySDM for Linux, OSX and Windows.
- Shima et al. 2009 Fig. 2
(Box model, coalescence only, test case employing Golovin analytical solution) - Berry 1967 Figs. 6, 8, 10
(Box model, coalescence only, test cases for realistic kernels) - Arabas & Shima 2017 Fig. 5
(Adiabatic parcel, monodisperse size spectrum activation/deactivation test case) - Yang et al. 2018 Fig. 2:
(Adiabatic parcel, polydisperse size spectrum activation/deactivation test case) - ICMW 2012 case 1
(2D prescripted flow stratocumulus-mimicking aerosol collisional processing test case)
In order to depict the PySDM API with a practical example, the following listings provide a sample code roughly reproducing the Figure 2 from Shima et al. 2009 paper. It is a coalescence-only set-up in which the initial particle size spectrum is exponential and is deterministically sampled to match the condition of each super-droplet having equal initial multiplicity:
from PySDM.physics import si
from PySDM.initialisation.spectral_sampling import constant_multiplicity
from PySDM.initialisation.spectra import Exponential
from PySDM.physics.formulae import volume
n_sd = 2**13
initial_spectrum = Exponential(norm_factor=8.39e12, scale=1.19e5 * si.um**3)
sampling_range = (volume(radius=10 * si.um), volume(radius=100 * si.um))
attributes = {}
attributes['volume'], attributes['n'] = constant_multiplicity(n_sd=n_sd, spectrum=initial_spectrum, range=sampling_range)
The key element of the PySDM interface is the Core
class which instances are used to manage the system state and control the simulation.
Instantiation of the Core
class is handled by the Builder
as exemplified below:
from PySDM import Builder
from PySDM.environments import Box
from PySDM.dynamics import Coalescence
from PySDM.dynamics.coalescence.kernels import Golovin
from PySDM.backends import Numba
from PySDM.state.products.particles_volume_spectrum import ParticlesVolumeSpectrum
builder = Builder(n_sd=n_sd, backend=Numba)
builder.set_environment(Box(dt=1 * si.s, dv=1e6 * si.m**3))
builder.add_dynamic(Coalescence(kernel=Golovin(b=1.5e3 / si.s)))
products = [ParticlesVolumeSpectrum()]
particles = builder.build(attributes, products)
The backend
argument may be set to Numba
or ThrustRTC
what translates to choosing the multi-threaded backend or the
GPU-resident computation mode, respectively.
The employed Box
environment corresponds to a zero-dimensional framework
(particle positions are not considered).
The vectors of particle multiplicities n
and particle volumes v
are
used to initialise super-droplet attributes.
The Coalescence
Monte-Carlo algorithm (Super Droplet Method) is registered as the only
dynamic in the system (other available dynamics representing
condensational growth and particle displacement).
Finally, the build()
method is used to obtain an instance
of Core
which can then be used to control time-stepping and
access simulation state.
The run(nt)
method advances the simulation by nt
timesteps.
In the listing below, its usage is interleaved with plotting logic
which displays a histogram of particle mass distribution
at selected timesteps:
from PySDM.physics.constants import rho_w
from matplotlib import pyplot
import numpy as np
radius_bins_edges = np.logspace(np.log10(10 * si.um), np.log10(5e3 * si.um), num=32)
for step in [0, 1200, 2400, 3600]:
particles.run(step - particles.n_steps)
pyplot.step(x=radius_bins_edges[:-1] / si.um,
y=particles.products['dv/dlnr'].get(radius_bins_edges) * rho_w / si.g,
where='post', label=f"t = {step}s")
pyplot.xscale('log')
pyplot.xlabel('particle radius [µm]')
pyplot.ylabel("dm/dlnr [g/m$^3$/(unit dr/r)]")
pyplot.legend()
pyplot.show()
The resultant plot looks as follows:
- backends:
- initialisation:
- multiplicities: integer-valued discretisation with sanity checks for errors due to type casting
- r_wet_init: kappa-Keohler-based equilibrium in unsaturated conditions (RH=1 used in root-finding above saturation)
- spatial_sampling: pseudorandom sampling using NumPy's default RNG
- spectra: Exponential and Lognormal classes
- spectral_sampling: linear, logarithmic and constant_multiplicity classes
- physics:
- environments:
- Box: bare zero-dimensional framework
- MoistLagrangianParcelAdiabatic: zero-dimensional adiabatic parcel framework
- MoistEulerian2DKinematic: two-dimensional prescribed-flow-coupled framework with Eulerian advection handled by MPyDATA
- dynamics:
- Coalescence
- Condensation
- solvers (working in arbitrary spectral coordinate specified through external class, defaults to logarithm of volume):
- Displacement: includes advection with the flow & sedimentation)
- EulerianAdvection
- Attributes (selected):
- Products (selected):
Development of PySDM is supported by the EU through a grant of the Foundation for Polish Science (POIR.04.04.00-00-5E1C/18).
copyright: Jagiellonian University
code licence: GPL v3
tutorials licence: CC-BY
- SCALE-SDM (Fortran):
https://github.com/Shima-Lab/SCALE-SDM_BOMEX_Sato2018/blob/master/contrib/SDM/sdm_coalescence.f90 - Pencil Code (Fortran):
https://github.com/pencil-code/pencil-code/blob/master/src/particles_coagulation.f90 - PALM LES (Fortran):
https://palm.muk.uni-hannover.de/trac/browser/palm/trunk/SOURCE/lagrangian_particle_model_mod.f90 - libcloudph++ (C++):
https://github.com/igfuw/libcloudphxx/blob/master/src/impl/particles_impl_coal.ipp - LCM1D (Python) https://github.com/SimonUnterstrasser/ColumnModel