PySDM 1.13

Pythonic particle-based (super-droplet) warm-rain/aqueous-chemistry cloud microphysics package with box, parcel & 1D/2D prescribed-flow examples in Python, Julia and Matlab

PySDM

PySDM package:

PySDM-examples package:

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 features 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. The Numba backend (aliased CPU) features multi-threaded parallelism for multi-core CPUs, it uses the just-in-time compilation technique based on the LLVM infrastructure. The ThrustRTC backend (aliased GPU) offers GPU-resident operation of PySDM leveraging the SIMT parallelisation model. Using the GPU backend requires nVidia hardware and CUDA driver.

For an overview paper on PySDM v1 (and the preferred item to cite if using PySDM), see Bartman et al. 2021 arXiv e-print (submitted to JOSS). For a list of talks and other materials on PySDM, see the project wiki.

A pdoc-generated documentation of PySDM public API is maintained at: https://atmos-cloud-sim-uj.github.io/PySDM

Dependencies and Installation

PySDM dependencies are: Numpy, Numba, SciPy, Pint, chempy, pyevtk, ThrustRTC and CURandRTC.

To install PySDM using pip, use: pip install PySDM (or pip install git+https://github.com/atmos-cloud-sim-uj/PySDM.git to get updates beyond the latest release).

Conda users may use pip as well, see the Installing non-conda packages section in the conda docs. Dependencies of PySDM are available at the following conda channels:

• numba: numba
• conda-forge: pyevtk, pint and
• fyplus: ThrustRTC, CURandRTC
• bjodah: chempy
• nvidia: cudatoolkit

For development purposes, we suggest cloning the repository and installing it using pip -e. Test-time dependencies are listed in the test-time-requirements.txt file.

PySDM examples are hosted in a separate repository and constitute the PySDM_examples package. The examples have additional dependencies listed in PySDM_examples package setup.py file. Running the examples requires the PySDM_examples package to be installed. Since the examples package includes Jupyter notebooks (and their execution requires write access), the suggested install and launch steps are:

git clone https://github.com/atmos-cloud-sim-uj/PySDM-examples.git
cd PySDM-examples
pip install -e .
jupyter-notebook

Alternatively, one can also install the examples package from pypi.org by using pip install PySDM-examples.

PySDM examples (Jupyter notebooks reproducing results from literature):

Examples are maintained at the PySDM-examples repository, see PySDM-examples README.md file for details.

Hello-world coalescence example in Python, Julia and Matlab

In order to depict the PySDM API with a practical example, the following listings provide sample code roughly reproducing the Figure 2 from Shima et al. 2009 paper using PySDM from Python, Julia and Matlab. 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:

Julia (click to expand)

using Pkg
Pkg.add("PyCall")
Pkg.add("Plots")
Pkg.add("PlotlyJS")

using PyCall
si = pyimport("PySDM.physics").si
ConstantMultiplicity = pyimport("PySDM.initialisation.spectral_sampling").ConstantMultiplicity
Exponential = pyimport("PySDM.physics.spectra").Exponential

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

Matlab (click to expand)

si = py.importlib.import_module('PySDM.physics').si;
ConstantMultiplicity = py.importlib.import_module('PySDM.initialisation.spectral_sampling').ConstantMultiplicity;
Exponential = py.importlib.import_module('PySDM.physics.spectra').Exponential;

n_sd = 2^15;
initial_spectrum = Exponential(pyargs(...
    'norm_factor', 8.39e12, ...
    'scale', 1.19e5 * si.um ^ 3 ...
));
tmp = ConstantMultiplicity(initial_spectrum).sample(int32(n_sd));
attributes = py.dict(pyargs('volume', tmp{1}, 'n', tmp{2}));

Python (click to expand)

from PySDM.physics import si
from PySDM.initialisation.spectral_sampling import ConstantMultiplicity
from PySDM.physics.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(initial_spectrum).sample(n_sd)

The key element of the PySDM interface is the Particulator class instances of which are used to manage the system state and control the simulation. Instantiation of the Particulator class is handled by the Builder as exemplified below:

Julia (click to expand)

Builder = pyimport("PySDM").Builder
Box = pyimport("PySDM.environments").Box
Coalescence = pyimport("PySDM.dynamics").Coalescence
Golovin = pyimport("PySDM.physics.coalescence_kernels").Golovin
CPU = pyimport("PySDM.backends").CPU
ParticlesVolumeSpectrum = pyimport("PySDM.products.state").ParticlesVolumeSpectrum

radius_bins_edges = 10 .^ range(log10(10*si.um), log10(5e3*si.um), length=32) 

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(radius_bins_edges)] 
particulator = builder.build(attributes, products)

