introducing PySDM tutorials with first two based on collision and con…

…densation classroom notebooks (#1164)

Co-authored-by: claresinger <[email protected]>
Co-authored-by: Agnieszka Makulska <[email protected]>
Co-authored-by: Oleksii Bulenok <[email protected]>
Co-authored-by: Sylwester Arabas <[email protected]>

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

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

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]")

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)

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}

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