kagrski/pysdm

Pythonic particle-based (super-droplet) cloud microphysics modelling with Jupyter examples

PySDM

PySDM is a package for simulating the dynamics of population of particles. It is intended to serve as a building block for simulation systems modelling fluid flows involving a dispersed phase, with PySDM being responsible for representation of the dispersed phase. Currently, the development is focused on atmospheric cloud physics applications, in particular on modelling the dynamics 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.

For a list of talks and other materials on PySDM, see the project wiki.

Dependencies

It is worth here to distinguish the dependencies of the PySDM core subpackage (named simply PySDM) indicated in the setup.py file vs. the dependencies of PySDM_examples and PySDM_tests subpackages indicated in the requirements.txt file.

PySDM core subpackage dependencies are: Numpy, Numba, SciPy, Pint and molmass.

The Numba backend named CPU 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 named GPU offers GPU-resident operation of PySDM leveraging the SIMT parallelisation model. Using the GPU backend requires nVidia hardware, CUDA driver and the ThrustRTC and CURandRTC Python packages.

PySDM ships with tutorial files depicting how to use the package from Julia and Matlab. These depend on Matlab-Python interface and on PyCall.jl, respectively.

Installation

To install PySDM, one may use: pip install --pre git+https://github.com/atmos-cloud-sim-uj/PySDM.git.

Demo notebooks reproducing results from literature:

0D box-model coalescence-only examples (work with both CPU and GPU backends):

• Shima et al. 2009 Fig. 2
  (Box model, coalescence only, test case employing Golovin analytical solution)
• Berry 1967 Figs. 5, 8 & 10
  (Box model, coalescence only, test cases for realistic kernels)

0D parcel-model condensation only examples (work with CPU backend only)

• 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)

2D kinematic-flow warm-rain example (works with CPU backend only)

• Arabas et al. 2015 Figs. 8 & 9:
  (2D prescribed-flow stratocumulus-mimicking aerosol collisional processing test case)

Hello-world example

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 ConstantMultiplicity
from PySDM.initialisation.spectra import Exponential

n_sd = 2**15
initial_spectrum = Exponential(norm_factor=8.39e12, scale=1.19e5 * si.um**3)
attributes = {}
attributes['volume'], attributes['n'] = ConstantMultiplicity(spectrum=initial_spectrum).sample(n_sd)

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 CPU
from PySDM.products.state import ParticlesVolumeSpectrum

builder = Builder(n_sd=n_sd, backend=CPU)
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 CPU or GPU 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.savefig('readme.svg')

The resultant plot looks as follows:

Package structure and API

• backends:
  • CPU=Numba: multi-threaded CPU backend using LLVM-powered just-in-time compilation
  • GPU=ThrustRTC: GPU-resident backend using NVRTC runtime compilation library for CUDA
• 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:
  • constants: physical constants partly imported from SciPy and molmass packages
  • dimensional_analysis: tool for enabling dimensional analysis of the code for unit tests (based on pint)
  • formulae: physical formulae partly imported from the Numba backend (e.g., for initialisation)
• environments:
  • Box: bare zero-dimensional framework
  • Parcel: zero-dimensional adiabatic parcel framework
  • Kinematic2D: two-dimensional prescribed-flow-coupled framework with Eulerian advection handled by PyMPDATA
• dynamics:
  • Coalescence
    • coalescence.kernels (selected)
      • Golovin
      • Geometric
      • Hydrodynamic
      • ...
  • Condensation
    • solvers (working in arbitrary spectral coordinate specified through external class, defaults to logarithm of volume):
      • default: bespoke solver with implicit-in-particle-size integration and adaptive timestepping (Numba only as of now, soon on all backends)
      • BDF: black-box SciPy-based solver for benchmarking (Numba backend only)
  • Displacement: includes advection with the flow & sedimentation)
  • EulerianAdvection
• Attributes (selected):
  • cell:
    • position_in_cell
    • cell_id
    • ...
  • droplet:
    • volume
    • multiplicities
    • critical_radius
    • ...
• Products (selected):
  • SuperDropletCount
  • ParticlesVolumeSpectrum
  • CondensationTimestep
  • ...

Credits:

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
licence: GPL v3

SDM patents (some expired, some withdrawn):

• https://patents.google.com/patent/US7756693B2
• https://patents.google.com/patent/EP1847939A3
• https://patents.google.com/patent/JP4742387B2
• https://patents.google.com/patent/CN101059821B

Other open-source SDM implementations:

• 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
• superdroplet (Cython/Numba/C++11/Fortran 2008/Julia)
  https://github.com/darothen/superdroplet