Matlab (click to expand)

Builder = py.importlib.import_module('PySDM').Builder;
Box = py.importlib.import_module('PySDM.environments').Box;
Coalescence = py.importlib.import_module('PySDM.dynamics').Coalescence;
Golovin = py.importlib.import_module('PySDM.physics.coalescence_kernels').Golovin;
CPU = py.importlib.import_module('PySDM.backends').CPU;
ParticlesVolumeSpectrum = py.importlib.import_module('PySDM.products.state').ParticlesVolumeSpectrum;

radius_bins_edges = logspace(log10(10 * si.um), log10(5e3 * si.um), 32);

builder = Builder(pyargs('n_sd', int32(n_sd), 'backend', CPU));
builder.set_environment(Box(pyargs('dt', 1 * si.s, 'dv', 1e6 * si.m ^ 3)));
builder.add_dynamic(Coalescence(pyargs('kernel', Golovin(1.5e3 / si.s))));
products = py.list({ ParticlesVolumeSpectrum(py.numpy.array(radius_bins_edges)) });
particulator = builder.build(attributes, products);

Python (click to expand)

import numpy as np
from PySDM import Builder
from PySDM.environments import Box
from PySDM.dynamics import Coalescence
from PySDM.physics.coalescence_kernels import Golovin
from PySDM.backends import CPU
from PySDM.products.state import ParticlesVolumeSpectrum

radius_bins_edges = np.logspace(np.log10(10 * si.um), np.log10(5e3 * si.um), num=32)

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(radius_bins_edges)]
particulator = 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. Finally, the build() method is used to obtain an instance of Particulator 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:

Julia (click to expand)

rho_w = pyimport("PySDM.physics.constants").rho_w
using Plots; plotlyjs()

for step = 0:1200:3600
    particulator.run(step - particulator.n_steps)
    plot!(
        radius_bins_edges[1:end-1] / si.um,
        particulator.products["dv/dlnr"].get()[:] * rho_w / si.g,
        linetype=:steppost,
        xaxis=:log,
        xlabel="particle radius [µm]",
        ylabel="dm/dlnr [g/m^3/(unit dr/r)]",
        label="t = $step s"
    )   
end
savefig("plot.svg")

Matlab (click to expand)

rho_w = py.importlib.import_module('PySDM.physics.constants').rho_w;

for step = 0:1200:3600
    particulator.run(int32(step - particulator.n_steps))
    x = radius_bins_edges / si.um;
    y = particulator.products{"dv/dlnr"}.get() * rho_w / si.g;
    stairs(...
        x(1:end-1), ... 
        double(py.array.array('d',py.numpy.nditer(y))), ...
        'DisplayName', sprintf("t = %d s", step) ...
    );
    hold on
end
hold off
set(gca,'XScale','log');
xlabel('particle radius [µm]')
ylabel("dm/dlnr [g/m^3/(unit dr/r)]")
legend()

Python (click to expand)

from PySDM.physics.constants import rho_w
from matplotlib import pyplot

