From cde2788008102a343f89cdb466706a506bda8ab3 Mon Sep 17 00:00:00 2001 From: Rafal Lukosz <93160829+lursz@users.noreply.github.com> Date: Tue, 10 Jun 2025 16:42:58 +0200 Subject: [PATCH 01/15] Alien attempt (#3) * simulation setup * uv lock * introduced all equations * oblate spheroid impl * AlienParcel class extending Parcel * PySDM package internal imports fix * Formulae bug fixes * Planet is now dataclass instead of dict * drop falling without pancake shape * cloudbase fix * figure settings * multicore support for simulation (#2) * parallel v1 * thread count fix * parallel * small parallel fixes * working multithread solution * Multi-thread + vector plots * black, black used everywhere --------- Co-authored-by: Piotr Kubala * pip install joblib * Tests (#1) * tests I guess * uv fix * planet impl * ground truth * gen plot function * unit tests, monoticity tests, sanity tests * tests are now in tests * added ground truth files * removed leftovers * npy files * tests tidyup * tests adapted for current structure but still useless * renamed test to comply with pytest convention * fix ground_truth generation script * refactor unittests * added accuracy tests * Added saturation test for cloud base in unit tests * Refactored imports and clean up unused code in simulation and parcel modules * Enhanced accuracy tests by adding assertions for output completeness and consistency in unit tests * Refactored calculation of simulated mass fraction evaporation point to initialize variable and improved error handling --------- Co-authored-by: Hevagog Co-authored-by: Mateusz Mazur * Notebook not running fix - added back necessary imports m8 u removed too much * black formatter * black for ipynb * Update PySDM/physics/particle_shape_and_density/oblate_spheroid.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> * Updated tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py Using dir shadowed the built-in dir() function. Renamed this variable (e.g., root_dir) to avoid confusion. Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> * Fixed pressure values in Planet subclasses to use correct scientific notation * Updated method signature for proper binding. Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> * Refactored Settings class to remove unused 'coord' parameter and updated diffusion_coordinate assignment * Removed unused Lofus et al. 2021 implementations from particle shape and density, terminal velocity, and ventilation modules * plots warnings fixes * removed files rogers and yau * notebook rerun * notebook fix * minor linter fixes --------- Co-authored-by: emmacware Co-authored-by: Piotr Kubala Co-authored-by: Hevagog Co-authored-by: Mateusz Mazur --- PySDM/physics/terminal_velocity/__init__.py | 4 +- .../Loftus_and_Wordsworth_2021/__init__.py | 8 + .../Loftus_and_Wordsworth_2021/figure_2.ipynb | 14723 ++++++++++++++++ .../Loftus_and_Wordsworth_2021/parcel.py | 62 + .../Loftus_and_Wordsworth_2021/planet.py | 118 + .../Loftus_and_Wordsworth_2021/settings.py | 46 + .../Loftus_and_Wordsworth_2021/simulation.py | 82 + .../Loftus_and_Wordsworth_2021/__init__.py | 0 .../accuracy_test.py | 231 + .../ground_truth/RHgrid.npy | Bin 0 -> 72128 bytes .../ground_truth/RHs.npy | Bin 0 -> 608 bytes .../ground_truth/gen_figure.py | 97 + .../ground_truth/m_frac_evap.npy | Bin 0 -> 72128 bytes .../ground_truth/r0grid.npy | Bin 0 -> 72128 bytes .../ground_truth/r_mins.npy | Bin 0 -> 608 bytes .../Loftus_and_Wordsworth_2021/unit_test.py | 238 + tests/examples_tests/conftest.py | 1 + 17 files changed, 15609 insertions(+), 1 deletion(-) create mode 100644 examples/PySDM_examples/Loftus_and_Wordsworth_2021/__init__.py create mode 100644 examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb create mode 100644 examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py create mode 100644 examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py create mode 100644 examples/PySDM_examples/Loftus_and_Wordsworth_2021/settings.py create mode 100644 examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py create mode 100644 tests/examples_tests/Loftus_and_Wordsworth_2021/__init__.py create mode 100644 tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py create mode 100644 tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/RHgrid.npy create mode 100644 tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/RHs.npy create mode 100644 tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py create mode 100644 tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/m_frac_evap.npy create mode 100644 tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/r0grid.npy create mode 100644 tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/r_mins.npy create mode 100644 tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py diff --git a/PySDM/physics/terminal_velocity/__init__.py b/PySDM/physics/terminal_velocity/__init__.py index 2e012dc474..3434f7ad3a 100644 --- a/PySDM/physics/terminal_velocity/__init__.py +++ b/PySDM/physics/terminal_velocity/__init__.py @@ -1,4 +1,6 @@ -"""terminal velocity formulae""" +""" +terminal velocity formulae +""" from .rogers_yau import RogersYau from .gunn_kinzer_1949 import GunnKinzer1949 diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/__init__.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/__init__.py new file mode 100644 index 0000000000..2d0f59fb30 --- /dev/null +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/__init__.py @@ -0,0 +1,8 @@ +""" +The Physics of Falling Raindrops in Diverse Planetary Atmospheres +[Loftus, K., & Wordsworth, R. D. 2021 (Journal of Geophysical Research: Planets, 126)](https://doi.org/10.1029/2020JE006653) + +""" + +from .settings import Settings +from .simulation import Simulation diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb new file mode 100644 index 0000000000..6a27e3ffa6 --- /dev/null +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb @@ -0,0 +1,14723 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "1e8d983b", + "metadata": {}, + "outputs": [], + "source": [ + "!pip install joblib" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "d8f644e9", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from numba import njit\n", + "from joblib import Parallel, delayed\n", + "from open_atmos_jupyter_utils import show_plot\n", + "\n", + "from PySDM import Formulae\n", + "from PySDM.physics import si\n", + "from scipy.optimize import fsolve\n", + "from PySDM_examples.Loftus_and_Wordsworth_2021 import Settings\n", + "from PySDM_examples.Loftus_and_Wordsworth_2021.planet import (\n", + " Planet,\n", + " EarthLike,\n", + " Earth,\n", + " EarlyMars,\n", + " Jupiter,\n", + " Saturn,\n", + " K2_18B,\n", + ")\n", + "from PySDM_examples.Loftus_and_Wordsworth_2021 import Simulation" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "92a3a574", + "metadata": {}, + "outputs": [], + "source": [ + "formulae = Formulae(\n", + " ventilation=\"PruppacherAndRasmussen1979\",\n", + " saturation_vapour_pressure=\"AugustRocheMagnus\",\n", + " diffusion_coordinate=\"WaterMassLogarithm\",\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "41f0ed6d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "300.0" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "new_Earth = EarthLike()\n", + "new_Earth.T_STP" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "d3a8f4b9", + "metadata": {}, + "outputs": [], + "source": [ + "radius_array = np.logspace(-4.5, -2.5, 50) * si.m\n", + "RH_array = np.linspace(0.25, 0.99, 50)\n", + "output_matrix = np.full((len(RH_array), len(radius_array)), np.nan)\n", + "const = formulae.constants\n", + "\n", + "\n", + "@njit()\n", + "def mix(dry, vap, ratio):\n", + " return (dry + ratio * vap) / (1 + ratio)\n", + "\n", + "\n", + "def compute_one_RH(i, RH):\n", + " \"\"\"\n", + " Compute one row of the output_matrix for a given RH.\n", + " Returns a 1D numpy array of length len(radius_array).\n", + " \"\"\"\n", + " new_Earth.RH_zref = RH\n", + "\n", + " pvs = formulae.saturation_vapour_pressure.pvs_water(new_Earth.T_STP)\n", + " initial_water_vapour_mixing_ratio = const.eps / (\n", + " new_Earth.p_STP / new_Earth.RH_zref / pvs - 1\n", + " )\n", + "\n", + " Rair = mix(const.Rd, const.Rv, initial_water_vapour_mixing_ratio)\n", + " c_p = mix(const.c_pd, const.c_pv, initial_water_vapour_mixing_ratio)\n", + "\n", + " def f(x):\n", + " return initial_water_vapour_mixing_ratio / (\n", + " initial_water_vapour_mixing_ratio + const.eps\n", + " ) * new_Earth.p_STP * (x / new_Earth.T_STP) ** (\n", + " c_p / Rair\n", + " ) - formulae.saturation_vapour_pressure.pvs_water(\n", + " x\n", + " )\n", + "\n", + " tdews = fsolve(f, [150, 300])\n", + " Tcloud = np.max(tdews)\n", + " Zcloud = (new_Earth.T_STP - Tcloud) * c_p / new_Earth.g_std\n", + " thstd = formulae.trivia.th_std(new_Earth.p_STP, new_Earth.T_STP)\n", + "\n", + " pcloud = formulae.hydrostatics.p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio(\n", + " new_Earth.p_STP, thstd, initial_water_vapour_mixing_ratio, Zcloud\n", + " )\n", + "\n", + " np.testing.assert_approx_equal(\n", + " actual=pcloud\n", + " * (\n", + " initial_water_vapour_mixing_ratio\n", + " / (initial_water_vapour_mixing_ratio + const.eps)\n", + " )\n", + " / formulae.saturation_vapour_pressure.pvs_water(Tcloud),\n", + " desired=1,\n", + " significant=4,\n", + " )\n", + "\n", + " output = None\n", + " row_data = np.full(len(radius_array), np.nan)\n", + " for j, r in enumerate(radius_array[::-1]):\n", + " settings = Settings(\n", + " planet=new_Earth,\n", + " r_wet=r,\n", + " mass_of_dry_air=1e5 * si.kg,\n", + " initial_water_vapour_mixing_ratio=initial_water_vapour_mixing_ratio,\n", + " pcloud=pcloud,\n", + " Zcloud=Zcloud,\n", + " Tcloud=Tcloud,\n", + " formulae=formulae,\n", + " )\n", + " simulation = Simulation(settings)\n", + " try:\n", + " output = simulation.run()\n", + " if output[\"z\"][-1] > 0:\n", + " row_data[j] = np.nan\n", + " break\n", + " else:\n", + " row_data[j] = 1 - (output[\"r\"][-1] / (r * 1e6))\n", + " except Exception as _:\n", + " break\n", + "\n", + " return i, row_data, output" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "4352de81", + "metadata": {}, + "outputs": [], + "source": [ + "all_rows = Parallel(n_jobs=os.cpu_count())(\n", + " delayed(compute_one_RH)(i, RH) for i, RH in enumerate(RH_array[::-1])\n", + ")\n", + "\n", + "last_output = None\n", + "for i, row_data, output in all_rows:\n", + " output_matrix[i] = row_data\n", + " last_output = output" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "8e6027d8", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2025-06-10T15:10:38.519047\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.10.3, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "53496ecd340048b79ff90faa9da71b6e", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(HTML(value=\".\\\\tmp_nham5pj.pdf
\"), HTML(val…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "row_data = output_matrix[18, ::-1] # Reverse the row for plotting\n", + "fig, ax = plt.subplots(\n", + " 2,\n", + " 1,\n", + " figsize=(6, 6),\n", + " sharex=True,\n", + " gridspec_kw={\"height_ratios\": [1, 3]},\n", + " constrained_layout=True,\n", + ")\n", + "ax[0].plot(radius_array, row_data, label=f\"Surface RH = {RH_array[-18]:.2f} %\")\n", + "ax[0].set_ylabel(\"Mass evaporated (%)\")\n", + "ax[0].legend()\n", + "\n", + "h = ax[1].contourf(\n", + " radius_array,\n", + " RH_array[::-1],\n", + " output_matrix[:, ::-1],\n", + " cmap=\"gray_r\",\n", + " levels=np.linspace(0, 1, 100),\n", + ")\n", + "ax[1].set_xscale(\"log\")\n", + "\n", + "# Add labels\n", + "ax[1].set_xlabel(\"Radius (µm)\")\n", + "ax[1].set_ylabel(\"Surface RH (%)\")\n", + "\n", + "cbar = fig.colorbar(h, ax=ax, shrink=0.4)\n", + "cbar.set_label(\"fraction Mass evaporated\")\n", + "contour_levels = [0.1] # Define the level for the contour\n", + "ax[1].contour(\n", + " radius_array,\n", + " RH_array[::-1],\n", + " output_matrix[:, ::-1],\n", + " levels=contour_levels,\n", + " colors=\"red\",\n", + " linewidths=1.5,\n", + ")\n", + "show_plot()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "369fa24a", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2025-06-10T15:10:39.361235\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.10.3, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "0baa5c73262744e4890c68dcfaa37b13", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(HTML(value=\".\\\\tmps5jrn8sn.pdf
\"), HTML(val…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, axs = plt.subplots(1, 3, figsize=(10, 8), sharey=True)\n", + "axs[0].plot(last_output[\"r\"], last_output[\"z\"], label=\"Radius\")\n", + "axs[0].set_ylabel(\"Height (m)\")\n", + "axs[0].set_xlabel(\"Radius (um)\")\n", + "axs[0].legend()\n", + "axs[1].plot(last_output[\"S\"], last_output[\"z\"], label=\"Supersaturation\")\n", + "axs[1].set_xlabel(\"Supersaturation\")\n", + "axs[1].set_ylabel(\"Height (m)\")\n", + "axs[1].legend()\n", + "axs[2].plot(last_output[\"t\"], last_output[\"z\"], label=\"Height\")\n", + "axs[2].set_ylabel(\"Height (m)\")\n", + "axs[2].set_xlabel(\"Time (s)\")\n", + "axs[2].legend()\n", + "show_plot()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "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.13.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py new file mode 100644 index 0000000000..7ac5eb0dac --- /dev/null +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py @@ -0,0 +1,62 @@ +from PySDM.environments.parcel import Parcel + + +class AlienParcel(Parcel): + def __init__( + self, + dt, + mass_of_dry_air: float, + pcloud: float, + initial_water_vapour_mixing_ratio: float, + Tcloud: float, + w: float = 0, + zcloud: float = 0, + mixed_phase=False, + ): + super().__init__( + dt=dt, + mass_of_dry_air=mass_of_dry_air, + p0=pcloud, + initial_water_vapour_mixing_ratio=initial_water_vapour_mixing_ratio, + T0=Tcloud, + w=w, + z0=zcloud, + mixed_phase=mixed_phase, + variables=None, + ) + + def advance_parcel_vars(self): + """ + Compute new values of displacement, dry-air density and volume, + and write them to self._tmp and self.mesh.dv + """ + dt = self.particulator.dt + formulae = self.particulator.formulae + T = self["T"][0] + p = self["p"][0] + + dz_dt = -self.particulator.attributes["terminal velocity"].to_ndarray()[0] + water_vapour_mixing_ratio = ( + self["water_vapour_mixing_ratio"][0] + - self.delta_liquid_water_mixing_ratio / 2 + ) + + drho_dz = formulae.hydrostatics.drho_dz( + p=p, + T=T, + water_vapour_mixing_ratio=water_vapour_mixing_ratio, + lv=formulae.latent_heat_vapourisation.lv(T), + d_liquid_water_mixing_ratio__dz=( + self.delta_liquid_water_mixing_ratio / dz_dt / dt + ), + ) + drhod_dz = drho_dz + + self.particulator.backend.explicit_euler(self._tmp["z"], dt, dz_dt) + self.particulator.backend.explicit_euler( + self._tmp["rhod"], dt, dz_dt * drhod_dz + ) + + self.mesh.dv = formulae.trivia.volume_of_density_mass( + (self._tmp["rhod"][0] + self["rhod"][0]) / 2, self.mass_of_dry_air + ) diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py new file mode 100644 index 0000000000..0e4211b397 --- /dev/null +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py @@ -0,0 +1,118 @@ +from PySDM.physics.constants import si +from dataclasses import dataclass +from typing import Optional, Dict, Any + + +@dataclass +class Planet: + g_std: float + T_STP: float + p_STP: float + RH_zref: float + dry_molar_conc_H2: float + dry_molar_conc_He: float + dry_molar_conc_N2: float + dry_molar_conc_O2: float + dry_molar_conc_CO2: float + H_LCL: Optional[float] = None + + def to_dict(self) -> Dict[str, Any]: + return vars(self) + + +@dataclass +class EarthLike(Planet): + g_std: float = 9.82 * si.metre / si.second**2 + T_STP: float = 300 * si.kelvin + p_STP: float = 1.01325 * 1e6 * si.pascal + RH_zref: float = 0.75 + dry_molar_conc_H2: float = 0 + dry_molar_conc_He: float = 0 + dry_molar_conc_N2: float = 1 + dry_molar_conc_O2: float = 0 + dry_molar_conc_CO2: float = 0 + H_LCL: float = 8.97 * si.kilometre + + +@dataclass +class Earth(Planet): + g_std: float = 9.82 * si.metre / si.second**2 + T_STP: float = 290 * si.kelvin + p_STP: float = 1.01325 * 1e6 * si.pascal + RH_zref: float = 0.75 + dry_molar_conc_H2: float = 0 + dry_molar_conc_He: float = 0 + dry_molar_conc_N2: float = 0.8 + dry_molar_conc_O2: float = 0.2 + dry_molar_conc_CO2: float = 0 + H_LCL: float = 8.41 * si.kilometre + + +@dataclass +class EarlyMars(Planet): + g_std: float = 3.71 * si.metre / si.second**2 + T_STP: float = 290 * si.kelvin + p_STP: float = 2 * 1e6 * si.pascal + RH_zref: float = 0.75 + dry_molar_conc_H2: float = 0 + dry_molar_conc_He: float = 0 + dry_molar_conc_N2: float = 0 + dry_molar_conc_O2: float = 0 + dry_molar_conc_CO2: float = 1 + H_LCL: float = 14.5 * si.kilometre + + +@dataclass +class Jupiter(Planet): + g_std: float = 24.84 * si.metre / si.second**2 + T_STP: float = 274 * si.kelvin + p_STP: float = 4.85 * 1e6 * si.pascal + RH_zref: float = 1 + dry_molar_conc_H2: float = 0.864 + dry_molar_conc_He: float = 0.136 + dry_molar_conc_N2: float = 0 + dry_molar_conc_O2: float = 0 + dry_molar_conc_CO2: float = 0 + H_LCL: float = 39.8 * si.kilometre + + +@dataclass +class Saturn(Planet): + g_std: float = 10.47 * si.metre / si.second**2 + T_STP: float = 284 * si.kelvin + p_STP: float = 10.4 * 1e6 * si.pascal + RH_zref: float = 1 + dry_molar_conc_H2: float = 0.88 + dry_molar_conc_He: float = 0.12 + dry_molar_conc_N2: float = 0 + dry_molar_conc_O2: float = 0 + dry_molar_conc_CO2: float = 0 + H_LCL: float = 99.2 * si.kilometre + + +@dataclass +class K2_18B(Planet): + g_std: float = 12.44 * si.metre / si.second**2 + T_STP: float = 275 * si.kelvin + p_STP: float = 0.1 * 1e6 * si.pascal + RH_zref: float = 1 + dry_molar_conc_H2: float = 0.9 + dry_molar_conc_He: float = 0.1 + dry_molar_conc_N2: float = 0 + dry_molar_conc_O2: float = 0 + dry_molar_conc_CO2: float = 0 + H_LCL: float = 56.6 * si.kilometre + + +@dataclass +class CompositeTest(Planet): + g_std: float = 9.82 * si.metre / si.second**2 + T_STP: float = 275 * si.kelvin + p_STP: float = 0.75 * 1e6 * si.pascal + RH_zref: float = 1 + dry_molar_conc_H2: float = 0.1 + dry_molar_conc_He: float = 0.1 + dry_molar_conc_N2: float = 0.1 + dry_molar_conc_O2: float = 0.1 + dry_molar_conc_CO2: float = 0.1 + H_LCL: Optional[float] = None diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/settings.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/settings.py new file mode 100644 index 0000000000..b0a81752a7 --- /dev/null +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/settings.py @@ -0,0 +1,46 @@ +# Planetary Properties, Loftus and Wordsworth 2021 Table 1 + +from pystrict import strict + +from PySDM import Formulae +from PySDM.dynamics import condensation +from PySDM.physics.constants import si +from PySDM_examples.Loftus_and_Wordsworth_2021.planet import Planet + + +@strict +class Settings: + def __init__( + self, + r_wet: float, + mass_of_dry_air: float, + planet: Planet, + initial_water_vapour_mixing_ratio: float, + pcloud: float, + Zcloud: float, + Tcloud: float, + formulae: Formulae = None, + ): + self.formulae = formulae or Formulae( + saturation_vapour_pressure="AugustRocheMagnus", + diffusion_coordinate="WaterMassLogarithm", + ) + + self.initial_water_vapour_mixing_ratio = initial_water_vapour_mixing_ratio + self.p0 = planet.p_STP + self.RH0 = planet.RH_zref + self.kappa = 0.2 + self.T0 = planet.T_STP + self.z_half = 150 * si.metres + self.dt = 1 * si.second + self.pcloud = pcloud + self.Zcloud = Zcloud + self.Tcloud = Tcloud + + self.r_wet = r_wet + self.mass_of_dry_air = mass_of_dry_air + self.n_output = 500 + + self.rtol_x = 0.5 * (condensation.DEFAULTS.rtol_x) + self.rtol_thd = condensation.DEFAULTS.rtol_thd + self.dt_cond_range = condensation.DEFAULTS.cond_range diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py new file mode 100644 index 0000000000..071dde8bf2 --- /dev/null +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py @@ -0,0 +1,82 @@ +import numpy as np + +import PySDM.products as PySDM_products +from PySDM.backends import CPU +from PySDM.builder import Builder +from PySDM.dynamics import AmbientThermodynamics, Condensation +from PySDM.physics import constants as const +from PySDM_examples.Loftus_and_Wordsworth_2021.parcel import AlienParcel + + +class Simulation: + def __init__(self, settings, backend=CPU): + builder = Builder( + backend=backend( + formulae=settings.formulae, + **( + {"override_jit_flags": {"parallel": False}} + if backend == CPU + else {} + ), + ), + n_sd=1, + environment=AlienParcel( + dt=settings.dt, + mass_of_dry_air=settings.mass_of_dry_air, + pcloud=settings.pcloud, + zcloud=settings.Zcloud, + initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio, + Tcloud=settings.Tcloud, + ), + ) + + builder.add_dynamic(AmbientThermodynamics()) + builder.add_dynamic( + Condensation( + rtol_x=settings.rtol_x, + rtol_thd=settings.rtol_thd, + dt_cond_range=settings.dt_cond_range, + ) + ) + builder.request_attribute("terminal velocity") + + attributes = {} + r_dry = 1e-10 + attributes["dry volume"] = settings.formulae.trivia.volume(radius=r_dry) + attributes["kappa times dry volume"] = attributes["dry volume"] * settings.kappa + attributes["multiplicity"] = np.array([1], dtype=np.int64) + r_wet = settings.r_wet + attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet) + products = [ + PySDM_products.MeanRadius(name="radius_m1", unit="um"), + PySDM_products.ParcelDisplacement(name="z"), + PySDM_products.AmbientRelativeHumidity(name="RH", unit="%"), + PySDM_products.Time(name="t"), + ] + + self.particulator = builder.build(attributes, products) + + def save(self, output): + cell_id = 0 + output["r"].append( + self.particulator.products["radius_m1"].get(unit=const.si.m)[cell_id] + ) + + output["z"].append(self.particulator.products["z"].get()[cell_id]) + output["S"].append(self.particulator.products["RH"].get()[cell_id] / 100 - 1) + output["t"].append(self.particulator.products["t"].get()) + + def run(self): + output = { + "r": [], + "S": [], + "z": [], + "t": [], + } + + self.save(output) + while self.particulator.environment["z"][0] > 0 and output["r"][-1] > 1e-6: + self.particulator.run(1) + self.save(output) + + return output diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/__init__.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py new file mode 100644 index 0000000000..75f57bb05d --- /dev/null +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py @@ -0,0 +1,231 @@ +import pytest +import os +import numpy as np +from scipy.optimize import fsolve + +from PySDM import Formulae +from PySDM.physics import si +from PySDM_examples.Loftus_and_Wordsworth_2021 import Settings, Simulation +from PySDM_examples.Loftus_and_Wordsworth_2021.planet import EarthLike + + +class GroundTruthLoader: + def __init__(self, groundtruth_dir_path, n_samples=2, random_seed=2137): + self.dir_path = groundtruth_dir_path + self.RHs = None + self.r0grid = None + self.m_frac_evap = None + self.n_samples = n_samples # Number of random samples to test + np.random.seed(random_seed) # reproducible random samples during debugging + + def __enter__(self): + try: + self.RHs = np.load(os.path.join(self.dir_path, "RHs.npy")) + self.r0grid = np.load(os.path.join(self.dir_path, "r0grid.npy")) + self.m_frac_evap = np.load(os.path.join(self.dir_path, "m_frac_evap.npy")) + return self + except FileNotFoundError as e: + print(f"Error loading ground truth files: {e}") + raise + except Exception as e: + print(f"An unexpected error occurred while loading ground truth data: {e}") + raise + + def __exit__(self, exc_type, exc_val, exc_tb): + pass + + +class TestNPYComparison: + @staticmethod + def _mix(dry_prop, vap_prop, ratio): + return (dry_prop + ratio * vap_prop) / (1 + ratio) + + def _calculate_cloud_properties( + self, planet: EarthLike, surface_RH: float, formulae_instance: Formulae + ): + const = formulae_instance.constants + + planet.RH_zref = surface_RH + + pvs_stp = formulae_instance.saturation_vapour_pressure.pvs_water(planet.T_STP) + + initial_water_vapour_mixing_ratio = const.eps / ( + planet.p_STP / planet.RH_zref / pvs_stp - 1 + ) + + R_air_mix = self._mix(const.Rd, const.Rv, initial_water_vapour_mixing_ratio) + cp_mix = self._mix(const.c_pd, const.c_pv, initial_water_vapour_mixing_ratio) + + def solve_Tcloud(T_candidate): + pv_ad = ( + initial_water_vapour_mixing_ratio + / (initial_water_vapour_mixing_ratio + const.eps) + * planet.p_STP + * (T_candidate / planet.T_STP) ** (cp_mix / R_air_mix) + ) + pvs_tc = formulae_instance.saturation_vapour_pressure.pvs_water(T_candidate) + return pv_ad - pvs_tc + + Tcloud_solutions = fsolve(solve_Tcloud, [150.0, 300.0]) + Tcloud = np.max(Tcloud_solutions) + + Zcloud = (planet.T_STP - Tcloud) * cp_mix / planet.g_std + + th_std = formulae_instance.trivia.th_std(planet.p_STP, planet.T_STP) + + pcloud = formulae_instance.hydrostatics.p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio( + planet.p_STP, th_std, initial_water_vapour_mixing_ratio, Zcloud + ) + return initial_water_vapour_mixing_ratio, Tcloud, Zcloud, pcloud + + def test_figure_2_replication_accuracy(self): + current_dir = os.path.dirname(os.path.abspath(__file__)) + groundtruth_dir = os.path.abspath(os.path.join(current_dir, "ground_truth")) + if not os.path.isdir(groundtruth_dir): + pytest.fail(f"Groundtruth directory not found at {groundtruth_dir}") + + formulae = Formulae( + ventilation="PruppacherAndRasmussen1979", + saturation_vapour_pressure="AugustRocheMagnus", + diffusion_coordinate="WaterMassLogarithm", + ) + + with GroundTruthLoader(groundtruth_dir) as gt: + if gt.RHs is None or gt.r0grid is None or gt.m_frac_evap is None: + pytest.fail("Ground truth data (.npy files) not loaded properly.") + + n_rh_values = len(gt.RHs) + n_radius_values = gt.r0grid.shape[1] + + total_points = n_rh_values * n_radius_values + n_samples = min(gt.n_samples, total_points) + + if n_samples == 0: + pytest.skip("No data points available to sample.") + + all_indices = np.array( + [(i, j) for i in range(n_rh_values) for j in range(n_radius_values)] + ) + + sampled_indices_flat = np.random.choice(len(all_indices), n_samples, replace=False) + sampled_ij_pairs = all_indices[sampled_indices_flat] + + for i_rh, j_r in sampled_ij_pairs: + current_planet_state = EarthLike() + + current_rh = gt.RHs[i_rh] + current_r_m = gt.r0grid[0, j_r] + expected_m_frac_evap = gt.m_frac_evap[i_rh, j_r] + + try: + iwvmr, Tcloud, Zcloud, pcloud = self._calculate_cloud_properties( + current_planet_state, current_rh, formulae + ) + simulated_m_frac_evap_point = self.calc_simulated_m_frac_evap_point( + current_planet_state, + formulae, + i_rh, + j_r, + current_rh, + current_r_m, + expected_m_frac_evap, + iwvmr, + Tcloud, + Zcloud, + pcloud, + ) + except Exception as e: + print( + f"Warning: Error in _calculate_cloud_properties for RH={current_rh} (sample idx {i_rh},{j_r}): {e}." + ) + + error_context = ( + f"Sample (RH_idx={i_rh}, R_idx={j_r}), " + f"RH={current_rh:.4f}, R_m={current_r_m:.3e}. " + f"Expected: {expected_m_frac_evap}, Got: {simulated_m_frac_evap_point}" + ) + + if np.isnan(expected_m_frac_evap): + assert np.isnan( + simulated_m_frac_evap_point + ), f"NaN Mismatch. {error_context} (Expected NaN, got non-NaN)" + else: + assert not np.isnan( + simulated_m_frac_evap_point + ), f"NaN Mismatch. {error_context} (Expected non-NaN, got NaN)" + np.testing.assert_allclose( + simulated_m_frac_evap_point, + expected_m_frac_evap, + rtol=1e-1, # Relative tolerance + atol=1e-1, # Absolute tolerance + err_msg=f"Value Mismatch. {error_context}", + ) + + def calc_simulated_m_frac_evap_point( + self, + current_planet_state, + formulae, + i_rh, + j_r, + current_rh, + current_r_m, + expected_m_frac_evap, + iwvmr, + Tcloud, + Zcloud, + pcloud, + ): + + simulated_m_frac_evap_point = np.nan + + if np.isnan(current_r_m) or current_r_m <= 0: + print(f"Warning: Invalid radius current_r_m={current_r_m} for sample idx {i_rh},{j_r}.") + else: + settings = Settings( + planet=current_planet_state, + r_wet=current_r_m, + mass_of_dry_air=1e5 * si.kg, + initial_water_vapour_mixing_ratio=iwvmr, + pcloud=pcloud, + Zcloud=Zcloud, + Tcloud=Tcloud, + formulae=formulae, + ) + simulation = Simulation(settings) + + try: + output = simulation.run() + if ( + output + and "r" in output + and len(output["r"]) > 0 + and "z" in output + and len(output["z"]) > 0 + ): + final_radius_um = output["r"][-1] + if np.isnan(final_radius_um) or final_radius_um < 0: + final_radius_m = final_radius_um * 1e-6 + if final_radius_m < 0: # Non-physical radius + simulated_m_frac_evap_point = 1.0 # 1.0 means fully evaporated + else: + simulated_m_frac_evap_point = np.nan + else: + final_radius_m = final_radius_um * 1e-6 + if current_r_m == 0: + frac_evap = 1.0 + else: + frac_evap = 1.0 - (final_radius_m / current_r_m) ** 3 + frac_evap = np.clip(frac_evap, 0.0, 1.0) + simulated_m_frac_evap_point = frac_evap + else: + simulated_m_frac_evap_point = np.nan + except Exception as e: + print( + f"Warning: Simulation run failed for RH={current_rh}, r={current_r_m} (sample idx {i_rh},{j_r}): {e}." + ) + if np.isclose(expected_m_frac_evap, 1.0, atol=1e-6): + simulated_m_frac_evap_point = 1.0 + else: + simulated_m_frac_evap_point = np.nan + + return simulated_m_frac_evap_point diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/RHgrid.npy b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/RHgrid.npy new file mode 100644 index 0000000000000000000000000000000000000000..e02761b427cb379b6be665a7fbf61a190db2633e GIT binary patch literal 72128 zcmeI5OG}ht7={(3B8m>U2#uOh8k1S1q#e+n2vLzrE2LaxOgWOGS))j3QP6=xjIa_4 zGIUWAAwq~sMPLw-7?qle<{=$Vqa)#>La4rXW`04dXK~HuCYKl2eK*d$@8{h^*$1<7 zLPN?!D(!i$($jAH9-BQWf0x~1v*(w%%iPYhCrjLUF7xwLXHluk`g!RY=Q)@4@0i5x z4qIHp_H7Q^CENeKR)+XrpM*37;D85spt%R_PuJB}Mtl|$4tPKhwC_M^W|ZUSjxR#O z0T1Yb_8qXss>7YxEhHT9fF5Yyfyuhmn`r^g#O#JT^7G)+;0&@PHm@-vQH7 zPu*7`;eZG9K>H4?x4s>8w)Y7M2Rxt$+IOHcxw?6%yI)8+-~m0*z5~`TgVrZl146<9 z59opR9oTE$R;(Bl5)ODk547)q|M*qZkdW3fG&taaIfjP^YCSOTDms(D2?+;0pa#&gKgFHCk0Ul`Xfw?u;;T{nZ4tPKhwC{lbxM=FAkmijh9Pj`SH1{BI-eNR( z{)&W$`N9|v@Ia>rjpngW&M_h3fCuzI`wj%ohm8f#i<0m#=VtK$4|ICq^gmGWyO8EV zt8l;rJkZ>Oz=ogVlfC)4Q>(p>z810LXk z<{kv@w@wD{JCpD*_kHmI4|ICq|MNJWNg>V6IXK_}9%$}C;6CrH?O3m_J+?02?so&2ikWa6;M}4tPKhwC_Mg)s{!C8)k)s10K)=?K_b8u|K{!XI4l!-~m0* zz5|PE@4Rn(I4dL^@PHm@-+|_%X}R@-vqHiF59opR9Vn{(v-n`6S4cSE0X@*Z12F?F zw`xv!g@gkh&;#u|;E62WTvg{45)ODk547(?+UA6W z10K)=?K==TQ8QZRoD&ibct8)d??B(yjPl}!IU(VI2lPPu4%{mUU6nsRCnOy3fF5Yy zfvo3k^|{gWLc#$L=z;bfSiMlP^O$R1NI2jDJV{IWgP>k|?Vct8)d??6I_J2T3$AS4{{fF5YyfgfS9)*~bqgoFbg&;#u| U(3Cb-5#F>QBpmR79%$czf9(F>CIA2c literal 0 HcmV?d00001 diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/RHs.npy b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/RHs.npy new file mode 100644 index 0000000000000000000000000000000000000000..73132c62032c60981b5d5716004d69315df29ade GIT binary patch literal 608 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I$d20EHL3bhL411<(IxM06?x!|;Fsf+f-{>nPvO)l93i9<#Hm+ha;bBlPF zdBuJ{h-$fN4?S?0#)3!w*zrGZrg*kE8npP@lM;` zwa)@8xMvSC{ZP?;d$2KX51{%V*slOtpZE}J&O`e|u!9~!&3kMQv<&FN$53;h*k1vf z+PmQi)cmLRU~5l5h1&DX9_-ww&!F}_w}<%ig?-Qtxz^W#FYNh$UzK_?=>^pO7xow#|jdroFV^>ygF2|N2Y&s&2)(+a+Gvo1MF1v?219{bzv* zCsxdQWxpZVJ#OLcSN1v6{(YV)^V(kH=GM+hv9Il)O6SUV&wFjZC{AT-^ZnQMKt*0P z3UBPCpGn% z*1oOeUXsW2xAsokr~YuzcxV6rr%#_{<~#d6HY_SetKQjH)bHM@_4=K?DKNH_b>7>5 TWYGdz{oa0^`@<@cb?@x~9yk!q literal 0 HcmV?d00001 diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py new file mode 100644 index 0000000000..ff8d780ac4 --- /dev/null +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py @@ -0,0 +1,97 @@ +################################################################ +# make LoWo21 Figure 2 +# r_min, fraction raindrop mass evaporated as functions of RH +################################################################ +import numpy as np +import matplotlib.pyplot as plt +import os + +# load results +root_dir = os.path.dirname(os.path.abspath(__file__)) +RHs = np.load(os.path.join(root_dir, "RHs.npy")) +r0grid = np.load(os.path.join(root_dir, "r0grid.npy")) +RHgrid = np.load(os.path.join(root_dir, "RHgrid.npy")) +m_frac_evap = np.load(os.path.join(root_dir, "m_frac_evap.npy")) +r_mins = np.load(os.path.join(root_dir, "r_mins.npy")) + +i_RH75 = 29 # index for RH=0.75 + +# make figure 2 +# only put colorbar on lower panel, still line up x-axes +f, axs = plt.subplots( + 2, + 2, + sharex="col", + figsize=(6, 7), + gridspec_kw={"height_ratios": [1, 3], "width_ratios": [20, 1]}, +) +plt.subplots_adjust(hspace=0.05) +axs[0, 0].set_xscale("log") +axs[1, 0].set_xlabel(r"$r_0$ [mm]") +axs[1, 0].set_ylabel("surface RH [ ]") +axs[0, 0].tick_params(right=True, which="both") +axs[1, 0].tick_params(right=True, which="both") +axs[1, 0].tick_params(top=True, which="both") +axs[0, 0].tick_params(top=True, which="both") +axs[0, 0].set_ylabel("fraction mass \n evaporated [ ]") +axs[0, 0].set_xlim(r0grid[0, 0] * 1e3, r0grid[0, -1] * 1e3) +axs[0, 0].set_ylim(-0.04, 1.04) +levels_smooth = np.linspace(0, 1, 250) +cmesh = axs[1, 0].contourf( + r0grid * 1e3, + RHgrid, + m_frac_evap, + cmap=plt.cm.binary, + vmin=0, + vmax=1, + levels=levels_smooth, +) + + +cb = f.colorbar(cmesh, cax=axs[1, 1]) +cb.solids.set_edgecolor("face") +axs[0, 1].axis("off") +cb.solids.set_edgecolor("face") +cb.solids.set_linewidth(1e-5) +cb.set_label("fraction mass evaporated [ ]") +cb.set_ticks([0, 0.1, 0.25, 0.5, 0.75, 1]) +axs[1, 0].axhline(0.75, lw=0.5, ls="--", c="plum") +axs[1, 0].plot(r_mins * 1e3, RHs, lw=3, c="darkviolet", zorder=10) +c_10 = axs[1, 0].contour( + r0grid * 1e3, + RHgrid, + m_frac_evap, + colors="indigo", + linewidths=1, + linestyles="--", + levels=[0.1], +) +cb.add_lines(c_10) +axs[1, 0].clabel(c_10, c_10.levels, fmt={0.1: "10% mass evaporated"}, fontsize="smaller") +axs[1, 0].fill_betweenx(RHs, 1e-2, (r_mins - 1e-6) * 1e3, edgecolor="k", facecolor="w", hatch="//") +axs[1, 0].annotate(text="TOTAL \\n EVAPORATION", xy=(0.04, 0.55), c="k", backgroundcolor="w") +axs[1, 0].annotate( + text=r"$r_\mathrm{min}$", + xy=(0.35, 0.27), + c="darkviolet", + backgroundcolor="w", + size=8, +) + +axs[0, 0].scatter(r_mins[i_RH75] * 1e3, 0.99, color="darkviolet", zorder=10) +axs[0, 0].axvline(r_mins[i_RH75] * 1e3, lw=0.5, c="darkviolet", ls="--") +axs[0, 0].plot(r0grid[0, :] * 1e3, m_frac_evap[i_RH75, :], lw=2.05, c="k") +axs[0, 0].axhline(1, c="w", lw=3) +axs[0, 0].plot([1e-2, r_mins[i_RH75] * 1e3], [1, 1], c="k", lw=2.05, ls="--") +axs[0, 0].annotate(text="surface RH = 0.75", xy=(0.8, 0.85), c="plum", size=8) +axs[0, 0].annotate(text=r"$r_\mathrm{min}$", xy=(0.13, 0.05), c="darkviolet", size=8) + +figs_path = os.path.join(dir, "figs") +os.mkdir(figs_path) if not os.path.exists(figs_path) else None +plt.savefig( + os.path.join(figs_path, "fig02.pdf"), + transparent=True, + bbox_inches="tight", + pad_inches=0.5, +) +plt.close() diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/m_frac_evap.npy b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/m_frac_evap.npy new file mode 100644 index 0000000000000000000000000000000000000000..657765da4b1600296772f69df82fe9e09157c51c GIT binary patch literal 72128 zcmeFaX*d;d)b?wPB1tL@A}N%jlF+cE5}|}jQidc-=CP7ZijZW^)HaWq=W&~7ws~fo z=b;P{)mhK`?R+`cx!!Z$bDsBk_~PQ@zUt>*_qx~p_quggLGb}Kl_k|{KAop#PfYnF zxcM&Xit!0@^XVFynwx4FY8aX7JpGT~OKTaJJ;mRf>1i21#h>{_g$21!pA!}kY<`5xTKPjpozv z;k){$!hKUfeS2q}LdgV7*RosPG98Bs+OEZ`f@8pZct*SV?C0SR0Yl!7p#t*s|4eV!h9Z0jw|(8<=`+q%6NCY&{@pKc|A>WcMOPJu!oS-vrhP|pXM zqMY*g@wq_JFCC?4&w&X_q}+S&EFjE1Rkk~k3C3HG_ZhLy zG9v{;=M_?5!jQbdc_kSntIU6YaZUt6exrd`d;(C>-J*4vjDv-@G=~Z$V?lUQ%(pBl z8WbIWb_$-0+T_&#_z$|tssFB=Vt+RE!~MlP?2zhJ6umSD$F7FT>CDeSsld5`9`9*z z_YIQBmW+8R*Ks(Zpcq;tH3nLmU&ANpM`22mEw{dH80J6f$YxF^w!8^09#2)xrqS~vv)&;A)rhh`7 zcS8SjwSdl_?cjFDxbU|}8xT(Bw*I=>3NF`1SV+hz!c7@HxsfLD=CzphY9s@{OH+@S zUn3A6&N2-q)&pxS3xn)PEeNl7oJc%Z1FJtSoYoGl0@As@=-Q)|&}p)?`n0$l>YlNX zNTy}L`!>SOPOKDy7@GLhca{KSSV|4|1_>yixO(o+7s2XVdbJk%LSU>d+?gw$4`S5= zLjoUj!B{Z)mDbN}aL#F!$T!P^)#Nqa6sAm|NSLmXv(sTBkcFYlCJhJz8U{BNQbG8m zmX`RL6!70~e&rQkG7$d$b7NCY1j#*fuO@60faLn#F|9BT$d22L8aZMi*+k`MnoBe! ze=vHz%!Zr-p{J6`5u3c)mmRe>?{_oJ`CaSZATij+wqB zZ5o!?P9)2!OaVJrS|07{1Ps-*=P-wl11GJg?;h1Lm~W{H?&lbV%r5B}yRl)Y&)>}w zP&@>`mEtItaf7I{Xm%bBr2t25w%NqZ0r*{zL3`J*4Yf>^B?Dc4 z)w`HajiCOp#`kGBa_Zn2XMv7dxR%E$E5KI+LECfBr+ujc@5f?S1$I|Lt4h#Knq=gZ zWaXn8%`#Bk`!MJ5u~J~Q<_^mHTMQiA0!K9`Nsw&DAlF`91bVUSK5mo(@Yr@gs-Ha{ zG%s;ekL%=uUb1OZe?vBKY-cV#cqa>}oW88EO=Z9Y&&Le~|8(%zmsl}XO9Kk?@{LPE zslfeh;I!TD6fj=bnK;Cp3=0lowkLQKft8B3yH7p=j9pj)+&todmD{%J^K1-ghBQ{I zX-5O`)a{41vyq@!?DHk`Qv`5-KQrodC47@to4opO#j7s)@$1iG=Hb*)7b}gob1;!Y zL(3yH3nHtp16iA(_4J;hSslpxR00(83tvUm44 z!7ZUJ<8u{c5E+nJ{chg~B=6S`VtnesQOUh+qOum;`0QiD_SAsetrcNzLKUz!mtHTR zsRWsAN^fU_%E4goTT?NGGDv^-B36Q}6c}YgYa^zLLDfd+k3<^@#8l{MqC$&+q_>Nb zm|Xy3jA732X7hlfvQSL?axO5&4R5E5%?9CrR=Y(7v%o_|_w1#L49JdX_X@R02j|AV zgY@EQK-w}U#I++8Jmeh*3RjXr^Sm6#=e{JMY@agInE4J1UEbVO90@@2j1P!-90!_} z9w)vP#sC$cDp6B38d&Roe=#qM1Y@_REe_9+Q{h{rqz{Dycc7z>&v@u2uQqx0--1^I zzes)(Wyq-krPF#ra}dgsaZXNt7Hp3UFyEV*hT}QCd+#_;fktccLqWbtDC=<1l&c>H zDml{J2iq|iyDXNQAwCMs#ICY+nh|)lCaFQyJp^ZN2e|Jj8Uzh`JArel6o_gM(ec$C zfWJjEPVA0-P{<&EhcmwiHbj^9HE!>Q%Yx58ov`kN4LPOGqrcnXZ;kh^SkE??I(VTl zNU{~qm>I$pO$(ffP8PE5Yy!o#Sq51W8LC)5zi~Hh1XY<>zZJ)NaE^z#;SAJMS48U^ zX=`96NcASCK^0W3NIcy#R{__Sv`-{Cl>@1s`*Yy=GMG8Q@lI!>1gx`4Ze_I=Lju`# zQl*Fl&gLo~cr1&+So!@*#HRuvG2c1HS(*n6vAW;CvF8Fu@(V{bLN+L>Xio*uWkJ=m zw|zC=GJt=D(rc)g4n&i4g}l4dfJonbH*_WiWIic8N^VOAV^?t_osc9DzRwtQF8Mnc zKh^H1^v8oF?`)noe;g1;NTEbmHetWThE;#o&A6M$Lb8iFk?2P@c3#|}hoXU4%qZx)ggMTNJo1lX6BDH508OWF9 zK6q#}!m71MNtAIt)HOM;y(ZQ|zRS)Z;f>UgFg8uDt@(pJX&N8kK`Q zZS{;JR~fK!4kw6>lz@oB)M=mmVqkyKY^xegg3e=cX3^3`kX>am+4>kcb(|*Yh)*6= zeH{!B9M6Hy;Xo@Rm28N)&Eq%Fl?h~#b$DQz0p#JL95w!Q(2HJNvYSflQ4vgAl`2@MT{f1=tiS zbaPDy;4xKj52tG%xVskFdldEnYuw;EFLO6YUwSv#Y}yH1To1f>{{#1ZM>pJgoZ2AS zLNNd1nO4xVFBlB^*$frOa@wv`G=X~l`{TL^WJrKqp5HebVX3b~ok_JGVy*~U?r^Gw z)xe~Of-TiRwS2}fSfUDIQq9@Y^D2OMXfm@yr5t8tw`9DhLrzVzDJxc&z_qL=F9t%3 z!8xn{=d2qExOXyf`mhzj#2uD_gcAk8N~4*__#_V|ru_aUW#*v2B*Iw3n+;+o{_PBi z&jb#&XZB9A84%<9ax2rHG~nktK0+!>g{l@QoqOL>z{70b!^0K*b76*aUs+L4bsfI{ zlm9yq9i_ZZKaL0D*gf`FiLt;soN#O}2Xe}`Fz@BpC?K{jy`(xF2_ET51_piMAiO4# z@+>$Ega?1$E;0%Q))!|}?$l52^bXqO)h4h0OY&+g;fS@!#2n;u z6+T_}nuX8&!lfoBXTY0=osPX^3JhL6y=(Vi65^zfebiqa2a*K;$(GPDIB4H76LEhO zq*A7vw;dP(?uZ=wgqa~QwY$*VQa1=xyH!QrBd5mRW-6XA8GxJg&&>6m`apV-<@9QH z53uw$d^5MG^>T3N^firr#>x(((L?!Wmt$B_Y@ji=U+HKc*q zhxQ~tzf@p+TB`8FECs~g`>if(A*W{e^EOC{ko;4@;b(6GXznuesox!sxrcSDB>hdzM4q!zdtMBX6l*!98E%!miM$a3J0(djCs53`mk^-}?dTDQ?3Btz$$W z==x{-NrV6;H(EVc8aY+SANbfPV3Svyy!tQ2tFe5C3r_2E(3nRd2t>|8d+--is*5vl zkuUM!M*S50Zg|u;s67eO9Rs&5{*8mBm}Xc()EMY7jIM;Mjl!KU7Z0U_BcN6_r1ugz zW#d)2;{xibe3^93{A`^2+7H`1>JPxS9A&O0+ddfCuvc67-UIxsR4?AHc7fx3ea&a} zPT+5{2l}pdSV-g0$}wnz=WP1hKJ00QzY+fLrP`W7{`cXHrsyWP=f8WM;}aRw#p1bG z85@DofnkgA-g?X{Z9Vn%I_46tMm0YBQVni5wzf5`R)TkyoZ#VC6=2}}>v`6$a;SR5 z&CFa@3eG-%E+u>{fel7SJBcU7u%I$!`0_RhtW%{vas?EEvrq%gpNxF4riggmS;z&# zvP>@LjU2F!S+{Y>&H{2dJ!_E^<`MKBe>gpo4wBKAHtqzbf$%oUQmSq$7!TEG{=J5J z>YJheFA?Mvf5G$V=eXw^V{oW)NdRHtR7Z!(IN(l@x4p;}3*5Z4R_8sULAd7O(^ifs zNN&kI{v|5{G;a*VYcD+7-oq7TW|Bbs|h{jxkg{{l|2cw|G9P4JWZv-^=A1oD_7=oqS zAD_3E55m3C8#?#UM;{rz97yvJ_0%nT4=T(hFh1EAQyJ6)OY+O_?~iqX;BJ;Nf7wou z_@eBfU)~OJ+Ieyzs%;>(ej<&DrWGuGD_fjN&A?+RJh9-}1Tw#m@d;RyAx5;)+D3;cl^rYeO%asS+T-83KgMNHg65_5soQdajB z19GbG5bf%#EKvP0zoEi_{?xYZ6Vq|&P$itVepV+9G>;yQD>$7B&Kg_;r`b|KG2_eY zJ9NoFrtM7~*`5g81EpDiPGfGWsE#;c8V3X^Q~9!n7$BWr^(0@7hJ|19?+wV2Kn{GY zz;1?nK31lS8pp$dRmpql&r&E5lg-kzn~0#Po@~(b_8VwM9(>I8Dj2A~ZK;s)3Iu{C zUDSuFuOKOJ@rR$oA0$s4unHmgZt`lASN~h{icxxHZE(*#Y(Kab`!RPG>VHit+>o0A z*QoQtejQU_d9za%EwdWsZkwOkrD2EP~f9lEGI3dR%VE-O1nfV_)}!F6B= zlq!ph&*l$;qU^T6zYY!_A0EFxQv=DV!}m4Bt3mkf@cj#6mEckHTJy-!3aEM|@X(x824$<# z$|0tuAX9q4;%Z1lG&-5gu;s&}r6bLE%Dw$|&qXNOcSRj=rMz)wvn`hd9qfJ#2!p zgF>jgVo%XU4sGcm3F#IeB zR!O<6D)%yoVWHg8lPU$3_UhP|OeLT?d}q>Xf&>#6Z2lT0V-9`ae&6ZT0+6w_n`Lj$ zgXEJ^3}t&Uk1*%TvSg7Boq9jbNpqQyY>>~VL&yNh?b_z&1=1m=rs7D_cq)*6KR$er zi@sFGxWDw8n!%w>)=F%Et@ z5l6l|j{#56Smq1SQAqfFx`AqO7#8kSw9FL_frs!K?W2G}xV(PW?eu#J$V=WmEcdq` zXdaLM3f<8MGi&bQFD~`K+xT(oPanGATF{BJy^|f#f6M(5moDa~o;!tc{zhMFBO8Jv zTA)n%Ik%opGqlg23%(=P1c@S|VQ0C>F!8S;X+o$0@A<08bXV)3^Mh4VmPHLX%Onr8 zR#btC(}Jz|nMx3hxn5fXuZh;va-h1o8lZGpLe;i^++~Q{(SB}{4Eolt)CpG7Rdn3 zSlaXd2JoIw@b+1c&{R;I<9;P&l>!{$Zy5_TkW+bP61HiHKv^&F7Og`6+|V+_Y8~gj zq(2|TF6Q9Fp_*_58U((Z|o>FpakxKByocuj8E4?q9yxQc|{}#Mjj0|(lJunaSp(@XLQ)XfC+P8yu zL}ySxB<((zGX)_*a(A_^Oo9aS(LeWE#v!w~Td%`p4EAuH=Kp+X6uOmCOijCoL3WOl zuoO82eRJWzgzN{wbD5Td+5mG2{ENR%)%64A_6Pl~{k@Ryqtfln+yfyj4@R8tbpfNu zvTJK@2UyO^4W1NhhwfOr1m*fx@OCndW3*|3spgG@of6I9?(#J4@(%Q;LP{zc<{E*@ z^bK?GP(A1v@{4|0tp&!$X>!rY8hEo}Yv+eKbWK}_36c2<7`fiN)L~T)yo*DLog8IA zq+y?qB9}n6+gXw7;9_9?>(AW#mISLVyv;nLg+M3{Nq#(^57}zPZw~Fx18Y%hTCpcN zAd_kEAf_q{tXs;!MJ5v%FNHql7)%GvE25?G9%RGEll>&?}CQo?n zPX^Y#=Fy*o@b0)@m;TMo1X$Q-E>!vsIpuMXM!hQr$o_Qd+Op9gY0MGG-i*1L||D4o7qJrDjhW*BJ)`cZJ@ucGo~4X8i3< zC`BDLDI{UH#~(;C#|45NaNg_w;ApYM2R|?SQZJ$xNM1T<5OUiSh@rI~RM@@%fyG#x zqw(`5uQqx0pPEZUici0k78gB=3JC@Z)`H5nJ9c$nc#pu=y|qtS)Kw_ z^9BT}&;76=aLMOiXfIT?WgKm4>V^-q{|-oV;=a!(?>P6n4oDEVIJ=Jq=RPj2+cW;H zAoimqRav?PF-KGY)v%IItWY^bJUDOS-ckxx9eF>FyO%(|=mQh_$Hmae@VqSd z3eJ81elD#U6@rX~h}==Td|(tXIvkJwxgzzNiwZm6_%8n6}zoDIM0hI2m^%u%eN9BwpPI)4yd}^4=f`GtYzgYMy0Ek=U6!pcv0^!sr#w;)Nr9x$k*BE?3vwXTI z)XN)4K|6CVo!!&qu)+}`Q%yU!!(m%l1u-$%VQi0Qbm=nog4#q z-v+6c&Ji$<@wqAOJq(%lysk4EL$F_{Vqe*XLAWD0GWeN~0u#5siL_ns2S=R>&seoy zha!M0h7lj3_w$&WS{+z^?8sLSGCJvOF4t-yqMc7@flv*`=W1+9x+{U#CvEbm zVgM7?vyo`T}fVkGvZZ}^5odPdx zb?EXzFK?*)g?KK+#5lScyJZ6@_)dlKdL~FRWZiXj%z%Zk&qnJ<(}APZZj!q{6*Q+` zd^Cy2975DBH9NOt;AZsnnf{jur0<1UN49sGK>+&m`rA~({gF%W+5#E;kYle&HjX|(tkIg+s^K7@S< zXU68apUeXHLJjlP?r8{#IT^03JO%p?Rhs6uPQZ2YZR3mD;~-F0D!cFRD3IQTiv-4w z0EPBWmEO}~a2IrP&N@8=>{Z0`rYwVSb3w`Q+sXjc{W(6GH{S=m>@wH+8G2zOo?}PK zg>JY<`!K6!NF_#HV?R? zvsX%Ua=`zmrC$RZ=FcPNIy{{+p;KSe`^=sUpq$DYoy*621%f>qc0VDMK<(l0DrEA*{uU1Jg4NL&V3nmD*Hy$^r#6j;b!#M$ z0!?HrjPZ_do6#c@PdE^$^h3KRL&0NInp?9B@Aw>88JM-d0jcs+-RVccKsxg8^@&hl)fHuVLH(fNX3zq~OAdPKm>+Y7n5W8kT@Cs5Vw&OOMC zzEsw|y5Che?^&`kR0e&*dqI|D$1kox)sn3*68RBx4x(qbwmWU|YLi$0DR`AFdtUKR z-W-shFfs2uHw#^cm+kIEPeaW8t^_5{DR^+Em$VQx0dbwdUH$^&fVaUr@05?iy-0=Y z=S@anjD~57d~z6;DsG8fS|5Z*Cw2O6cTu2jd<8T<4*=t_0G%Sdmzi09%=$H|2YfA= zHDmg_;7!DZdtc9VLe+Uk8!-DA;NnGE=Jl8=)KL=$tC=goz=srE>{<>~ zFLUg8c`-L7ZDuCdig(9+%ONoV#V}HGfVIMc1lcQQAs?%7-n*wImsyQEs{i7Fn3-J6 zYd+d~Pzd|*?7yuox@Uoy(EVPKZJEIT^|dSAw{)=nuCNp>mjv`aJAdjtif60so#iPd5G-0x|g>~tOS3<8p}d*sm^ zocFl4M_6q51Gn=orfD@l;0~FyC~o!vQm@viiZ;%Bvmb99X7mE$!$5rIUw z9zb|~@xp^4)KO;dhjhwu-V3$wAXd2mq1s#W@RT!1jxP!B<97l=!HLdZ!uw5LZSv~> zS6-=)e_4wEI0vWNavrw)o&j;L`57g>X&}a+r z*zYh3dfh6A$4`xb`r9)$&6q>CZaYxjn>Prl3pFRJ0@0V6o7_<37=SIwKgGs5`ha;< z`o~_e9mlEx(2Fmr7OW$?t7`f%fBvnPd|KSR*D&`y4WWz^3-r>kx8k=Owd>rc_w<%rNN- zNCDx_L;h;k$v{e*w@7YE1P^-MRgQMlQ9i}~_PF2k2r?Upg;=1Fxy!7G(LlNv;_&N0 z6c|5n7`ux-c#4a(@gooMK4z(kvV$`WJTzXgX#K=F?|#jd)pYdDpGw}m5QBNnf&u+A z`KY6q_8d=|Lmj21K9q0?ee-zDf6~#o-|PD0Koal)?xn3GuS&guvXTF0*%ar!p*o{l zQS{9h7d6wiWiNrfnd8UA4 zZFK}HyYd~zO7Av#waKghPkCkLcj@Ft(>yo{_PMyr&VfA1%wZ;V1_+%YOlt?GAvw+I zPrutFd_JC6*>gGo@#}Qrh}lM z^GR9yE(HkPMa|kC{ZRMdE=An87dm%NHWlV}!_NTr2)B(+IHb;Lbw;fN-o%C}LUS7w zddwY_e$oobSCR+&(3j$-;b8q<(*&y@Q@?28J)hvyQB6~(M$oKv`uF#59S|=4{>-_% z7W7676o_)w5HmQtU`eb5)jK5S#GMteRP(vrFt`kQXDUtpTq%XlL0XgSKgFQ8rtG=5 zkp#j8Q_gpB-kX>+d2?E=09IFE?t)<+sP=mO+7q4w&eP`uU03ivW^Lu0lP=!JteW3C z@FN{q^F>mvywiYOC3$=17WUw&heV7X!u}SLp5QbM&U@C6PD?DJZ{8<*`Z;X^5UO;P zX+`iphTqGf-T`w`M`#7JX7Mg&I*{06iaDt-lJ2UknCIj%HVZ`Ge4!^}4~skA?`?N8 zlzM@=^Hvp(;taeq4zCZmJXA*UvBuc!9l z$halu52*jLJUDQvGo8zra0AJ!nR}gPT!1vf z39Bq0fpjUW{@5)ikkpzSRgOU&rTcKl-~mV6BYb3(k9Gh8O}g?~S^G_1ZSv~>6JD`B zy7J(S`#i*)@KGsPn}u^j8plf?&cIwszBh%0zLd#wEQ9bQ_Dsr$ANoEH`8H9}k0r*y z`=x8mGWzGHjgh8!QDmYCEwPPg){stQ_z8k}1@6%lF$rzI1M8d9DTctwwgM1~-FPNR?c^ zRue>|JxptqAVXP@DO1p`1`rFrCPyb<2gV)}xo$6NKu7j|Dr0;VMow-ha}Maqk{y-GSLN-CEsQKdmu-GEG43Ff4j ztM)4d;9blEb+zSfN#J2TTsA@f9Td|aMK18*e(&hQ*9?tVFn%GPx-}PjrT*$Vj0i=6 zCbMC%c2NZS*FpzNbkR2#xY)13i@9^jnaIo)%tv|i_hx5*14$nCH`J)3DE+|HnvHic zztv97PkaU9D$~$8QOrr@uji$PVE*7!!jum;&U@6mHHC7~k5bJ@oYVKj9O|-)<|*Wq zo_FerEjaI$yjJ0%d;(I#$vAQw`cWO)ITxpqPl~^K^$s}$)s}r<-n>K|bwMw$xZM#* zKHoXRR8dF8d$)!EKt4r=mcAqsfH2DG^P=+YCa*Sm^}m!?p?o&4wu#Tfth5o8p#Lnq zD^S`Mus99gl*Kx?yHlXtUH%}ibOLB>9hkJPjKi(uCWdR|QSgwo`*8_ya+xqwFCIC$|a zwOu{-4@ld5#r_t(x0VIfDb=9)`gY#4^-74bDpc3gt^mEd{jL2YW#B$*-Th2)0< z^V`nhI}P2M*UbN6AHtHKi*W(=;QcY$LJ2DXRo0FS--JA1q}C*Hbz*M|)v~Y<<~f~b z-}`)X%>*ifxQy|R4B*~%yN#5P26~|drp;QZpm+9@JKf0?kQ{bQXu>{WlGF6f;D6sS zH$hjRyfpz7ndo;ZqK=|4-IAAi8v~RT;ir2?qkxnrmRo|oQdGH@#EhvUfG8>T_*q&Q zkn`KjSYD$arBR$UaSu6F&2xqi{V2&XSy@+C-0NMO{1%viI*ORC&ozT{UbO}tgNz># zFZbTjFGN1wWLLFSK;Qgfh2XbECg5k9{s2f zEF7lksH3h1M#}X%13@ds^A|Vz=FFXJG3ZATCgoXD*4_c>PILHs7t~Ra1d*)M$f>{k zX?l9>fa?6Uh3?n3KzLj}fB4GlOVFBZR<3KxzG|3*Y;(1SRGwMjNWG;SZaocu ztr!c8Iwv9e;}5z>=?TcLa~o+X9|M=rr>)nq$0hyj7lPaRF!<*Lgyp0Tf&4Y^{Vs0? zLDi-6rK2hZma?|J*9hitGkcYYnp>luo$$G5ofacn1*6{b%PB z+JHUzblD4uR!Aj#a-5-H?)>OrT76Ixlytm%)M<$KG4*`ntFaA`-jk_29#{v?7n>Y* zr_=!ZM3VG{@hS+?pT2tLVkI!s6E+sVqmKGuS?w=T1_}R4zTKnX{odkWt{vt%XDSNs z*&AbDXQY;OSyCY+4~sTCL|{(pe7#w0elA%1_8o0p&IUp-t4_CSgV_|@_Zb(Fyyo$Wer zfuNEf;&kW@P|Y$&b#K_zaac^VY`7vm2ERsCuRm@W0Yb{VxOK~6AW73*SU58T9P7q^ zny3eX-)cI)O%_iFVN*b=2!y z+pqj|fgZIw4F5){4Bj#ML-y>>4 z@Iyt+cx5%1H|C3IFjWCvqNn7dK?PVhidlzG;d@sWJ#Fp|rJ(u9>%z}7CD7}-vW4$A z_MT^Ly?pgM_H-Wbnj64AVP}1*o5o&w*t<}%{4x)9RAkX3TiR@}4o~5Fp`Qtww0}Pf zOs4}krT%lnhcw`DxBl!QkqXvlne58xQt)n0@;<3D35bI%YrP4$*JBF%aSHc(!j(FE z&kf@HR*E4Xzn;Ln`OIajQ(jR(P-LT<-;FuWv-?7SrG$g#w>xh>J`Dp+F#~E^A?!g2 zp%bMr6M;}4=WDVa3`AZZ$}6fMAclz-y*!3JEl;n+F}&~x3h$DvP@OMGp1&IT1XmQ8^^%wWXo<;&tbyq6X@}rK@ z(*AX&8~OBi-e1*nqu``ee$Pc@ z1U&eSYhs3mz`;jHpAaz!I`z+mY-}i?ofhhwQrr*E-OYNV(|RFFvpbTr67!tALb>;T zcY?}}JDlEf9k8wJ8R2kF8w6FRCM#cQg`n!w-`$3qVd-)y)0My`IK7ysaJJrF6aKl`)e+}sUU0@8U#XjL{7CB_jN?77PY<}lcIljM_+r)==F>k1+ zcD~Lp0l|XfpFON`uP0V}uI47@%^g&`WuFv6r_pa-!~6Nrxj3}^$S4;&>s#o2vsIlRX7q6p^GWyi0X=MRady?$EE6lK8 zu!Zx)Y12(!ZSv}WidPJ5kEp&C&cO;LmrrYH2JVep8#;M zEz2^aj^cmw=cHH;>XNPJPdj2iXT-AlzPSvLDQ4hoF+)DpR(&$sl?K@?)m4XUQ-E+o zHS*$Dobw*&9@bh-goTighjxvk4<+Ou)w7H~)IJ$DaenMSII=EPXM^urvD_UpoIw4g z=03L96z^b|Z?%rHhXcvztni-3P@w!OZDag~`OVfKce5?nch~Y*zbJ!rp^&vDT1%sH4_TRRlL8pH8l)Qvz_G zXOitPwdf86KkH704$PPLI1C5&x?sOcI`KXY<~LifZQU#F1fbk@*;5@7`d z@t;c6t;i=g!lEg+Idbit?`vn|lfwZ+!Zs7sX|mK?ZW?X!YLi$0w|JF*#dl$A!#q5G zl6W$9Y8FO#hoyP22fbm3=H!J3Q_yH|(CcaZ1Vpi4*R^LG2jZmqi-(^_VdGZ&$-O)y z5cWQqRikDI-16y5J)H(2`0OX{IrO6nCBCU~zV3&qEk;12-wW5c)VjlLFmJB^yli`X zCkO*m*4=gNmAY-X7^;Ww^3^K|^8ILmpFurNC%u}1J+b%xPK72|pqQyhbCNNa&`2FY z*8r?79eH8P*n3A7o}1ZL10ICSRxDDu*E`a2A=keG|K3}fRWV~eN{uis5>yI0d}la0 zua|(jxs6>RJ?0O{7Y|Hk7a^DA_zgc7fbqz~abdSSkf9q~a7xGl)nmHOtF!oCjG=+w z*;|>A{YRC|R+|otwux&#+G((`A)!TKO9c+EGnS&=$>5rbK&Wz#31Zisy1i?uVVeGjA$VXt_uBA*5+_tIJ%Fqh@0VE-ER(`g;iZQ{tM z>5tZnl(b_`cP- zcv0Amdp#|lPfg3HpW?aYbSLsavPZe_^}ih4*PD#mi)CZKrtFJ7Uov5}Z^Mt80q=^x zi`$&_PXp(XztqK2c&}H&MjAvP%0r>Yr#mMJSQ&}?Q$6uc?^Cae!IyX-$I?n_reR<3 zkb?cL73>SP^P(2Nk9T@hS;rDvagV2U>yo(v`cNtJVrd7%fcxBlTiP(rcPTgQa_lg- z`JCtJ6AkpqZ=1i5c@hYcE{2bOeZ;;4uAf5lji{f#I%z*Vi$2s(#+xTokxy}P_O>_C zCvRYrbNY??iJ|&d6bbz&yOOFHQM@D0e!YB)4SjNp>(bwLVs7(u_Uq$kP(S4_xFwpP zPriTl)vgBAPn;aAm@i(uR9H0^De_#XTO!wI7KW(1+UobKky$=#y8iQ~1YU z1L0BmKaO_PPu+6=!pD$LXX`l<_aL9l3$|a?LH)Ge+!Iuf`l(?)eN_?pwAk`7ybJlH zy5X)5$S3E}0XZx=oMxq#H*Repu0ZV^Dx+{L427u3mbQj z9k^3J4YfB!-WmOw1RW>9Vos1j<=7<{ecbuaT_AK~M&Uf_79}Pd=au@Sa z6&!mBY=cm-C)f8}%>YDgjVm|v>jP`$!FYR@9#~aXEaDIC0+n3{9!Yj}0C#YfQi58ykE2X}cnmSN1WApH< zslh%gwzz9^Rghe#&!l#?68w$!X*&3o!-QqE?`QO(7{6wGwn;C6s$COi5A}<|LyYO7 z={cP9qWlF!P8EWt+Lzz{`|u8iB9nIUIQE;LSh)D+5$2-kgieekXF`>+mi0*y6YK%PvG{sM4|7}Uap1kY+e$mJ$i{C!6|M89$rDCUT5$^G- zLmXZI;(S+l{rDSoynh*3ojKeQ0m+s1R2Il54;hI8@BNr77d!B>liJ6WF%P98%xP|oJ;4WVRxlQze)4eA*X6=odFWApt6<#Yy>;8IfISDq zt*?|1jD7*heLqu0<8hDIvBG+c8Tq6?AYhMqO@g=9iPj19pA4&dBDW!*v@iF>$swP* z4!taiMLzBLXT8XSdpz^BBK-{X$LIYd%`|Ynd;X@@mJ9V0yS%O71kQIEz87kcgoH4FLVLL6V< zMLrco3ab<#pAIMPJMk9zWckWt?^EQ{?dPMzMmp&87wvx=q`k?jO-R9+)ACEa z?>_|JIW#=3Dhz@InPF9dg97AA$NldX`ap%cU9`EQ2lGcG#`{LPfcuHL@b+DuAmjG- z=u7Q(sFDpj!-2V|f2Xae&tt#2 z$G%5#6EzTHd08WM3-+5E$f@>7x{!{)LxkD_$AbF_hPV~+|p!C^)uR{MxQuRp<+gr>T=)H@U8p8eE%-Mv?c$Y_% zc%8j#1M>xLh6X>sV-D&UN0^```r`?-g1S-Yj}H*Ck2#?~?qEeb?TY(5)vq}sspyYy zjbBUqh5I||+LEhkn1doHERA)%1A@L9q1+Jl)4_!X0ZyFniiLwVO3)u?8VvC@M}K^0 z#j2S+&UaEqyuWWCpI9D<_v;{^ZVET75Rp$Jw-v7bHUp}={sR++sGsVX-V6OP2CCvm z$$zDhPye`<@<;}l%L=*Kse^n9KhsbljC}HF`OU|JeCl|VUM`A!x=mf;Zh(9$FU$2U z)�BgRckl9&hsMzd5h`q^m~0Bd0R$+21Ot&OjZ@oVOvq$6h^uYCs$JdO~J)!Cebu zU>IDGcS8;P&WrAlYWjwOE_nxy{Hq}_=CAtYdukA(M((qg{Tu-4+-oLHsrU}x(!p9% z0N(3;JiT*IOcxY3-ymyrc0g|u<5p#Xc9`Rj^Uuc~sTlX{FlXKt*e3iqs}%2GTxzIA z0zQ)A;6H(G6;MW+o_dnIwi5Vy--p``qlUG@$T{#51p@!0IX;9LvMrC<{kgdcGte zxs+7yqW+G)GdanTIv($I+HM9M#dr4h40+T%M^1T-HRU!$f~MtKo5XGO$=COr!c-XE z=RH0?85{}=*CcGL^-xDqIc^EUJsuV776-Bw-n}&Y$CbRp9u^NRw)Paf9}YIQj{S!{ zQTZ7W8_#i%ck0mg{2w^q?b~tn_ZQ3=kVY@gV$Oj0`%ClmUYzrK?3&w`(0^je_g=@l zJgRqRuZPT`|Kx8VCxdx%0$XN$vcU%+%`J9ici>&#E7EKIN64u>M!JP{=V{w9D&&++8=ZB!KF)XaUK{Su@V=b8@}T=u z%xAsVk4t)jpQrAQpdVU5)t~3Xtc0BEcQ<_&@(8GAieA<=sBQA*aPGOk zGc~dX>Xg-yu;uu$Se;zav5=1-O^L`g-X*_HYJ8*1Rjq0m{e~!r`SXNG5JQa$PnPJZ@hj`s8D8 zRJUs9DE1wQ1^j(>dKuq^Iaw69RGCqR}@V&(tZ71lE=3&*QE{fXC~U1-D+j&+F2x5b(f#9%Ep2p&II`Ho0gO z3(S>IP?DR}gFy4^lV?ZukyA10lIeb!huZo!k!=h)<=)!)=Z+6(3Kea6)`dR#Ju~x% zuTW1t8Is>|_6tz)U9IS##(8gA`}3WePdHaywmh1JoEqU{wXZ}@<*q%OUd4Wc3h|r! zWIuppd-k@e4CIvgPe+|g$f<+QO6h&bsk+Z+`FxO5AKu8`Ekhsb_nJ$<2jmo;=U#O~ z%tKX&{pPkoPAOek^NYYdltDpQ;2iqo<~I9;?xLQG{Kx#d9XVyKJy>ImoElQ0X1|V{ z68DOBXnleEgr_9}pO8}K zduS3$WbQ9}sE>n6-Up}2tWi*YNbSZYGy+3o8?Ak$Avm_9?t9LQLFjo)aJ_np0<5+> zkL+UUhf|g>H5Dg&z#xB%`}^rG5V>G4zl*IC;=<0Y+pD#MXXsV>HzlpWNENhRaiayA zZ(ozXH-Y)h_Wy&uH~;2B{oB4jlqij&q=X1XA%zO3B6C88NBOa6$rQrwPxbEG5ax&ny8TYXU@wb|zz`=j?k@C|l$w2wBRy&y zz1P2@@4&e6n%r_UG28g$jin57^U-wT6OH+oj!&Cz?0ZWb9DRQdVvd-?t)#~N0GwhO z6)}*(J{G;(tN$pVFDh$yY5O7cMGbsE=_iFf<{bO?dYYjx%2Y=4QUUVvW9OYe9)Rz& zoTz=44o)TN^-pL)PiZ+Vtl5fOYKxcgjWXm>4Nfm>U##eep{4#U1F2a*A+M^)|PF)emx=FyPgva@tZi7=M@y5Gy!6`qLW^W$osY>aZ*Kx?BI+Mm+6u~KZ z@kXO%QE*CIfOp^|IK?9Jd4nl9<;cH*E)$#@WOA6~ z0H-GJzWQ+*oZ=ssqI(8TF;~p~T+&2+aM7pb_4$9i`X9n8{^~uv;?7gVXv53+0O<)L z$v{)TPi%}tYUmxQR~aD=EJ|al(TZ1L=WaLJ_)e#bdfx{;;HRRU&-XNl&7`7+sL--7Y3fPwvyV*SLw5nPZ=1f z+I%l>B!^pr-cBkuko1CWJD-fyk)E@!JxjxCNoaZVY2W^8620v_O@6``Vv&A_uJr9^ zvcPq!BR(JdS&AQ)rY&J!FFE~!tTy^MchWL0wHA|{m;DDK%!|l5s?q88TWFdXng-A-=z z=@ICtTjf*?mM`$Rzjfo54C;Da%T%qR=sU2#_$)#NdTPk(^LbbF9puE3hz|4}3|;+m z;|w_U)xtBS1G)L*^@`C>D|Nj*FF%6Zoa|s~;RL6IeDjQ_k((#2cIE6tT~F5C`RP3R zHgl7rRrtWE@h=}`?N;h~S&a85kW1b3D_?C1PW_sXlh^@Hy_5~z9tlp#xwlT81gEy= za4D|^r{)U=qJAKkdXwrAO#@C@Y|G;m1*f==vu<$%r=$+|?3}m)9p%3vjtZQzdY1S| z5u94DQRrC$PKmkNe!B%u)!*DzA$9g2ul~b%bwu!iVbkCg>EAxUcGZ4@l*$RJg&2*I zZEr)qH(QJlQ*m8|=Mh7s#BPeMZRQ*4?p+F7t2018^m_K24)&3%*(ReaPkM@CI1(yK4nD|IK(wb~!^dJ?qTjV*9@VQ}XwQ!DJ@>|mo+4vi-4mkWZgjYN^2C$CvA zsDV>5?dIc6xL;Km%d+7f=3s7#u{ZID5q8;4&CB!Hzv;|F-I#~`T#TNHH9iRWw#cS+ zMd*)`NHY5S6LZ9cjNH#Lr$>33tzcb%{QSA7)dv-CLQ%1v5SvFXwcR6H4SUQBWcGvx zF@jTeH+^J=P~ZDiLtEJGj`|y2x$z%2Qt;U@Jy#T*TCO_g;)VL2OG%mETJ#^R=__0c zbRv`*Rh~f^=&Agl6CpIhF;L%atv1#Pzb29Fg^5E3|?!ZtRaEjM%hk39z=9|@X3#TvO`wlEP z+ykfj^fRh>z^Tt&+hBc;Iaz)3pIokhIa^$!H%Dpzh~6>E26r-K*}uC*4ETj~|Qe z=ItWM?lXs(e|{y(Pc(}iC)s!fq^DqB1r<;jiLwfy{ghuk&X?N3I20IFC)#lA=Y(siJ%I#tKO5XG2qRJD>2T^`~DxoJ*w2=j%T2 z$tG@V!#`PxBA*I5{=?(;M|`f{ZM#*GMiS$!Bj1x$GX0@#oW23`ddACl>>QE^y8~ks zvn=$~##V1WnRrq=ut_FrZyc$m)hP-(fP0{O9??i#h$apTL*Hg2ac9p$|MA5wk+^sI z*5(-Ynb+DKy5o2peFz6Zn$J-a^<7An-Ij)U9gr{3RT~mXes@+e2>IV9x z?lwD4yhWXFV*UD`e|=F0yKJ==d!yJ}+z(vC{!QJ|ZlL!STy=;VA=oq<2 zmwfVV&Q#@7w~&RdZI0Fln#rMbvxjOg8j1EX7pote z8VECG`w;VoI^yE=h$g@QeF>$Uk)y`Azc=u_Chy1>VtnxWTT%YcsB>j|ODa?lW1^Sc z^Q4R@({?DCkCnhj&}_YTvzXjcaastaMZItHtE5e<< z+LF>)72MftwqVlyfxRr_6@@gZ*l!ofooLH|e$G4khf4>MOUeHn@9sq|wORJsYAWov zyLsv5>l5ggN_Vf34a2>Kqxq*RSTPSHeJEyYniuAlQ(B}9Uy`11c9p06o`f8_vjMz*V${%Kaiiti9ZYW z(SQ%!l;Ls(oVss58E_VyO8333V+u}vp||~e5uEyRaPMB6n_F=y=8rC(O$P)BNSxvmO7T1$9oy|t#9stK=Scx>ypN}~JNUfZ_&6Io;y z$W2@>$30HNV^5{a$YHB5;!4paBsD`i>keNrnM;nln*6ztn4O=s7k*_v<{5*p&m+>}g)Sz!_B&;GxFpXGweS0v7zRt~Q-wgjK zqHk)?;=NDHiP(erO0G5lb20TlJ$Q9K;Qmd&#wBI+b^g@SH_*lYgHzJSZ2jL6-H%0E z7^iSo-EF9&!T^1pid8)~=EJe)>9XBqSQzpr3lj^&Q0R)xDmy#uGe0i(_}V4xJuv*x z`&I}0qg*8gEFU31PYxMXsX{IlRi0NRj60wiU%CsCpR<>DeAsi&3wq;8@P>8hJGVU@({Wqp&#rQYT4Grx>ns%o2< z;t>2^@t!0XcW_Ew?Aygt4utZ%OkI^9xzy6^{Fiz7z4gxD&kw=xt!`^m{{tQ6pYr;w zAUIVVpU&(8zsKtT!($HpoC2qzS>c;smD}4yLu|!SB`Zni}dMmkPW1Et(&k z%J1@7wHur|Zo*9?3{IVMx4&^6oZ{_U%aRYjCy-msb_kr>b5c!nrOu~OX2o2OT*}I& z>*85(YE3mgA2m3o9_my26@E{ruEDDle($6RBPR_wC3sdwBUBkV{ds%p2MUBz8#a8{ zQ}!RP{zG_m<2vPHMc)MZ%YX8~6?#c?&StORrEXHY=H%wsz)oT?^HTcwwXekGf)RCrb{pYV ze@pk>x`nvtOVpO;Hj(l50jrGmG!o_KRyW0?>Iw7EuKcgEb=Vu}@R6RnmK@gHp=KRh zMKVg?Wp8t+Bre?RnC{qpBJxfhDGCwgWFcrzNz7O&QT5f^(tNIjI9yyhFW*r_D75Sf z@%HGOs=NNXUmW*AkM7vj%b7=*n-t?{exo0OZIiaiL|j1pZ8(h0l0 z>GN*P@5x{1vY%5$-;{`Ufa^>W$+&Px#w`|pPo|iM=XX2_U8f~chr7iNTJ@o|`N+?8 zg&FsE;O;`;Xw=S~?}(~ipui4m>_3=%-tlQ3b99{G^Ze8o^xgBj&_2naHK~e0z|TfVpCw*S~&!f{x12I~hU)PVLMr zI;jgz%|89|aiy=bPgZ~aIymJiF3rk>egxK-s~z$7#Ni!7r=Je|UaV@_4H5Xgs2@Tv zgyHwj?n|^di(Km58{2em=qSD}=iT3-qc)J!p|`*(uOM0tD(ER|!8?Dn;rGs;mYEoU zjtXh_wDpJI+f;kH*W4I&JZFlNIsD#)Ncb;r_&rnp!i_D^QPKPziDKZ?9d(uw6Xa4I z-~2g4;rF=o7IIgEQxCUSgap9vS(hdEXd#!<@b9TQ2~Krn6ns4kPMzCqtFrYBzHdIW zuUH9lIbnukY4U_p@z+tc{}i|-AZDN_{f}4w&3U!zk9b5^sIwk zXM`lbQf`~MJwyyH%$rZX`$m30c`|FdYJjL&hWvbL-$!`uwJ#*>=p|+^BUF|{x`||( z9Y?f8CkZs0okdj{#e2zDph?v2u4_6wHOZkQ_xn8U% z>x*s*I{vLCjgpQwTdZq{ytaPlR^BSIT63C7XjLV#(0o((jrJ2+*ew5S828q9NlD8~ zVSjKVt+oRHbTKJO657Izegr8qk*t=zg`{@x%^$~V^YQt(M7!jS{>}|Xx6fTgf2Xq3 z%i$ZDq=1{dHrgNk2uE(7%p1hMJO11a&(34cINEgNLoep{eg~W$bwe)2BJD1skVqU} zo*Fhdf&ARQ>4oG@>_OoBLic$W_O$#ub%H?}`MJOY$;(fp$n>n(KHgc>_a0OK)}asK?6A@MlbDa$H*kUz{pY&3vuPZtkW2kk&1~O+{_`~cQI~U= z-_ta{yea{G2pVQP*6|=emsl#a%7WiBGa{Nsm@gjN9$U|cT;ls&3@OwT2yPu{*NAaFWpj&|7lTVOsD}~=v{qW{my&QO@eW-E=IQ8Y1wRXTs ze81uPNJXiCy!vm-E8ChWpNOCdvS;S>dc~7tMCG}Od+(kR^6)qBFu%|c(f#~ByW8j+ zX&hU7hOM}tTwN_En?w4@<)rS{GwnUtzovf7j;oslFDP(lmvxXhyAd7hfOfL+uc-ew z-&Qi2@bsWfb~EXcaZaOJ(?mRu2vj*^Z;P05c8xsdWthwQ#%iN#Ne&-J(d!1>-%~iO zrsj@ZDu4FUR}z+>^NG~k@>4f)&x%)y!E!Jq3wtoFOK8e6h;GUI!+DG7Ltr}~DtR7$ z?{|FR)ED%1)_t)s$KK$8#8s!fHYY*1Y-|Wl#GH(7^ifX}_`UPVlHTWV$0{~R$yOhG z&3O-Z8~EVP-VNT{fkTnREY>wW;sWl&Jl(G=Ivh?mduOm7eTh9%j)ipTxL-wSt&c1^ zk2|2oA0F%x48|SL??F~V=EX%bDsP7rlvoXDf zj#{(bDS9LP-d`W(xFF=_rU!}ON$la=GjPV30#2E3qdOOa`krOgQiC;g)UnGO?%sip zy4vmM{tSL^-v;lupP-|7>2w}&fm6E6T;oB|Q3f21v_+S;P>+F zx#Pbf^Da zd*S!AUoKa>t>`Gh5+X}t9t>8G2o9lYI7?i4fn8&9Hd zYE_8r@25RPt6A#f*wZeObL9TeqHqUk-&9mqN8e8RE+pq$ZfYgBW?#~sJKjv$dP{b@vXufUs*wG(q8DOHM!l_t|K?M@x3@v2^}TQz`>{lzc*FS*DwGb zRes7v-roVfk$Lb&A#~LKf$E=j$jw!_<;4t;o7b{mG%!VOo}w@}90(n?o2OfU2z^r7 z6|q5T&{4nXY+PHAM;!{JNLWHgsqk6-x(pq4Qn6J2H*{1%@UVUXbd*8;w9N=M(vzH(5RRpbw$(=9M!P?7bVz zK7X1v6LTpx-%lO>NT$Esawn+oX`Pbibeu~eT5XB)u1?8hbAWL?g&+G`O16K9=uE)A z-Y)u$*w`q>h@34o!IgEA=I;w{8%#t4Vn$!68F;3tv zOx&mpn|B!cCBx`;Ku0;yM8%(0#r>+qEf=|XFvpj&F5P1me%5HhE415~e7Ouzrf)s2JI{9=Ac~p;DUrsCX zsGi(m#xd*(=Hs^+W_yf2%vVbJ7U+|D_$}wm5cst1v!=~GCveH|i4vkefTnssj6 zV6#3obd<27NxL3&RD6iekqP8c=U!h6Lw+&_7=N7v&Yt^_olTnWxD1 z$Fb@s)W(US`3vs6xe;RQ;2mt+G(;{sMWwDD{ziuVV_d5D4v^-F)CcnJeWW3Oyu)ux zFZn#B`*Y5$o22gI>)gK3LAV3emY<}w6GlgcyU`)7Mg|^p`EKdQ0{kFAa{WX)|1D4f9hX1-WC(jq++3Rv=ng;z+^{rQvHdT-X zam{y~Qe|ZLt-|SpekIUZn&L_8i^=qT&#Osc$fMS+-jJqVK+J>`8D@^B(@OAt4II2kOmvcYt zVUK+XWzmr?J1}qjI6=^|DwI%$a(#LZq5pvMA*)Q4`_aLmO}PUb#w$--T|$BHI4tI#J!x3!G>rVD&np3~VW@Tp?|(S>KI>zPIV z5n8G19Xa3Q#ErUMdXL=Z3GivDsB`^S_&k|6`Uw-@Q(>L9;U3iW1|Mft-iObN@REod z2A@v0i~ha`9i@Egam_C1sK*D<3J0IJQ8Dn0E_Bo#A>oa}&{33P2WfujsM%6WEd}VP z0=n;mkD;SFvVW(21E1)+)Y5fz;U7mMmL}lyn#D_fPD4ixM%xW_fKRG1zXt=Lqb8H% zvh1LvT<83=pF&3ozM8w83qBQGHrmJu9p#f=7)K8sm6)C%BM2Q8`fauFbLc2B;jAQf z=qP!e@-G8Y`2P6@w|MBNf+Q`TQRpbE=tl+jp`+x>KTfF#gKrKVK_-IWUWh|~5Z^yu z{TJhvxkoZ=+$2@AbD{%VZGMePpglSVG(bNJ5~=8yf(Ia%?@d?T)t%GVa%`pKa+ z+bpx@dI@zSy{o=+57E(oV##+EeVrChP0JRrAK~Y3P4%KS62F&YVpnzx_8WBW{od6? z{@$^BN4vj~thdpW84IW@|ORB|7vH`VZncxV^SwFXvXzWV)Z2-&?O89Zn4$^+K2Wl<+$u#Z;EE%{7uF z-oLo*-U58ycAcq($LQlUxzcz25a#z(jjgT6;PdX*$Y?>|T`o?uEXh1K%or(I3U?ZoXF?I;vZM`-kFl z^bOn>JiZJ2g4s@IQZqqEy&bt}atJ<;o>13sxf2wlI@ zAJtD!^<^G>y132eq6%~rYfI;=V(`h0c7{6!K2J+eZ|pDlW{s=x#+gLC{OcVWZ zJb_QPK}X4GGO_N0jym^SwpkWBD(-zo5-;j{R%aZ$S3^gAvADQ(H*}Pl2j?v_1v-k|-$SHb6uhktnOuw*pI+!Ch3TNskNc&q}R#zl0JI(qv1K-MA3c?cL^Wv@M(&&ST(ehyoltP z>+!f7Q@Ad^>tiz!Kdw?2FxE(RU&yiibP9JuyF0t7it9+uHTQ*my0{Pe$bCcKVf0B| zNInz6*aQO9>tKMO`at++i5rcx^7Th-f_zUl(wu zkR0Z?u~0#UI}923lFtO=-kvgt7n5N&N!+(KdHO;o3EiE=CS&=LSlwq({*;+UTEAHQ zu*O|=&kfC10pZD{mOq~JwKVQniFVZA|GVPzay(_P$Ko{e&z#a^hXJQ zNO?XUgnh%JoD3u2QJH%>N&%nZeK{Af?~cMCP-51AdY(#d zw1PY4_PS{nWs{IYeJyc)`~W%BRl6w88}5V>rKrX006ty%I`TCGIn+LBlgusXGndGH ze#;Gfifk>eV!=MC{$pRd%8-*ge);vB2|ka3`%UC0)bqUmYPUs!PrL8b9eWQx&3(Ph zRCk|H?B{O`GQ;PkS`}Zk2cNEZu2-Ulj@oJJn(+pF+LF4(YU|eb) zd>*UxyZCwV>6FbSPCEF!PoIZ(h{NaYxLGak2|m@?zS5_Lj@mWL#KZxgCz&anY!08N zHN(R{3qCEr7*q`hpSX>v_B{ljF4`D{QR;r5`5~QG&ppDPwdZj%BO-)b9@|E z%E2d}(|f;eh0kL=cW{#*d|n1+&kIHPJde7q%m<;Pbd3*i+lWJ_9zP>C4xhLAoc~PB z5xiVQ`x6d-M1!p-3n~;fm=u!WA~4&-$jHXjL&v4*c#Vdm6iLQ`z=l zuejptoReL|vtjPus!v}D|G9H)Wv*=``E^)h{QVXZdXVi2wNDedV<3IcxxIlX)@jR| z$YP&xv$~;XV=d7gGd!zDT|;CPa$lNazf^Ror+mu=+>3ekNpZVM1rap`hTdJsqgEejKlZMW%&nm;nVcyg(G&)cVdUk4V|q%fx-kzU)pdBO6!XPL z5_Q(Jd_?_wa#zI0bV8R;+ZTI1m6#2$<3IW>nJnxZ7l^V;B3kDqgq{c|V&A06)cIBM zMCxpb)`sd>5;_$*?pln!oGEmNb$Z@mzHZ@(B}Wt~IGs?w=#0JtT7^;-8uXca&@IqK z;J)7R;adyX(^9ZhVJe^jKGjn5F>?lE-)?27)DG11tUGo*Is%_3PtQPo4}DR4*Hmq6 z!amOTc}8Mq;qz+dE=Wv(Pv37{Uc1s4Rgt5AwGDjo-g&h@82U**`OcXq$jRp|{NDy) zkJOimS(ZBJr;)%l96ZROUaeQt3_=byV|&?y4|Bzijx3o?&`*!m*8LGcU)0&BM~?S` zPlDe!|IG%Uyd-!Oa>1u-H)}8Vflq~rc8R>;)0#IU-=2X_CVJOvHi1w5<*$w6z$d!k z)R1EEDJ&|?R~vl#t3PVX4nC>f1%?ti z@$<~Rj(%r@j~uzIbZiIsXn*w_|JHxJ`u~_$0&)+#{Szlh5kKeMkoHmHTAsAm13JoZ zlY?5O>>%-Cdr$R93;mo2_k7F?>?1jX(RWwV_L67^wp0_LZgTExcX*V42N~cfa{6?% zotOp3x7BTF#XMND(S4R?vKf;tbY~igPNC`3%dz$3Vty4@CSM)VvRjVzeN{uO7HK9p zwpS4@(O(gX*dttVSc&Ue0_I^ZdkuFLlo3yrwV|AUO2{qduurcq79$^%`MjkA`iX&# z{?Wq%BERY8i8puhh}i7f)oa;vh^p&Uf7fCr5!<%MVK+?%`ZArA_Z&$l1(g~s0Z&p% zhEuijgXLr*#jhW~$|s2oQyP|&WfF;QWyOX2yTPaF9sAex#}dy?Nr|K(h9s`0IubSp zekO{D7*o8^<+k(mUz%_&|bH!eGU4FjrazYLO)5Y^E)Pu{VWboc20J@M1S?0?f_5p zLsgf?cxF6DKew!aKJHZ&w0_mmybgaSv|ISYKr*;^3|E+0ajGI@W9ag})o;I=M9&d^#s>f8hn{cX#Qu5C z+Y)>#_)w-wK@MfHcPt=(LQDMH}Q@H5A)2cWl(RX6(aLN_@*3Uuazjv93osb|O6F0()9S@7q~mQkBw1iwU_ z9X4#j&$}}lqCxYISN|{b>K47qIvtHE;I-IRM&rQVc*qH!u`|K8}9a!i>p|;9frF}N794F8$!6t@cnaRZdN-?2~eAiWBT9BLRH6--Q`U(Z|_*H|It%k&5#xxja}%g4H6Ne>x&3FN@E! zxtB*m|18g(+Lc2V8kKcYW;4k!>x*EYpC1YH%dwtaT(}S8^h2ukUMfBp##9j_xC=A? z^h?-7>@{!A`*@5W_p9WZsCSQlAgZtT4UQ&YPOqLWcJxgQ{1jKUa|Zg#cOLM#HIID< zmmaZ3XuKtY4Y#I^o5P9fJwKD9ccGt@`?$;Zha&fo)ExhTK7;o~tO>>F*SzLCqL>YT zcXi#cU7J7AUAI42Whd%*A5+Hu*rKo8k2aHP*$Z_*uEmYP&`;u8D%59@KjpD!S4x9V zOu^PI8&Jn15zI{#@QE~4rS1TqMzti~o(G@gsPxw+yTC`8OxW@}lL9C2jg>{n$L03< zeY@j;zNvzq#(UsXuyeBLA>>a>B8stwr{SKZ2;FFHU4%QF4=#%nfm+?UUr2F8A%M9|Ti%P2=yMj;cRVF*+!6%tm zei>o#DUhm3_7wOuzeUO027KBj^V{nk{GDfX%{71cyNaEkor)CDXK?A+hP~iZMY)Yh zwhZzo6+LAq@X6|mip_l~{PzlI{M-jV>2l>z(!r?hyQt)Q{ycqh# zWLkWkaVh-5Ur7pOd;EqolXzEtO#kE7|HJgAu6eeN44YE(A0h07y^D`Lp^hx-tz zclkzH*8N5#kLb%9-R&p5+jg=SX7my-|Hf3+{XOKncTF2%n~8k*c6l}IKY0FQW?Jz_Jz=kX;4-C;x#H~XnY7^=a`^d3!U5$fqV1vh zcl#*zoc9rjiJS_;|L)@N4D3B`{EU?RsoUX zj4-^yo`-tciF+Mo*(>{Qw- zp8@IKF*_JR%s9UX$lSs`nAcwl>^6px0%pdmA%4imPd~KaKY_WI{H-=6%*dab!V8Q#mOeO9Yl%T(W*qDo$ zVXvLYrff+$)bB!kB|T&@Km681lIbS;H05>PZAgZ%+n7Gk%Z&M9p=|y`uhFMTzfPoD z+!1+pruEe^_&U`QjlWOPUqJcE`&Se7JCjq%BFYc2ADE>j(-8A6Z!SjFhutF-zT4H5 zzu;2Cti_xqeBJd5Bliv9)3a>y^1*HF2d3SfP76M9-!b4Qg0EwLT#^)NjJjZCregwp zU2+vwe}@5l6zAQGeBjfiL9XBc_`2bp{g3y6Pr-_xl%8K+;m}68wcyj{6prv@E%>qR z3F416z%dyPA(wNwzsy-#6NUPn-^crbf5s5Rw(Z`7Vs&0-$LS#6ZreVtIczPfFEx(JT=9^s}~PusKF<(1v$STN8l5MDr5J9 zPtkdkZ(j?7R}|h`KJa(76z6a)@TrzF=9B~YquT`HAw)DS3w)w9r0~*j!z?c!cZCk!qpf5_>ng4-658)Lz zKmFLQlL(TQ*b{^8gl_6_{Ngj*FBalk?lWyB%smoZ4Ni@uY+s<<-I97T8E9m)RRH}^ z!gJ@@Y*EKs-PZe{3wwj9hdw3Q;LhINRa+O3KPAppHSaoGMk+L-e%Uye5T8iZroOL5 zM0wI~XPQKf#d+M5lX3JlqVeOq&P=}ZQRzEnJ9_!qt|v}=O* za2lyimHBx4MG8qg_q6=`dhm(rp;=}O`ZV+GGdRm|_jAj~k*y}!!y?nZyI%S|N&Lg| zRer~&^f!rzmI@{|6%(*Pq6k76n5zy0E1TK* z&2J+Yw~aKpb_aV6%Cnk33VR@j7QUdy2HoWFzI9&64Z6gPKEW9KHnR>E=VH!1z$Jrz|gj$ zn^ZsX9fWRjNM&%>U(rpQgV~NkH#zj3d>Fc-o1B$Xbe6!SHuiy&@1UEWwTQ;pKsUXY zmGiKKZi-OJVGn?AdKTDO+YT;q>063PKsN<4GjrsEO9ownONr1;8=0zFSm5irQbW(A zUO->y+tIz&&`tU@+@-jrKQ|ih^%) za-6+~;rnbOQ;#DbXTR2eIUMpXV|@7hWZcj=$^-p0AuZ>}3GQwE)1t_Z@0))} z9|nI{kU+V768ecx0F`RT8xtW9DeI7iiZfNXkP=j z`8T39_qs6Bu%C!ivD>Zj?Imw%_R0H?cN6yQ$A6!s=_JF|pR^^e<1XlHw)7PAmCIWF zHu9KmBGYR%Eb?d?NxcL0rYMDaQr7xPM&)BIF=+iTTs|jq?G*QY zGtB=M(_hy^-jsLJZUYVSrgV3y?gnt_*@2~!&)`z@YPa4&a7l)G|BJoQP07F2F1dh9 z!tJ5dRC?geh-dpD(OHo%72l+0duY1EWu?cXgBY6GsR&YsW<74JgaEX8C zCjoZoCa;}Wwpgj+`yBMI606T5y7zHUq&r{0aPNA(MxneqEBmdkr?#e1Ve@%TpYiSv<+ z^cuXjQHhH1cMiX24!)TImlny)lQz85@1C9i0v-iaQa>HW+fBW3WESsvI>+a~^C%Ac zbA7of`Tx^6^S}G&eaWQldGPWCvE=$#V~xH|2H94ld&h?eRl>>QQ-cHK&HPkMc6%Sm z^FEOIoW7UXb-S(p+0{j&4_B4wEq^7cJ2%z_-D)GBuXlFQv$l}*9gI6`>6*yJ!L6Zj zyx12d@>uD-WgWRy$v&pnT|;<3mA!exQbnSbgG{{xK9li3GV@yM6~xLeTDau|=3E-q zJ&M#TA>tJ?R6pX-KYs7BoGn`+iFOa&(|QB@Hs9-C32x3I%ta?F)gNXN)*0@)d z_}cdrt!_FAh!cGO#wV482G?6DPbU+>Bh!zqZzmCUo1QH*8`1Yvoz3Ct4=$YwFv*vV zC4xT^vy!-DFb~42v;+N46ee5QoyOQ>Kor>z7DeFR=BG+_N$fE&lg!`H2E8=eJru?M z20qT$aAz@c@NOS|HecwaEoU$0ctI~QrD$qqAO~kGI~qleyDwL7td23o++8$tdiXE+ zxROn@TUY9I)I2ZsN--DQq7wDV3j68`$%wm_JLY}_!?^U&N3c0lMel_x=6o;BoaupH z%CDZ0P;&-1ltYucp_ir`MGc=oFFj*Pb_sx9I^)-yt^5%6ccJBBK5!{T-Ms1$xU^Zb z`<4#$(jF`-h=X3DZFo$N(h0augAEB2v-e6!^ zxuZaQcbm9_OG}lTQsTg+-{TxIGgmSH)_*JRJh;RkUGuI5dWoC%n)oR6(&00rQd-~= zixO@559p=9nfs3U;8Mn1k!~`$RLyJXSqUzAw6flz1DEWkTwh%SmjYIA@hk_I_?wvP z|3EMCrnK$41}=@4inde3&xuhT&#nfS+#2RX%fKb$fFCDnz@^#+DOw$HDZnfH!D0Bh zfVwZWr$x}mVDIJ;2`&l#EV?Z!1inPh{`di>u=f^Zhjjvg_9!dNHkDNHfevjd0 zI3~khido=s=r*PbG>Ake1F(hIJ|3+T$uiwZPL6$OaN@XnxtE1N))A zddV0|;U0{-LRrSEDx%EUUNMaMmes~tSKkCzkYGy9<8@ZxQngm~-pBBBv$2^WCD2O_ zLuR>rg(Oi_kUiWC{g_FM`4Uw*#7(*E&g*MggttTa+^7ipF^#w+D#S1!9IZraWt2)1 zRi5Q!B`QKJnBk$MhoH z^?d6$_V5?xf$NJDX3nENcl%0OL1P$ERjt@GatHHv_bi?DE`<>G9S2J!X@Zg0h^tR_ zy&}3N)(PnhU|+#Zx$D<0&`VoWYrZEqm|o^w8HwgX(^yZJJ= z zivCNQzm`_5_UKEZFkQc63!RpFt?4N8CAa)eaX!qi?0L;%ERB5WjDnxO{XKln&h}b# zKqoDm-C)zYgM55%KGQUK#C-Igg$(kg$a%rTZO}4+j^PBCB+oVtuH;MW^}3_>gG85p(H@ZVLv&cX8xejVbIWr0AN=_382JD+M{cSdH)h!c#W82;C%o^*-lH?qoqz z!`_RRa)Euq4sVA4#Qsij$jpA(ZN=Z+rFm|!!l&5&HDrZP0zbL0E#l2_N*rIt>o$GX ze1%iKM!K=!ltb+N(b!e_|LNTMzx(sCPZxU%Moy4Qm7t4?c_V~o%f_9rO9u%{?&?z` z$^(QeRdBP}fj(07XKRYJdJh>lq+t!b(usRgJDe%j?L_w6O^dF<7SgzP<2Bu&Cel~4 zQ`Gr+0|_kH?HiL(N9Kdxl`U+oC05lIJAVmR;U1^3XvWjeM7i?WAzSVWvg^*Nn|4d3 z#A4rSbMbX0WLK2pkH4}-&ZYF#j^QS}gV-NC3#iu>7wqLV}>h)QD78~?8IaWRP zQ~!wEJ7Qi#8$OOsTko283Yk7i^fhOplWx9`KWv5lnV)BU2gT!wprzu}mr~T_J~>`6 ze}wx}`{AKE$E~*UT?!_!52VP(JV8DEEmL_Kc=TXjc;oOR)X6WFD-!T1+A~lo zePu5}K&Y4!bW*_0?h8MWck7E)O}QiQ{>xJ*cp5x96vAPC3_QBA7<%@MB|gW8+KSx3 zqcpmR(SGnKB*f^D4tPXoHZJ-DJo=s^>U$bIDp$3s?gfvk9e9mb=II`5$Nqi}9hs`Dfz6l0;87Oi%kz)HBgV5&x6Xn`!iw{3vV#>mY2YRV4-Bgnc|9zJJ<5%GwYPdtdCCbbB z8+t2T8rwjt(vMfQ-@!N^ul$D=wuU?k#ne;p{5QO{Yn?r&@TPAP;-AACbuw=1C*Bjn z3+#XJ@@NcIFF{A;v^-u#0hjo!rbh8nDz%{8!Wwx1gth*ufHLj6lOuVoQ!=Ablr#fzr)>lwk`%*Ebt=_1tQE~|MaIesK=j?>PMEwO)mmg?E@kQ9#9pT-xi40!K_B_`-U#gvA`ZSS>hIoT&-jAxtY#SW5zV1t<6_j`UQUYq zS?_}!iio_h#$BmM+l1eAFHxVP5iYp#0(V+2`<(vmfzMszs^M(p+lE=U6=L0>A5L(1 zmbjvCw72E(YV1kLzd221gnZj4`2qc?Bknl~?uvGYU;FH8-Q6{d~4-O?(JOvH~rWjwaM!wB{L+%gH6iCB%Mx(N zOM5xJ5q_;z=Y{ZD_%+qki|;AOk-9w3#kIk&t^0j#<5zGfbS89=3i?Rw;FyaN{90UJ zpJ|#Zbdx9TM{5=Qd;Z%xMM~g8*}?A;3do7Z#o9+?!PCnLrA0FEb07L5N>73}+cHC@ zq@aJoJNG|?J~IBRC2$$~XxenJ&=&g0>Ro@5NN8-oI8x)`AyP~fySi3Icb?E4e+_l0l)4h+mR^IUMiO5n9$OPXM8wos*!6r^lh?$%8mCY8B$XU3LrbJ-BTzUPGi%aNJqTJmv3=Z`O#|!I7;2viu zgJDHp5cUdeSi@_D{g=15&>b=WhxWJB?6vbHp@(KR<`sGqiqbU^PCjsGJEvDS@@y%( zJjVCM&_fUM42;ylp}}J2enxPpW#pdg8stVLscZN+F@N$>fAc0I=%HIWrz@-B&l+0F z6g1(_7&NqO{=lD!i;!OyLUXm=ho`#2_E>fND0Yp+UT1MysT%I z3O%&+Q~t&Qa40nV$dwZC$F^1c^MomKqCXFxKY||WyL5?K8G0z!Y1BdjdZ_X!&m)a% zgmU9Z8Hj$$l%X={`5x)a7e7dy2}n68ooT+S|S0?+!@*SCttxU$23-!3p z9WfOv9P*HJc()t>pGL1O3&XED2-X}vxD8yYbX<7GfUh4-vH5Sr&okY9E}srOdRy?( z75d1*E8VJK3BD~mo!|b4B!zM#t>0iA?;W=14BdFwuTfZEh}ZL0r(#(ig|e$JL#P9< zc=Z0eD|%^sm2>`zpDVeT>M@P?g!iVZ?|84e7{smc=>Uz)k$Jpn1}f$&oKie5{eM+= z=I>OlYaCw<8l<8`GAk8jW~5woJ)b zLWPncLx#lpK7Ygc;aumZ&t+NH^SOxbIUZ-1@#iitxbd8XBRWmf%0+`zHpt?Fj$Be6^Rs`BR_zN-fNUr7n03 z;#?Euqda~esr!@8r;K>|DXW6|pft{2YWwO6gNEuWKS zJf@Pp(gzjBu;0(Kmz`z>{}gB_#th;;T?wxHVJ7g;pmF=IW%O^WvZi(VqVb-{R)bia zb5hbAc~3(wf-1S$Q!o>Q_{YRIu1gI4+s>uJ-xGsT?~Sp)DuepapC`OiXA%EM6(0_s z!2bRf5B0QdzUZgselcW&_-9JGgwY1?*zc))c%Dwt#)BOga$cB=;@T)x5AWpo4NgbH zJN#RQV#AP!nK{i0<-$Am%5EdST=2XaT`jPOcSiCz2G7DfmU15*4kHiiT^#ga!8FjrNBpDPv(e`&yrU}X`1uXIV=-EH$^_o2 zkdpc}0PoDDBwaH?{4?j0Auq$pBCVu?9sZRDK#3ce(BX;q`uu@ z$eGv&`W!uVsR5*q7saJ8X*BhliD}JXWZ9u5;v*qmy-E@%{cz@1BlYj4lRC@9ODP^7 z1V}!oI~LGJ{M3C!vzx?En*@Z?h^GX)$2CYD%A4cEGs0k1*Lgy=v)LlVUx8K2l7wP+ zZ3PHTe!DP7{;$ZbTt--V{ecIe=lva$gbR<}H;{E}OgiicEi((}$@QYt*ZUFvZ+(@g zFPHLk+c&DB;sESI~ zj0yeqrJO2x6BXUkP)czI6}Be7L%q69B{MY$ecQf^$FviXhkf_k*teWVO_+zIpUTAf za4Kh}de9GPJ<^+1RgSsj2^}Tns0XRXc2V(q>6kkuJ=6X-gEIVGx7G1PD)sE1TG$o2 zCscC$R`wyi6!b6KY8b>}AAhR6pfGv{LI3rPOFOr=H}=yFD_%3<8HdVw*&KLgd^+HjIP$Fhnlt>ZsP~LCmeX<& z-)Ii0CX`&qb3dvRKM&9NC2!|9hiCd#O_F-xnd5ER-+v;$xmwkfdJ~=*SwE*J4bO~D zv8K4-ncL<$8KUqEU8O37%x~hkmar?&2L5=rN>34LG1~nK z`IY_?$&3)hH{6yn??VvZjD0W~$U}T%MsM+41Lfff zyw7Q&Qy2Ny-SuMO@Qh(uMgMhpMzOH9<_n0a3LPLkhk<_Kgb9& zU-=1?U!WuhK?_#O&=4plFQxqjB@79qn(EtjlcVl1(pu(LePh4^szJ*$RX z-=6;hA7P5as`KQ&dgG^aNF4Wi&T5I=Z-D+GUBcnrW--L4|DCT2`eikhhlZ)n6czn$4*BfeD_NcX}Lksq2 z-WbR`Vh;Iikn;3{I?7+T*UZzjnrc(&bsg}kpg62##hWnC$?VzP>NVBH6z8bmnO7V| z6ul)hZ4(Q5RjT)-*6d44hb?b$v*!!S@L+X;q!Hp7w>Q$lomuFUp32qqMxT&XADeXq zgX&}cUH0)O@+jul<4M_%DUpLsjjm5IM_?vAN-!PsEw#LLhn*9tJ|E4M3K7VkEOwaw zYjq4I5_FbZz6gDa#y6&}=b#UW?HvEQzA(IBA7=ZzC*F%?KD%fe9*B8F=eAc)qW>q} zF`gBReCfjH&6A!uPvk9|lWh$0B|Z0xrv0ebm=qR@XQE!?^38erlRNgo2Qt3=iFhTc z_jb7e;+65Qvz(jZ5tl3pJ_Y1So1Jb|-f>2MVf=4`^N3SC52$s>UqV05k)^wt@CUQx z_YxcUqjT_Qr9S*YiTk?R!5?cPM-HaLAE%t7CwbtHH4^$C$a}B;X`f2}MtvqAv_Dn@ z{&*dgxqcJ;F=72klO6sD5HkNJ2Y;wt%M-Y51P|Gr(VI5Jb9{LqgcJT?RPC|6rjIzn z@>cC8_=B-SYe$bZybyh9?OQGIv~V^#M})!V7pMt;JTWP|n5~Y_^L@(Q)R15A-tvH@ zg6G-5IFi(jJBvrOw!t6U3j*?%;g1XKfk(I%;G6Tl$p-L8^N4exr3|hgZarwZAHFEs zrk8+tMP2ycVqD0V!e(<-3K6d~+m)3>iIM#(ZECXqo5=A*6u-&FU+-w#5K{zFN)j2#qkyOb>KRU_u?-> znxdgX0+U9Q9bkXS1Q}JGUbP_2z-46%X!d1)F!4@C@%?KZVC3hv`Y!OA`{gG+p!$r; zCbHl5?&W>nUhq@6Ig2p4K~&`fn7dEifl$BtvmGJ#u>-#mik!aKN4{sm!Ezu}^ZqV&hJlq z+Q9=XzO{rV25L%g!M(ZNcI0=chlVD;0Uf2^wYP#6f9&KT^l{2QO@6P=uDmBC-gD#6 ZA8sQ2aWSG1bV`dp#sWKy4sB>){s(E&21)<` literal 0 HcmV?d00001 diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/r0grid.npy b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/r0grid.npy new file mode 100644 index 0000000000000000000000000000000000000000..5f4e75e31d867f6c5ca2998a9315b8516510993c GIT binary patch literal 72128 zcmeI*`&*287zS`fbRvyHIv5pGrWn&fWtv+ll}bhFAd(J*RZX!DrLfkjq*!(-5^}ys zOCz;*9gA(TwpD63nVdyMN)Bu7zp&S}zuoiOynn#=`o7P5Kc9Qs^br5RaJ4?m`>Ynl z#-+q4L@q*+EZ#|EBNW9aDpD2EOQRAMv2p)CpBlX+CGPX{l!ehrai9NNIZAAV6C5Pg zHo|n_fBR{Nd)k^?S$4~l>$4x0zUr1|{i&Mxb4`zYLyDbW;la1^nX14Cn{#^QQPouo zTUNc3H#9}7sgmEzjba^6wkCg&8=a|+cV7Kbp0uG>?2`Mpe7ybv+sVg1$s--VR0JpV zfn%i6&PD2K=x=`eyg{)Vx->M~zV=kd+!JdC^|Y&F*Mj~!yZ82m&#qfH%R>7h&^2*E zs-XsI(l_V|8#PcZkjT4B`eXkK@gL4B2cXpPqLCz66J=L38-)|K@Fef*`dQ}Ms5v8c zG8?Lcb(4R6S}o9pwDGIdgTn_xU)Qig&uS3v^v~nbmi6Z4|9j0Cvbhxi6I8K=X^cs%?$(8*4t$0;JP(Hx}-^dDrkxLh*D+VmqqO$}z~IsReA&bK2m-Sq6;gnS{MPceA>Mm!3^CDXQ9 z)Q-aG%KA8w!W=n8j<+V7SzwE@5}R*X;Pz`n$^1>0a6XtOf9*XQquujdGxWyb!Nx8h zufN7XyZPw7+%ge58tgg>wpd}iG^4jX)f!f3PnnA6kA!?#RVkW|O zhv~sIS0_}QODx$r(iykQHSWCDa)Eu=8bQMc7wGoPQQz@)5^k2++Gl>0VS9M4Q=5() zwoP8=rIx@;*OT?37L(C*s6ceTVlqrRpB;Y`;)dJ%&9)ZwxZz5{r)K*ickCXnp7qpk z3dCA|S+Dwdz<6cMv#d%FjCIRyYESb-*37zV^E{`*+xc12NJB5|-CEq4*ye>{#{-^x z|BE-IX|L9|miZvKShA;J-!!a%bDds+FH9b0@69jtLtb0yvSo+-vB>24T|w1!xEDTO zf4MOLHNShyy4D7w{oXx;4&xxCxmx*tQyB!o>0J(;^JXA0^`RhIcP6UL%8dt~pNU83 zthYNY4@Th8C#^#rLeSRKRrK^-2m;Sbj&xSf!s6p4GDTr1_TR|UKD8nYYmT}4cSnZf zr$=X6a(p6S{n~%_Mfq$9FAo(xmPI0=URF5WV-8Np=9mP8%*EFHyW`xGq7Y^_BH-e* zdFY%1-%8ayoSBxaa?OlJ#3k#hbc^}unW~H~t(_0;zO?}Z(_+x7U8;FwLM$5MOhw~f z#-c^5HhJ->IBXQ3ZojfR9+evhe$4P+fT}D(bn}>nc(6P-)=D!07dl*@9_mcMk^-%D zM_LzQn9@8at$8ude$?_wQ7yp&uT8zG_AjvK#^GLz0g31;s(a~ekc92p{cAjrCPAhX zakj!U8Ra9sJy-uI89c9eUh%x*dByXJ_bcA7c)#NPiuWtN5AuDG?}L0FCpddT#U z=^@iYrr%D#oqjw0cKYq?n_%Ar`zF{o!G0X}qUvVV~M zgX|w<{~-HP*_X<`RQ9E^FO~hy?006rGy9#{N6$Wb_R+JCo_+KN#cJr%px;ivoqjw0 zcKYq~+v&H{Z>Qf*zny+N{dW58^xNsT({HEWPQRUgJN9^Bwr{7M$oqjw0cKYq~+v&H{Z>Qf*zny+N{dW58^xNsT({HEWPQRUgJN9^Bwr{7M$oqjw0cKYq~+v&H{Z>Qf*zny+N{dW58 z^xNsT({HEWPQRUgJN9^Bwr{7M$oqjw0cKYq~+v&H{ zZ>Qf*zny+N{dW58^xNsT({HEWPQRUgJN9^Bw J|G)kAe*hZO^Fsgt literal 0 HcmV?d00001 diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/r_mins.npy b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/r_mins.npy new file mode 100644 index 0000000000000000000000000000000000000000..dba178b0df375d983ec6a4987c4693411186ffba GIT binary patch literal 608 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I$d20EHL3bhL41Fk9o-GG@}TSF)&j{oXUhn?-6 zGC~$as5{#qFgQ11Uzd};N8#T$b3ZxS-}xuJ^;?{yee*i!mnY6U*n5WV>lATyu746TiZXdUijB5(8~S_$Cfp_ z!Yu7&(waiFVl3<%OQf&&WSQH?-g{*7sL{;+&8}D7<}*$0_i!5D>DgdnpKxlo>Z9Yv z_Q_UO-dZn=?A>~;vx+#3?B^WX-?U4^(Eh=p`{iE)4eYaZyUr$j)3eX6zCZ7phpzpa zX*;xx)@s`y|LW7UaFdq(^+W&U553j2clwsobizv0UZnbu)XP#0`+K{)SIw=|x5RFGW`t9pzE9cThX}{B)#(z0uJNhi+)f+uyzU-0A6IS^EonUuvsa z$=Ej>Tg9KgLCXHvXRUY3)g 0 + assert planet.T_STP > 0 + assert planet.p_STP > 0 + assert planet.RH_zref >= 0 + assert planet.RH_zref <= 1 + + # atmospheric composition sums to 1 or less + total_conc = ( + planet.dry_molar_conc_H2 + + planet.dry_molar_conc_He + + planet.dry_molar_conc_N2 + + planet.dry_molar_conc_O2 + + planet.dry_molar_conc_CO2 + ) + assert ( + total_conc <= 1.01 + ), f"Total molar concentration {total_conc} exceeds 1.01 for {planet.__class__.__name__}" + + def test_water_vapour_mixing_ratio_calculation(self): + """Test water vapour mixing ratio calculation.""" + with self._get_test_resources() as (formulae, earth_like): + const = formulae.constants + planet = earth_like + + pvs = formulae.saturation_vapour_pressure.pvs_water(planet.T_STP) + initial_water_vapour_mixing_ratio = const.eps / ( + planet.p_STP / planet.RH_zref / pvs - 1 + ) + + assert initial_water_vapour_mixing_ratio > 0 + assert initial_water_vapour_mixing_ratio < 0.1 # Should be less than 10% + + def test_alien_parcel_initialization(self): + """Test AlienParcel class initialization.""" + parcel = AlienParcel( + dt=1.0 * si.second, + mass_of_dry_air=1e5 * si.kg, + pcloud=90000 * si.pascal, + initial_water_vapour_mixing_ratio=0.01, + Tcloud=280 * si.kelvin, + w=0, + zcloud=1000 * si.m, + ) + assert parcel is not None + + @pytest.mark.parametrize( + "r_wet_val, mass_of_dry_air_val, iwvmr_val, pcloud_val, Zcloud_val, Tcloud_val", + [ + (1e-4, 1e5, 0.01, 90000, 1000, 280), + (1e-5, 1e4, 0.005, 80000, 500, 270), + (2e-4, 2e5, 0.02, 95000, 1500, 290), + ], + ) + def test_simulation_class( + self, + r_wet_val, + mass_of_dry_air_val, + iwvmr_val, + pcloud_val, + Zcloud_val, + Tcloud_val, + ): + """Test Simulation class initialization and basic functionality with parametrized settings.""" + with self._get_test_resources() as (formulae, earth_like): + planet = earth_like + + settings = Settings( + planet=planet, + r_wet=r_wet_val * si.m, + mass_of_dry_air=mass_of_dry_air_val * si.kg, + initial_water_vapour_mixing_ratio=iwvmr_val, + pcloud=pcloud_val * si.pascal, + Zcloud=Zcloud_val * si.m, + Tcloud=Tcloud_val * si.kelvin, + formulae=formulae, + ) + + simulation = Simulation(settings) + + assert hasattr(simulation, "particulator") + assert hasattr(simulation, "run") + assert hasattr(simulation, "save") + + products = simulation.particulator.products + required_products = ["radius_m1", "z", "RH", "t"] + for product in required_products: + assert product in products + + @pytest.mark.parametrize( + "r_wet_val, mass_of_dry_air_val, iwvmr_val, pcloud_val, Zcloud_val, Tcloud_val", + [ + (1e-5, 1e5, 0.01, 90000, 100, 280), + (1e-5, 1e4, 0.005, 80000, 500, 270), + (2e-4, 2e5, 0.02, 95000, 1500, 290), + ], + ) + def test_simulation_run_basic( + self, + r_wet_val, + mass_of_dry_air_val, + iwvmr_val, + pcloud_val, + Zcloud_val, + Tcloud_val, + ): + """Test basic simulation run functionality.""" + with self._get_test_resources() as (formulae, earth_like): + planet = earth_like + + settings = Settings( + planet=planet, + r_wet=r_wet_val * si.m, + mass_of_dry_air=mass_of_dry_air_val * si.kg, + initial_water_vapour_mixing_ratio=iwvmr_val, + pcloud=pcloud_val * si.pascal, + Zcloud=Zcloud_val * si.m, + Tcloud=Tcloud_val * si.kelvin, + formulae=formulae, + ) + + simulation = Simulation(settings) + output = simulation.run() + + assert "r" in output + assert "S" in output + assert "z" in output + assert "t" in output + + assert output["r"] is not None + assert output["S"] is not None + assert output["z"] is not None + assert output["t"] is not None + + assert len(output["r"]) > 0, "Output array 'r' is empty" + assert len(output["S"]) > 0, "Output array 'S' is empty" + assert len(output["z"]) > 0, "Output array 'z' is empty" + assert len(output["t"]) > 0, "Output array 't' is empty" + + lengths = [len(output[key]) for key in output.keys()] + assert all( + l == lengths[0] for l in lengths + ), "Not all output arrays have the same length" + + def test_saturation_at_cloud_base(self): + formulae = Formulae( + ventilation="PruppacherAndRasmussen1979", + saturation_vapour_pressure="AugustRocheMagnus", + diffusion_coordinate="WaterMassLogarithm", + ) + + new_Earth = EarthLike() + new_Earth.T_STP + + RH_array = np.linspace(0.25, 0.99, 50) + const = formulae.constants + + def mix(dry, vap, ratio): + return (dry + ratio * vap) / (1 + ratio) + + for i, RH in enumerate(RH_array[::-1]): + new_Earth.RH_zref = RH + + pvs = formulae.saturation_vapour_pressure.pvs_water(new_Earth.T_STP) + initial_water_vapour_mixing_ratio = const.eps / ( + new_Earth.p_STP / new_Earth.RH_zref / pvs - 1 + ) + + Rair = mix(const.Rd, const.Rv, initial_water_vapour_mixing_ratio) + c_p = mix(const.c_pd, const.c_pv, initial_water_vapour_mixing_ratio) + + def f(x): + return initial_water_vapour_mixing_ratio / ( + initial_water_vapour_mixing_ratio + const.eps + ) * new_Earth.p_STP * (x / new_Earth.T_STP) ** ( + c_p / Rair + ) - formulae.saturation_vapour_pressure.pvs_water( + x + ) + + tdews = fsolve(f, [150, 300]) + Tcloud = np.max(tdews) + Zcloud = (new_Earth.T_STP - Tcloud) * c_p / new_Earth.g_std + thstd = formulae.trivia.th_std(new_Earth.p_STP, new_Earth.T_STP) + + pcloud = formulae.hydrostatics.p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio( + new_Earth.p_STP, thstd, initial_water_vapour_mixing_ratio, Zcloud + ) + + np.testing.assert_approx_equal( + actual=pcloud + * ( + initial_water_vapour_mixing_ratio + / (initial_water_vapour_mixing_ratio + const.eps) + ) + / formulae.saturation_vapour_pressure.pvs_water(Tcloud), + desired=1, + significant=4, + ) diff --git a/tests/examples_tests/conftest.py b/tests/examples_tests/conftest.py index 2b57a0daa5..0b18e0c3ea 100644 --- a/tests/examples_tests/conftest.py +++ b/tests/examples_tests/conftest.py @@ -84,6 +84,7 @@ def findfiles(path, regex): "utils", "Zaba_et_al_2025", "Gonfiantini_1986", + "Loftus_and_Wordsworth_2021", ], } From 29168e0e374a3bfb834bd2cb7a8707b2d3245a37 Mon Sep 17 00:00:00 2001 From: Hevagog Date: Tue, 17 Jun 2025 23:40:55 +0200 Subject: [PATCH 02/15] fix exception messages and more fixtures in tests --- .../accuracy_test.py | 295 +++++++++--------- 1 file changed, 152 insertions(+), 143 deletions(-) diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py index 75f57bb05d..6f107c8374 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py @@ -1,7 +1,11 @@ +from __future__ import annotations + import pytest import os import numpy as np from scipy.optimize import fsolve +import warnings +from typing import Tuple, List, Generator from PySDM import Formulae from PySDM.physics import si @@ -10,13 +14,15 @@ class GroundTruthLoader: - def __init__(self, groundtruth_dir_path, n_samples=2, random_seed=2137): + def __init__( + self, groundtruth_dir_path: str, n_samples: int = 2, random_seed: int = 2137 + ): self.dir_path = groundtruth_dir_path self.RHs = None self.r0grid = None self.m_frac_evap = None - self.n_samples = n_samples # Number of random samples to test - np.random.seed(random_seed) # reproducible random samples during debugging + self.n_samples = n_samples + np.random.seed(random_seed) def __enter__(self): try: @@ -25,16 +31,64 @@ def __enter__(self): self.m_frac_evap = np.load(os.path.join(self.dir_path, "m_frac_evap.npy")) return self except FileNotFoundError as e: - print(f"Error loading ground truth files: {e}") - raise + pytest.fail(f"Error loading ground truth files: {e}") except Exception as e: - print(f"An unexpected error occurred while loading ground truth data: {e}") - raise + pytest.fail(f"Unexpected error loading ground truth data: {e}") def __exit__(self, exc_type, exc_val, exc_tb): pass +@pytest.fixture(scope="module") +def ground_truth_data(request) -> Generator[GroundTruthLoader, None, None]: + current_dir = os.path.dirname(os.path.abspath(request.fspath)) + groundtruth_dir = os.path.abspath(os.path.join(current_dir, "ground_truth")) + if not os.path.isdir(groundtruth_dir): + pytest.fail(f"Groundtruth directory not found at {groundtruth_dir}") + with GroundTruthLoader(groundtruth_dir) as gt: + yield gt + + +@pytest.fixture +def ground_truth_sample(ground_truth_data: GroundTruthLoader) -> List[dict]: + gt = ground_truth_data + n_rh_values = len(gt.RHs) + n_radius_values = gt.r0grid.shape[1] + total_points = n_rh_values * n_radius_values + n_samples = min(gt.n_samples, total_points) + if n_samples == 0: + pytest.skip("No data points available to sample.") + all_indices = np.array( + [(i, j) for i in range(n_rh_values) for j in range(n_radius_values)] + ) + sampled_indices_flat = np.random.choice(len(all_indices), n_samples, replace=False) + sampled_ij_pairs = all_indices[sampled_indices_flat] + return [ + { + "rh": gt.RHs[i_rh], + "r_m": gt.r0grid[0, j_r], + "expected_m_frac_evap": gt.m_frac_evap[i_rh, j_r], + "i_rh": i_rh, + "j_r": j_r, + } + for i_rh, j_r in sampled_ij_pairs + ] + + +@pytest.fixture(scope="module") +def static_arrays() -> Tuple[np.ndarray, np.ndarray, np.ndarray, object]: + formulae = Formulae( + ventilation="PruppacherAndRasmussen1979", + saturation_vapour_pressure="AugustRocheMagnus", + diffusion_coordinate="WaterMassLogarithm", + ) + radius_array = np.logspace(-4.5, -2.5, 50) * si.m + RH_array = np.linspace(0.25, 0.99, 50) + output_matrix = np.full((len(RH_array), len(radius_array)), np.nan) + const = formulae.constants + return radius_array, RH_array, output_matrix, const + + class TestNPYComparison: @staticmethod def _mix(dry_prop, vap_prop, ratio): @@ -44,15 +98,11 @@ def _calculate_cloud_properties( self, planet: EarthLike, surface_RH: float, formulae_instance: Formulae ): const = formulae_instance.constants - planet.RH_zref = surface_RH - pvs_stp = formulae_instance.saturation_vapour_pressure.pvs_water(planet.T_STP) - initial_water_vapour_mixing_ratio = const.eps / ( planet.p_STP / planet.RH_zref / pvs_stp - 1 ) - R_air_mix = self._mix(const.Rd, const.Rv, initial_water_vapour_mixing_ratio) cp_mix = self._mix(const.c_pd, const.c_pv, initial_water_vapour_mixing_ratio) @@ -78,154 +128,113 @@ def solve_Tcloud(T_candidate): ) return initial_water_vapour_mixing_ratio, Tcloud, Zcloud, pcloud - def test_figure_2_replication_accuracy(self): - current_dir = os.path.dirname(os.path.abspath(__file__)) - groundtruth_dir = os.path.abspath(os.path.join(current_dir, "ground_truth")) - if not os.path.isdir(groundtruth_dir): - pytest.fail(f"Groundtruth directory not found at {groundtruth_dir}") - + def test_figure_2_replication_accuracy(self, ground_truth_sample, static_arrays): formulae = Formulae( ventilation="PruppacherAndRasmussen1979", saturation_vapour_pressure="AugustRocheMagnus", diffusion_coordinate="WaterMassLogarithm", ) - - with GroundTruthLoader(groundtruth_dir) as gt: - if gt.RHs is None or gt.r0grid is None or gt.m_frac_evap is None: - pytest.fail("Ground truth data (.npy files) not loaded properly.") - - n_rh_values = len(gt.RHs) - n_radius_values = gt.r0grid.shape[1] - - total_points = n_rh_values * n_radius_values - n_samples = min(gt.n_samples, total_points) - - if n_samples == 0: - pytest.skip("No data points available to sample.") - - all_indices = np.array( - [(i, j) for i in range(n_rh_values) for j in range(n_radius_values)] + for sample in ground_truth_sample: + planet = EarthLike() + rh = sample["rh"] + r_m = sample["r_m"] + expected = sample["expected_m_frac_evap"] + i_rh = sample["i_rh"] + j_r = sample["j_r"] + try: + iwvmr, Tcloud, Zcloud, pcloud = self._calculate_cloud_properties( + planet, rh, formulae + ) + simulated = self.calc_simulated_m_frac_evap_point( + planet, + formulae, + i_rh, + j_r, + rh, + r_m, + expected, + iwvmr, + Tcloud, + Zcloud, + pcloud, + ) + except Exception as e: + pytest.fail( + f"Error in _calculate_cloud_properties for RH={rh} (sample idx {i_rh},{j_r}): {e}." + ) + error_context = ( + f"Sample (RH_idx={i_rh}, R_idx={j_r}), RH={rh:.4f}, R_m={r_m:.3e}. " + f"Expected: {expected}, Got: {simulated}" ) - - sampled_indices_flat = np.random.choice(len(all_indices), n_samples, replace=False) - sampled_ij_pairs = all_indices[sampled_indices_flat] - - for i_rh, j_r in sampled_ij_pairs: - current_planet_state = EarthLike() - - current_rh = gt.RHs[i_rh] - current_r_m = gt.r0grid[0, j_r] - expected_m_frac_evap = gt.m_frac_evap[i_rh, j_r] - - try: - iwvmr, Tcloud, Zcloud, pcloud = self._calculate_cloud_properties( - current_planet_state, current_rh, formulae - ) - simulated_m_frac_evap_point = self.calc_simulated_m_frac_evap_point( - current_planet_state, - formulae, - i_rh, - j_r, - current_rh, - current_r_m, - expected_m_frac_evap, - iwvmr, - Tcloud, - Zcloud, - pcloud, - ) - except Exception as e: - print( - f"Warning: Error in _calculate_cloud_properties for RH={current_rh} (sample idx {i_rh},{j_r}): {e}." - ) - - error_context = ( - f"Sample (RH_idx={i_rh}, R_idx={j_r}), " - f"RH={current_rh:.4f}, R_m={current_r_m:.3e}. " - f"Expected: {expected_m_frac_evap}, Got: {simulated_m_frac_evap_point}" + if np.isnan(expected): + assert np.isnan( + simulated + ), f"NaN Mismatch. {error_context} (Expected NaN, got non-NaN)" + else: + assert not np.isnan( + simulated + ), f"NaN Mismatch. {error_context} (Expected non-NaN, got NaN)" + np.testing.assert_allclose( + simulated, + expected, + rtol=1e-1, + atol=1e-1, + err_msg=f"Value Mismatch. {error_context}", ) - if np.isnan(expected_m_frac_evap): - assert np.isnan( - simulated_m_frac_evap_point - ), f"NaN Mismatch. {error_context} (Expected NaN, got non-NaN)" - else: - assert not np.isnan( - simulated_m_frac_evap_point - ), f"NaN Mismatch. {error_context} (Expected non-NaN, got NaN)" - np.testing.assert_allclose( - simulated_m_frac_evap_point, - expected_m_frac_evap, - rtol=1e-1, # Relative tolerance - atol=1e-1, # Absolute tolerance - err_msg=f"Value Mismatch. {error_context}", - ) - def calc_simulated_m_frac_evap_point( self, - current_planet_state, + planet, formulae, i_rh, j_r, - current_rh, - current_r_m, - expected_m_frac_evap, + rh, + r_m, + expected, iwvmr, Tcloud, Zcloud, pcloud, ): - - simulated_m_frac_evap_point = np.nan - - if np.isnan(current_r_m) or current_r_m <= 0: - print(f"Warning: Invalid radius current_r_m={current_r_m} for sample idx {i_rh},{j_r}.") - else: - settings = Settings( - planet=current_planet_state, - r_wet=current_r_m, - mass_of_dry_air=1e5 * si.kg, - initial_water_vapour_mixing_ratio=iwvmr, - pcloud=pcloud, - Zcloud=Zcloud, - Tcloud=Tcloud, - formulae=formulae, - ) - simulation = Simulation(settings) - - try: - output = simulation.run() - if ( - output - and "r" in output - and len(output["r"]) > 0 - and "z" in output - and len(output["z"]) > 0 - ): - final_radius_um = output["r"][-1] - if np.isnan(final_radius_um) or final_radius_um < 0: - final_radius_m = final_radius_um * 1e-6 - if final_radius_m < 0: # Non-physical radius - simulated_m_frac_evap_point = 1.0 # 1.0 means fully evaporated - else: - simulated_m_frac_evap_point = np.nan - else: - final_radius_m = final_radius_um * 1e-6 - if current_r_m == 0: - frac_evap = 1.0 - else: - frac_evap = 1.0 - (final_radius_m / current_r_m) ** 3 - frac_evap = np.clip(frac_evap, 0.0, 1.0) - simulated_m_frac_evap_point = frac_evap - else: - simulated_m_frac_evap_point = np.nan - except Exception as e: - print( - f"Warning: Simulation run failed for RH={current_rh}, r={current_r_m} (sample idx {i_rh},{j_r}): {e}." - ) - if np.isclose(expected_m_frac_evap, 1.0, atol=1e-6): - simulated_m_frac_evap_point = 1.0 + if np.isnan(r_m) or r_m <= 0: + pytest.fail(f"Invalid radius r_m={r_m} for sample idx {i_rh},{j_r}.") + settings = Settings( + planet=planet, + r_wet=r_m, + mass_of_dry_air=1e5 * si.kg, + initial_water_vapour_mixing_ratio=iwvmr, + pcloud=pcloud, + Zcloud=Zcloud, + Tcloud=Tcloud, + formulae=formulae, + ) + simulation = Simulation(settings) + try: + output = simulation.run() + if ( + output + and "r" in output + and len(output["r"]) > 0 + and "z" in output + and len(output["z"]) > 0 + ): + final_radius_um = output["r"][-1] + if np.isnan(final_radius_um) or final_radius_um < 0: + final_radius_m = final_radius_um * 1e-6 + if final_radius_m < 0: + return 1.0 + return np.nan + final_radius_m = final_radius_um * 1e-6 + if r_m == 0: + frac_evap = 1.0 else: - simulated_m_frac_evap_point = np.nan - - return simulated_m_frac_evap_point + frac_evap = 1.0 - (final_radius_m / r_m) ** 3 + return np.clip(frac_evap, 0.0, 1.0) + return np.nan + except Exception as e: + warnings.warn( + f"Simulation run failed for RH={rh:.4f}, r={r_m:.3e} (sample idx {i_rh},{j_r}): {type(e).__name__}: {e}" + ) + if np.isclose(expected, 1.0, atol=1e-6): + return 1.0 + return np.nan From f55d91ae0e97238920ab402061a9f6e9e2c1466c Mon Sep 17 00:00:00 2001 From: Mateusz Mazur Date: Thu, 19 Jun 2025 22:23:36 +0200 Subject: [PATCH 03/15] refactor: streamlined dz_dt computation in Parcel class and updated planetary pressure units --- PySDM/environments/parcel.py | 5 ++- .../Loftus_and_Wordsworth_2021/parcel.py | 37 +------------------ .../Loftus_and_Wordsworth_2021/planet.py | 21 +++++++---- 3 files changed, 20 insertions(+), 43 deletions(-) diff --git a/PySDM/environments/parcel.py b/PySDM/environments/parcel.py index 59e65c54bc..a6adc96e04 100644 --- a/PySDM/environments/parcel.py +++ b/PySDM/environments/parcel.py @@ -106,7 +106,7 @@ def advance_parcel_vars(self): T = self["T"][0] p = self["p"][0] - dz_dt = self.w((self.particulator.n_steps + 1 / 2) * dt) # "mid-point" + dz_dt = self._compute_dz_dt(dt) water_vapour_mixing_ratio = ( self["water_vapour_mixing_ratio"][0] - self.delta_liquid_water_mixing_ratio / 2 @@ -133,6 +133,9 @@ def advance_parcel_vars(self): (self._tmp["rhod"][0] + self["rhod"][0]) / 2, self.mass_of_dry_air ) + def _compute_dz_dt(self, dt): + self.w((self.particulator.n_steps + 1 / 2) * dt) # "mid-point" + def get_thd(self): return self["thd"] diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py index 7ac5eb0dac..3223e23772 100644 --- a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py @@ -25,38 +25,5 @@ def __init__( variables=None, ) - def advance_parcel_vars(self): - """ - Compute new values of displacement, dry-air density and volume, - and write them to self._tmp and self.mesh.dv - """ - dt = self.particulator.dt - formulae = self.particulator.formulae - T = self["T"][0] - p = self["p"][0] - - dz_dt = -self.particulator.attributes["terminal velocity"].to_ndarray()[0] - water_vapour_mixing_ratio = ( - self["water_vapour_mixing_ratio"][0] - - self.delta_liquid_water_mixing_ratio / 2 - ) - - drho_dz = formulae.hydrostatics.drho_dz( - p=p, - T=T, - water_vapour_mixing_ratio=water_vapour_mixing_ratio, - lv=formulae.latent_heat_vapourisation.lv(T), - d_liquid_water_mixing_ratio__dz=( - self.delta_liquid_water_mixing_ratio / dz_dt / dt - ), - ) - drhod_dz = drho_dz - - self.particulator.backend.explicit_euler(self._tmp["z"], dt, dz_dt) - self.particulator.backend.explicit_euler( - self._tmp["rhod"], dt, dz_dt * drhod_dz - ) - - self.mesh.dv = formulae.trivia.volume_of_density_mass( - (self._tmp["rhod"][0] + self["rhod"][0]) / 2, self.mass_of_dry_air - ) + def _compute_dz_dt(self, dt): + return -self.particulator.attributes["terminal velocity"].to_ndarray()[0] \ No newline at end of file diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py index 0e4211b397..fbffa399e8 100644 --- a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py @@ -1,3 +1,10 @@ +""" +Planetary properties used in calculations. + +Values are primarily taken from Table 1 of Loftus & Wordsworth (2021), unless otherwise noted. +Each variable represents a physical property or atmospheric composition relevant for cloud microphysics modeling. +""" + from PySDM.physics.constants import si from dataclasses import dataclass from typing import Optional, Dict, Any @@ -24,7 +31,7 @@ def to_dict(self) -> Dict[str, Any]: class EarthLike(Planet): g_std: float = 9.82 * si.metre / si.second**2 T_STP: float = 300 * si.kelvin - p_STP: float = 1.01325 * 1e6 * si.pascal + p_STP: float = 1.01325 * 1e5 * si.Pa RH_zref: float = 0.75 dry_molar_conc_H2: float = 0 dry_molar_conc_He: float = 0 @@ -38,7 +45,7 @@ class EarthLike(Planet): class Earth(Planet): g_std: float = 9.82 * si.metre / si.second**2 T_STP: float = 290 * si.kelvin - p_STP: float = 1.01325 * 1e6 * si.pascal + p_STP: float = 1.01325 * 1e5 * si.Pa RH_zref: float = 0.75 dry_molar_conc_H2: float = 0 dry_molar_conc_He: float = 0 @@ -52,7 +59,7 @@ class Earth(Planet): class EarlyMars(Planet): g_std: float = 3.71 * si.metre / si.second**2 T_STP: float = 290 * si.kelvin - p_STP: float = 2 * 1e6 * si.pascal + p_STP: float = 2 * 1e5 * si.Pa RH_zref: float = 0.75 dry_molar_conc_H2: float = 0 dry_molar_conc_He: float = 0 @@ -66,7 +73,7 @@ class EarlyMars(Planet): class Jupiter(Planet): g_std: float = 24.84 * si.metre / si.second**2 T_STP: float = 274 * si.kelvin - p_STP: float = 4.85 * 1e6 * si.pascal + p_STP: float = 4.85 * 1e5 * si.Pa RH_zref: float = 1 dry_molar_conc_H2: float = 0.864 dry_molar_conc_He: float = 0.136 @@ -80,7 +87,7 @@ class Jupiter(Planet): class Saturn(Planet): g_std: float = 10.47 * si.metre / si.second**2 T_STP: float = 284 * si.kelvin - p_STP: float = 10.4 * 1e6 * si.pascal + p_STP: float = 10.4 * 1e5 * si.Pa RH_zref: float = 1 dry_molar_conc_H2: float = 0.88 dry_molar_conc_He: float = 0.12 @@ -94,7 +101,7 @@ class Saturn(Planet): class K2_18B(Planet): g_std: float = 12.44 * si.metre / si.second**2 T_STP: float = 275 * si.kelvin - p_STP: float = 0.1 * 1e6 * si.pascal + p_STP: float = 0.1 * 1e5 * si.Pa RH_zref: float = 1 dry_molar_conc_H2: float = 0.9 dry_molar_conc_He: float = 0.1 @@ -108,7 +115,7 @@ class K2_18B(Planet): class CompositeTest(Planet): g_std: float = 9.82 * si.metre / si.second**2 T_STP: float = 275 * si.kelvin - p_STP: float = 0.75 * 1e6 * si.pascal + p_STP: float = 0.75 * 1e5 * si.Pa RH_zref: float = 1 dry_molar_conc_H2: float = 0.1 dry_molar_conc_He: float = 0.1 From bb15488e5504aa40779d462029ec045417b91113 Mon Sep 17 00:00:00 2001 From: Mateusz Mazur Date: Thu, 19 Jun 2025 22:57:13 +0200 Subject: [PATCH 04/15] style: improved code formatting and readability --- .../Loftus_and_Wordsworth_2021/parcel.py | 2 +- .../ground_truth/gen_figure.py | 12 +++++++++--- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py index 3223e23772..56380c5750 100644 --- a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py @@ -26,4 +26,4 @@ def __init__( ) def _compute_dz_dt(self, dt): - return -self.particulator.attributes["terminal velocity"].to_ndarray()[0] \ No newline at end of file + return -self.particulator.attributes["terminal velocity"].to_ndarray()[0] diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py index ff8d780ac4..67626bb2e9 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py @@ -67,9 +67,15 @@ levels=[0.1], ) cb.add_lines(c_10) -axs[1, 0].clabel(c_10, c_10.levels, fmt={0.1: "10% mass evaporated"}, fontsize="smaller") -axs[1, 0].fill_betweenx(RHs, 1e-2, (r_mins - 1e-6) * 1e3, edgecolor="k", facecolor="w", hatch="//") -axs[1, 0].annotate(text="TOTAL \\n EVAPORATION", xy=(0.04, 0.55), c="k", backgroundcolor="w") +axs[1, 0].clabel( + c_10, c_10.levels, fmt={0.1: "10% mass evaporated"}, fontsize="smaller" +) +axs[1, 0].fill_betweenx( + RHs, 1e-2, (r_mins - 1e-6) * 1e3, edgecolor="k", facecolor="w", hatch="//" +) +axs[1, 0].annotate( + text="TOTAL \\n EVAPORATION", xy=(0.04, 0.55), c="k", backgroundcolor="w" +) axs[1, 0].annotate( text=r"$r_\mathrm{min}$", xy=(0.35, 0.27), From 3f4f9a873322eed60957b5a8ccfa409ac2df3e92 Mon Sep 17 00:00:00 2001 From: lursz Date: Fri, 20 Jun 2025 09:46:40 +0200 Subject: [PATCH 05/15] =?UTF-8?q?=F0=9F=94=A7=20PR=20fixes?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Introduced const for MIN DROPLET RADIUS --- .../Loftus_and_Wordsworth_2021/simulation.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py index 071dde8bf2..11062d54b2 100644 --- a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py @@ -7,6 +7,8 @@ from PySDM.physics import constants as const from PySDM_examples.Loftus_and_Wordsworth_2021.parcel import AlienParcel +MIN_DROPLET_RADIUS = 1e-6 + class Simulation: def __init__(self, settings, backend=CPU): @@ -75,7 +77,10 @@ def run(self): } self.save(output) - while self.particulator.environment["z"][0] > 0 and output["r"][-1] > 1e-6: + while ( + self.particulator.environment["z"][0] > 0 + and output["r"][-1] > MIN_DROPLET_RADIUS + ): self.particulator.run(1) self.save(output) From 7354f937390f6d38c14a2d675089899634c574b4 Mon Sep 17 00:00:00 2001 From: Piotr Kubala Date: Mon, 23 Jun 2025 22:06:40 +0200 Subject: [PATCH 06/15] some linter fixes --- .../accuracy_test.py | 28 +++++++---- .../ground_truth/gen_figure.py | 4 +- .../Loftus_and_Wordsworth_2021/unit_test.py | 47 ++++++++++++------- 3 files changed, 50 insertions(+), 29 deletions(-) diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py index 6f107c8374..24e3a3c1b3 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py @@ -1,12 +1,14 @@ +# pylint: disable=missing-module-docstring from __future__ import annotations -import pytest import os -import numpy as np -from scipy.optimize import fsolve import warnings from typing import Tuple, List, Generator +import pytest +import numpy as np +from scipy.optimize import fsolve + from PySDM import Formulae from PySDM.physics import si from PySDM_examples.Loftus_and_Wordsworth_2021 import Settings, Simulation @@ -32,8 +34,9 @@ def __enter__(self): return self except FileNotFoundError as e: pytest.fail(f"Error loading ground truth files: {e}") - except Exception as e: - pytest.fail(f"Unexpected error loading ground truth data: {e}") + pytest.fail("Ground truth data not loaded successfully.") + + return self def __exit__(self, exc_type, exc_val, exc_tb): pass @@ -123,12 +126,14 @@ def solve_Tcloud(T_candidate): th_std = formulae_instance.trivia.th_std(planet.p_STP, planet.T_STP) - pcloud = formulae_instance.hydrostatics.p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio( + pcloud = \ + formulae_instance.hydrostatics\ + .p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio( planet.p_STP, th_std, initial_water_vapour_mixing_ratio, Zcloud ) return initial_water_vapour_mixing_ratio, Tcloud, Zcloud, pcloud - def test_figure_2_replication_accuracy(self, ground_truth_sample, static_arrays): + def test_figure_2_replication_accuracy(self, ground_truth_sample): formulae = Formulae( ventilation="PruppacherAndRasmussen1979", saturation_vapour_pressure="AugustRocheMagnus", @@ -142,7 +147,8 @@ def test_figure_2_replication_accuracy(self, ground_truth_sample, static_arrays) i_rh = sample["i_rh"] j_r = sample["j_r"] try: - iwvmr, Tcloud, Zcloud, pcloud = self._calculate_cloud_properties( + iwvmr, Tcloud, Zcloud, pcloud = \ + self._calculate_cloud_properties( planet, rh, formulae ) simulated = self.calc_simulated_m_frac_evap_point( @@ -160,7 +166,8 @@ def test_figure_2_replication_accuracy(self, ground_truth_sample, static_arrays) ) except Exception as e: pytest.fail( - f"Error in _calculate_cloud_properties for RH={rh} (sample idx {i_rh},{j_r}): {e}." + f"Error in _calculate_cloud_properties for RH={rh} " + + f"(sample idx {i_rh},{j_r}): {e}." ) error_context = ( f"Sample (RH_idx={i_rh}, R_idx={j_r}), RH={rh:.4f}, R_m={r_m:.3e}. " @@ -233,7 +240,8 @@ def calc_simulated_m_frac_evap_point( return np.nan except Exception as e: warnings.warn( - f"Simulation run failed for RH={rh:.4f}, r={r_m:.3e} (sample idx {i_rh},{j_r}): {type(e).__name__}: {e}" + f"Simulation run failed for RH={rh:.4f}, r={r_m:.3e} " +\ + f"(sample idx {i_rh},{j_r}): {type(e).__name__}: {e}" ) if np.isclose(expected, 1.0, atol=1e-6): return 1.0 diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py index 67626bb2e9..cebc9301e4 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py @@ -67,9 +67,7 @@ levels=[0.1], ) cb.add_lines(c_10) -axs[1, 0].clabel( - c_10, c_10.levels, fmt={0.1: "10% mass evaporated"}, fontsize="smaller" -) +axs[1, 0].clabel(c_10, c_10.levels, fmt={0.1: "10% mass evaporated"}, fontsize="smaller") axs[1, 0].fill_betweenx( RHs, 1e-2, (r_mins - 1e-6) * 1e3, edgecolor="k", facecolor="w", hatch="//" ) diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py index d5b6425e39..b19a0b912b 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py @@ -1,3 +1,5 @@ +# pylint: disable=missing-module-docstring +from functools import partial from contextlib import contextmanager from PySDM import Formulae from PySDM.physics import si @@ -54,7 +56,8 @@ def test_planet_classes(self): ) assert ( total_conc <= 1.01 - ), f"Total molar concentration {total_conc} exceeds 1.01 for {planet.__class__.__name__}" + ), \ + f"Total molar concentration {total_conc} exceeds 1.01 for {planet.__class__.__name__}" def test_water_vapour_mixing_ratio_calculation(self): """Test water vapour mixing ratio calculation.""" @@ -100,7 +103,9 @@ def test_simulation_class( Zcloud_val, Tcloud_val, ): - """Test Simulation class initialization and basic functionality with parametrized settings.""" + """ + Test Simulation class initialization and basic functionality with parametrized settings. + """ with self._get_test_resources() as (formulae, earth_like): planet = earth_like @@ -178,7 +183,7 @@ def test_simulation_run_basic( lengths = [len(output[key]) for key in output.keys()] assert all( - l == lengths[0] for l in lengths + length == lengths[0] for length in lengths ), "Not all output arrays have the same length" def test_saturation_at_cloud_base(self): @@ -189,7 +194,6 @@ def test_saturation_at_cloud_base(self): ) new_Earth = EarthLike() - new_Earth.T_STP RH_array = np.linspace(0.25, 0.99, 50) const = formulae.constants @@ -197,7 +201,16 @@ def test_saturation_at_cloud_base(self): def mix(dry, vap, ratio): return (dry + ratio * vap) / (1 + ratio) - for i, RH in enumerate(RH_array[::-1]): + def f(x, water_mixing_ratio, p_stp, t_stp, c_p, Rair): + return water_mixing_ratio / ( + water_mixing_ratio + const.eps + ) * p_stp * (x / t_stp) ** ( + c_p / Rair + ) - formulae.saturation_vapour_pressure.pvs_water( + x + ) + + for RH in RH_array[::-1]: new_Earth.RH_zref = RH pvs = formulae.saturation_vapour_pressure.pvs_water(new_Earth.T_STP) @@ -208,21 +221,23 @@ def mix(dry, vap, ratio): Rair = mix(const.Rd, const.Rv, initial_water_vapour_mixing_ratio) c_p = mix(const.c_pd, const.c_pv, initial_water_vapour_mixing_ratio) - def f(x): - return initial_water_vapour_mixing_ratio / ( - initial_water_vapour_mixing_ratio + const.eps - ) * new_Earth.p_STP * (x / new_Earth.T_STP) ** ( - c_p / Rair - ) - formulae.saturation_vapour_pressure.pvs_water( - x - ) - - tdews = fsolve(f, [150, 300]) + tdews = fsolve( + partial( + f, + water_mixing_ratio=initial_water_vapour_mixing_ratio, + p_stp=new_Earth.p_STP, + t_stp=new_Earth.T_STP, + c_p=c_p, + Rair=Rair, + ), + [150, 300] + ) Tcloud = np.max(tdews) Zcloud = (new_Earth.T_STP - Tcloud) * c_p / new_Earth.g_std thstd = formulae.trivia.th_std(new_Earth.p_STP, new_Earth.T_STP) - pcloud = formulae.hydrostatics.p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio( + pcloud = formulae.hydrostatics\ + .p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio( new_Earth.p_STP, thstd, initial_water_vapour_mixing_ratio, Zcloud ) From 1d29d0c4f7954d94269c9686b53cf610595c7e8e Mon Sep 17 00:00:00 2001 From: Piotr Kubala Date: Mon, 23 Jun 2025 22:19:23 +0200 Subject: [PATCH 07/15] accuracy test linter --- .../Loftus_and_Wordsworth_2021/accuracy_test.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py index 24e3a3c1b3..ecfe58535b 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py @@ -164,7 +164,7 @@ def test_figure_2_replication_accuracy(self, ground_truth_sample): Zcloud, pcloud, ) - except Exception as e: + except Exception as e: # pylint: disable=broad-except pytest.fail( f"Error in _calculate_cloud_properties for RH={rh} " + f"(sample idx {i_rh},{j_r}): {e}." @@ -238,7 +238,7 @@ def calc_simulated_m_frac_evap_point( frac_evap = 1.0 - (final_radius_m / r_m) ** 3 return np.clip(frac_evap, 0.0, 1.0) return np.nan - except Exception as e: + except Exception as e: # pylint: disable=broad-except warnings.warn( f"Simulation run failed for RH={rh:.4f}, r={r_m:.3e} " +\ f"(sample idx {i_rh},{j_r}): {type(e).__name__}: {e}" From 70464d613883e2197f18ac65465a01b7ba6af836 Mon Sep 17 00:00:00 2001 From: Piotr Kubala Date: Mon, 23 Jun 2025 22:59:22 +0200 Subject: [PATCH 08/15] accuracy linter fixed --- .../accuracy_test.py | 123 +++++++++--------- 1 file changed, 61 insertions(+), 62 deletions(-) diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py index ecfe58535b..1824a04002 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py @@ -42,6 +42,7 @@ def __exit__(self, exc_type, exc_val, exc_tb): pass +# pylint: disable=redefined-outer-name @pytest.fixture(scope="module") def ground_truth_data(request) -> Generator[GroundTruthLoader, None, None]: current_dir = os.path.dirname(os.path.abspath(request.fspath)) @@ -52,6 +53,7 @@ def ground_truth_data(request) -> Generator[GroundTruthLoader, None, None]: yield gt +# pylint: disable=redefined-outer-name @pytest.fixture def ground_truth_sample(ground_truth_data: GroundTruthLoader) -> List[dict]: gt = ground_truth_data @@ -78,6 +80,7 @@ def ground_truth_sample(ground_truth_data: GroundTruthLoader) -> List[dict]: ] +# pylint: disable=redefined-outer-name @pytest.fixture(scope="module") def static_arrays() -> Tuple[np.ndarray, np.ndarray, np.ndarray, object]: formulae = Formulae( @@ -141,80 +144,76 @@ def test_figure_2_replication_accuracy(self, ground_truth_sample): ) for sample in ground_truth_sample: planet = EarthLike() - rh = sample["rh"] - r_m = sample["r_m"] - expected = sample["expected_m_frac_evap"] - i_rh = sample["i_rh"] - j_r = sample["j_r"] try: iwvmr, Tcloud, Zcloud, pcloud = \ self._calculate_cloud_properties( - planet, rh, formulae + planet, sample["rh"], formulae ) - simulated = self.calc_simulated_m_frac_evap_point( - planet, - formulae, - i_rh, - j_r, - rh, - r_m, - expected, - iwvmr, - Tcloud, - Zcloud, - pcloud, + settings = Settings( + planet=planet, + r_wet=sample["r_m"], + mass_of_dry_air=1e5 * si.kg, + initial_water_vapour_mixing_ratio=iwvmr, + pcloud=pcloud, + Zcloud=Zcloud, + Tcloud=Tcloud, + formulae=formulae, ) + simulated = TestNPYComparison.calc_simulated_m_frac_evap_point( + sample["i_rh"], + sample["j_r"], + sample["rh"], + sample["expected_m_frac_evap"], + settings + ) + expected = sample["expected_m_frac_evap"] + error_context = ( + f"Sample (RH_idx={sample['i_rh']}, R_idx={sample['j_r']}), " + + f"RH={sample['rh']:.4f}, R_m={sample['r_m']:.3e}. " + f"Expected: {expected}, Got: {simulated}" + ) + if np.isnan(expected): + assert np.isnan( + simulated + ), f"NaN Mismatch. {error_context} (Expected NaN, got non-NaN)" + else: + assert not np.isnan( + simulated + ), f"NaN Mismatch. {error_context} (Expected non-NaN, got NaN)" + np.testing.assert_allclose( + simulated, + expected, + rtol=1e-1, + atol=1e-1, + err_msg=f"Value Mismatch. {error_context}", + ) except Exception as e: # pylint: disable=broad-except pytest.fail( - f"Error in _calculate_cloud_properties for RH={rh} " + - f"(sample idx {i_rh},{j_r}): {e}." - ) - error_context = ( - f"Sample (RH_idx={i_rh}, R_idx={j_r}), RH={rh:.4f}, R_m={r_m:.3e}. " - f"Expected: {expected}, Got: {simulated}" - ) - if np.isnan(expected): - assert np.isnan( - simulated - ), f"NaN Mismatch. {error_context} (Expected NaN, got non-NaN)" - else: - assert not np.isnan( - simulated - ), f"NaN Mismatch. {error_context} (Expected non-NaN, got NaN)" - np.testing.assert_allclose( - simulated, - expected, - rtol=1e-1, - atol=1e-1, - err_msg=f"Value Mismatch. {error_context}", + f"Error in _calculate_cloud_properties for RH={sample['rh']} " + + f"(sample idx {sample['i_rh']},{sample['j_r']}): {e}." ) + @staticmethod def calc_simulated_m_frac_evap_point( - self, - planet, - formulae, i_rh, j_r, rh, - r_m, expected, - iwvmr, - Tcloud, - Zcloud, - pcloud, + settings ): - if np.isnan(r_m) or r_m <= 0: - pytest.fail(f"Invalid radius r_m={r_m} for sample idx {i_rh},{j_r}.") - settings = Settings( - planet=planet, - r_wet=r_m, - mass_of_dry_air=1e5 * si.kg, - initial_water_vapour_mixing_ratio=iwvmr, - pcloud=pcloud, - Zcloud=Zcloud, - Tcloud=Tcloud, - formulae=formulae, - ) + # settings = Settings( + # planet=planet, + # r_wet=sample["r_m"], + # mass_of_dry_air=1e5 * si.kg, + # initial_water_vapour_mixing_ratio=iwvmr, + # pcloud=pcloud, + # Zcloud=Zcloud, + # Tcloud=Tcloud, + # formulae=formulae, + # ) + + if np.isnan(settings.r_wet) or settings.r_wet <= 0: + pytest.fail(f"Invalid radius r_m={settings.r_wet} for sample idx {i_rh},{j_r}.") simulation = Simulation(settings) try: output = simulation.run() @@ -232,15 +231,15 @@ def calc_simulated_m_frac_evap_point( return 1.0 return np.nan final_radius_m = final_radius_um * 1e-6 - if r_m == 0: + if settings.r_wet == 0: frac_evap = 1.0 else: - frac_evap = 1.0 - (final_radius_m / r_m) ** 3 + frac_evap = 1.0 - (final_radius_m / settings.r_wet) ** 3 return np.clip(frac_evap, 0.0, 1.0) return np.nan except Exception as e: # pylint: disable=broad-except warnings.warn( - f"Simulation run failed for RH={rh:.4f}, r={r_m:.3e} " +\ + f"Simulation run failed for RH={rh:.4f}, r={settings.r_wet:.3e} " +\ f"(sample idx {i_rh},{j_r}): {type(e).__name__}: {e}" ) if np.isclose(expected, 1.0, atol=1e-6): From 736038cc37d8f8856faba90e0090e3c4974edc9b Mon Sep 17 00:00:00 2001 From: Piotr Kubala Date: Mon, 23 Jun 2025 23:19:47 +0200 Subject: [PATCH 09/15] linter fixed --- .../accuracy_test.py | 13 +------ .../Loftus_and_Wordsworth_2021/unit_test.py | 37 ++++++++++++------- 2 files changed, 24 insertions(+), 26 deletions(-) diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py index 1824a04002..d498b264b7 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py @@ -187,7 +187,7 @@ def test_figure_2_replication_accuracy(self, ground_truth_sample): atol=1e-1, err_msg=f"Value Mismatch. {error_context}", ) - except Exception as e: # pylint: disable=broad-except + except ValueError as e: pytest.fail( f"Error in _calculate_cloud_properties for RH={sample['rh']} " + f"(sample idx {sample['i_rh']},{sample['j_r']}): {e}." @@ -201,17 +201,6 @@ def calc_simulated_m_frac_evap_point( expected, settings ): - # settings = Settings( - # planet=planet, - # r_wet=sample["r_m"], - # mass_of_dry_air=1e5 * si.kg, - # initial_water_vapour_mixing_ratio=iwvmr, - # pcloud=pcloud, - # Zcloud=Zcloud, - # Tcloud=Tcloud, - # formulae=formulae, - # ) - if np.isnan(settings.r_wet) or settings.r_wet <= 0: pytest.fail(f"Invalid radius r_m={settings.r_wet} for sample idx {i_rh},{j_r}.") simulation = Simulation(settings) diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py index b19a0b912b..8d611f96ea 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py @@ -1,4 +1,6 @@ # pylint: disable=missing-module-docstring + +from collections import namedtuple from functools import partial from contextlib import contextmanager from PySDM import Formulae @@ -23,7 +25,8 @@ class TestLoftusWordsworth2021: @contextmanager - def _get_test_resources(self): + @staticmethod + def _get_test_resources(): formulae = Formulae( ventilation="PruppacherAndRasmussen1979", saturation_vapour_pressure="AugustRocheMagnus", @@ -61,7 +64,7 @@ def test_planet_classes(self): def test_water_vapour_mixing_ratio_calculation(self): """Test water vapour mixing ratio calculation.""" - with self._get_test_resources() as (formulae, earth_like): + with TestLoftusWordsworth2021._get_test_resources() as (formulae, earth_like): const = formulae.constants planet = earth_like @@ -86,6 +89,7 @@ def test_alien_parcel_initialization(self): ) assert parcel is not None + # pylint: disable=too-many-arguments,too-many-positional-arguments @pytest.mark.parametrize( "r_wet_val, mass_of_dry_air_val, iwvmr_val, pcloud_val, Zcloud_val, Tcloud_val", [ @@ -106,7 +110,7 @@ def test_simulation_class( """ Test Simulation class initialization and basic functionality with parametrized settings. """ - with self._get_test_resources() as (formulae, earth_like): + with TestLoftusWordsworth2021._get_test_resources() as (formulae, earth_like): planet = earth_like settings = Settings( @@ -131,6 +135,7 @@ def test_simulation_class( for product in required_products: assert product in products + # pylint: disable=too-many-arguments,too-many-positional-arguments @pytest.mark.parametrize( "r_wet_val, mass_of_dry_air_val, iwvmr_val, pcloud_val, Zcloud_val, Tcloud_val", [ @@ -149,7 +154,7 @@ def test_simulation_run_basic( Tcloud_val, ): """Test basic simulation run functionality.""" - with self._get_test_resources() as (formulae, earth_like): + with TestLoftusWordsworth2021._get_test_resources() as (formulae, earth_like): planet = earth_like settings = Settings( @@ -201,11 +206,11 @@ def test_saturation_at_cloud_base(self): def mix(dry, vap, ratio): return (dry + ratio * vap) / (1 + ratio) - def f(x, water_mixing_ratio, p_stp, t_stp, c_p, Rair): + def f(x, water_mixing_ratio, params): return water_mixing_ratio / ( water_mixing_ratio + const.eps - ) * p_stp * (x / t_stp) ** ( - c_p / Rair + ) * params.p_stp * (x / params.t_stp) ** ( + params.c_p / params.Rair ) - formulae.saturation_vapour_pressure.pvs_water( x ) @@ -213,22 +218,26 @@ def f(x, water_mixing_ratio, p_stp, t_stp, c_p, Rair): for RH in RH_array[::-1]: new_Earth.RH_zref = RH - pvs = formulae.saturation_vapour_pressure.pvs_water(new_Earth.T_STP) initial_water_vapour_mixing_ratio = const.eps / ( - new_Earth.p_STP / new_Earth.RH_zref / pvs - 1 + new_Earth.p_STP / new_Earth.RH_zref / \ + formulae.saturation_vapour_pressure.pvs_water(new_Earth.T_STP) - 1 ) - Rair = mix(const.Rd, const.Rv, initial_water_vapour_mixing_ratio) c_p = mix(const.c_pd, const.c_pv, initial_water_vapour_mixing_ratio) tdews = fsolve( partial( f, water_mixing_ratio=initial_water_vapour_mixing_ratio, - p_stp=new_Earth.p_STP, - t_stp=new_Earth.T_STP, - c_p=c_p, - Rair=Rair, + params=namedtuple( + "params", + ["p_stp", "t_stp", "c_p", "Rair"], + )( + p_stp=new_Earth.p_STP, + t_stp=new_Earth.T_STP, + c_p=c_p, + Rair=mix(const.Rd, const.Rv, initial_water_vapour_mixing_ratio), + ) ), [150, 300] ) From 217919ede4e1bbdf0abf523c4ea37d33f8d901c9 Mon Sep 17 00:00:00 2001 From: Piotr Kubala Date: Mon, 23 Jun 2025 23:23:59 +0200 Subject: [PATCH 10/15] pre-commit run --- .../accuracy_test.py | 35 ++++++++----------- .../ground_truth/gen_figure.py | 4 ++- .../Loftus_and_Wordsworth_2021/unit_test.py | 16 ++++----- 3 files changed, 25 insertions(+), 30 deletions(-) diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py index d498b264b7..fc0cb54817 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py @@ -129,9 +129,7 @@ def solve_Tcloud(T_candidate): th_std = formulae_instance.trivia.th_std(planet.p_STP, planet.T_STP) - pcloud = \ - formulae_instance.hydrostatics\ - .p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio( + pcloud = formulae_instance.hydrostatics.p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio( planet.p_STP, th_std, initial_water_vapour_mixing_ratio, Zcloud ) return initial_water_vapour_mixing_ratio, Tcloud, Zcloud, pcloud @@ -145,8 +143,7 @@ def test_figure_2_replication_accuracy(self, ground_truth_sample): for sample in ground_truth_sample: planet = EarthLike() try: - iwvmr, Tcloud, Zcloud, pcloud = \ - self._calculate_cloud_properties( + iwvmr, Tcloud, Zcloud, pcloud = self._calculate_cloud_properties( planet, sample["rh"], formulae ) settings = Settings( @@ -164,12 +161,12 @@ def test_figure_2_replication_accuracy(self, ground_truth_sample): sample["j_r"], sample["rh"], sample["expected_m_frac_evap"], - settings + settings, ) expected = sample["expected_m_frac_evap"] error_context = ( - f"Sample (RH_idx={sample['i_rh']}, R_idx={sample['j_r']}), " + - f"RH={sample['rh']:.4f}, R_m={sample['r_m']:.3e}. " + f"Sample (RH_idx={sample['i_rh']}, R_idx={sample['j_r']}), " + + f"RH={sample['rh']:.4f}, R_m={sample['r_m']:.3e}. " f"Expected: {expected}, Got: {simulated}" ) if np.isnan(expected): @@ -189,20 +186,16 @@ def test_figure_2_replication_accuracy(self, ground_truth_sample): ) except ValueError as e: pytest.fail( - f"Error in _calculate_cloud_properties for RH={sample['rh']} " + - f"(sample idx {sample['i_rh']},{sample['j_r']}): {e}." + f"Error in _calculate_cloud_properties for RH={sample['rh']} " + + f"(sample idx {sample['i_rh']},{sample['j_r']}): {e}." ) @staticmethod - def calc_simulated_m_frac_evap_point( - i_rh, - j_r, - rh, - expected, - settings - ): + def calc_simulated_m_frac_evap_point(i_rh, j_r, rh, expected, settings): if np.isnan(settings.r_wet) or settings.r_wet <= 0: - pytest.fail(f"Invalid radius r_m={settings.r_wet} for sample idx {i_rh},{j_r}.") + pytest.fail( + f"Invalid radius r_m={settings.r_wet} for sample idx {i_rh},{j_r}." + ) simulation = Simulation(settings) try: output = simulation.run() @@ -226,10 +219,10 @@ def calc_simulated_m_frac_evap_point( frac_evap = 1.0 - (final_radius_m / settings.r_wet) ** 3 return np.clip(frac_evap, 0.0, 1.0) return np.nan - except Exception as e: # pylint: disable=broad-except + except Exception as e: # pylint: disable=broad-except warnings.warn( - f"Simulation run failed for RH={rh:.4f}, r={settings.r_wet:.3e} " +\ - f"(sample idx {i_rh},{j_r}): {type(e).__name__}: {e}" + f"Simulation run failed for RH={rh:.4f}, r={settings.r_wet:.3e} " + + f"(sample idx {i_rh},{j_r}): {type(e).__name__}: {e}" ) if np.isclose(expected, 1.0, atol=1e-6): return 1.0 diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py index cebc9301e4..67626bb2e9 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py @@ -67,7 +67,9 @@ levels=[0.1], ) cb.add_lines(c_10) -axs[1, 0].clabel(c_10, c_10.levels, fmt={0.1: "10% mass evaporated"}, fontsize="smaller") +axs[1, 0].clabel( + c_10, c_10.levels, fmt={0.1: "10% mass evaporated"}, fontsize="smaller" +) axs[1, 0].fill_betweenx( RHs, 1e-2, (r_mins - 1e-6) * 1e3, edgecolor="k", facecolor="w", hatch="//" ) diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py index 8d611f96ea..911122af28 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py @@ -59,8 +59,7 @@ def test_planet_classes(self): ) assert ( total_conc <= 1.01 - ), \ - f"Total molar concentration {total_conc} exceeds 1.01 for {planet.__class__.__name__}" + ), f"Total molar concentration {total_conc} exceeds 1.01 for {planet.__class__.__name__}" def test_water_vapour_mixing_ratio_calculation(self): """Test water vapour mixing ratio calculation.""" @@ -219,8 +218,10 @@ def f(x, water_mixing_ratio, params): new_Earth.RH_zref = RH initial_water_vapour_mixing_ratio = const.eps / ( - new_Earth.p_STP / new_Earth.RH_zref / \ - formulae.saturation_vapour_pressure.pvs_water(new_Earth.T_STP) - 1 + new_Earth.p_STP + / new_Earth.RH_zref + / formulae.saturation_vapour_pressure.pvs_water(new_Earth.T_STP) + - 1 ) c_p = mix(const.c_pd, const.c_pv, initial_water_vapour_mixing_ratio) @@ -237,16 +238,15 @@ def f(x, water_mixing_ratio, params): t_stp=new_Earth.T_STP, c_p=c_p, Rair=mix(const.Rd, const.Rv, initial_water_vapour_mixing_ratio), - ) + ), ), - [150, 300] + [150, 300], ) Tcloud = np.max(tdews) Zcloud = (new_Earth.T_STP - Tcloud) * c_p / new_Earth.g_std thstd = formulae.trivia.th_std(new_Earth.p_STP, new_Earth.T_STP) - pcloud = formulae.hydrostatics\ - .p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio( + pcloud = formulae.hydrostatics.p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio( new_Earth.p_STP, thstd, initial_water_vapour_mixing_ratio, Zcloud ) From b8f819d996314304b5faea6c7478f84bbfa43bc1 Mon Sep 17 00:00:00 2001 From: Piotr Kubala Date: Tue, 24 Jun 2025 01:48:07 +0200 Subject: [PATCH 11/15] linter other fixes --- docs/bibliography.json | 7 +++++ .../accuracy_test.py | 11 +++---- .../ground_truth/gen_figure.py | 4 ++- .../Loftus_and_Wordsworth_2021/unit_test.py | 29 ++++++++++++------- 4 files changed, 34 insertions(+), 17 deletions(-) diff --git a/docs/bibliography.json b/docs/bibliography.json index 44b4be2186..baf49dae51 100644 --- a/docs/bibliography.json +++ b/docs/bibliography.json @@ -925,5 +925,12 @@ ], "title": "A physically constrained classical description of the homogeneous nucleation of ice in water", "label": "Koop and Murray 2016 (J. Chem. Phys. 145)" + }, + "https://doi.org/10.1029/2020JE006653": { + "usages": [ + "PySDM/examples/PySDM_examples/Loftus_and_Wordsworth_2021/__init__.py" + ], + "title": "The Physics of Falling Raindrops in Diverse Planetary Atmospheres", + "label": "Kaitlyn Loftus, Robin D. Wordsworth (JGR Planets 2021)" } } diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py index fc0cb54817..f23f2fc4f5 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py @@ -9,11 +9,12 @@ import numpy as np from scipy.optimize import fsolve -from PySDM import Formulae -from PySDM.physics import si from PySDM_examples.Loftus_and_Wordsworth_2021 import Settings, Simulation from PySDM_examples.Loftus_and_Wordsworth_2021.planet import EarthLike +from PySDM import Formulae +from PySDM.physics import si + class GroundTruthLoader: def __init__( @@ -122,14 +123,14 @@ def solve_Tcloud(T_candidate): pvs_tc = formulae_instance.saturation_vapour_pressure.pvs_water(T_candidate) return pv_ad - pvs_tc - Tcloud_solutions = fsolve(solve_Tcloud, [150.0, 300.0]) - Tcloud = np.max(Tcloud_solutions) + Tcloud = np.max(fsolve(solve_Tcloud, [150.0, 300.0])) Zcloud = (planet.T_STP - Tcloud) * cp_mix / planet.g_std th_std = formulae_instance.trivia.th_std(planet.p_STP, planet.T_STP) - pcloud = formulae_instance.hydrostatics.p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio( + hydro = formulae_instance.hydrostatics + pcloud = hydro.p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio( planet.p_STP, th_std, initial_water_vapour_mixing_ratio, Zcloud ) return initial_water_vapour_mixing_ratio, Tcloud, Zcloud, pcloud diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py index 67626bb2e9..fd63165b9d 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py @@ -1,10 +1,12 @@ +# pylint: disable=missing-module-docstring + ################################################################ # make LoWo21 Figure 2 # r_min, fraction raindrop mass evaporated as functions of RH ################################################################ +import os import numpy as np import matplotlib.pyplot as plt -import os # load results root_dir = os.path.dirname(os.path.abspath(__file__)) diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py index 911122af28..a15d1374b5 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py @@ -3,8 +3,6 @@ from collections import namedtuple from functools import partial from contextlib import contextmanager -from PySDM import Formulae -from PySDM.physics import si import pytest import numpy as np from scipy.optimize import fsolve @@ -21,6 +19,9 @@ from PySDM_examples.Loftus_and_Wordsworth_2021.parcel import AlienParcel from PySDM_examples.Loftus_and_Wordsworth_2021 import Settings +from PySDM import Formulae +from PySDM.physics import si + class TestLoftusWordsworth2021: @@ -57,9 +58,10 @@ def test_planet_classes(self): + planet.dry_molar_conc_O2 + planet.dry_molar_conc_CO2 ) - assert ( - total_conc <= 1.01 - ), f"Total molar concentration {total_conc} exceeds 1.01 for {planet.__class__.__name__}" + assert total_conc <= 1.01, ( + f"Total molar concentration {total_conc} " + + f"exceeds 1.01 for {planet.__class__.__name__}" + ) def test_water_vapour_mixing_ratio_calculation(self): """Test water vapour mixing ratio calculation.""" @@ -88,7 +90,7 @@ def test_alien_parcel_initialization(self): ) assert parcel is not None - # pylint: disable=too-many-arguments,too-many-positional-arguments + # pylint: disable=too-many-arguments @pytest.mark.parametrize( "r_wet_val, mass_of_dry_air_val, iwvmr_val, pcloud_val, Zcloud_val, Tcloud_val", [ @@ -134,7 +136,7 @@ def test_simulation_class( for product in required_products: assert product in products - # pylint: disable=too-many-arguments,too-many-positional-arguments + # pylint: disable=too-many-arguments @pytest.mark.parametrize( "r_wet_val, mass_of_dry_air_val, iwvmr_val, pcloud_val, Zcloud_val, Tcloud_val", [ @@ -185,7 +187,7 @@ def test_simulation_run_basic( assert len(output["z"]) > 0, "Output array 'z' is empty" assert len(output["t"]) > 0, "Output array 't' is empty" - lengths = [len(output[key]) for key in output.keys()] + lengths = [len(one_output) for one_output in output.values()] assert all( length == lengths[0] for length in lengths ), "Not all output arrays have the same length" @@ -243,11 +245,16 @@ def f(x, water_mixing_ratio, params): [150, 300], ) Tcloud = np.max(tdews) - Zcloud = (new_Earth.T_STP - Tcloud) * c_p / new_Earth.g_std thstd = formulae.trivia.th_std(new_Earth.p_STP, new_Earth.T_STP) - pcloud = formulae.hydrostatics.p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio( - new_Earth.p_STP, thstd, initial_water_vapour_mixing_ratio, Zcloud + hydro = formulae.hydrostatics + pcloud = ( + hydro.p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio( + new_Earth.p_STP, + thstd, + initial_water_vapour_mixing_ratio, + (new_Earth.T_STP - Tcloud) * c_p / new_Earth.g_std, + ) ) np.testing.assert_approx_equal( From e0874d5680186ba06eb28647f69c7ac82b94b41d Mon Sep 17 00:00:00 2001 From: Piotr Kubala Date: Tue, 24 Jun 2025 01:58:07 +0200 Subject: [PATCH 12/15] bibliography fix --- docs/bibliography.json | 2 +- .../Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/docs/bibliography.json b/docs/bibliography.json index baf49dae51..8d15c97406 100644 --- a/docs/bibliography.json +++ b/docs/bibliography.json @@ -928,7 +928,7 @@ }, "https://doi.org/10.1029/2020JE006653": { "usages": [ - "PySDM/examples/PySDM_examples/Loftus_and_Wordsworth_2021/__init__.py" + "examples/PySDM_examples/Loftus_and_Wordsworth_2021/__init__.py" ], "title": "The Physics of Falling Raindrops in Diverse Planetary Atmospheres", "label": "Kaitlyn Loftus, Robin D. Wordsworth (JGR Planets 2021)" diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py index fd63165b9d..89fec0af21 100644 --- a/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py +++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py @@ -43,7 +43,7 @@ r0grid * 1e3, RHgrid, m_frac_evap, - cmap=plt.cm.binary, + cmap=plt.cm.binary, # pylint: disable=no-member vmin=0, vmax=1, levels=levels_smooth, @@ -95,7 +95,10 @@ axs[0, 0].annotate(text=r"$r_\mathrm{min}$", xy=(0.13, 0.05), c="darkviolet", size=8) figs_path = os.path.join(dir, "figs") -os.mkdir(figs_path) if not os.path.exists(figs_path) else None + +if not os.path.exists(figs_path): + os.mkdir(figs_path) + plt.savefig( os.path.join(figs_path, "fig02.pdf"), transparent=True, From 4d0299aa4edbf711cbf28c6d29596f720da0cf60 Mon Sep 17 00:00:00 2001 From: Piotr Kubala Date: Tue, 24 Jun 2025 02:11:27 +0200 Subject: [PATCH 13/15] parcel small fix --- PySDM/environments/parcel.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PySDM/environments/parcel.py b/PySDM/environments/parcel.py index a6adc96e04..fdff10de91 100644 --- a/PySDM/environments/parcel.py +++ b/PySDM/environments/parcel.py @@ -134,7 +134,7 @@ def advance_parcel_vars(self): ) def _compute_dz_dt(self, dt): - self.w((self.particulator.n_steps + 1 / 2) * dt) # "mid-point" + return self.w((self.particulator.n_steps + 1 / 2) * dt) # "mid-point" def get_thd(self): return self["thd"] From 7f114a7008e199eed368db4f3c1e2e8901deccad Mon Sep 17 00:00:00 2001 From: Piotr Kubala Date: Tue, 24 Jun 2025 02:29:15 +0200 Subject: [PATCH 14/15] linter fixes --- .../Loftus_and_Wordsworth_2021/__init__.py | 1 + .../Loftus_and_Wordsworth_2021/figure_2.ipynb | 4 +--- .../PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py | 5 +++++ .../PySDM_examples/Loftus_and_Wordsworth_2021/planet.py | 9 ++++++--- .../Loftus_and_Wordsworth_2021/settings.py | 8 ++++++-- .../Loftus_and_Wordsworth_2021/simulation.py | 7 ++++++- 6 files changed, 25 insertions(+), 9 deletions(-) diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/__init__.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/__init__.py index 2d0f59fb30..2b6bbaec44 100644 --- a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/__init__.py +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/__init__.py @@ -1,3 +1,4 @@ +# pylint: disable=line-too-long """ The Physics of Falling Raindrops in Diverse Planetary Atmospheres [Loftus, K., & Wordsworth, R. D. 2021 (Journal of Geophysical Research: Planets, 126)](https://doi.org/10.1029/2020JE006653) diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb index 6a27e3ffa6..3630950b03 100644 --- a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb @@ -114,9 +114,7 @@ " initial_water_vapour_mixing_ratio + const.eps\n", " ) * new_Earth.p_STP * (x / new_Earth.T_STP) ** (\n", " c_p / Rair\n", - " ) - formulae.saturation_vapour_pressure.pvs_water(\n", - " x\n", - " )\n", + " ) - formulae.saturation_vapour_pressure.pvs_water(x)\n", "\n", " tdews = fsolve(f, [150, 300])\n", " Tcloud = np.max(tdews)\n", diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py index 56380c5750..6b1b2e30fb 100644 --- a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py @@ -1,3 +1,7 @@ +""" +Modified Parcel class for the Loftus and Wordsworth 2021 example. +""" + from PySDM.environments.parcel import Parcel @@ -25,5 +29,6 @@ def __init__( variables=None, ) + # pylint: disable=unused-argument def _compute_dz_dt(self, dt): return -self.particulator.attributes["terminal velocity"].to_ndarray()[0] diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py index fbffa399e8..d360626619 100644 --- a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py @@ -1,14 +1,17 @@ """ Planetary properties used in calculations. -Values are primarily taken from Table 1 of Loftus & Wordsworth (2021), unless otherwise noted. -Each variable represents a physical property or atmospheric composition relevant for cloud microphysics modeling. +Values are primarily taken from Table 1 of Loftus & Wordsworth (2021), +unless otherwise noted. +Each variable represents a physical property or atmospheric +composition relevant for cloud microphysics modeling. """ -from PySDM.physics.constants import si from dataclasses import dataclass from typing import Optional, Dict, Any +from PySDM.physics.constants import si + @dataclass class Planet: diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/settings.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/settings.py index b0a81752a7..8f5d6f95ad 100644 --- a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/settings.py +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/settings.py @@ -1,15 +1,19 @@ -# Planetary Properties, Loftus and Wordsworth 2021 Table 1 +""" +Planetary Properties, Loftus and Wordsworth 2021 Table 1 +""" from pystrict import strict +from PySDM_examples.Loftus_and_Wordsworth_2021.planet import Planet + from PySDM import Formulae from PySDM.dynamics import condensation from PySDM.physics.constants import si -from PySDM_examples.Loftus_and_Wordsworth_2021.planet import Planet @strict class Settings: + # pylint: disable=too-many-arguments def __init__( self, r_wet: float, diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py index 11062d54b2..a182885e99 100644 --- a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py @@ -1,11 +1,16 @@ +""" +Simulation for the Loftus and Wordsworth (2021) example. +""" + import numpy as np +from PySDM_examples.Loftus_and_Wordsworth_2021.parcel import AlienParcel + import PySDM.products as PySDM_products from PySDM.backends import CPU from PySDM.builder import Builder from PySDM.dynamics import AmbientThermodynamics, Condensation from PySDM.physics import constants as const -from PySDM_examples.Loftus_and_Wordsworth_2021.parcel import AlienParcel MIN_DROPLET_RADIUS = 1e-6 From dba682dc335fcc133c2bfd90e98e596b9e6e3c9a Mon Sep 17 00:00:00 2001 From: Piotr Kubala Date: Tue, 24 Jun 2025 02:45:55 +0200 Subject: [PATCH 15/15] last linter --- .../Loftus_and_Wordsworth_2021/figure_2.ipynb | 38 ++++++++----------- 1 file changed, 15 insertions(+), 23 deletions(-) diff --git a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb index 3630950b03..7fa0adb654 100644 --- a/examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb +++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb @@ -12,7 +12,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "id": "d8f644e9", "metadata": {}, "outputs": [], @@ -28,15 +28,7 @@ "from PySDM.physics import si\n", "from scipy.optimize import fsolve\n", "from PySDM_examples.Loftus_and_Wordsworth_2021 import Settings\n", - "from PySDM_examples.Loftus_and_Wordsworth_2021.planet import (\n", - " Planet,\n", - " EarthLike,\n", - " Earth,\n", - " EarlyMars,\n", - " Jupiter,\n", - " Saturn,\n", - " K2_18B,\n", - ")\n", + "from PySDM_examples.Loftus_and_Wordsworth_2021.planet import EarthLike\n", "from PySDM_examples.Loftus_and_Wordsworth_2021 import Simulation" ] }, @@ -56,7 +48,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "id": "41f0ed6d", "metadata": {}, "outputs": [ @@ -73,12 +65,12 @@ ], "source": [ "new_Earth = EarthLike()\n", - "new_Earth.T_STP" + "print(new_Earth.T_STP)" ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "id": "d3a8f4b9", "metadata": {}, "outputs": [], @@ -94,12 +86,13 @@ " return (dry + ratio * vap) / (1 + ratio)\n", "\n", "\n", - "def compute_one_RH(i, RH):\n", + "# pylint: disable=redefined-outer-name\n", + "def compute_one_RH(index, RH_value):\n", " \"\"\"\n", " Compute one row of the output_matrix for a given RH.\n", " Returns a 1D numpy array of length len(radius_array).\n", " \"\"\"\n", - " new_Earth.RH_zref = RH\n", + " new_Earth.RH_zref = RH_value\n", "\n", " pvs = formulae.saturation_vapour_pressure.pvs_water(new_Earth.T_STP)\n", " initial_water_vapour_mixing_ratio = const.eps / (\n", @@ -155,28 +148,27 @@ " if output[\"z\"][-1] > 0:\n", " row_data[j] = np.nan\n", " break\n", - " else:\n", - " row_data[j] = 1 - (output[\"r\"][-1] / (r * 1e6))\n", - " except Exception as _:\n", + " row_data[j] = 1 - (output[\"r\"][-1] / (r * 1e6))\n", + " except Exception as _: # pylint: disable=broad-exception-caught\n", " break\n", "\n", - " return i, row_data, output" + " return index, row_data, output" ] }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "id": "4352de81", "metadata": {}, "outputs": [], "source": [ "all_rows = Parallel(n_jobs=os.cpu_count())(\n", - " delayed(compute_one_RH)(i, RH) for i, RH in enumerate(RH_array[::-1])\n", + " delayed(compute_one_RH)(index, RH_value) for index, RH_value in enumerate(RH_array[::-1])\n", ")\n", "\n", "last_output = None\n", - "for i, row_data, output in all_rows:\n", - " output_matrix[i] = row_data\n", + "for index, row_data, output in all_rows:\n", + " output_matrix[index] = row_data\n", " last_output = output" ] },