diff --git a/examples/PySDM_examples/Graf_et_al_2019/figure_4.ipynb b/examples/PySDM_examples/Graf_et_al_2019/figure_4.ipynb
index 6342f5a51a..2c490575c5 100644
--- a/examples/PySDM_examples/Graf_et_al_2019/figure_4.ipynb
+++ b/examples/PySDM_examples/Graf_et_al_2019/figure_4.ipynb
@@ -18,21 +18,17 @@
"metadata": {
"collapsed": false
},
- "source": [
- "based on [Graf et al. 2019](https://doi.org/10.5194/acp-19-747-2019)"
- ]
+ "source": "based on [Graf et al. 2019](https://doi.org/10.5194/acp-19-747-2019)"
},
{
"cell_type": "code",
- "execution_count": 1,
"id": "b2c87898f6737935",
"metadata": {
"ExecuteTime": {
- "end_time": "2024-12-06T13:01:19.373289Z",
- "start_time": "2024-12-06T13:01:19.360208Z"
+ "end_time": "2026-04-29T15:42:23.464727Z",
+ "start_time": "2026-04-29T15:42:23.461432Z"
}
},
- "outputs": [],
"source": [
"import os, sys\n",
"os.environ['NUMBA_THREADING_LAYER'] = 'workqueue' # PySDM & PyMPDATA don't work with TBB; OpenMP has extra dependencies on macOS\n",
@@ -40,25 +36,28 @@
" !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', 'PySDM')"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 1
},
{
"cell_type": "code",
- "execution_count": 13,
"id": "2432d43d15e97625",
"metadata": {
+ "collapsed": false,
"ExecuteTime": {
- "end_time": "2024-01-11T08:43:41.429069Z",
- "start_time": "2024-01-11T08:43:41.423952Z"
- },
- "collapsed": false
+ "end_time": "2026-04-29T15:42:26.604258Z",
+ "start_time": "2026-04-29T15:42:23.468474Z"
+ }
},
- "outputs": [],
"source": [
"import math\n",
"import numpy as np\n",
+ "from matplotlib import pyplot\n",
+ "from scipy.interpolate import CubicSpline\n",
+ "\n",
"from open_atmos_jupyter_utils import show_plot\n",
- "from PySDM.physics import si\n",
+ "from PySDM.physics import si, in_unit\n",
"from PySDM.environments import Parcel\n",
"from PySDM.backends import CPU\n",
"from PySDM import Builder, Formulae\n",
@@ -68,28 +67,38 @@
"from PySDM.initialisation.hygroscopic_equilibrium import equilibrate_wet_radii\n",
"from PySDM.initialisation import discretise_multiplicities\n",
"from PySDM import products"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 2
},
{
"cell_type": "code",
- "execution_count": 2,
"id": "81d815991be0cc9",
"metadata": {
+ "collapsed": false,
"ExecuteTime": {
- "end_time": "2024-01-11T08:42:35.791569Z",
- "start_time": "2024-01-11T08:42:35.352112Z"
- },
- "collapsed": false
+ "end_time": "2026-04-29T15:42:28.893670Z",
+ "start_time": "2026-04-29T15:42:26.610236Z"
+ }
},
- "outputs": [],
"source": [
- "formulae = Formulae()\n",
+ "formulae = Formulae(\n",
+ " isotope_equilibrium_fractionation_factors=\"Majoube1971+MerlivatAndNief1967+Majoube1970\",\n",
+ " isotope_kinetic_fractionation_factors=\"JouzelAndMerlivat1984\",\n",
+ " isotope_diffusivity_ratios=\"Stewart1975\",\n",
+ " isotope_meteoric_water_line=\"Dansgaard1964\"\n",
+ ")\n",
"const = formulae.constants\n",
"pvs_water = formulae.saturation_vapour_pressure.pvs_water\n",
"\n",
"RH0 = .75\n",
"p0 = 950 * si.hPa\n",
"T0 = const.T0 + 12 * si.K\n",
+ "\n",
+ "delta_2H = -150\n",
+ "R0_2H = formulae.trivia.isotopic_delta_2_ratio(delta_2H * const.PER_MILLE, const.VSMOW_R_2H)\n",
+ "R0_18O = formulae.trivia.isotopic_delta_2_ratio(formulae.isotope_meteoric_water_line.d18O_of_d2H(delta_2H * const.PER_MILLE), const.VSMOW_R_18O)\n",
+ "\n",
"alt_initial = 500 * si.m \n",
"alt_final = 3500 * si.m\n",
"\n",
@@ -99,20 +108,20 @@
"n_sd = 1\n",
"dt = 10 * si.s\n",
"vertical_velocity = 2.5 * si.m / si.s"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 3
},
{
"cell_type": "code",
- "execution_count": 3,
"id": "dada2ae5b4083f0e",
"metadata": {
+ "collapsed": false,
"ExecuteTime": {
- "end_time": "2024-01-11T08:42:35.796706Z",
- "start_time": "2024-01-11T08:42:35.793581Z"
- },
- "collapsed": false
+ "end_time": "2026-04-29T15:42:28.914742Z",
+ "start_time": "2026-04-29T15:42:28.912209Z"
+ }
},
- "outputs": [],
"source": [
"env = Parcel(\n",
" # given in the paper\n",
@@ -124,106 +133,115 @@
" dt=dt,\n",
" mass_of_dry_air=1e3 * si.kg,\n",
")"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 4
},
{
"cell_type": "code",
- "execution_count": 4,
"id": "47f3f883e4e67ed",
"metadata": {
+ "collapsed": false,
"ExecuteTime": {
- "end_time": "2024-01-11T08:42:36.874577Z",
- "start_time": "2024-01-11T08:42:35.804074Z"
- },
- "collapsed": false
+ "end_time": "2026-04-29T15:42:31.407724Z",
+ "start_time": "2026-04-29T15:42:28.925970Z"
+ }
},
- "outputs": [],
"source": [
"builder = Builder(backend=CPU(formulae=formulae), n_sd=n_sd, environment=env)\n",
"builder.add_dynamic(AmbientThermodynamics())\n",
"builder.add_dynamic(Condensation())"
- ]
+ ],
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/agnieszkazaba/PycharmProjects/PySDM/PySDM/backends/numba.py:57: UserWarning: Disabling Numba threading due to ARM64 CPU (atomics do not work yet)\n",
+ " warnings.warn(\n"
+ ]
+ }
+ ],
+ "execution_count": 5
},
{
"cell_type": "code",
- "execution_count": 5,
"id": "f98070a3d550d65",
"metadata": {
+ "collapsed": false,
"ExecuteTime": {
- "end_time": "2024-01-11T08:42:36.878312Z",
- "start_time": "2024-01-11T08:42:36.876718Z"
- },
- "collapsed": false
+ "end_time": "2026-04-29T15:42:31.417746Z",
+ "start_time": "2026-04-29T15:42:31.415357Z"
+ }
},
- "outputs": [],
"source": [
"spectrum = Lognormal(norm_factor=1e4 / si.mg, m_mode=50 * si.nm, s_geom=1.5)\n",
"kappa = .5 * si.dimensionless\n",
"cloud_range = (.5 * si.um, 25 * si.um)\n",
"output_interval = 4\n",
"output_points = 40"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 6
},
{
"cell_type": "code",
- "execution_count": 6,
"id": "5ec5eabc3c1d3872",
"metadata": {
+ "collapsed": false,
"ExecuteTime": {
- "end_time": "2024-01-11T08:42:36.883Z",
- "start_time": "2024-01-11T08:42:36.880940Z"
- },
- "collapsed": false
+ "end_time": "2026-04-29T15:42:31.423965Z",
+ "start_time": "2026-04-29T15:42:31.420878Z"
+ }
},
- "outputs": [],
"source": [
"r_dry, specific_concentration = spectral_sampling.Logarithmic(spectrum).sample_deterministic(n_sd)"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 7
},
{
"cell_type": "code",
- "execution_count": 7,
"id": "ce639f01f48ff348",
"metadata": {
+ "collapsed": false,
"ExecuteTime": {
- "end_time": "2024-01-11T08:42:37.072808Z",
- "start_time": "2024-01-11T08:42:36.886365Z"
- },
- "collapsed": false
+ "end_time": "2026-04-29T15:42:31.743451Z",
+ "start_time": "2026-04-29T15:42:31.426548Z"
+ }
},
- "outputs": [],
"source": [
"v_dry = formulae.trivia.volume(radius=r_dry)"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 8
},
{
"cell_type": "code",
- "execution_count": 8,
"id": "28897f499330f699",
"metadata": {
+ "collapsed": false,
"ExecuteTime": {
- "end_time": "2024-01-11T08:42:38.909344Z",
- "start_time": "2024-01-11T08:42:37.088206Z"
- },
- "collapsed": false
+ "end_time": "2026-04-29T15:42:34.538618Z",
+ "start_time": "2026-04-29T15:42:31.747969Z"
+ }
},
- "outputs": [],
"source": [
"r_wet = equilibrate_wet_radii(r_dry=r_dry, environment=builder.particulator.environment, kappa_times_dry_volume=kappa * v_dry)"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 9
},
{
"cell_type": "code",
- "execution_count": 9,
"id": "8cf40dfb18e7e082",
"metadata": {
+ "collapsed": false,
"ExecuteTime": {
- "end_time": "2024-01-11T08:42:38.913147Z",
- "start_time": "2024-01-11T08:42:38.909929Z"
- },
- "collapsed": false
+ "end_time": "2026-04-29T15:42:34.544696Z",
+ "start_time": "2026-04-29T15:42:34.542450Z"
+ }
},
- "outputs": [],
"source": [
"attributes = {\n",
" 'multiplicity': discretise_multiplicities(specific_concentration * builder.particulator.environment.mass_of_dry_air),\n",
@@ -231,43 +249,45 @@
" 'kappa times dry volume': kappa * v_dry,\n",
" 'volume': formulae.trivia.volume(radius=r_wet)\n",
"}"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 10
},
{
"cell_type": "code",
- "execution_count": 10,
"id": "f36faa7368db791e",
"metadata": {
+ "collapsed": false,
"ExecuteTime": {
- "end_time": "2024-01-11T08:42:39.067069Z",
- "start_time": "2024-01-11T08:42:38.920056Z"
- },
- "collapsed": false
+ "end_time": "2026-04-29T15:42:34.715111Z",
+ "start_time": "2026-04-29T15:42:34.547868Z"
+ }
},
- "outputs": [],
"source": [
"particulator = builder.build(\n",
" attributes=attributes,\n",
" products=(\n",
" products.PeakSaturation(name='S_max'),\n",
" products.ParcelDisplacement(name=\"z\"),\n",
- " products.AmbientTemperature(name=\"T\")\n",
+ " products.AmbientTemperature(name=\"T\"),\n",
+ " products.AmbientWaterVapourMixingRatio(name=\"w\", var=\"water_vapour_mixing_ratio\"),\n",
+ "\n",
" )\n",
")"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 11
},
{
"cell_type": "code",
- "execution_count": 11,
"id": "843a224f900fb94d",
"metadata": {
+ "collapsed": false,
"ExecuteTime": {
- "end_time": "2024-01-11T08:42:45.925004Z",
- "start_time": "2024-01-11T08:42:39.077019Z"
- },
- "collapsed": false
+ "end_time": "2026-04-29T15:42:53.624386Z",
+ "start_time": "2026-04-29T15:42:34.718804Z"
+ }
},
- "outputs": [],
"source": [
"n_steps = total_displacement / (dt * vertical_velocity)\n",
"\n",
@@ -290,1199 +310,167 @@
" \n",
" if levels['0C'] is None and output[\"T\"][-1] < const.T0:\n",
" levels['0C'] = .5 * (output['z'][-1] + output[\"z\"][-2])"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 12
},
{
"cell_type": "code",
- "execution_count": 12,
"id": "a363f02a36a714d",
"metadata": {
"ExecuteTime": {
- "end_time": "2024-01-11T08:42:46.314435Z",
- "start_time": "2024-01-11T08:42:45.927080Z"
- },
- "collapsed": false
+ "end_time": "2026-04-29T15:42:56.585776Z",
+ "start_time": "2026-04-29T15:42:53.636881Z"
+ }
+ },
+ "source": [
+ "z = np.asarray(output['z']) + alt_initial\n",
+ "q = formulae.trivia.mixing_ratio_to_specific_content(np.asarray(output['w']))\n",
+ "R_Rayleigh_eql = np.full_like(q, np.nan)\n",
+ "R_Rayleigh_eql[0] = R0_2H\n",
+ "R_Rayleigh_eql_18O = np.copy(R_Rayleigh_eql)\n",
+ "R_Rayleigh_eql_18O[0] = R0_18O\n",
+ "R_Rayleigh_kin = np.copy(R_Rayleigh_eql)\n",
+ "\n",
+ "\n",
+ "T = np.asarray(output['T'])\n",
+ "S_max = np.asarray(output['S_max'])\n",
+ "\n",
+ "alpha_eql_l = formulae.isotope_equilibrium_fractionation_factors.alpha_l_2H\n",
+ "alpha_eql_i = formulae.isotope_equilibrium_fractionation_factors.alpha_i_2H\n",
+ "alpha_eql_18O_l = formulae.isotope_equilibrium_fractionation_factors.alpha_l_18O\n",
+ "alpha_eql_18O_i = formulae.isotope_equilibrium_fractionation_factors.alpha_i_18O\n",
+ "alpha_kin = formulae.isotope_kinetic_fractionation_factors.alpha_kinetic\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "def fc(T):\n",
+ " T0 = 273.15 * si.K\n",
+ " T_min = 250.15 * si.K\n",
+ " if T <= T_min:\n",
+ " return 0.0\n",
+ " elif T >= T0:\n",
+ " return 1.0\n",
+ " else:\n",
+ " return CubicSpline((T_min, T0), (0.0, 1.0), bc_type=((1, 0.0), (1, 0.0)))(T)\n",
+ "R2d = formulae.trivia.isotopic_ratio_2_delta\n",
+ "for i in range(1, len(q)):\n",
+ " T_mid = (T[i]+T[i-1])/2\n",
+ " alpha_eql = fc(T_mid) * alpha_eql_l(T_mid) + (1 - fc(T_mid)) * alpha_eql_i(T_mid)\n",
+ " alpha_eql_18O = fc(T_mid) * alpha_eql_18O_l(T_mid) + (1 - fc(T_mid)) * alpha_eql_18O_i(T_mid)\n",
+ " q_ratio = q[i] / q[i-1]\n",
+ " R_Rayleigh_eql[i] = R_Rayleigh_eql[i-1] * q_ratio**(alpha_eql - 1)\n",
+ " alpha_kinetic = alpha_kin(alpha_eql, S_max[i], formulae.isotope_diffusivity_ratios.ratio_2H_heavy_to_light(T_mid))\n",
+ " R_Rayleigh_kin[i] = R_Rayleigh_kin[i-1] * q_ratio**(alpha_eql * alpha_kinetic - 1)\n",
+ " R_Rayleigh_eql_18O[i] = R_Rayleigh_eql_18O[i-1] * q_ratio**(alpha_eql_18O - 1)\n",
+ " \n",
+ "delta_2H_Rayleigh_eql = R2d(R_Rayleigh_eql, const.VSMOW_R_2H)\n",
+ "delta_2H_Rayleigh_kin = R2d(R_Rayleigh_kin, const.VSMOW_R_2H)\n",
+ "\n",
+ "d_excess = formulae.isotope_meteoric_water_line.excess_d(delta_2H_Rayleigh_eql, R2d(R_Rayleigh_eql_18O, const.VSMOW_R_18O))"
+ ],
+ "outputs": [],
+ "execution_count": 13
+ },
+ {
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2026-04-29T15:42:57.405968Z",
+ "start_time": "2026-04-29T15:42:56.592250Z"
+ }
},
+ "cell_type": "code",
+ "source": [
+ "fig, axs = pyplot.subplots(2, 2, squeeze=False, figsize=(8, 8), sharey=True)\n",
+ "xya1 = axs[0,0]\n",
+ "xya2 = xya1.twiny()\n",
+ "\n",
+ "xya1.plot(\n",
+ " in_unit(S_max, const.PER_CENT),\n",
+ " z\n",
+ ")\n",
+ "xya2.plot(T - const.T0, z)\n",
+ "\n",
+ "xya1.set_xlim(75, 107)\n",
+ "xya2.set_xlim(-10, 16)\n",
+ "xya1.set_ylim(alt_initial, alt_final)\n",
+ "xya1.set_xlabel(\"RH / %\")\n",
+ "xya2.set_xlabel(\"T / C\")\n",
+ "\n",
+ "\n",
+ "xyc = axs[1,0]\n",
+ "xlabel = \"$\\\\delta^2\\\\mathrm{H}$\"\n",
+ "for delta_2H, label, ls in zip(\n",
+ " (delta_2H_Rayleigh_eql, delta_2H_Rayleigh_kin),\n",
+ " (f\"{xlabel} eql\", f\"{xlabel} kin\"),\n",
+ " (\"-\", \"--\")\n",
+ "):\n",
+ " xyc.plot(in_unit(delta_2H, const.PER_MILLE), z, label=label, ls=ls)\n",
+ "\n",
+ "xyc.set_xlabel(f\"{xlabel} [‰]\")\n",
+ "xyc.set_xlim(-250, -50)\n",
+ "xyc.legend()\n",
+ "\n",
+ "xyd = axs[1,1]\n",
+ "xyd.plot(in_unit(d_excess, const.PER_MILLE), z)\n",
+ "xyd.set_xlabel(\"d [‰]\")\n",
+ "xyd.set_xlim(0, 20)\n",
+ "\n",
+ "for ax in axs.flat:\n",
+ " ax.set_ylabel(\"z / m\")\n",
+ " for level in levels.values():\n",
+ " ax.axhline(y=level + alt_initial, ls='--', color='black', lw=0.8)\n",
+ "\n",
+ "show_plot('fig_4.pdf')"
+ ],
+ "id": "6a4bfd9a8f959d0c",
"outputs": [
{
"data": {
- "image/svg+xml": [
- "\n",
- "\n",
- "\n"
- ],
"text/plain": [
- ""
- ]
+ ""
+ ],
+ "image/svg+xml": "\n\n\n"
},
"metadata": {},
- "output_type": "display_data"
+ "output_type": "display_data",
+ "jetTransient": {
+ "display_id": null
+ }
},
{
"data": {
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./fig_4.pdf
\"), HTML(value=\"./fig_4a.pdf
\")"
- ]
+ "version_minor": 0,
+ "model_id": "954f218c5aa1416e8fea48a7d71ef888"
+ }
},
"metadata": {},
- "output_type": "display_data"
+ "output_type": "display_data",
+ "jetTransient": {
+ "display_id": null
+ }
}
],
- "source": [
- "from matplotlib import pyplot\n",
- "\n",
- "fig, axs = pyplot.subplots(1, 1, squeeze=False)\n",
- "xy1 = axs[0,0]\n",
- "xy2 = xy1.twiny()\n",
- "\n",
- "xy1.plot(\n",
- " np.asarray(output['S_max']) * 100,\n",
- " np.asarray(output['z']) + alt_initial\n",
- ")\n",
- "xy2.plot(\n",
- " np.asarray(output['T']) - const.T0,\n",
- " np.asarray(output['z']) + alt_initial\n",
- ")\n",
- "for level in levels.values():\n",
- " xy1.axhline(y=level + alt_initial, linestyle='--', color='black')\n",
- "\n",
- "xy1.set_xlim(75, 107)\n",
- "xy2.set_xlim(-10, 16)\n",
- "xy1.set_ylim(alt_initial, alt_final)\n",
- "xy1.set_xlabel(\"RH / %\")\n",
- "xy2.set_xlabel(\"T / C\")\n",
- "xy1.set_ylabel(\"z / m\")\n",
- "show_plot('fig_4a.pdf')"
- ]
+ "execution_count": 14
},
{
- "cell_type": "code",
- "execution_count": 12,
- "id": "5acab115cd92fd9b",
"metadata": {
"ExecuteTime": {
- "end_time": "2024-01-11T08:42:46.317574Z",
- "start_time": "2024-01-11T08:42:46.313381Z"
- },
- "collapsed": false
- },
- "outputs": [],
- "source": []
- },
- {
- "cell_type": "code",
- "execution_count": 12,
- "id": "d0173c4dbd123f7b",
- "metadata": {
- "ExecuteTime": {
- "end_time": "2024-01-11T08:42:46.323642Z",
- "start_time": "2024-01-11T08:42:46.319974Z"
- },
- "collapsed": false
+ "end_time": "2026-04-29T15:42:57.417229Z",
+ "start_time": "2026-04-29T15:42:57.415402Z"
+ }
},
- "outputs": [],
- "source": []
- },
- {
"cell_type": "code",
- "execution_count": 12,
- "id": "430e8a3643291372",
- "metadata": {
- "ExecuteTime": {
- "end_time": "2024-01-11T08:42:46.326817Z",
- "start_time": "2024-01-11T08:42:46.322560Z"
- },
- "collapsed": false
- },
+ "source": "",
+ "id": "652a377137f61167",
"outputs": [],
- "source": []
+ "execution_count": null
}
],
"metadata": {
diff --git a/tests/smoke_tests/parcel_d/graf_et_al_2019/test_fig_4.py b/tests/smoke_tests/parcel_d/graf_et_al_2019/test_fig_4.py
index e8a0dfed32..a9bb85238d 100644
--- a/tests/smoke_tests/parcel_d/graf_et_al_2019/test_fig_4.py
+++ b/tests/smoke_tests/parcel_d/graf_et_al_2019/test_fig_4.py
@@ -3,11 +3,13 @@
from pathlib import Path
import pytest
+import numpy as np
from open_atmos_jupyter_utils import notebook_vars
from PySDM_examples import Graf_et_al_2019
from PySDM.physics import si
+from PySDM.physics.constants import PER_MILLE
PLOT = False
@@ -19,15 +21,46 @@ def variables_fixture():
)
-@pytest.mark.parametrize(
- "level_name, expected_height", (("CB", 1 * si.km), ("0C", 2.24 * si.km))
-)
-def test_fig_4(variables, level_name, expected_height):
- # arrange
- tolerance = 50 * si.m
+class TestFig4:
+ @staticmethod
+ @pytest.mark.parametrize(
+ "level_name, expected_height", (("CB", 1 * si.km), ("0C", 2.24 * si.km))
+ )
+ def test_fig_4(variables, level_name, expected_height):
+ # arrange
+ tolerance = 50 * si.m
+
+ # act
+ actual = variables["levels"][level_name] + variables["alt_initial"]
+
+ # assert
+ assert abs(expected_height - actual) < tolerance
+
+ @staticmethod
+ @pytest.mark.parametrize("R_name", ("R_Rayleigh_eql", "R_Rayleigh_kin"))
+ def test_fig_4c(variables, R_name):
+ # arrange
+ below_cloud = (
+ variables["z"] < variables["levels"]["CB"] + variables["alt_initial"]
+ )
+
+ # act
+ actual = variables[R_name][below_cloud]
+ desired = variables["R0_2H"]
+
+ # assert
+ np.testing.assert_array_almost_equal(
+ actual,
+ desired,
+ )
+
+ @staticmethod
+ def test_fig_4d(variables):
+ # arrange
+ desired = 10
- # act
- actual = variables["levels"][level_name] + variables["alt_initial"]
+ # act
+ actual = variables["d_excess"] / PER_MILLE
- # assert
- assert abs(expected_height - actual) < tolerance
+ # assert
+ np.testing.assert_allclose(actual, desired, rtol=4e-2)