for step in [0, 1200, 2400, 3600]:
    particulator.run(step - particulator.n_steps)
    pyplot.step(x=radius_bins_edges[:-1] / si.um,
                y=particulator.products['dv/dlnr'].get()[0] * 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 (generated with the Python code) looks as follows:

Hello-world condensation example in Python, Julia and Matlab

In the following example, a condensation-only setup is used with the adiabatic Parcel environment. An initial Lognormal spectrum of dry aerosol particles is first initialised to equilibrium wet size for the given initial humidity. Subsequent particle growth due to Condensation of water vapour (coupled with the release of latent heat) causes a subset of particles to activate into cloud droplets. Results of the simulation are plotted against vertical ParcelDisplacement and depict the evolution of Supersaturation, CloudDropletEffectiveRadius, CloudDropletConcentration and the WaterMixingRatio .

Julia (click to expand)

using PyCall
using Plots; plotlyjs()
si = pyimport("PySDM.physics").si
spectral_sampling = pyimport("PySDM.initialisation").spectral_sampling
multiplicities = pyimport("PySDM.initialisation").multiplicities
spectra = pyimport("PySDM.physics").spectra
r_wet_init = pyimport("PySDM.initialisation").r_wet_init
CPU = pyimport("PySDM.backends").CPU
AmbientThermodynamics = pyimport("PySDM.dynamics").AmbientThermodynamics
Condensation = pyimport("PySDM.dynamics").Condensation
Parcel = pyimport("PySDM.environments").Parcel
Builder = pyimport("PySDM").Builder
products = pyimport("PySDM.products")

env = Parcel(
    dt=.25 * si.s,
    mass_of_dry_air=1e3 * si.kg,
    p0=1122 * si.hPa,
    q0=20 * si.g / si.kg,
    T0=300 * si.K,
    w= 2.5 * si.m / si.s
)
spectrum=spectra.Lognormal(norm_factor=1e4/si.mg, m_mode=50*si.nm, s_geom=1.4)
kappa = .5 * si.dimensionless
cloud_range = (.5 * si.um, 25 * si.um)
output_interval = 4
output_points = 40
n_sd = 256

builder = Builder(backend=CPU, n_sd=n_sd)
builder.set_environment(env)
builder.add_dynamic(AmbientThermodynamics())
builder.add_dynamic(Condensation())

r_dry, specific_concentration = spectral_sampling.Logarithmic(spectrum).sample(n_sd)
v_dry = builder.formulae.trivia.volume(radius=r_dry)
r_wet = r_wet_init(r_dry, env, kappa * v_dry)


attributes = Dict()
attributes["n"] = multiplicities.discretise_n(specific_concentration * env.mass_of_dry_air)
attributes["dry volume"] = v_dry
attributes["kappa times dry volume"] = kappa * v_dry
attributes["volume"] = builder.formulae.trivia.volume(radius=r_wet) 

particulator = builder.build(attributes, products=[
    products.PeakSupersaturation(),
    products.CloudDropletEffectiveRadius(radius_range=cloud_range),
    products.CloudDropletConcentration(radius_range=cloud_range),
    products.WaterMixingRatio(radius_range=cloud_range),
    products.ParcelDisplacement()
])
    
cell_id=1
output = Dict()
for (_, product) in particulator.products
    output[product.name] = Array{Float32}(undef, output_points+1)
    output[product.name][1] = product.get()[cell_id]
end 
    
for step = 2:output_points+1
    particulator.run(steps=output_interval)
    for (_, product) in particulator.products
        output[product.name][step] = product.get()[cell_id]
    end 
end 

plots = []
ylbl = particulator.products["z"].unit
for (_, product) in particulator.products
    if product.name != "z"
        append!(plots, [plot(output[product.name], output["z"], ylabel=ylbl, xlabel=product.unit, title=product.name)])
    end
    global ylbl = ""
end
plot(plots..., layout=(1, length(output)-1))
savefig("parcel.svg")

Matlab (click to expand)

si = py.importlib.import_module('PySDM.physics').si;
spectral_sampling = py.importlib.import_module('PySDM.initialisation').spectral_sampling;
multiplicities = py.importlib.import_module('PySDM.initialisation').multiplicities;
spectra = py.importlib.import_module('PySDM.physics').spectra;
r_wet_init = py.importlib.import_module('PySDM.initialisation').r_wet_init;
CPU = py.importlib.import_module('PySDM.backends').CPU;
AmbientThermodynamics = py.importlib.import_module('PySDM.dynamics').AmbientThermodynamics;
Condensation = py.importlib.import_module('PySDM.dynamics').Condensation;
Parcel = py.importlib.import_module('PySDM.environments').Parcel;
Builder = py.importlib.import_module('PySDM').Builder;
products = py.importlib.import_module('PySDM.products');

env = Parcel(pyargs( ...
    'dt', .25 * si.s, ...
    'mass_of_dry_air', 1e3 * si.kg, ...
    'p0', 1122 * si.hPa, ...
    'q0', 20 * si.g / si.kg, ...
    'T0', 300 * si.K, ...
    'w', 2.5 * si.m / si.s ...
));
spectrum = spectra.Lognormal(pyargs('norm_factor', 1e4/si.mg, 'm_mode', 50 * si.nm, 's_geom', 1.4));
kappa = .5;
cloud_range = py.tuple({.5 * si.um, 25 * si.um});
output_interval = 4;
output_points = 40;
n_sd = 256;

builder = Builder(pyargs('backend', CPU, 'n_sd', int32(n_sd)));
builder.set_environment(env);
builder.add_dynamic(AmbientThermodynamics())
builder.add_dynamic(Condensation())

tmp = spectral_sampling.Logarithmic(spectrum).sample(int32(n_sd));
r_dry = tmp{1};
v_dry = builder.formulae.trivia.volume(r_dry);
specific_concentration = tmp{2};
r_wet = r_wet_init(r_dry, env, kappa * v_dry);

attributes = py.dict(pyargs( ...
    'n', multiplicities.discretise_n(specific_concentration * env.mass_of_dry_air), ...
    'dry volume', v_dry, ...
    'kappa times dry volume', kappa * v_dry, ... 
    'volume', builder.formulae.trivia.volume(r_wet) ...
));

particulator = builder.build(attributes, py.list({ ...
    products.PeakSupersaturation(), ...
    products.CloudDropletEffectiveRadius(pyargs('radius_range', cloud_range)), ...
    products.CloudDropletConcentration(pyargs('radius_range', cloud_range)), ...
    products.WaterMixingRatio(pyargs('radius_range', cloud_range)) ...
    products.ParcelDisplacement() ...
}));

cell_id = int32(0);
output_size = [output_points+1, length(py.list(particulator.products.keys()))];
output_types = repelem({'double'}, output_size(2));
output_names = [cellfun(@string, cell(py.list(particulator.products.keys())))];
output = table(...
    'Size', output_size, ...
    'VariableTypes', output_types, ...
    'VariableNames', output_names ...
);
for pykey = py.list(keys(particulator.products))
    get = py.getattr(particulator.products{pykey{1}}.get(), '__getitem__');
    key = string(pykey{1});
    output{1, key} = get(cell_id);
end

for i=2:output_points+1
    particulator.run(pyargs('steps', int32(output_interval)));
    for pykey = py.list(keys(particulator.products))
        get = py.getattr(particulator.products{pykey{1}}.get(), '__getitem__');
        key = string(pykey{1});
        output{i, key} = get(cell_id);
    end
end

i=1;
for pykey = py.list(keys(particulator.products))
    product = particulator.products{pykey{1}};
    if string(product.name) ~= "z"
        subplot(1, width(output)-1, i);
        plot(output{:, string(pykey{1})}, output.z, '-o');
        title(string(product.name), 'Interpreter', 'none');
        xlabel(string(product.unit));
    end
    if i == 1
        ylabel(string(particulator.products{"z"}.unit));
    end
    i=i+1;
end
saveas(gcf, "parcel.svg")

Python (click to expand)

import numpy as np
from matplotlib import pyplot
from PySDM.physics import si, spectra
from PySDM.initialisation import spectral_sampling, multiplicities, r_wet_init
from PySDM.backends import CPU
from PySDM.dynamics import AmbientThermodynamics, Condensation
from PySDM.environments import Parcel
from PySDM import Builder, products

env = Parcel(
    dt=.25 * si.s,
    mass_of_dry_air=1e3 * si.kg,
    p0=1122 * si.hPa,
    q0=20 * si.g / si.kg,
    T0=300 * si.K,
    w=2.5 * si.m / si.s
)
spectrum = spectra.Lognormal(norm_factor=1e4 / si.mg, m_mode=50 * si.nm, s_geom=1.5)
kappa = .5 * si.dimensionless
cloud_range = (.5 * si.um, 25 * si.um)
output_interval = 4
output_points = 40
n_sd = 256

builder = Builder(backend=CPU, n_sd=n_sd)
builder.set_environment(env)
builder.add_dynamic(AmbientThermodynamics())
builder.add_dynamic(Condensation())

r_dry, specific_concentration = spectral_sampling.Logarithmic(spectrum).sample(n_sd)
v_dry = builder.formulae.trivia.volume(radius=r_dry)
r_wet = r_wet_init(r_dry, env, kappa * v_dry)

attributes = {
    'n': multiplicities.discretise_n(specific_concentration * env.mass_of_dry_air),
    'dry volume': v_dry,
    'kappa times dry volume': kappa * v_dry,
    'volume': builder.formulae.trivia.volume(radius=r_wet)
}

particulator = builder.build(attributes, products=[
    products.PeakSupersaturation(),
    products.CloudDropletEffectiveRadius(radius_range=cloud_range),
    products.CloudDropletConcentration(radius_range=cloud_range),
    products.WaterMixingRatio(radius_range=cloud_range),
    products.ParcelDisplacement()
])

cell_id = 0
output = {product.name: [product.get()[cell_id]] for product in particulator.products.values()}

for step in range(output_points):
    particulator.run(steps=output_interval)
    for product in particulator.products.values():
        output[product.name].append(product.get()[cell_id])

fig, axs = pyplot.subplots(1, len(particulator.products) - 1, sharey="all")
for i, (key, product) in enumerate(particulator.products.items()):
    if key != 'z':
        axs[i].plot(output[key], output['z'], marker='.')
        axs[i].set_title(product.name)
        axs[i].set_xlabel(product.unit)
        axs[i].grid()
axs[0].set_ylabel(particulator.products['z'].unit)
pyplot.savefig('parcel.svg')

The resultant plot (generated with the Matlab code) looks as follows:

Contributing, reporting issues, seeking support

Submitting new code to the project, please preferably use GitHub pull requests (or the PySDM-examples PR site if working on examples) - it helps to keep record of code authorship, track and archive the code review workflow and allows to benefit from the continuous integration setup which automates execution of tests with the newly added code.

As of now, the copyright to the entire PySDM codebase is with the Jagiellonian University, and code contributions are assumed to imply transfer of copyright. Should there be a need to make an exception, please indicate it when creating a pull request or contributing code in any other way. In any case, the license of the contributed code must be compatible with GPL v3.

Developing the code, we follow The Way of Python and the KISS principle. The codebase has greatly benefited from PyCharm code inspections.

Issues regarding any incorrect, unintuitive or undocumented bahaviour of PySDM are best to be reported on the GitHub issue tracker. Feature requests are recorded in the "Ideas..." PySDM wiki page.

We encourage to use the GitHub Discussions feature (rather than the issue tracker) for seeking support in understanding, using and extending PySDM code.

Please use the PySDM issue-tracking and dicsussion infrastructure for PySDM-examples as well. We look forward to your contributions and feedback.

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

Related resources and open-source projects

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 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
• NTLP (FORTRAN)
  https://github.com/Folca/NTLP/blob/SuperDroplet/les.F

non-SDM probabilistic particle-based coagulation solvers

• PartMC (Fortran):
  https://github.com/compdyn/partmc

Python models with discrete-particle (moving-sectional) representation of particle size spectrum

• pyrcel: https://github.com/darothen/pyrcel
• PyBox: https://github.com/loftytopping/PyBox
• py-cloud-parcel-model: https://github.com/emmasimp/py-cloud-parcel-model