diff --git a/PySDM/environments/parcel.py b/PySDM/environments/parcel.py
index 59e65c54bc..fdff10de91 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):
+ return self.w((self.particulator.n_steps + 1 / 2) * dt) # "mid-point"
+
def get_thd(self):
return self["thd"]
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/docs/bibliography.json b/docs/bibliography.json
index 44b4be2186..8d15c97406 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": [
+ "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/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..2b6bbaec44
--- /dev/null
+++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/__init__.py
@@ -0,0 +1,9 @@
+# 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)
+
+"""
+
+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..7fa0adb654
--- /dev/null
+++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb
@@ -0,0 +1,14713 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1e8d983b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "!pip install joblib"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "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 EarthLike\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": null,
+ "id": "41f0ed6d",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "300.0"
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "new_Earth = EarthLike()\n",
+ "print(new_Earth.T_STP)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "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",
+ "# 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_value\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(x)\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",
+ " row_data[j] = 1 - (output[\"r\"][-1] / (r * 1e6))\n",
+ " except Exception as _: # pylint: disable=broad-exception-caught\n",
+ " break\n",
+ "\n",
+ " return index, row_data, output"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4352de81",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "all_rows = Parallel(n_jobs=os.cpu_count())(\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 index, row_data, output in all_rows:\n",
+ " output_matrix[index] = row_data\n",
+ " last_output = output"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "8e6027d8",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\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"
+ ],
+ "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..6b1b2e30fb
--- /dev/null
+++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py
@@ -0,0 +1,34 @@
+"""
+Modified Parcel class for the Loftus and Wordsworth 2021 example.
+"""
+
+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,
+ )
+
+ # 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
new file mode 100644
index 0000000000..d360626619
--- /dev/null
+++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py
@@ -0,0 +1,128 @@
+"""
+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 dataclasses import dataclass
+from typing import Optional, Dict, Any
+
+from PySDM.physics.constants import si
+
+
+@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 * 1e5 * si.Pa
+ 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 * 1e5 * si.Pa
+ 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 * 1e5 * si.Pa
+ 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 * 1e5 * si.Pa
+ 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 * 1e5 * si.Pa
+ 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 * 1e5 * si.Pa
+ 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 * 1e5 * si.Pa
+ 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..8f5d6f95ad
--- /dev/null
+++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/settings.py
@@ -0,0 +1,50 @@
+"""
+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
+
+
+@strict
+class Settings:
+ # pylint: disable=too-many-arguments
+ 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..a182885e99
--- /dev/null
+++ b/examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py
@@ -0,0 +1,92 @@
+"""
+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
+
+MIN_DROPLET_RADIUS = 1e-6
+
+
+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] > MIN_DROPLET_RADIUS
+ ):
+ 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..f23f2fc4f5
--- /dev/null
+++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py
@@ -0,0 +1,230 @@
+# pylint: disable=missing-module-docstring
+from __future__ import annotations
+
+import os
+import warnings
+from typing import Tuple, List, Generator
+
+import pytest
+import numpy as np
+from scipy.optimize import fsolve
+
+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__(
+ 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
+ np.random.seed(random_seed)
+
+ 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:
+ pytest.fail(f"Error loading ground truth files: {e}")
+ pytest.fail("Ground truth data not loaded successfully.")
+
+ return self
+
+ 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))
+ 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
+
+
+# pylint: disable=redefined-outer-name
+@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
+ ]
+
+
+# pylint: disable=redefined-outer-name
+@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):
+ 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 = 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)
+
+ 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
+
+ def test_figure_2_replication_accuracy(self, ground_truth_sample):
+ formulae = Formulae(
+ ventilation="PruppacherAndRasmussen1979",
+ saturation_vapour_pressure="AugustRocheMagnus",
+ diffusion_coordinate="WaterMassLogarithm",
+ )
+ for sample in ground_truth_sample:
+ planet = EarthLike()
+ try:
+ iwvmr, Tcloud, Zcloud, pcloud = self._calculate_cloud_properties(
+ planet, sample["rh"], 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,
+ )
+ 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 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}."
+ )
+
+ @staticmethod
+ 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}."
+ )
+ 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 settings.r_wet == 0:
+ frac_evap = 1.0
+ else:
+ 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={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
+ return np.nan
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 0000000000..e02761b427
Binary files /dev/null and b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/RHgrid.npy differ
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 0000000000..73132c6203
Binary files /dev/null and b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/RHs.npy differ
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..89fec0af21
--- /dev/null
+++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py
@@ -0,0 +1,108 @@
+# 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
+
+# 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, # pylint: disable=no-member
+ 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")
+
+if not os.path.exists(figs_path):
+ os.mkdir(figs_path)
+
+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 0000000000..657765da4b
Binary files /dev/null and b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/m_frac_evap.npy differ
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 0000000000..5f4e75e31d
Binary files /dev/null and b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/r0grid.npy differ
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 0000000000..dba178b0df
Binary files /dev/null and b/tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/r_mins.npy differ
diff --git a/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py b/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py
new file mode 100644
index 0000000000..a15d1374b5
--- /dev/null
+++ b/tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py
@@ -0,0 +1,269 @@
+# pylint: disable=missing-module-docstring
+
+from collections import namedtuple
+from functools import partial
+from contextlib import contextmanager
+import pytest
+import numpy as np
+from scipy.optimize import fsolve
+
+from PySDM_examples.Loftus_and_Wordsworth_2021.planet import (
+ EarthLike,
+ Earth,
+ EarlyMars,
+ Jupiter,
+ Saturn,
+ K2_18B,
+)
+from PySDM_examples.Loftus_and_Wordsworth_2021.simulation import Simulation
+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:
+
+ @contextmanager
+ @staticmethod
+ def _get_test_resources():
+ formulae = Formulae(
+ ventilation="PruppacherAndRasmussen1979",
+ saturation_vapour_pressure="AugustRocheMagnus",
+ diffusion_coordinate="WaterMassLogarithm",
+ )
+ earth_like = EarthLike()
+ try:
+ yield formulae, earth_like
+ finally:
+ pass
+
+ def test_planet_classes(self):
+ """Test planet class instantiation and basic properties."""
+ planets = [EarthLike(), Earth(), EarlyMars(), Jupiter(), Saturn(), K2_18B()]
+
+ for planet in planets:
+ assert planet.g_std > 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} "
+ + f"exceeds 1.01 for {planet.__class__.__name__}"
+ )
+
+ def test_water_vapour_mixing_ratio_calculation(self):
+ """Test water vapour mixing ratio calculation."""
+ with TestLoftusWordsworth2021._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
+
+ # pylint: disable=too-many-arguments
+ @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 TestLoftusWordsworth2021._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
+
+ # pylint: disable=too-many-arguments
+ @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 TestLoftusWordsworth2021._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(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"
+
+ def test_saturation_at_cloud_base(self):
+ formulae = Formulae(
+ ventilation="PruppacherAndRasmussen1979",
+ saturation_vapour_pressure="AugustRocheMagnus",
+ diffusion_coordinate="WaterMassLogarithm",
+ )
+
+ new_Earth = EarthLike()
+
+ RH_array = np.linspace(0.25, 0.99, 50)
+ const = formulae.constants
+
+ def mix(dry, vap, ratio):
+ return (dry + ratio * vap) / (1 + ratio)
+
+ def f(x, water_mixing_ratio, params):
+ return water_mixing_ratio / (
+ water_mixing_ratio + const.eps
+ ) * params.p_stp * (x / params.t_stp) ** (
+ params.c_p / params.Rair
+ ) - formulae.saturation_vapour_pressure.pvs_water(
+ x
+ )
+
+ for RH in RH_array[::-1]:
+ 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
+ )
+
+ 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,
+ 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],
+ )
+ Tcloud = np.max(tdews)
+ thstd = formulae.trivia.th_std(new_Earth.p_STP, new_Earth.T_STP)
+
+ 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(
+ 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",
],
}