diff --git a/.github/workflows/tests+artifacts+pypi.yml b/.github/workflows/tests+artifacts+pypi.yml
index 094824944..74ea7e90f 100644
--- a/.github/workflows/tests+artifacts+pypi.yml
+++ b/.github/workflows/tests+artifacts+pypi.yml
@@ -102,10 +102,11 @@ jobs:
pip install -e .[tests]
- run: |
# TODO #682
- pylint --unsafe-load-any-extension=y --disable=fixme,invalid-name,missing-function-docstring,missing-class-docstring,protected-access,duplicate-code $(git ls-files '*.py' | grep -v ^examples)
+ pylint --unsafe-load-any-extension=y --disable=fixme,invalid-name,missing-function-docstring,missing-class-docstring,protected-access,duplicate-code $(git ls-files '*.py' | grep -v -e ^examples -e ^tutorials)
- run: |
# TODO #682
pylint --max-module-lines=550 --unsafe-load-any-extension=y --disable=fixme,too-many-function-args,unsubscriptable-object,consider-using-with,protected-access,too-many-statements,too-many-public-methods,too-many-branches,duplicate-code,invalid-name,missing-function-docstring,missing-module-docstring,missing-class-docstring,too-many-locals,too-many-instance-attributes,too-few-public-methods,too-many-arguments,c-extension-no-member $(git ls-files '*.py' | grep ^examples)
+ pylint --max-module-lines=550 --unsafe-load-any-extension=y --disable=fixme,too-many-function-args,unsubscriptable-object,consider-using-with,protected-access,too-many-statements,too-many-public-methods,too-many-branches,duplicate-code,invalid-name,missing-function-docstring,missing-module-docstring,missing-class-docstring,too-many-locals,too-many-instance-attributes,too-few-public-methods,too-many-arguments,c-extension-no-member $(git ls-files '*.py' | grep ^tutorials)
- run: |
# TODO #682
nbqa pylint --unsafe-load-any-extension=y --disable=fixme,duplicate-code,invalid-name,trailing-whitespace,line-too-long,missing-function-docstring,wrong-import-position,missing-module-docstring,wrong-import-order,ungrouped-imports,no-member,too-many-locals,unnecessary-lambda-assignment $(git ls-files '*.ipynb')
@@ -125,7 +126,7 @@ jobs:
matrix:
platform: [ubuntu-latest, macos-12, windows-latest]
python-version: ["3.8", "3.10"]
- test-suite: ["unit_tests", "smoke_tests/no_env", "smoke_tests/box", "smoke_tests/parcel", "smoke_tests/kinematic_1d", "smoke_tests/kinematic_2d"]
+ test-suite: ["unit_tests", "smoke_tests/no_env", "smoke_tests/box", "smoke_tests/parcel", "smoke_tests/kinematic_1d", "smoke_tests/kinematic_2d", "tutorials_tests"]
exclude:
- test-suite: "devops_tests"
python-version: "3.8"
@@ -160,6 +161,10 @@ jobs:
- if: startsWith(matrix.platform, 'ubuntu-')
run: echo NUMBA_THREADING_LAYER=omp >> $GITHUB_ENV
+ # install devops_tests for tutorials_tests
+ - if: matrix.test-suite == 'tutorials_tests'
+ run: pip install -r tests/devops_tests/requirements.txt
+
- env:
GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
run: pytest --durations=10 --timeout=900 --timeout_method=thread -p no:unraisableexception -We tests/${{ matrix.test-suite }}
diff --git a/README.md b/README.md
index 58ebafef2..9762bbcbe 100644
--- a/README.md
+++ b/README.md
@@ -37,10 +37,17 @@ The package features a Pythonic high-performance implementation of the
There is a growing set of example Jupyter notebooks exemplifying how to perform
various types of calculations and simulations using PySDM.
-Most of the example notebooks reproduce resutls and plot from literature, see below for
+Most of the example notebooks reproduce results and plot from literature, see below for
a list of examples and links to the notebooks (which can be either executed or viewed
"in the cloud").
+There are also a growing set of tutorials, also in the form of Jupyter notebooks.
+These tutorials are intended for teaching purposes and include short explanations of cloud microphysical
+ concepts paired with widgets for running interactive simulations using PySDM.
+Each tutorial also comes with a set of questions at the end that can be used as homework problems.
+Like the examples, these tutorials can be executed or viewed "in the cloud" making it an especially
+ easy way for students to get started.
+
PySDM has two alternative parallel number-crunching backends
available: multi-threaded CPU backend based on [Numba](http://numba.pydata.org/)
and GPU-resident backend built on top of [ThrustRTC](https://pypi.org/project/ThrustRTC/).
diff --git a/examples/PySDM_examples/Pyrcel/profile_plotter.py b/examples/PySDM_examples/Pyrcel/profile_plotter.py
new file mode 100644
index 000000000..1be4a431f
--- /dev/null
+++ b/examples/PySDM_examples/Pyrcel/profile_plotter.py
@@ -0,0 +1,67 @@
+import numpy as np
+from matplotlib import pyplot
+from open_atmos_jupyter_utils import show_plot
+
+from PySDM.physics.constants import si
+
+
+class ProfilePlotter:
+ def __init__(self, settings, legend=True, log_base=10):
+ self.settings = settings
+ self.format = "pdf"
+ self.legend = legend
+ self.log_base = log_base
+ self.ax = pyplot
+ self.fig = pyplot
+
+ def show(self):
+ pyplot.tight_layout()
+ show_plot()
+
+ def save(self, file):
+ # self.finish()
+ pyplot.savefig(file, format=self.format)
+
+ def plot(self, output):
+ self.plot_data(self.settings, output)
+
+ def plot_data(self, settings, output):
+ _, axs = pyplot.subplots(1, 2, sharey=True, figsize=(10, 5))
+ axS = axs[0]
+ axS.plot(
+ np.asarray(output["products"]["S_max"]) - 100,
+ output["products"]["z"],
+ color="black",
+ )
+ axS.set_ylabel("Displacement [m]")
+ axS.set_xlabel("Supersaturation [%]")
+ axS.set_xlim(0, 0.7)
+ axS.set_ylim(0, 250)
+ axS.text(0.3, 52, f"max S = {np.nanmax(output['products']['S_max'])-100:.2f}%")
+ axS.grid()
+
+ axT = axS.twiny()
+ axT.xaxis.label.set_color("red")
+ axT.tick_params(axis="x", colors="red")
+ axT.plot(output["products"]["T"], output["products"]["z"], color="red")
+ rng = (272, 274)
+ axT.set_xlim(*rng)
+ axT.set_xticks(np.linspace(*rng, num=5))
+ axT.set_xlabel("Temperature [K]")
+
+ axR = axs[1]
+ axR.set_xscale("log")
+ axR.set_xlim(1e-2, 1e2)
+ for drop_id, volume in enumerate(output["attributes"]["volume"]):
+ axR.plot(
+ settings.formulae.trivia.radius(volume=np.asarray(volume)) / si.um,
+ output["products"]["z"],
+ color="magenta" if drop_id < settings.n_sd_per_mode[0] else "blue",
+ label="mode 1"
+ if drop_id == 0
+ else "mode 2"
+ if drop_id == settings.n_sd_per_mode[0]
+ else "",
+ )
+ axR.legend(loc="upper right")
+ axR.set_xlabel("Droplet radius [μm]")
diff --git a/examples/PySDM_examples/Pyrcel/tutorial_settings.py b/examples/PySDM_examples/Pyrcel/tutorial_settings.py
new file mode 100644
index 000000000..12838a2e3
--- /dev/null
+++ b/examples/PySDM_examples/Pyrcel/tutorial_settings.py
@@ -0,0 +1,68 @@
+from typing import Dict
+
+import numpy as np
+from pystrict import strict
+
+from PySDM import Formulae
+from PySDM.initialisation.impl.spectrum import Spectrum
+
+
+@strict
+class Settings:
+ def __init__(
+ self,
+ dz: float,
+ n_sd_per_mode: tuple,
+ aerosol_modes_by_kappa: Dict[float, Spectrum],
+ vertical_velocity: float,
+ initial_temperature: float,
+ initial_pressure: float,
+ initial_relative_humidity: float,
+ displacement: float,
+ formulae: Formulae,
+ ):
+ self.formulae = formulae
+ self.n_sd_per_mode = n_sd_per_mode
+ self.aerosol_modes_by_kappa = aerosol_modes_by_kappa
+
+ const = self.formulae.constants
+ self.vertical_velocity = vertical_velocity
+ self.initial_pressure = initial_pressure
+ self.initial_temperature = initial_temperature
+ pv0 = (
+ initial_relative_humidity
+ * formulae.saturation_vapour_pressure.pvs_Celsius(
+ initial_temperature - const.T0
+ )
+ )
+ self.initial_vapour_mixing_ratio = const.eps * pv0 / (initial_pressure - pv0)
+ self.t_max = displacement / vertical_velocity
+ self.timestep = dz / vertical_velocity
+ self.output_interval = self.timestep
+
+ @property
+ def initial_air_density(self):
+ const = self.formulae.constants
+ dry_air_density = (
+ self.formulae.trivia.p_d(
+ self.initial_pressure, self.initial_vapour_mixing_ratio
+ )
+ / self.initial_temperature
+ / const.Rd
+ )
+ return dry_air_density * (1 + self.initial_vapour_mixing_ratio)
+
+ @property
+ def nt(self) -> int:
+ nt = self.t_max / self.timestep
+ nt_int = round(nt)
+ np.testing.assert_almost_equal(nt, nt_int)
+ return nt_int
+
+ @property
+ def steps_per_output_interval(self) -> int:
+ return int(self.output_interval / self.timestep)
+
+ @property
+ def output_steps(self) -> np.ndarray:
+ return np.arange(0, self.nt + 1, self.steps_per_output_interval)
diff --git a/examples/PySDM_examples/Pyrcel/tutorial_simulation.py b/examples/PySDM_examples/Pyrcel/tutorial_simulation.py
new file mode 100644
index 000000000..feb17c72b
--- /dev/null
+++ b/examples/PySDM_examples/Pyrcel/tutorial_simulation.py
@@ -0,0 +1,90 @@
+import numpy as np
+from PySDM_examples.utils import BasicSimulation
+
+from PySDM import Builder
+from PySDM.backends import CPU
+from PySDM.backends.impl_numba.test_helpers import scipy_ode_condensation_solver
+from PySDM.dynamics import AmbientThermodynamics, Condensation
+from PySDM.environments import Parcel
+from PySDM.initialisation import equilibrate_wet_radii
+from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity
+from PySDM.physics import si
+
+
+class Simulation(BasicSimulation):
+ def __init__(
+ self, settings, products=None, scipy_solver=False, rtol_thd=1e-10, rtol_x=1e-10
+ ):
+ env = Parcel(
+ dt=settings.timestep,
+ p0=settings.initial_pressure,
+ initial_water_vapour_mixing_ratio=settings.initial_vapour_mixing_ratio,
+ T0=settings.initial_temperature,
+ w=settings.vertical_velocity,
+ mass_of_dry_air=44 * si.kg,
+ )
+ n_sd = sum(settings.n_sd_per_mode)
+ builder = Builder(n_sd=n_sd, backend=CPU(formulae=settings.formulae))
+ builder.set_environment(env)
+ builder.add_dynamic(AmbientThermodynamics())
+ builder.add_dynamic(Condensation(rtol_thd=rtol_thd, rtol_x=rtol_x))
+
+ volume = env.mass_of_dry_air / settings.initial_air_density
+ attributes = {
+ k: np.empty(0) for k in ("dry volume", "kappa times dry volume", "n")
+ }
+ for i, (kappa, spectrum) in enumerate(settings.aerosol_modes_by_kappa.items()):
+ sampling = ConstantMultiplicity(spectrum)
+ r_dry, n_per_volume = sampling.sample(settings.n_sd_per_mode[i])
+ v_dry = settings.formulae.trivia.volume(radius=r_dry)
+ attributes["n"] = np.append(attributes["n"], n_per_volume * volume)
+ attributes["dry volume"] = np.append(attributes["dry volume"], v_dry)
+ attributes["kappa times dry volume"] = np.append(
+ attributes["kappa times dry volume"], v_dry * kappa
+ )
+ r_wet = equilibrate_wet_radii(
+ r_dry=settings.formulae.trivia.radius(volume=attributes["dry volume"]),
+ environment=env,
+ kappa_times_dry_volume=attributes["kappa times dry volume"],
+ )
+ attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet)
+
+ super().__init__(
+ particulator=builder.build(attributes=attributes, products=products)
+ )
+ if scipy_solver:
+ scipy_ode_condensation_solver.patch_particulator(self.particulator)
+
+ self.output_attributes = {
+ "volume": tuple([] for _ in range(self.particulator.n_sd))
+ }
+ self.settings = settings
+
+ self.__sanity_checks(attributes, volume)
+
+ def __sanity_checks(self, attributes, volume):
+ for attribute in attributes.values():
+ assert attribute.shape[0] == self.particulator.n_sd
+ np.testing.assert_approx_equal(
+ sum(attributes["multiplicity"]) / volume,
+ sum(
+ mode.norm_factor
+ for mode in self.settings.aerosol_modes_by_kappa.values()
+ ),
+ significant=4,
+ )
+
+ def _save(self, output):
+ for key, attr in self.output_attributes.items():
+ attr_data = self.particulator.attributes[key].to_ndarray()
+ for drop_id in range(self.particulator.n_sd):
+ attr[drop_id].append(attr_data[drop_id])
+ super()._save(output)
+
+ def run(self, observers=()):
+ for observer in observers:
+ self.particulator.observers.append(observer)
+ output_products = super()._run(
+ self.settings.nt, self.settings.steps_per_output_interval
+ )
+ return {"products": output_products, "attributes": self.output_attributes}
diff --git a/examples/PySDM_examples/Shima_et_al_2009/tutorial_example.py b/examples/PySDM_examples/Shima_et_al_2009/tutorial_example.py
new file mode 100644
index 000000000..123cc676b
--- /dev/null
+++ b/examples/PySDM_examples/Shima_et_al_2009/tutorial_example.py
@@ -0,0 +1,42 @@
+from PySDM.backends import CPU
+from PySDM.builder import Builder
+from PySDM.dynamics import Coalescence
+from PySDM.environments import Box
+from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity
+from PySDM.products import ParticleVolumeVersusRadiusLogarithmSpectrum, WallTime
+
+
+def run(settings, observers=()):
+ builder = Builder(n_sd=settings.n_sd, backend=CPU(formulae=settings.formulae))
+ builder.set_environment(Box(dv=settings.dv, dt=settings.dt))
+ attributes = {}
+ sampling = ConstantMultiplicity(settings.spectrum)
+ attributes["volume"], attributes["n"] = sampling.sample(settings.n_sd)
+ coalescence = Coalescence(
+ collision_kernel=settings.kernel, adaptive=settings.adaptive
+ )
+ builder.add_dynamic(coalescence)
+ products = (
+ ParticleVolumeVersusRadiusLogarithmSpectrum(
+ settings.radius_bins_edges, name="dv/dlnr"
+ ),
+ WallTime(),
+ )
+ particulator = builder.build(attributes, products)
+ if hasattr(settings, "u_term") and "terminal velocity" in particulator.attributes:
+ particulator.attributes["terminal velocity"].approximation = settings.u_term(
+ particulator
+ )
+
+ for observer in observers:
+ particulator.observers.append(observer)
+
+ vals = {}
+ particulator.products["wall time"].reset()
+ for step in settings.output_steps:
+ particulator.run(step - particulator.n_steps)
+ vals[step] = particulator.products["dv/dlnr"].get()[0]
+ vals[step][:] *= settings.rho
+
+ exec_time = particulator.products["wall time"].get()
+ return vals, exec_time
diff --git a/examples/PySDM_examples/Shima_et_al_2009/tutorial_plotter.py b/examples/PySDM_examples/Shima_et_al_2009/tutorial_plotter.py
new file mode 100644
index 000000000..ec4f239ab
--- /dev/null
+++ b/examples/PySDM_examples/Shima_et_al_2009/tutorial_plotter.py
@@ -0,0 +1,129 @@
+import numpy as np
+from matplotlib import pyplot
+from open_atmos_jupyter_utils import show_plot
+
+from PySDM.physics.constants import si
+
+
+class SpectrumColors:
+ def __init__(self, begining="#2cbdfe", end="#b317b1"):
+ self.b = begining
+ self.e = end
+
+ def __call__(self, value: float):
+ bR, bG, bB = int(self.b[1:3], 16), int(self.b[3:5], 16), int(self.b[5:7], 16)
+ eR, eG, eB = int(self.e[1:3], 16), int(self.e[3:5], 16), int(self.e[5:7], 16)
+ R = bR + int((eR - bR) * value)
+ G = bG + int((eG - bG) * value)
+ B = bB + int((eB - bB) * value)
+ result = f"#{hex(R)[2:4]}{hex(G)[2:4]}{hex(B)[2:4]}"
+ return result
+
+
+class SpectrumPlotter:
+ def __init__(self, settings, grid=True, legend=True, log_base=10):
+ self.settings = settings
+ self.format = "pdf"
+ self.colors = SpectrumColors()
+ self.smooth = False
+ self.smooth_scope = 2
+ self.legend = legend
+ self.grid = grid
+ self.xlabel = "particle radius [µm]"
+ self.ylabel = "dm/dlnr [g/m$^3$]"
+ self.log_base = log_base
+ self.ax = pyplot
+ self.fig = pyplot
+ self.finished = False
+
+ def finish(self):
+ if self.finished:
+ return
+ self.finished = True
+ if self.grid:
+ self.ax.grid()
+
+ self.ax.xscale("log")
+ self.ax.xlabel(self.xlabel)
+ self.ax.ylabel(self.ylabel)
+ if self.legend:
+ self.ax.legend()
+
+ def show(self):
+ self.finish()
+ pyplot.tight_layout()
+ show_plot()
+
+ def save(self, file):
+ self.finish()
+ pyplot.savefig(file, format=self.format)
+
+ def plot(self, spectrum, t):
+ self.plot_analytic_solution(self.settings, t)
+ self.plot_data(self.settings, t, spectrum)
+
+ def plot_analytic_solution(self, settings, t):
+ def analytic_solution(x):
+ return settings.norm_factor * settings.kernel.analytic_solution(
+ x=x, t=t, x_0=settings.X0, N_0=settings.n_part
+ )
+
+ if t == 0:
+ analytic_solution = settings.spectrum.size_distribution
+
+ volume_bins_edges = self.settings.formulae.trivia.volume(
+ settings.radius_bins_edges
+ )
+ dm = np.diff(volume_bins_edges)
+ dr = np.diff(settings.radius_bins_edges)
+
+ pdf_m_x = volume_bins_edges[:-1] + dm / 2
+ pdf_m_y = analytic_solution(pdf_m_x)
+
+ pdf_r_x = settings.radius_bins_edges[:-1] + dr / 2
+ pdf_r_y = pdf_m_y * dm / dr * pdf_r_x
+
+ x = pdf_r_x * si.metres / si.micrometres
+ y_true = (
+ pdf_r_y
+ * self.settings.formulae.trivia.volume(radius=pdf_r_x)
+ * settings.rho
+ / settings.dv
+ * si.kilograms
+ / si.grams
+ )
+
+ self.ax.plot(x, y_true, color="black")
+
+ def plot_data(self, settings, t, spectrum):
+ if self.smooth:
+ scope = self.smooth_scope
+ if t != 0:
+ new = np.copy(spectrum)
+ for _ in range(2):
+ for i in range(scope, len(spectrum) - scope):
+ new[i] = np.mean(spectrum[i - scope : i + scope + 1])
+ scope = 1
+ for i in range(scope, len(spectrum) - scope):
+ spectrum[i] = np.mean(new[i - scope : i + scope + 1])
+
+ x = settings.radius_bins_edges[:-scope]
+ dx = np.diff(x)
+ self.ax.plot(
+ (x[:-1] + dx / 2) * si.metres / si.micrometres,
+ spectrum[:-scope] * si.kilograms / si.grams,
+ label=f"t = {t}s",
+ color=self.colors(
+ t / (self.settings.output_steps[-1] * self.settings.dt)
+ ),
+ )
+ else:
+ self.ax.step(
+ settings.radius_bins_edges[:-1] * si.metres / si.micrometres,
+ spectrum * si.kilograms / si.grams,
+ where="post",
+ label=f"t = {t}s",
+ color=self.colors(
+ t / (self.settings.output_steps[-1] * self.settings.dt)
+ ),
+ )
diff --git a/examples/PySDM_examples/Shima_et_al_2009/tutorial_settings.py b/examples/PySDM_examples/Shima_et_al_2009/tutorial_settings.py
new file mode 100644
index 000000000..6343e518d
--- /dev/null
+++ b/examples/PySDM_examples/Shima_et_al_2009/tutorial_settings.py
@@ -0,0 +1,35 @@
+from typing import Optional
+
+import numpy as np
+from pystrict import strict
+
+from PySDM import Formulae
+from PySDM.dynamics.collisions.collision_kernels import Golovin
+from PySDM.initialisation import spectra
+from PySDM.physics import si
+
+
+@strict
+class Settings:
+ def __init__(self, steps: Optional[list] = None):
+ steps = steps or [0, 1200, 2400, 3600]
+ self.formulae = Formulae()
+ self.n_sd = 2**13
+ self.n_part = 2**23 / si.metre**3
+ self.X0 = self.formulae.trivia.volume(radius=30.531 * si.micrometres)
+ self.dv = 1e6 * si.metres**3
+ self.norm_factor = self.n_part * self.dv
+ self.rho = 1000 * si.kilogram / si.metre**3
+ self.dt = 1 * si.seconds
+ self.adaptive = False
+ self.seed = 44
+ self.steps = steps
+ self.kernel = Golovin(b=1.5e3 / si.second)
+ self.spectrum = spectra.Exponential(norm_factor=self.norm_factor, scale=self.X0)
+ self.radius_bins_edges = np.logspace(
+ np.log10(10 * si.um), np.log10(5e4 * si.um), num=256, endpoint=True
+ )
+
+ @property
+ def output_steps(self):
+ return [int(step / self.dt) for step in self.steps]
diff --git a/tests/tutorials_tests/__init__.py b/tests/tutorials_tests/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/tutorials_tests/conftest.py b/tests/tutorials_tests/conftest.py
new file mode 100644
index 000000000..f6c0d7a1c
--- /dev/null
+++ b/tests/tutorials_tests/conftest.py
@@ -0,0 +1,17 @@
+# pylint: disable=missing-module-docstring
+import pathlib
+
+import pytest
+
+from ..examples_tests.conftest import findfiles
+
+PYSDM_TUTORIALS_ABS_PATH = (
+ pathlib.Path(__file__).parent.parent.parent.absolute().joinpath("tutorials")
+)
+
+
+@pytest.fixture(
+ params=(path for path in findfiles(PYSDM_TUTORIALS_ABS_PATH, r".*\.(ipynb)$")),
+)
+def notebook_filename(request):
+ return request.param
diff --git a/tests/tutorials_tests/test_run_notebooks.py b/tests/tutorials_tests/test_run_notebooks.py
new file mode 100644
index 000000000..026173cff
--- /dev/null
+++ b/tests/tutorials_tests/test_run_notebooks.py
@@ -0,0 +1,6 @@
+# pylint: disable=missing-module-docstring
+from ..devops_tests.test_notebooks import test_run_notebooks as _impl
+
+
+def test_run_notebooks(notebook_filename, tmp_path):
+ _impl(notebook_filename, tmp_path)
diff --git a/tutorials/collisions/__init__.py b/tutorials/collisions/__init__.py
new file mode 100644
index 000000000..44b071d42
--- /dev/null
+++ b/tutorials/collisions/__init__.py
@@ -0,0 +1,5 @@
+# pylint: disable=invalid-name
+"""
+collision-coalescence tutorial following setup from
+[Shima et al. 2009](https://doi.org/10.1002/qj.441)
+"""
diff --git a/tutorials/collisions/collection_droplet.svg b/tutorials/collisions/collection_droplet.svg
new file mode 100644
index 000000000..d8bac55fc
--- /dev/null
+++ b/tutorials/collisions/collection_droplet.svg
@@ -0,0 +1,169 @@
+
diff --git a/tutorials/collisions/collisions_playground.ipynb b/tutorials/collisions/collisions_playground.ipynb
new file mode 100644
index 000000000..177b09f45
--- /dev/null
+++ b/tutorials/collisions/collisions_playground.ipynb
@@ -0,0 +1,241 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "[](https://github.com/open-atmos/PySDM/tree/main/tutorials/collisions/collisions_playground.ipynb) \n",
+ "[](https://mybinder.org/v2/gh/open-atmos/PySDM/tree/main/tutorials/collisions/collisions_playground.ipynb) \n",
+ "[](https://colab.research.google.com/github/open-atmos/PySDM/tree/main/tutorials/collisions/collisions_playground.ipynb)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Cloud Microphysics: Part 2\n",
+ "- Collisions and coalescence of cloud droplets\n",
+ "\n",
+ "Based on Fig. 2 from Shima et al. 2009 (Q. J. R. Meteorol. Soc. 135) \"_The super‐droplet method for the numerical simulation of clouds and precipitation: a particle‐based and probabilistic microphysics model coupled with a non‐hydrostatic model_.\" \n",
+ "https://doi.org/10.1002/qj.441"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Collision/coalescence \n",
+ "\n",
+ "The process of droplets colliding and coalescence, together refered to as collection, is the mechanism by which cloud droplets grow and eventually grow large enough to precipitate.\n",
+ "The collection process depends on these two processes, first two droplets colliding, and second that collision resulting in the coalescence of a new larger droplet.\n",
+ "\n",
+ "In models we parameterize this collection process stochastically by solving what is known as the SCE: Stochastic Collection Equation.\n",
+ "And we write the probability that two droplets collide (collision rate) in terms of a \"kernel\": $K(x,y)$, where $x$ and $y$ are the sizes of the two droplets.\n",
+ "\n",
+ "In this example, we consider the most basic kernel called the Golovin kernel, which is a linear kernel of the form $K(x,y) = b(x+y)$.\n",
+ "\n",
+ "Below is a drawing from Lamb and Verlinde's \"_The Physics and Chemistry of Clouds_\" illustrating the geometry of droplet collisions.\n",
+ "\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## PySDM box model widget\n",
+ "\n",
+ "In this homework assignment, and with this `PySDM` example notebook, you have the chance to explore how particle size, collision kernel, and spectral resolution (number of superdroplets used to represent the droplet size distribution) influence the growth of a population of cloud droplets as they undergo collision and coalescence. \n",
+ "\n",
+ "In this box setup, we can focus on only the collision/coalescence process while ignoring the hygroscopic growth of particles and activation of aerosols considered in Part 1, as well as fluid flow and mixing from a 2D or 3D simulation."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2023-10-16T10:50:29.614596281Z",
+ "start_time": "2023-10-16T10:50:29.474658084Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "import sys\n",
+ "if 'google.colab' in sys.modules:\n",
+ " !pip --quiet install \"open_atmos_jupyter_utils\"\n",
+ " from open_atmos_jupyter_utils import pip_install_on_colab\n",
+ " pip_install_on_colab('PySDM')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2023-10-16T10:50:34.479445450Z",
+ "start_time": "2023-10-16T10:50:29.475070307Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "from numpy import errstate\n",
+ "import os\n",
+ "\n",
+ "from PySDM import Formulae\n",
+ "from PySDM.dynamics.collisions.collision_kernels import Golovin\n",
+ "from PySDM.initialisation import spectra\n",
+ "from PySDM.physics import si\n",
+ "\n",
+ "from PySDM_examples.utils import widgets\n",
+ "\n",
+ "from PySDM_examples.Shima_et_al_2009.tutorial_plotter import SpectrumPlotter\n",
+ "from PySDM_examples.Shima_et_al_2009.tutorial_settings import Settings\n",
+ "from PySDM_examples.Shima_et_al_2009.tutorial_example import run"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2023-10-16T10:50:34.608979949Z",
+ "start_time": "2023-10-16T10:50:34.508186707Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "def demo(*, _freezer, _n, _b, _r, _smooth):\n",
+ " frm = Formulae()\n",
+ " with _freezer:\n",
+ " with errstate(all='raise'):\n",
+ " n_step = 3600\n",
+ " n_plot = 3\n",
+ " settings = Settings(steps=[i * (n_step // n_plot) for i in range(n_plot + 1)])\n",
+ " settings.n_sd = 2 ** _n\n",
+ " settings.adaptive = True\n",
+ " settings.dt = 10\n",
+ " settings.kernel = Golovin(b=_b / si.second)\n",
+ " settings.X0 = frm.trivia.volume(radius=_r * si.micrometres)\n",
+ " settings.spectrum = spectra.Exponential(\n",
+ " norm_factor=settings.norm_factor, scale=settings.X0\n",
+ " )\n",
+ " states, _ = run(settings, (widgets.ProgbarUpdater(progbar, settings.output_steps[-1]),))\n",
+ "\n",
+ " with errstate(invalid='ignore'):\n",
+ " plotter = SpectrumPlotter(settings)\n",
+ " plotter.smooth = _smooth\n",
+ " for step, state in states.items():\n",
+ " plotter.plot(state, step * settings.dt)\n",
+ " plotter.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Widget\n",
+ "\n",
+ "Play around with the widget and change the number of superdroplets ($n_{SD} = 2^X$), slope parameter ($b$) in the Golovin collision kernel, and the scale parameter in the droplet size distribution ($r$).\n",
+ "\n",
+ "
\n",
+ "Note: Running the box model takes a few seconds, so be patient after you move one of the sliders.
\n",
+ "\n",
+ "The plot generated shows the evolution of the droplet size distribution at various time points (in color). Plotted underneath (in black) is the analytical solution, which exists for this simple Golovin collision kernel."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2023-10-16T10:50:54.426815690Z",
+ "start_time": "2023-10-16T10:50:34.522224837Z"
+ },
+ "pycharm": {
+ "is_executing": true
+ }
+ },
+ "outputs": [],
+ "source": [
+ "style = {'description_width': 'initial'}\n",
+ "n_SD = widgets.IntSlider(value=12, min=6, max=18, step=1, \n",
+ " description='log2(nSD)', continuous_update=False, style=style)\n",
+ "b = widgets.FloatSlider(value=1.5e3, min=1e3, max=2e3, step=1e2, \n",
+ " description='b (s-1)', continuous_update=False,\n",
+ " readout_format='.1e', style=style)\n",
+ "r = widgets.IntSlider(value=30, min=25, max=35, step=1, \n",
+ " description='r (um)', continuous_update=False, style=style)\n",
+ "sliders = widgets.HBox([n_SD, b, r])\n",
+ "\n",
+ "smooth = widgets.Checkbox(value=True, description='smooth plot')\n",
+ "options = [smooth]\n",
+ "boxes = widgets.HBox(options)\n",
+ "freezer = widgets.Freezer([n_SD, b, r])\n",
+ "inputs = {'_freezer': freezer, '_n': n_SD, '_b': b, '_r': r, '_smooth': smooth}\n",
+ "progbar = widgets.IntProgress(min=0, max=100, description='%')\n",
+ "\n",
+ "if 'CI' not in os.environ:\n",
+ " widgets.display(sliders, boxes, progbar, widgets.interactive_output(demo, inputs))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "pycharm": {
+ "is_executing": true
+ }
+ },
+ "source": [
+ "## Questions\n",
+ "\n",
+ "1. How does the shape of the droplet size distribution change over time as particles collide and coalesce?\n",
+ "\n",
+ "2. For the Golovin collision kernel there is an analytical solution, plotted in the black curve. \n",
+ "How many superdroplets are needed to robustly simulate droplet collection?\n",
+ "\n",
+ "3. What does the `b` parameter in the collision kernel control?\n",
+ "\n",
+ "4. How does the mean radius of the droplets affect the collision rate? "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2023-10-16T10:50:54.440313749Z",
+ "start_time": "2023-10-16T10:50:54.426253375Z"
+ }
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3.10.9 ('pysdm')",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.9"
+ },
+ "vscode": {
+ "interpreter": {
+ "hash": "b14f34a08619f4a218d80d7380beed3f0c712c89ff93e7183219752d640ed427"
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/tutorials/condensation/__init__.py b/tutorials/condensation/__init__.py
new file mode 100644
index 000000000..b393f540a
--- /dev/null
+++ b/tutorials/condensation/__init__.py
@@ -0,0 +1,5 @@
+# pylint: disable=invalid-name
+"""
+condensation tutorial based on the test case from
+[Pyrcel package docs](https://pyrcel.readthedocs.io/)
+"""
diff --git a/tutorials/condensation/condensation_playground.ipynb b/tutorials/condensation/condensation_playground.ipynb
new file mode 100644
index 000000000..9e098cf17
--- /dev/null
+++ b/tutorials/condensation/condensation_playground.ipynb
@@ -0,0 +1,333 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "[](https://github.com/open-atmos/PySDM/tree/main/tutorials/condensation/condensation_playground.ipynb) \n",
+ "[](https://mybinder.org/v2/gh/open-atmos/PySDM/tree/main/tutorials/condensation/condensation_playground.ipynb) \n",
+ "[](https://colab.research.google.com/github/open-atmos/PySDM/tree/main/tutorials/condensation/condensation_playground.ipynb)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Cloud Microphysics: Part 1\n",
+ "- Activation of aerosol particles into cloud droplets\n",
+ "- Exploring how size/composition affect condensational growth\n",
+ "\n",
+ "Based on Example Figure from Pyrcel code documentation https://pyrcel.readthedocs.io/en/latest/examples/basic_run.html"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Droplet activation \n",
+ "(for more info read Ch. 6 of Lohmann's _An Introduction to Clouds_)\n",
+ "\n",
+ "#### Köhler curve:\n",
+ "- Curvature effect (Kelvin equation), describes the increase in vapor pressure over a curved surface\n",
+ "compared to a flat surface and how this depends on the surface tension and radius of the droplet. \n",
+ "$e_s$ is the saturation vapor pressure over a surface of radius $r$, so $e_s(\\infty)$ is the \n",
+ "saturation vapor pressure over a flat surface. $\\sigma$ is the surface tension, $\\rho$ is the density\n",
+ "of the solution, $R_v$ is the gas constant for water vapor, and $T$ is the temperature.\n",
+ "\n",
+ "$e_s(r) = e_s(\\infty) \\exp \\left( \\frac{2 \\sigma}{r \\rho R_v T} \\right)$\n",
+ "\n",
+ "\n",
+ "
\n",
+ "Fun fact: Based on the curvature considerations alone, saturation ratio in the atmosphere would need to be 5-10 for water to condense homogeneously, aka it would be extremely humid! Fortunately, we have aerosols that can serve as nuclei for water vapor to condense onto, and supersaturations in Earth's atmosphere rarely exceed 1%.\n",
+ "
\n",
+ "\n",
+ "- Solute effect (Raoult's law), describes the reduction of vapor pressure over a flat surface due\n",
+ "to the presence of soluble material, aka aerosol.\n",
+ "$\\kappa$ is refered to as the hygroscopicity, defined as the inverse of the water activity ($a_w$).\n",
+ "Again, $e_s$ is the saturation vapor pressure of pure water, and now $e_*$ is the vapor pressure \n",
+ "of the solution with $n_s$ moles of solute and $n_w$ moles of water.\n",
+ "\n",
+ "$\\kappa = \\frac{1}{a_w} = \\frac{e_s(\\infty)}{e_*(\\infty)} = \\frac{n_s + n_w}{n_w}$\n",
+ "\n",
+ "The hygroscopicity (inverse of water activity) is defined as the ratio of the total number of \n",
+ "moles of solute plus water to the number of moles of water.\n",
+ "\n",
+ "- Putting it together, the Köhler curve, or $\\kappa$-Köhler curve, describes the hygroscopic \n",
+ "growth of particles, and the maximum of this curve, describes the point of activation from\n",
+ "an aerosol into a cloud droplet.\n",
+ "$S$ is the saturation ratio, which is usually linerarized as follows:\n",
+ "\n",
+ "$S(r) = \\frac{e_*(r)}{e_s(\\infty)} \\approx 1 + \\frac{a}{r} - \\frac{b}{r^3}$\n",
+ "\n",
+ "\n",
+ "Fig 6.11 from Lohmann. You can see a characteristic Köhler curve with the critical radius ($r_{act}$) and supersaturation ($S_{act}$) which separate the stable (aerosol or \"solution droplet\") and unstable (cloud droplet) regimes labeled.\n",
+ "\n",
+ "\n",
+ "
\n",
+ "Other considerations: Surface tension: The surface tension $\\sigma$ in the Kelvin equation is usually assumed as constant $\\sigma = \\sigma_w = 72$ mN, but complex chemistry of the aerosol can sometimes actually\n",
+ "modify the effective surface tension of the growing cloud droplet.\n",
+ "
\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## PySDM parcel model widget\n",
+ "\n",
+ "In this homework assignment, and with this `PySDM` example notebook, you have the chance to explore how particle size, number concentration, and chemical composition, influence the bulk properties of a cloud by using a parcel model.\n",
+ "\n",
+ "A parcel model takes a parcel of air and lifts it adiabatically, greatly simplifying the dynamics taking place in a real cloud, but resolving the microphysical processes we are interested in exploring here.\n",
+ "\n",
+ "We initialize this parcel with a bimodal aerosol composed of two lognormal modes. The first mode is fixed, while the widget will let you play with the properties of the second mode. The default configuration represents a typical case in a marine environment. The first mode is smaller, more numerous sulfate aerosol, and the second mode is larger radii, less numerous, highly hygroscopic sea salt aerosol. \n",
+ "\n",
+ "You can play around with the widget at the bottom to change the initial aerosol properties, while keeping the dynamics fixed (i.e. updraft velocity `w = 1 * si.m / si.s` or temperature `T0 = 274 * si.K`). "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2023-10-16T10:48:42.215442231Z",
+ "start_time": "2023-10-16T10:48:42.165173760Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "# import PySDM library on google colab\n",
+ "import sys\n",
+ "if 'google.colab' in sys.modules:\n",
+ " !pip --quiet install \"open-atmos-jupyter-utils\"\n",
+ " from open_atmos_jupyter_utils import pip_install_on_colab\n",
+ " pip_install_on_colab('PySDM-examples')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2023-10-16T10:48:46.446157420Z",
+ "start_time": "2023-10-16T10:48:42.168839122Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "# import functions for creating interactive widget\n",
+ "from PySDM_examples.utils import widgets\n",
+ "import numpy as np\n",
+ "from numpy import errstate\n",
+ "import os\n",
+ "\n",
+ "# import PySDM tools for initializing and running a cloud parcel model\n",
+ "from PySDM import Formulae\n",
+ "from PySDM.physics import si\n",
+ "from PySDM.initialisation.spectra import Lognormal\n",
+ "from PySDM.products import (\n",
+ " ParcelDisplacement, AmbientTemperature, AmbientRelativeHumidity,\n",
+ " ParticleSizeSpectrumPerVolume, ParticleVolumeVersusRadiusLogarithmSpectrum\n",
+ ")\n",
+ "\n",
+ "# import tools for running and plotting this tutorial\n",
+ "from PySDM_examples.Pyrcel.tutorial_settings import Settings\n",
+ "from PySDM_examples.Pyrcel.tutorial_simulation import Simulation\n",
+ "from PySDM_examples.Pyrcel.profile_plotter import ProfilePlotter"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2023-10-16T10:48:46.455145740Z",
+ "start_time": "2023-10-16T10:48:46.451616736Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "# create progress bar for widget\n",
+ "progbar = widgets.IntProgress(min=0, max=100, description='%')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2023-10-16T10:48:46.504641095Z",
+ "start_time": "2023-10-16T10:48:46.458931372Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "# create initial aerosol distribution\n",
+ "# run cloud parcel model\n",
+ "# save and plot results\n",
+ "\n",
+ "# k2, N2, and r2 are the hygroscopicity, number concentration, and mean radius\n",
+ "# of the second Lognormal aerosol mode \n",
+ "def demo(*, _freezer, _k2, _N2, _r2):\n",
+ " with _freezer:\n",
+ " with errstate(all='raise'):\n",
+ " settings = Settings(\n",
+ " dz = 1 * si.m,\n",
+ " n_sd_per_mode = (10, 10),\n",
+ " aerosol_modes_by_kappa = {\n",
+ " .54: Lognormal(\n",
+ " norm_factor=850 / si.cm ** 3,\n",
+ " m_mode=15 * si.nm,\n",
+ " s_geom=1.6\n",
+ " ),\n",
+ " _k2: Lognormal(\n",
+ " norm_factor=_N2 / si.cm ** 3,\n",
+ " m_mode=_r2 * si.nm,\n",
+ " s_geom=1.2\n",
+ " )\n",
+ " },\n",
+ " vertical_velocity = 1.0 * si.m / si.s,\n",
+ " initial_pressure = 775 * si.mbar,\n",
+ " initial_temperature = 274 * si.K,\n",
+ " initial_relative_humidity = 0.98,\n",
+ " displacement = 250 * si.m,\n",
+ " formulae = Formulae(constants={'MAC': .3})\n",
+ " )\n",
+ " dry_radius_bin_edges = np.logspace(\n",
+ " np.log10(1e-3 * si.um),\n",
+ " np.log10(5e0 * si.um),\n",
+ " 33, endpoint=False\n",
+ " )\n",
+ " simulation = Simulation(\n",
+ " settings,\n",
+ " products=(\n",
+ " ParcelDisplacement(\n",
+ " name='z'),\n",
+ " AmbientRelativeHumidity(\n",
+ " name='S_max', unit='%', var='RH'),\n",
+ " AmbientTemperature(\n",
+ " name='T'),\n",
+ " ParticleSizeSpectrumPerVolume(\n",
+ " name='dry:dN/dR', radius_bins_edges=dry_radius_bin_edges, dry=True),\n",
+ " ParticleVolumeVersusRadiusLogarithmSpectrum(\n",
+ " name='dry:dV/dlnR', radius_bins_edges=dry_radius_bin_edges, dry=True),\n",
+ " ),\n",
+ " scipy_solver=False,\n",
+ " )\n",
+ " output = simulation.run((widgets.ProgbarUpdater(progbar, settings.output_steps[-1]),))\n",
+ "\n",
+ " with errstate(invalid='ignore'):\n",
+ " plotter = ProfilePlotter(settings)\n",
+ " plotter.plot(output)\n",
+ " plotter.show()\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Widget\n",
+ "\n",
+ "Play around with the widget and change the hygroscopicity ($\\kappa_2$), number concentration ($N_2$), and mean radius ($r_2$) of the second (\"sea salt\") mode. \n",
+ "\n",
+ "
\n",
+ "Note: Running the parcel model takes a few seconds, so be patient after you move one of the sliders.
\n",
+ "\n",
+ "The plots generated show (on the left) the profile of supersaturation ($S-1$, black) and temperature ($T$, red) and (on the right) profiles of droplet radius for each super particle. In pink are particles from the first mode (sulfate) and in blue are particles from the second mode (in the default case, sea salt)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2023-10-16T10:49:25.298263790Z",
+ "start_time": "2023-10-16T10:48:46.504445432Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "# create widget\n",
+ "# use to explore how the hygroscopicity, number concentration, and mean radius\n",
+ "# of the initial aerosol distribution affect the bulk parcel properties\n",
+ "# like maximum supersaturation and temperature profile\n",
+ "\n",
+ "style = {'description_width': 'initial'}\n",
+ "k2 = widgets.FloatSlider(value=1.2, min=0.2, max=1.4, step=0.1, description='κ2',\n",
+ " continuous_update=False, readout_format='.1f', style=style)\n",
+ "N2 = widgets.IntSlider(value=10, min=5, max=50, step=5, description='N2 (cm-3)',\n",
+ " continuous_update=False, style=style)\n",
+ "r2 = widgets.IntSlider(value=850, min=200, max=1000, step=50, description='r2 (nm)',\n",
+ " continuous_update=False, style=style)\n",
+ "sliders = widgets.HBox([k2, N2, r2])\n",
+ "freezer = widgets.Freezer([k2, N2, r2])\n",
+ "inputs = {'_freezer': freezer, '_k2': k2, '_N2': N2, '_r2': r2}\n",
+ "\n",
+ "if 'CI' not in os.environ:\n",
+ " widgets.display(sliders, progbar, widgets.interactive_output(demo, inputs))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Questions\n",
+ "\n",
+ "1. Extremes: Which combination of (changeable) parameters leads to the largest maximum supersaturation? Which to the smallest? Why?\n",
+ "\n",
+ "2. Sensitivity: Is the cloud more sensitive to changes in aerosol size, number, or composition? Explain how you are measuring this.\n",
+ "\n",
+ "3. Albedo: The albedo of a cloud is very dependent on the size of the individual droplets it is composed of. \n",
+ "We can express the cloud albedo ($\\alpha$) in terms of the cloud optical thickness ($\\tau$) and a dimensionless asymmetry parameter ($g$)\n",
+ "that describes the relative amount of radiation that is forward-scattered vs. backscattered. \n",
+ "$$\\alpha = \\frac{(1-g) \\tau}{2 + (1-g)\\tau}$$\n",
+ "Cloud droplets (order 1-10$\\mu$m) tend to be strongly forward-scattering with an asymmetry parameter around $g=0.85$.\n",
+ "The cloud optical thickness can be written in terms of the liquid water path through the cloud (LWP) and effective radius of the droplets ($r_e$).\n",
+ "$$\\tau = \\frac{3 LWP}{2 \\rho_w r_e}$$\n",
+ "\n",
+ "
Write down an expression for the cloud albedo. Assuming a fixed liquid water path, what is the sensitivity of albedo to droplet effective radius? This sensitivity is known as the \"Twomey effect.\"
\n",
+ "
Describe how the albedo would change given changes to the initial aerosol size distribution.
\n",
+ "\n",
+ "\n",
+ "4. Real-world pollution regulations: How would you expect temperatures in Los Angeles to change if a policy was implemented that cut pollution in half. You can assume that this policy will also reduce the number of aerosols that can serve as cloud nuclei in half. Qualitative answers are good."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2023-10-16T10:49:25.298929952Z",
+ "start_time": "2023-10-16T10:49:25.296187349Z"
+ }
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3.10.9 ('pysdm')",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.9"
+ },
+ "vscode": {
+ "interpreter": {
+ "hash": "b14f34a08619f4a218d80d7380beed3f0c712c89ff93e7183219752d640ed427"
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/tutorials/condensation/kohler_curve.svg b/tutorials/condensation/kohler_curve.svg
new file mode 100644
index 000000000..4791fd6b0
--- /dev/null
+++ b/tutorials/condensation/kohler_curve.svg
@@ -0,0 +1,242 @@
+
![](\"./collection_droplet.svg\")
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## PySDM box model widget\n",
+ "\n",
+ "In this homework assignment, and with this `PySDM` example notebook, you have the chance to explore how particle size, collision kernel, and spectral resolution (number of superdroplets used to represent the droplet size distribution) influence the growth of a population of cloud droplets as they undergo collision and coalescence. \n",
+ "\n",
+ "In this box setup, we can focus on only the collision/coalescence process while ignoring the hygroscopic growth of particles and activation of aerosols considered in Part 1, as well as fluid flow and mixing from a 2D or 3D simulation."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2023-10-16T10:50:29.614596281Z",
+ "start_time": "2023-10-16T10:50:29.474658084Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "import sys\n",
+ "if 'google.colab' in sys.modules:\n",
+ " !pip --quiet install \"open_atmos_jupyter_utils\"\n",
+ " from open_atmos_jupyter_utils import pip_install_on_colab\n",
+ " pip_install_on_colab('PySDM')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2023-10-16T10:50:34.479445450Z",
+ "start_time": "2023-10-16T10:50:29.475070307Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "from numpy import errstate\n",
+ "import os\n",
+ "\n",
+ "from PySDM import Formulae\n",
+ "from PySDM.dynamics.collisions.collision_kernels import Golovin\n",
+ "from PySDM.initialisation import spectra\n",
+ "from PySDM.physics import si\n",
+ "\n",
+ "from PySDM_examples.utils import widgets\n",
+ "\n",
+ "from PySDM_examples.Shima_et_al_2009.tutorial_plotter import SpectrumPlotter\n",
+ "from PySDM_examples.Shima_et_al_2009.tutorial_settings import Settings\n",
+ "from PySDM_examples.Shima_et_al_2009.tutorial_example import run"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2023-10-16T10:50:34.608979949Z",
+ "start_time": "2023-10-16T10:50:34.508186707Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "def demo(*, _freezer, _n, _b, _r, _smooth):\n",
+ " frm = Formulae()\n",
+ " with _freezer:\n",
+ " with errstate(all='raise'):\n",
+ " n_step = 3600\n",
+ " n_plot = 3\n",
+ " settings = Settings(steps=[i * (n_step // n_plot) for i in range(n_plot + 1)])\n",
+ " settings.n_sd = 2 ** _n\n",
+ " settings.adaptive = True\n",
+ " settings.dt = 10\n",
+ " settings.kernel = Golovin(b=_b / si.second)\n",
+ " settings.X0 = frm.trivia.volume(radius=_r * si.micrometres)\n",
+ " settings.spectrum = spectra.Exponential(\n",
+ " norm_factor=settings.norm_factor, scale=settings.X0\n",
+ " )\n",
+ " states, _ = run(settings, (widgets.ProgbarUpdater(progbar, settings.output_steps[-1]),))\n",
+ "\n",
+ " with errstate(invalid='ignore'):\n",
+ " plotter = SpectrumPlotter(settings)\n",
+ " plotter.smooth = _smooth\n",
+ " for step, state in states.items():\n",
+ " plotter.plot(state, step * settings.dt)\n",
+ " plotter.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Widget\n",
+ "\n",
+ "Play around with the widget and change the number of superdroplets ($n_{SD} = 2^X$), slope parameter ($b$) in the Golovin collision kernel, and the scale parameter in the droplet size distribution ($r$).\n",
+ "\n",
+ "![](\"./kohler_curve.svg\")
\n",
+ "Fig 6.11 from Lohmann. You can see a characteristic Köhler curve with the critical radius ($r_{act}$) and supersaturation ($S_{act}$) which separate the stable (aerosol or \"solution droplet\") and unstable (cloud droplet) regimes labeled.\n",
+ "\n",
+ "\n",
+ "