From 11be0004ce9867645561aeb9be6cb2e3db0313ea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Tue, 27 May 2025 11:22:04 +0200 Subject: [PATCH 01/23] add isotopes to simulation kinematic 2d --- .../PySDM_examples/Strzabala_2026/__init__.py | 0 .../kinematic_2d_isotopes.ipynb | 56 +++++++++++++++++++ .../utils/kinematic_2d/simulation.py | 9 ++- 3 files changed, 64 insertions(+), 1 deletion(-) create mode 100644 examples/PySDM_examples/Strzabala_2026/__init__.py create mode 100644 examples/PySDM_examples/Strzabala_2026/kinematic_2d_isotopes.ipynb diff --git a/examples/PySDM_examples/Strzabala_2026/__init__.py b/examples/PySDM_examples/Strzabala_2026/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/examples/PySDM_examples/Strzabala_2026/kinematic_2d_isotopes.ipynb b/examples/PySDM_examples/Strzabala_2026/kinematic_2d_isotopes.ipynb new file mode 100644 index 0000000000..2eaf4b0108 --- /dev/null +++ b/examples/PySDM_examples/Strzabala_2026/kinematic_2d_isotopes.ipynb @@ -0,0 +1,56 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "initial_id", + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "from PySDM_examples.Arabas_et_al_2015 import Settings, SpinUp\n", + "from PySDM_examples.utils.kinematic_2d import Simulation" + ] + }, + { + "metadata": {}, + "cell_type": "code", + "outputs": [], + "execution_count": null, + "source": [ + "settings = Settings()\n", + "storage = Storage()\n", + "simulation = Simulation(settings, storage, SpinUp=SpinUp)\n", + "advectees_init_profiles = {\n", + " \"HDO vapour mixing ratio\": np.full(shape=settings.grid[1], fill_value=1e-6),\n", + " \"H18O vapour mixing ratio\": np.full(shape=settings.grid[1], fill_value=1e-6),\n", + " \"H17O vapour mixing ratio\": np.full(shape=settings.grid[1], fill_value=1e-6),\n", + "}\n", + "simulation.reinit(additional_advectees_initial_profiles=advectees_init_profiles)" + ], + "id": "c340b54ecdfefc85" + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/PySDM_examples/utils/kinematic_2d/simulation.py b/examples/PySDM_examples/utils/kinematic_2d/simulation.py index cbe0c85e12..247b033245 100644 --- a/examples/PySDM_examples/utils/kinematic_2d/simulation.py +++ b/examples/PySDM_examples/utils/kinematic_2d/simulation.py @@ -1,4 +1,6 @@ import numpy as np +from numpy import ndarray + from PySDM_examples.utils.kinematic_2d.make_default_product_collection import ( make_default_product_collection, ) @@ -32,7 +34,11 @@ def __init__(self, settings, storage, SpinUp, backend_class=CPU): def products(self): return self.particulator.products - def reinit(self, products=None): + def reinit( + self, + products=None, + additional_advectees_initial_profiles: Dict[str, ndarray] = None, + ): formulae = self.settings.formulae backend = self.backend_class(formulae=formulae) environment = Kinematic2D( @@ -70,6 +76,7 @@ def reinit(self, products=None): initial_profiles = { "th": self.settings.initial_dry_potential_temperature_profile, "water_vapour_mixing_ratio": self.settings.initial_vapour_mixing_ratio_profile, + **(additional_advectees_initial_profiles or {}), } advectees = dict( ( From 678e3e5bc888301fdbfc9cb302610de90b5338b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Tue, 27 May 2025 16:38:55 +0200 Subject: [PATCH 02/23] add run simulation; fix numpy import and dict declaration --- .../kinematic_2d_isotopes.ipynb | 146 ++++++++++++++++-- .../utils/kinematic_2d/simulation.py | 5 +- 2 files changed, 139 insertions(+), 12 deletions(-) diff --git a/examples/PySDM_examples/Strzabala_2026/kinematic_2d_isotopes.ipynb b/examples/PySDM_examples/Strzabala_2026/kinematic_2d_isotopes.ipynb index 2eaf4b0108..345f4560b3 100644 --- a/examples/PySDM_examples/Strzabala_2026/kinematic_2d_isotopes.ipynb +++ b/examples/PySDM_examples/Strzabala_2026/kinematic_2d_isotopes.ipynb @@ -1,23 +1,55 @@ { "cells": [ + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-05-27T14:35:58.795167Z", + "start_time": "2025-05-27T14:35:58.791816Z" + } + }, + "cell_type": "code", + "source": [ + "import os\n", + "import subprocess\n", + "SUBPROCESS_ENV = os.environ.copy()" + ], + "id": "5dd595c4c0edac06", + "outputs": [], + "execution_count": 1 + }, { "cell_type": "code", - "execution_count": null, "id": "initial_id", "metadata": { - "collapsed": true + "collapsed": true, + "ExecuteTime": { + "end_time": "2025-05-27T14:36:02.233992Z", + "start_time": "2025-05-27T14:35:58.798955Z" + } }, - "outputs": [], "source": [ + "import numpy as np\n", + "\n", "from PySDM_examples.Arabas_et_al_2015 import Settings, SpinUp\n", - "from PySDM_examples.utils.kinematic_2d import Simulation" - ] + "from PySDM_examples.utils.kinematic_2d import Simulation, Storage\n", + "from PySDM.exporters import VTKExporter\n", + "from PySDM_examples.utils import ProgBarController\n", + "import PySDM_examples\n", + "import glob\n", + "import platform\n", + "import pathlib" + ], + "outputs": [], + "execution_count": 2 }, { - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2025-05-27T14:36:10.341002Z", + "start_time": "2025-05-27T14:36:02.342762Z" + } + }, "cell_type": "code", - "outputs": [], - "execution_count": null, "source": [ "settings = Settings()\n", "storage = Storage()\n", @@ -29,7 +61,103 @@ "}\n", "simulation.reinit(additional_advectees_initial_profiles=advectees_init_profiles)" ], - "id": "c340b54ecdfefc85" + "id": "c340b54ecdfefc85", + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/agnieszkazaba/PycharmProjects/PySDM/PySDM/backends/numba.py:48: UserWarning: Disabling Numba threading due to ARM64 CPU (atomics do not work yet)\n", + " warnings.warn(\n" + ] + } + ], + "execution_count": 3 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-05-27T14:37:26.546427Z", + "start_time": "2025-05-27T14:36:10.357658Z" + } + }, + "cell_type": "code", + "source": [ + "vtk_exporter = VTKExporter(path='.') \n", + "\n", + "simulation.run(ProgBarController(\"progress:\"), vtk_exporter=vtk_exporter)\n", + "vtk_exporter.write_pvd()" + ], + "id": "ab898240155a0e78", + "outputs": [ + { + "data": { + "text/plain": [ + "FloatProgress(value=0.0, description='progress:', max=1.0)" + ], + "application/vnd.jupyter.widget-view+json": { + "version_major": 2, + "version_minor": 0, + "model_id": "00c2aab707a84b0db24eab9669549106" + } + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "execution_count": 4 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-05-27T14:37:58.157595Z", + "start_time": "2025-05-27T14:37:26.568105Z" + } + }, + "cell_type": "code", + "source": [ + "pvanim = pathlib.Path(PySDM_examples.__file__).parent / \"utils\" / \"pvanim.py\"\n", + "\n", + "product = pathlib.Path(\"./output/sd_products.pvd\").absolute()\n", + "attributes = pathlib.Path(\"./output/sd_attributes.pvd\").absolute()\n", + "\n", + "try:\n", + " result = subprocess.run(\n", + " [\n", + " \"pvpython\",\n", + " \"--force-offscreen-rendering\",\n", + " str(pvanim),\n", + " str(product),\n", + " str(attributes),\n", + " str(pathlib.Path('./output').absolute()),\n", + " \"--animationname\", \"docs_intro_animation.ogv\",\n", + " \"--animationframename\", \"last_animation_frame.pdf\"\n", + " ],\n", + " check=platform.system() != \"Windows\",\n", + " capture_output=True,\n", + " text=True,\n", + " env=SUBPROCESS_ENV,\n", + " )\n", + "except subprocess.CalledProcessError as e:\n", + " print(e.stderr)\n", + " assert False" + ], + "id": "ffc5ec194b89182d", + "outputs": [], + "execution_count": 5 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-05-27T14:37:58.175984Z", + "start_time": "2025-05-27T14:37:58.172186Z" + } + }, + "cell_type": "code", + "source": "", + "id": "c9130b52152a4dea", + "outputs": [], + "execution_count": null } ], "metadata": { diff --git a/examples/PySDM_examples/utils/kinematic_2d/simulation.py b/examples/PySDM_examples/utils/kinematic_2d/simulation.py index 247b033245..23e0a92673 100644 --- a/examples/PySDM_examples/utils/kinematic_2d/simulation.py +++ b/examples/PySDM_examples/utils/kinematic_2d/simulation.py @@ -1,5 +1,4 @@ import numpy as np -from numpy import ndarray from PySDM_examples.utils.kinematic_2d.make_default_product_collection import ( make_default_product_collection, @@ -34,10 +33,10 @@ def __init__(self, settings, storage, SpinUp, backend_class=CPU): def products(self): return self.particulator.products - def reinit( + def reinit( # pylint: disable=too-many-locals, too-many-branches, too-many-statements self, products=None, - additional_advectees_initial_profiles: Dict[str, ndarray] = None, + additional_advectees_initial_profiles: dict[str, np.ndarray] = None, ): formulae = self.settings.formulae backend = self.backend_class(formulae=formulae) From 1e80a351d175a39d95529d092679f2556566593a Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 25 Jun 2025 16:53:30 +0200 Subject: [PATCH 03/23] draft logic for backend isotope_fractionation routine; placeholder file for timescale attributes --- PySDM/attributes/isotopes/timescales.py | 4 ++ .../impl_numba/methods/isotope_methods.py | 53 ++++++++++++++++++- PySDM/particulator.py | 18 ++++++- 3 files changed, 71 insertions(+), 4 deletions(-) create mode 100644 PySDM/attributes/isotopes/timescales.py diff --git a/PySDM/attributes/isotopes/timescales.py b/PySDM/attributes/isotopes/timescales.py new file mode 100644 index 0000000000..6224ac50e8 --- /dev/null +++ b/PySDM/attributes/isotopes/timescales.py @@ -0,0 +1,4 @@ +# TODO: register one class per isotope, +# and in each class call the relevant physics formula for a given isotope + +# TODO: code analogous to that in the delta.py file here diff --git a/PySDM/backends/impl_numba/methods/isotope_methods.py b/PySDM/backends/impl_numba/methods/isotope_methods.py index f33e6761dc..9fb57cb3b7 100644 --- a/PySDM/backends/impl_numba/methods/isotope_methods.py +++ b/PySDM/backends/impl_numba/methods/isotope_methods.py @@ -24,5 +24,54 @@ def body(output, ratio, reference_ratio): def isotopic_delta(self, output, ratio, reference_ratio): self._isotopic_delta_body(output.data, ratio.data, reference_ratio) - def isotopic_fractionation(self): - pass + @cached_property + def _isotopic_fractionation_body(self): + @numba.njit(**self.default_jit_flags) + def body( + *, + multiplicity, + signed_timescale, + moles, + dt, + ambient_isotope_mixing_ratio, + cell_id, + cell_volume, + dry_air_density, + molar_mass, + ): + for sd_id in numba.prange(multiplicity.shape[0]): + dn = dt / signed_timescale[sd_id] * moles[sd_id] + moles[sd_id] += dn + mass_of_dry_air = ( + dry_air_density[cell_id[sd_id]] * cell_volume[cell_id[sd_id]] + ) + ambient_isotope_mixing_ratio[cell_id[sd_id]] -= ( + dn * multiplicity[sd_id] * molar_mass / mass_of_dry_air + ) + + return body + + def isotopic_fractionation( + self, + *, + isotope, + multiplicity, + signed_timescale, + moles, + dt, + ambient_isotope_mixing_ratio, + cell_id, + cell_volume, + dry_air_density, + ): + self._isotopic_fractionation_body( + multiplicity=multiplicity.data, + signed_timescale=signed_timescale.data, + moles=moles.data, + dt=dt, + ambient_isotope_mixing_ratio=ambient_isotope_mixing_ratio.data, + cell_id=cell_id.data, + cell_volume=cell_volume.data, + dry_air_density=dry_air_density.data, + molar_mass=getattr(self.formulae.constants, f"M_{isotope}"), + ) diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 32ae0697d5..3a74a85893 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -439,9 +439,23 @@ def calculate_displacement( n_substeps=n_substeps, ) - def isotopic_fractionation(self, heavy_isotopes: tuple): - self.backend.isotopic_fractionation() + def isotopic_fractionation( + self, heavy_isotopes: tuple, ambient_isotope_mixing_ratios + ): for isotope in heavy_isotopes: + self.backend.isotopic_fractionation( + isotope=isotope, + multiplicity=self.attributes["multiplicity"], + signed_timescale=self.attributes[f"isotope_timescale_{isotope}"], + moles=self.attributes[f"moles_{isotope}"], + dt=self.dt, + dry_air_density=self.environment[f"dry_air_density"], + cell_volume=self.environment.mesh.dv, + cell_id=self.attributes["cell id"], + ambient_isotope_mixing_ratio=self.environment[ + f"mixing_ratio_{isotope}" + ], + ) self.attributes.mark_updated(f"moles_{isotope}") def seeding( From 4d76b220b0f237b6d229612661909fb3021deb4d Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 27 Jun 2025 17:46:14 +0200 Subject: [PATCH 04/23] Bolin's number! --- PySDM/attributes/isotopes/__init__.py | 2 +- PySDM/attributes/isotopes/bolins_numbers.py | 41 +++++++++++ PySDM/attributes/isotopes/timescales.py | 4 -- .../impl_numba/methods/isotope_methods.py | 70 ++++++++++++++++--- PySDM/particulator.py | 11 ++- .../miyake_et_al_1968.py | 4 ++ tests/unit_tests/attributes/test_isotopes.py | 49 ++++++++++++- .../backends/test_isotope_methods.py | 15 +++- .../dynamics/test_isotopic_fractionation.py | 1 - 9 files changed, 176 insertions(+), 21 deletions(-) create mode 100644 PySDM/attributes/isotopes/bolins_numbers.py delete mode 100644 PySDM/attributes/isotopes/timescales.py diff --git a/PySDM/attributes/isotopes/__init__.py b/PySDM/attributes/isotopes/__init__.py index 26875f6fe3..849fed6a34 100644 --- a/PySDM/attributes/isotopes/__init__.py +++ b/PySDM/attributes/isotopes/__init__.py @@ -3,4 +3,4 @@ """ from .moles import Moles1H, Moles16O, MolesLightWater -from ..isotopes import delta +from ..isotopes import delta, bolins_numbers diff --git a/PySDM/attributes/isotopes/bolins_numbers.py b/PySDM/attributes/isotopes/bolins_numbers.py new file mode 100644 index 0000000000..ca74bfd93f --- /dev/null +++ b/PySDM/attributes/isotopes/bolins_numbers.py @@ -0,0 +1,41 @@ +""" +# TODO: define Bolin's number +# TODO: consider positive/negative values? +# TODO: comment on total vs. light approximation +""" + +from PySDM.attributes.impl import DerivedAttribute, register_attribute +from PySDM.dynamics.isotopic_fractionation import HEAVY_ISOTOPES + + +class BolinsNumberImpl(DerivedAttribute): + def __init__(self, builder, *, heavy_isotope: str): + self.moles_heavy_isotope = builder.get_attribute("moles_" + heavy_isotope) + self.molar_mass = getattr( + builder.particulator.formulae.constants, f"M_{heavy_isotope}" + ) + super().__init__( + builder, + name="Bolin's number for " + heavy_isotope, + dependencies=(self.moles_heavy_isotope,), + ) + + def recalculate(self): + self.particulator.bolins_number( + output=self.data, + molar_mass=self.molar_mass, + moles_heavy_isotope=self.moles_heavy_isotope.data, + ) + + +def make_bolins_number_factory(heavy_isotope: str): + def _factory(builder): + return BolinsNumberImpl(builder, heavy_isotope=heavy_isotope) + + return _factory + + +for heavy_isotope in HEAVY_ISOTOPES: + register_attribute(name=f"Bolin's number for {heavy_isotope}")( + make_bolins_number_factory(heavy_isotope) + ) diff --git a/PySDM/attributes/isotopes/timescales.py b/PySDM/attributes/isotopes/timescales.py deleted file mode 100644 index 6224ac50e8..0000000000 --- a/PySDM/attributes/isotopes/timescales.py +++ /dev/null @@ -1,4 +0,0 @@ -# TODO: register one class per isotope, -# and in each class call the relevant physics formula for a given isotope - -# TODO: code analogous to that in the delta.py file here diff --git a/PySDM/backends/impl_numba/methods/isotope_methods.py b/PySDM/backends/impl_numba/methods/isotope_methods.py index 9fb57cb3b7..d623a52484 100644 --- a/PySDM/backends/impl_numba/methods/isotope_methods.py +++ b/PySDM/backends/impl_numba/methods/isotope_methods.py @@ -26,27 +26,49 @@ def isotopic_delta(self, output, ratio, reference_ratio): @cached_property def _isotopic_fractionation_body(self): - @numba.njit(**self.default_jit_flags) + @numba.njit(**{**self.default_jit_flags, "parallel": False}) def body( *, multiplicity, - signed_timescale, - moles, - dt, + dm_total, + bolins_number, + moles_light, + moles_heavy, + molar_mass_light, + molar_mass_heavy, ambient_isotope_mixing_ratio, cell_id, cell_volume, dry_air_density, - molar_mass, ): - for sd_id in numba.prange(multiplicity.shape[0]): - dn = dt / signed_timescale[sd_id] * moles[sd_id] - moles[sd_id] += dn + # input: + # tau_heavy (ignoring population/curvature effects) + # tau_light (ditto!) + # dm_total (actual - incl. population/curvature effects) + # output: + # dm_heavy = dm_total * tau_light / tau_heavy + + # tau' = m'/dm' * dt + # tau = m/dm * dt + # dm'/dm = tau/tau' * m'/m + + # dm' = dm * tau/tau' * m'/m + # ^^^^^^^ + # Bolin's 1/c1 + + for sd_id in range(multiplicity.shape[0]): + mass_ratio_heavy_to_light = (moles_heavy[sd_id] * molar_mass_heavy) / ( + moles_light[sd_id] * molar_mass_light + ) + dm_heavy_approx = ( + dm_total[sd_id] / bolins_number * mass_ratio_heavy_to_light + ) + moles_heavy[sd_id] += dm_heavy_approx / molar_mass_heavy mass_of_dry_air = ( dry_air_density[cell_id[sd_id]] * cell_volume[cell_id[sd_id]] ) ambient_isotope_mixing_ratio[cell_id[sd_id]] -= ( - dn * multiplicity[sd_id] * molar_mass / mass_of_dry_air + dm_heavy_approx * multiplicity[sd_id] / mass_of_dry_air ) return body @@ -54,7 +76,7 @@ def body( def isotopic_fractionation( self, *, - isotope, + molar_mass, multiplicity, signed_timescale, moles, @@ -73,5 +95,31 @@ def isotopic_fractionation( cell_id=cell_id.data, cell_volume=cell_volume.data, dry_air_density=dry_air_density.data, - molar_mass=getattr(self.formulae.constants, f"M_{isotope}"), + molar_mass=molar_mass, + ) + + @cached_property + def _bolins_number_body(self): + ff = self.formulae_flattened + + @numba.njit(**self.default_jit_flags) + def body(output, molar_mass, cell_id, moles_heavy_isotope, relative_humidity): + for i in numba.prange(output.shape[0]): # pylint: disable=not-an-iterable + output[i] = ff.isotope_relaxation_timescale__bolins_number( + moles_heavy_isotope=moles_heavy_isotope[i], + relative_humidity=relative_humidity[cell_id[i]], + molar_mass=molar_mass, + ) + + return body + + def bolins_number( + self, *, output, molar_mass, cell_id, moles_heavy_isotope, relative_humidity + ): + self._bolins_number_body( + output=output.data, + molar_mass=molar_mass, + cell_id=cell_id.data, + moles_heavy_isotope=moles_heavy_isotope.data, + relative_humidity=relative_humidity.data, ) diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 3a74a85893..bf2ff5692a 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -444,7 +444,7 @@ def isotopic_fractionation( ): for isotope in heavy_isotopes: self.backend.isotopic_fractionation( - isotope=isotope, + molar_mass=getattr(self.formulae.constants, f"M_{isotope}"), multiplicity=self.attributes["multiplicity"], signed_timescale=self.attributes[f"isotope_timescale_{isotope}"], moles=self.attributes[f"moles_{isotope}"], @@ -458,6 +458,15 @@ def isotopic_fractionation( ) self.attributes.mark_updated(f"moles_{isotope}") + def bolins_number(self, *, output, molar_mass, moles_heavy_isotope): + self.backend.bolins_number( + output=output, + molar_mass=molar_mass, + cell_id=self.attributes["cell id"], + moles_heavy_isotope=moles_heavy_isotope, + relative_humidity=self.environment["RH"], + ) + def seeding( self, *, diff --git a/PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py b/PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py index 5514e55b67..cc6645e462 100644 --- a/PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py +++ b/PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py @@ -11,3 +11,7 @@ def tau(const, e_s, D, M, vent_coeff, radius, alpha, temperature): return (radius**2 * alpha * const.rho_w * const.R_str * temperature) / ( 3 * e_s * D * M * vent_coeff ) + + @staticmethod + def bolins_number(const, moles_heavy_isotope, relative_humidity, molar_mass): + return 44.0 diff --git a/tests/unit_tests/attributes/test_isotopes.py b/tests/unit_tests/attributes/test_isotopes.py index 6cb3c2dd66..004190e6c0 100644 --- a/tests/unit_tests/attributes/test_isotopes.py +++ b/tests/unit_tests/attributes/test_isotopes.py @@ -5,7 +5,7 @@ import numpy as np import pytest -from PySDM import Builder +from PySDM import Builder, Formulae from PySDM.dynamics.isotopic_fractionation import HEAVY_ISOTOPES from PySDM.environments import Box from PySDM.physics import constants_defaults, si @@ -109,3 +109,50 @@ def test_delta_attribute(backend_class, isotope, heavier_water_specific_content) ), significant=5, ) + + @staticmethod + @pytest.mark.parametrize("heavy_isotope", HEAVY_ISOTOPES) + @pytest.mark.parametrize( + "moles_heavy_isotope, relative_humidity, expected_tau", + ( + (0, 0.99, 44), + (0, 1.01, 44), + ), + ) # TODO: check sign! + @pytest.mark.parametrize("variant", ("MiyakeEtAl1968",)) + def test_bolins_number_attribute( + backend_class, + heavy_isotope: str, + moles_heavy_isotope: float, + relative_humidity: float, + expected_tau: float, + variant: str, + ): + if backend_class.__name__ != "Numba": + pytest.skip("# TODO - isotopes on GPU") + + # arrange + n_sd = 1 + builder = Builder( + n_sd=n_sd, + backend=backend_class( + formulae=Formulae(isotope_relaxation_timescale=variant) + ), + environment=Box(dt=np.nan, dv=np.nan), + ) + attribute_name = f"Bolin's number for {heavy_isotope}" + builder.request_attribute(attribute_name) + particulator = builder.build( + attributes={ + "multiplicity": np.ones(n_sd), + "signed water mass": np.ones(n_sd) * si.ng, + f"moles_{heavy_isotope}": np.ones(n_sd) * moles_heavy_isotope, + } + ) + particulator.environment["RH"] = relative_humidity + + # act + value = particulator.attributes[attribute_name].to_ndarray() + + # assert + assert value == expected_tau diff --git a/tests/unit_tests/backends/test_isotope_methods.py b/tests/unit_tests/backends/test_isotope_methods.py index 05643be667..74f1cb05bb 100644 --- a/tests/unit_tests/backends/test_isotope_methods.py +++ b/tests/unit_tests/backends/test_isotope_methods.py @@ -8,11 +8,22 @@ class TestIsotopeMethods: @staticmethod def test_isotopic_fractionation(backend_instance): + """checks if for a given dm_dt of total liquid water, the changes + in the amounts of moles of heavy isotopes in droplets, + and in the ambient air are OK""" # arrange - backend = backend_instance + # tau = + # dm_dt = + # ambient_isotope_mixing_ratio= # act - backend.isotopic_fractionation() + # backend_instance.isotopic_fractionation( + # + # ) + + # assert + # assert dm_iso_dt == ... + # assert ambient_isotope_mixing_ratio == ... @staticmethod def test_isotopic_delta(backend_instance): diff --git a/tests/unit_tests/dynamics/test_isotopic_fractionation.py b/tests/unit_tests/dynamics/test_isotopic_fractionation.py index ffe2d15841..f37f23c6be 100644 --- a/tests/unit_tests/dynamics/test_isotopic_fractionation.py +++ b/tests/unit_tests/dynamics/test_isotopic_fractionation.py @@ -80,7 +80,6 @@ def test_call_marks_all_isotopes_as_updated(backend_instance): particulator.dynamics["IsotopicFractionation"]() # assert - for isotope in HEAVY_ISOTOPES: # pylint:disable=protected-access assert ( From 8ff4a0b627faa9587c065a2353f5ea493cac5466 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Mon, 30 Jun 2025 18:45:34 +0200 Subject: [PATCH 05/23] work in progress on isotopic_fractionation method and test --- PySDM/attributes/isotopes/__init__.py | 2 +- .../{bolins_numbers.py => bolin_numbers.py} | 6 +-- PySDM/attributes/physics/__init__.py | 1 + .../physics/diffusional_growth_mass_change.py | 13 +++++ .../impl_numba/methods/isotope_methods.py | 34 +++++++------- PySDM/dynamics/isotopic_fractionation.py | 2 + PySDM/particulator.py | 23 ++++++--- .../dynamics/test_isotopic_fractionation.py | 47 ++++++++++++++++++- 8 files changed, 100 insertions(+), 28 deletions(-) rename PySDM/attributes/isotopes/{bolins_numbers.py => bolin_numbers.py} (89%) create mode 100644 PySDM/attributes/physics/diffusional_growth_mass_change.py diff --git a/PySDM/attributes/isotopes/__init__.py b/PySDM/attributes/isotopes/__init__.py index 849fed6a34..a0c57cde30 100644 --- a/PySDM/attributes/isotopes/__init__.py +++ b/PySDM/attributes/isotopes/__init__.py @@ -3,4 +3,4 @@ """ from .moles import Moles1H, Moles16O, MolesLightWater -from ..isotopes import delta, bolins_numbers +from ..isotopes import delta, bolin_numbers diff --git a/PySDM/attributes/isotopes/bolins_numbers.py b/PySDM/attributes/isotopes/bolin_numbers.py similarity index 89% rename from PySDM/attributes/isotopes/bolins_numbers.py rename to PySDM/attributes/isotopes/bolin_numbers.py index ca74bfd93f..18c19ab7f6 100644 --- a/PySDM/attributes/isotopes/bolins_numbers.py +++ b/PySDM/attributes/isotopes/bolin_numbers.py @@ -8,7 +8,7 @@ from PySDM.dynamics.isotopic_fractionation import HEAVY_ISOTOPES -class BolinsNumberImpl(DerivedAttribute): +class BolinNumberImpl(DerivedAttribute): def __init__(self, builder, *, heavy_isotope: str): self.moles_heavy_isotope = builder.get_attribute("moles_" + heavy_isotope) self.molar_mass = getattr( @@ -16,12 +16,12 @@ def __init__(self, builder, *, heavy_isotope: str): ) super().__init__( builder, - name="Bolin's number for " + heavy_isotope, + name="Bolin number for " + heavy_isotope, dependencies=(self.moles_heavy_isotope,), ) def recalculate(self): - self.particulator.bolins_number( + self.particulator.bolin_number( output=self.data, molar_mass=self.molar_mass, moles_heavy_isotope=self.moles_heavy_isotope.data, diff --git a/PySDM/attributes/physics/__init__.py b/PySDM/attributes/physics/__init__.py index 44416bd9d9..699deab9c3 100644 --- a/PySDM/attributes/physics/__init__.py +++ b/PySDM/attributes/physics/__init__.py @@ -18,3 +18,4 @@ from .terminal_velocity import TerminalVelocity from .volume import Volume from .reynolds_number import ReynoldsNumber +from .diffusional_growth_mass_change import DiffusionalGrowthMassChange diff --git a/PySDM/attributes/physics/diffusional_growth_mass_change.py b/PySDM/attributes/physics/diffusional_growth_mass_change.py new file mode 100644 index 0000000000..4c38b1e55d --- /dev/null +++ b/PySDM/attributes/physics/diffusional_growth_mass_change.py @@ -0,0 +1,13 @@ +from PySDM.attributes.impl import register_attribute, DerivedAttribute + + +@register_attribute() +class DiffusionalGrowthMassChange(DerivedAttribute): + def __init__(self, builder): + attr = builder.get_attribute("water mass") + super().__init__( + builder, name="diffusional growth mass change", dependencies=(attr,) + ) + + def recalculate(self): + pass # TODO diff --git a/PySDM/backends/impl_numba/methods/isotope_methods.py b/PySDM/backends/impl_numba/methods/isotope_methods.py index d623a52484..97d9818eab 100644 --- a/PySDM/backends/impl_numba/methods/isotope_methods.py +++ b/PySDM/backends/impl_numba/methods/isotope_methods.py @@ -26,6 +26,8 @@ def isotopic_delta(self, output, ratio, reference_ratio): @cached_property def _isotopic_fractionation_body(self): + molar_mass_light = self.formulae.constants.M_1H2_16O + @numba.njit(**{**self.default_jit_flags, "parallel": False}) def body( *, @@ -34,7 +36,6 @@ def body( bolins_number, moles_light, moles_heavy, - molar_mass_light, molar_mass_heavy, ambient_isotope_mixing_ratio, cell_id, @@ -57,16 +58,15 @@ def body( # Bolin's 1/c1 for sd_id in range(multiplicity.shape[0]): + moles_light = signed_water_mass[sd_id] / molar_mass_light mass_ratio_heavy_to_light = (moles_heavy[sd_id] * molar_mass_heavy) / ( - moles_light[sd_id] * molar_mass_light + moles_light * molar_mass_light ) dm_heavy_approx = ( dm_total[sd_id] / bolins_number * mass_ratio_heavy_to_light ) moles_heavy[sd_id] += dm_heavy_approx / molar_mass_heavy - mass_of_dry_air = ( - dry_air_density[cell_id[sd_id]] * cell_volume[cell_id[sd_id]] - ) + mass_of_dry_air = dry_air_density[cell_id[sd_id]] * cell_volume ambient_isotope_mixing_ratio[cell_id[sd_id]] -= ( dm_heavy_approx * multiplicity[sd_id] / mass_of_dry_air ) @@ -76,11 +76,12 @@ def body( def isotopic_fractionation( self, *, - molar_mass, multiplicity, - signed_timescale, - moles, - dt, + dm_total, + bolin_number, + signed_water_mass, + moles_heavy, + molar_mass_heavy, ambient_isotope_mixing_ratio, cell_id, cell_volume, @@ -88,14 +89,15 @@ def isotopic_fractionation( ): self._isotopic_fractionation_body( multiplicity=multiplicity.data, - signed_timescale=signed_timescale.data, - moles=moles.data, - dt=dt, + dm_total=dm_total.data, + bolin_number=bolin_number.data, + signed_water_mass=signed_water_mass.data, + moles_heavy=moles_heavy.data, + molar_mass_heavy=molar_mass_heavy, ambient_isotope_mixing_ratio=ambient_isotope_mixing_ratio.data, cell_id=cell_id.data, - cell_volume=cell_volume.data, + cell_volume=cell_volume, dry_air_density=dry_air_density.data, - molar_mass=molar_mass, ) @cached_property @@ -113,10 +115,10 @@ def body(output, molar_mass, cell_id, moles_heavy_isotope, relative_humidity): return body - def bolins_number( + def bolin_number( self, *, output, molar_mass, cell_id, moles_heavy_isotope, relative_humidity ): - self._bolins_number_body( + self._bolin_number_body( output=output.data, molar_mass=molar_mass, cell_id=cell_id.data, diff --git a/PySDM/dynamics/isotopic_fractionation.py b/PySDM/dynamics/isotopic_fractionation.py index dffeea5160..1ca10df28a 100644 --- a/PySDM/dynamics/isotopic_fractionation.py +++ b/PySDM/dynamics/isotopic_fractionation.py @@ -33,8 +33,10 @@ def register(self, builder): f"{Condensation.__name__} needs to be registered to run prior to {self.__class__}" ) + builder.request_attribute("diffusional growth mass change") for isotope in self.isotopes: builder.request_attribute(f"moles_{isotope}") + builder.request_attribute(f"Bolin number for {isotope}") def __call__(self): self.particulator.isotopic_fractionation(self.isotopes) diff --git a/PySDM/particulator.py b/PySDM/particulator.py index bf2ff5692a..0fab24d9bd 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -439,16 +439,25 @@ def calculate_displacement( n_substeps=n_substeps, ) - def isotopic_fractionation( - self, heavy_isotopes: tuple, ambient_isotope_mixing_ratios - ): + def isotopic_fractionation(self, heavy_isotopes: tuple): for isotope in heavy_isotopes: self.backend.isotopic_fractionation( - molar_mass=getattr(self.formulae.constants, f"M_{isotope}"), + molar_mass_heavy=getattr( + self.formulae.constants, + { + "2H": "M_2H_1H_16O", + "3H": "M_3H_1H_16O", + "17O": "M_1H2_17O", + "18O": "M_1H2_18O", + }[isotope], + ), + dm_total=self.attributes[ + f"diffusional growth mass change" + ], # TODO: total vs. heavy + bolin_number=self.attributes[f"Bolin number for {isotope}"], + signed_water_mass=self.attributes["signed water mass"], multiplicity=self.attributes["multiplicity"], - signed_timescale=self.attributes[f"isotope_timescale_{isotope}"], - moles=self.attributes[f"moles_{isotope}"], - dt=self.dt, + moles_heavy=self.attributes[f"moles_{isotope}"], dry_air_density=self.environment[f"dry_air_density"], cell_volume=self.environment.mesh.dv, cell_id=self.attributes["cell id"], diff --git a/tests/unit_tests/dynamics/test_isotopic_fractionation.py b/tests/unit_tests/dynamics/test_isotopic_fractionation.py index f37f23c6be..914ce24a09 100644 --- a/tests/unit_tests/dynamics/test_isotopic_fractionation.py +++ b/tests/unit_tests/dynamics/test_isotopic_fractionation.py @@ -7,7 +7,7 @@ import numpy as np import pytest -from PySDM import Builder +from PySDM import Builder, Formulae from PySDM.dynamics import Condensation, IsotopicFractionation from PySDM.dynamics.isotopic_fractionation import HEAVY_ISOTOPES from PySDM.environments import Box @@ -90,3 +90,48 @@ def test_call_marks_all_isotopes_as_updated(backend_instance): "multiplicity" ].timestamp ) + + @staticmethod + def test_no_isotope_fractionation_if_droplet_size_unchanged(backend_class): + """neither a bug nor a feature :) - just a simplification (?)""" + # arrange + ambient_initial_isotope_mixing_ratio = 44.0 + particle_initial_isotope_content = 666.0 * si.moles + cell_volume = 1 * si.m**3 + + attributes = DUMMY_ATTRIBUTES.copy() + attributes["moles_2H"] = particle_initial_isotope_content + attributes["water mass"] = 1 * si.ng + attributes["multiplicity"] = np.ones(1) + + builder = Builder( + n_sd=1, + backend=backend_class( + formulae=Formulae(isotope_relaxation_timescale="MiyakeEtAl1968"), + ), + environment=Box(dv=cell_volume, dt=-1 * si.s), + ) + builder.add_dynamic(Condensation()) + builder.add_dynamic(IsotopicFractionation(isotopes=("2H",))) + particulator = builder.build(attributes=attributes, products=()) + particulator.environment["RH"] = 1 + particulator.environment["dry_air_density"] = 1 * si.kg / si.m**3 + particulator.environment["mixing_ratio_2H"] = ( + ambient_initial_isotope_mixing_ratio + ) + assert ( + particulator.attributes["moles_2H"][0] == particle_initial_isotope_content + ) + + # act + particulator.attributes["diffusional growth mass change"].data[:] = 0 + particulator.dynamics["IsotopicFractionation"]() + + # assert + assert ( + particulator.environment["mixing_ratio_2H"][0] + == ambient_initial_isotope_mixing_ratio + ) + assert ( + particulator.attributes["moles_2H"][0] == particle_initial_isotope_content + ) From 2f00403b12ecc16ab22dfea5d1795c3eedcd5ae0 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 27 Jun 2025 17:46:14 +0200 Subject: [PATCH 06/23] Bolin's number! # Conflicts: # PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py --- PySDM/attributes/isotopes/bolins_numbers.py | 41 +++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 PySDM/attributes/isotopes/bolins_numbers.py diff --git a/PySDM/attributes/isotopes/bolins_numbers.py b/PySDM/attributes/isotopes/bolins_numbers.py new file mode 100644 index 0000000000..ca74bfd93f --- /dev/null +++ b/PySDM/attributes/isotopes/bolins_numbers.py @@ -0,0 +1,41 @@ +""" +# TODO: define Bolin's number +# TODO: consider positive/negative values? +# TODO: comment on total vs. light approximation +""" + +from PySDM.attributes.impl import DerivedAttribute, register_attribute +from PySDM.dynamics.isotopic_fractionation import HEAVY_ISOTOPES + + +class BolinsNumberImpl(DerivedAttribute): + def __init__(self, builder, *, heavy_isotope: str): + self.moles_heavy_isotope = builder.get_attribute("moles_" + heavy_isotope) + self.molar_mass = getattr( + builder.particulator.formulae.constants, f"M_{heavy_isotope}" + ) + super().__init__( + builder, + name="Bolin's number for " + heavy_isotope, + dependencies=(self.moles_heavy_isotope,), + ) + + def recalculate(self): + self.particulator.bolins_number( + output=self.data, + molar_mass=self.molar_mass, + moles_heavy_isotope=self.moles_heavy_isotope.data, + ) + + +def make_bolins_number_factory(heavy_isotope: str): + def _factory(builder): + return BolinsNumberImpl(builder, heavy_isotope=heavy_isotope) + + return _factory + + +for heavy_isotope in HEAVY_ISOTOPES: + register_attribute(name=f"Bolin's number for {heavy_isotope}")( + make_bolins_number_factory(heavy_isotope) + ) From ac9c7436c40de75a05d21fb191f3e7f0091f2f51 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Fri, 4 Jul 2025 13:23:21 +0200 Subject: [PATCH 07/23] change bolins/Bolin's to bolin/Bolin --- PySDM/attributes/isotopes/bolin_numbers.py | 10 +++++----- PySDM/attributes/isotopes/bolins_numbers.py | 16 ++++++++-------- .../impl_numba/methods/isotope_methods.py | 8 ++++---- PySDM/particulator.py | 4 ++-- .../miyake_et_al_1968.py | 7 +++---- tests/unit_tests/attributes/test_isotopes.py | 4 ++-- 6 files changed, 24 insertions(+), 25 deletions(-) diff --git a/PySDM/attributes/isotopes/bolin_numbers.py b/PySDM/attributes/isotopes/bolin_numbers.py index 18c19ab7f6..9246052751 100644 --- a/PySDM/attributes/isotopes/bolin_numbers.py +++ b/PySDM/attributes/isotopes/bolin_numbers.py @@ -1,5 +1,5 @@ """ -# TODO: define Bolin's number +# TODO: define Bolin number # TODO: consider positive/negative values? # TODO: comment on total vs. light approximation """ @@ -28,14 +28,14 @@ def recalculate(self): ) -def make_bolins_number_factory(heavy_isotope: str): +def make_bolin_number_factory(heavy_isotope: str): def _factory(builder): - return BolinsNumberImpl(builder, heavy_isotope=heavy_isotope) + return BolinNumberImpl(builder, heavy_isotope=heavy_isotope) return _factory for heavy_isotope in HEAVY_ISOTOPES: - register_attribute(name=f"Bolin's number for {heavy_isotope}")( - make_bolins_number_factory(heavy_isotope) + register_attribute(name=f"Bolin number for {heavy_isotope}")( + make_bolin_number_factory(heavy_isotope) ) diff --git a/PySDM/attributes/isotopes/bolins_numbers.py b/PySDM/attributes/isotopes/bolins_numbers.py index ca74bfd93f..9246052751 100644 --- a/PySDM/attributes/isotopes/bolins_numbers.py +++ b/PySDM/attributes/isotopes/bolins_numbers.py @@ -1,5 +1,5 @@ """ -# TODO: define Bolin's number +# TODO: define Bolin number # TODO: consider positive/negative values? # TODO: comment on total vs. light approximation """ @@ -8,7 +8,7 @@ from PySDM.dynamics.isotopic_fractionation import HEAVY_ISOTOPES -class BolinsNumberImpl(DerivedAttribute): +class BolinNumberImpl(DerivedAttribute): def __init__(self, builder, *, heavy_isotope: str): self.moles_heavy_isotope = builder.get_attribute("moles_" + heavy_isotope) self.molar_mass = getattr( @@ -16,26 +16,26 @@ def __init__(self, builder, *, heavy_isotope: str): ) super().__init__( builder, - name="Bolin's number for " + heavy_isotope, + name="Bolin number for " + heavy_isotope, dependencies=(self.moles_heavy_isotope,), ) def recalculate(self): - self.particulator.bolins_number( + self.particulator.bolin_number( output=self.data, molar_mass=self.molar_mass, moles_heavy_isotope=self.moles_heavy_isotope.data, ) -def make_bolins_number_factory(heavy_isotope: str): +def make_bolin_number_factory(heavy_isotope: str): def _factory(builder): - return BolinsNumberImpl(builder, heavy_isotope=heavy_isotope) + return BolinNumberImpl(builder, heavy_isotope=heavy_isotope) return _factory for heavy_isotope in HEAVY_ISOTOPES: - register_attribute(name=f"Bolin's number for {heavy_isotope}")( - make_bolins_number_factory(heavy_isotope) + register_attribute(name=f"Bolin number for {heavy_isotope}")( + make_bolin_number_factory(heavy_isotope) ) diff --git a/PySDM/backends/impl_numba/methods/isotope_methods.py b/PySDM/backends/impl_numba/methods/isotope_methods.py index 97d9818eab..b34f43499f 100644 --- a/PySDM/backends/impl_numba/methods/isotope_methods.py +++ b/PySDM/backends/impl_numba/methods/isotope_methods.py @@ -33,7 +33,7 @@ def body( *, multiplicity, dm_total, - bolins_number, + bolin_number, moles_light, moles_heavy, molar_mass_heavy, @@ -63,7 +63,7 @@ def body( moles_light * molar_mass_light ) dm_heavy_approx = ( - dm_total[sd_id] / bolins_number * mass_ratio_heavy_to_light + dm_total[sd_id] / bolin_number * mass_ratio_heavy_to_light ) moles_heavy[sd_id] += dm_heavy_approx / molar_mass_heavy mass_of_dry_air = dry_air_density[cell_id[sd_id]] * cell_volume @@ -101,13 +101,13 @@ def isotopic_fractionation( ) @cached_property - def _bolins_number_body(self): + def _bolin_number_body(self): ff = self.formulae_flattened @numba.njit(**self.default_jit_flags) def body(output, molar_mass, cell_id, moles_heavy_isotope, relative_humidity): for i in numba.prange(output.shape[0]): # pylint: disable=not-an-iterable - output[i] = ff.isotope_relaxation_timescale__bolins_number( + output[i] = ff.isotope_relaxation_timescale__bolin_number( moles_heavy_isotope=moles_heavy_isotope[i], relative_humidity=relative_humidity[cell_id[i]], molar_mass=molar_mass, diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 0fab24d9bd..5b724ca0ec 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -467,8 +467,8 @@ def isotopic_fractionation(self, heavy_isotopes: tuple): ) self.attributes.mark_updated(f"moles_{isotope}") - def bolins_number(self, *, output, molar_mass, moles_heavy_isotope): - self.backend.bolins_number( + def bolin_number(self, *, output, molar_mass, moles_heavy_isotope): + self.backend.bolin_number( output=output, molar_mass=molar_mass, cell_id=self.attributes["cell id"], diff --git a/PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py b/PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py index 62971323a4..625b598c82 100644 --- a/PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py +++ b/PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py @@ -6,7 +6,7 @@ class MiyakeEtAl1968: # pylint:disable=too-few-public-methods def __init__(self, _): pass - + @staticmethod def tau( const, rho_s, radius, D_iso, D, S, R_liq, alpha, R_vap, Fk @@ -22,8 +22,7 @@ def tau( diffusivity * theta, where theta from eq. (25) is ventilation_coefficient """ return (radius**2 * alpha * const.rho_w) / (3 * rho_s * D) - - + @staticmethod - def bolins_number(const, moles_heavy_isotope, relative_humidity, molar_mass): + def bolin_number(const, moles_heavy_isotope, relative_humidity, molar_mass): return 44.0 diff --git a/tests/unit_tests/attributes/test_isotopes.py b/tests/unit_tests/attributes/test_isotopes.py index 004190e6c0..5c98b58ef8 100644 --- a/tests/unit_tests/attributes/test_isotopes.py +++ b/tests/unit_tests/attributes/test_isotopes.py @@ -120,7 +120,7 @@ def test_delta_attribute(backend_class, isotope, heavier_water_specific_content) ), ) # TODO: check sign! @pytest.mark.parametrize("variant", ("MiyakeEtAl1968",)) - def test_bolins_number_attribute( + def test_bolin_number_attribute( backend_class, heavy_isotope: str, moles_heavy_isotope: float, @@ -140,7 +140,7 @@ def test_bolins_number_attribute( ), environment=Box(dt=np.nan, dv=np.nan), ) - attribute_name = f"Bolin's number for {heavy_isotope}" + attribute_name = f"Bolin number for {heavy_isotope}" builder.request_attribute(attribute_name) particulator = builder.build( attributes={ From a5a37ba0f4d9d51ad0810d2b6e4e2663595c24d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Fri, 4 Jul 2025 13:26:40 +0200 Subject: [PATCH 08/23] pylint hints addressed --- PySDM/particulator.py | 2 +- .../physics/isotope_relaxation_timescale/miyake_et_al_1968.py | 4 +++- tests/unit_tests/attributes/test_isotopes.py | 2 +- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 5b724ca0ec..4faa2c3b34 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -458,7 +458,7 @@ def isotopic_fractionation(self, heavy_isotopes: tuple): signed_water_mass=self.attributes["signed water mass"], multiplicity=self.attributes["multiplicity"], moles_heavy=self.attributes[f"moles_{isotope}"], - dry_air_density=self.environment[f"dry_air_density"], + dry_air_density=self.environment["dry_air_density"], cell_volume=self.environment.mesh.dv, cell_id=self.attributes["cell id"], ambient_isotope_mixing_ratio=self.environment[ diff --git a/PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py b/PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py index 625b598c82..e164884184 100644 --- a/PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py +++ b/PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py @@ -24,5 +24,7 @@ def tau( return (radius**2 * alpha * const.rho_w) / (3 * rho_s * D) @staticmethod - def bolin_number(const, moles_heavy_isotope, relative_humidity, molar_mass): + def bolin_number( + const, moles_heavy_isotope, relative_humidity, molar_mass + ): # pylint: disable=unused-argument return 44.0 diff --git a/tests/unit_tests/attributes/test_isotopes.py b/tests/unit_tests/attributes/test_isotopes.py index 5c98b58ef8..aaafe3c73d 100644 --- a/tests/unit_tests/attributes/test_isotopes.py +++ b/tests/unit_tests/attributes/test_isotopes.py @@ -127,7 +127,7 @@ def test_bolin_number_attribute( relative_humidity: float, expected_tau: float, variant: str, - ): + ): # pylint: disable=too-many-arguments if backend_class.__name__ != "Numba": pytest.skip("# TODO - isotopes on GPU") From 707dec337daadc7604bdb4bbb19784a6d809c41d Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 4 Jul 2025 16:50:22 +0200 Subject: [PATCH 09/23] implement diffusional growth mass change; draft a test for it; other code cleanups --- PySDM/attributes/isotopes/moles.py | 19 +++------ .../physics/diffusional_growth_mass_change.py | 14 +++++-- .../impl_numba/methods/isotope_methods.py | 40 ++++++++---------- PySDM/particulator.py | 4 +- .../test_diffusional_growth_mass_change.py | 41 +++++++++++++++++++ .../dynamics/test_isotopic_fractionation.py | 4 +- .../dynamics/test_vapour_deposition_on_ice.py | 1 + 7 files changed, 79 insertions(+), 44 deletions(-) create mode 100644 tests/unit_tests/attributes/test_diffusional_growth_mass_change.py diff --git a/PySDM/attributes/isotopes/moles.py b/PySDM/attributes/isotopes/moles.py index 9f4b08b003..4bf71cb938 100644 --- a/PySDM/attributes/isotopes/moles.py +++ b/PySDM/attributes/isotopes/moles.py @@ -45,24 +45,15 @@ def recalculate(self): class MolesLightWater(Helper): def __init__(self, builder): const = builder.formulae.constants - M_H2O = 2 * const.M_1H + const.M_16O super().__init__( builder=builder, name="moles light water", attrs_to_multiplier={ - builder.get_attribute("moles_2H"): -( - const.M_1H * const.M_2H + const.M_16O - ) - / M_H2O, - builder.get_attribute("moles_3H"): -( - const.M_1H * const.M_3H + const.M_16O - ) - / M_H2O, - builder.get_attribute("moles_17O"): -(2 * const.M_1H + const.M_17O) - / M_H2O, - builder.get_attribute("moles_18O"): -(2 * const.M_1H + const.M_18O) - / M_H2O, - builder.get_attribute("signed water mass"): 1 / M_H2O, + builder.get_attribute("moles_2H"): -const.M_2H_1H_16O / const.M_1H2_16O, + builder.get_attribute("moles_3H"): -const.M_3H_1H_16O / const.M_1H2_16O, + builder.get_attribute("moles_17O"): -const.M_1H2_17O / const.M_1H2_16O, + builder.get_attribute("moles_18O"): -const.M_1H2_18O / const.M_1H2_16O, + builder.get_attribute("signed water mass"): 1 / const.M_1H2_16O, }, ) diff --git a/PySDM/attributes/physics/diffusional_growth_mass_change.py b/PySDM/attributes/physics/diffusional_growth_mass_change.py index 4c38b1e55d..63100a4b26 100644 --- a/PySDM/attributes/physics/diffusional_growth_mass_change.py +++ b/PySDM/attributes/physics/diffusional_growth_mass_change.py @@ -4,10 +4,18 @@ @register_attribute() class DiffusionalGrowthMassChange(DerivedAttribute): def __init__(self, builder): - attr = builder.get_attribute("water mass") + self.water_mass = builder.get_attribute("water mass") super().__init__( - builder, name="diffusional growth mass change", dependencies=(attr,) + builder, + name="diffusional growth mass change", + dependencies=(self.water_mass,), ) + builder.particulator.observers.append(self) + assert "Collision" not in builder.particulator.dynamics + self.notify() + + def notify(self): + self.data[:] = -self.water_mass.data def recalculate(self): - pass # TODO + self.data[:] += self.water_mass.data diff --git a/PySDM/backends/impl_numba/methods/isotope_methods.py b/PySDM/backends/impl_numba/methods/isotope_methods.py index b34f43499f..27981151be 100644 --- a/PySDM/backends/impl_numba/methods/isotope_methods.py +++ b/PySDM/backends/impl_numba/methods/isotope_methods.py @@ -26,7 +26,6 @@ def isotopic_delta(self, output, ratio, reference_ratio): @cached_property def _isotopic_fractionation_body(self): - molar_mass_light = self.formulae.constants.M_1H2_16O @numba.njit(**{**self.default_jit_flags, "parallel": False}) def body( @@ -34,7 +33,7 @@ def body( multiplicity, dm_total, bolin_number, - moles_light, + signed_water_mass, moles_heavy, molar_mass_heavy, ambient_isotope_mixing_ratio, @@ -42,33 +41,28 @@ def body( cell_volume, dry_air_density, ): - # input: - # tau_heavy (ignoring population/curvature effects) - # tau_light (ditto!) + # assume Bo = tau' / tau_total + # = (m'/dm') / (m/dm) + # = m'/m * dm/dm' + # input (per droplet: + # Bo (discarding curvature/population effects) # dm_total (actual - incl. population/curvature effects) # output: - # dm_heavy = dm_total * tau_light / tau_heavy - - # tau' = m'/dm' * dt - # tau = m/dm * dt - # dm'/dm = tau/tau' * m'/m - - # dm' = dm * tau/tau' * m'/m - # ^^^^^^^ - # Bolin's 1/c1 + # dm_heavy = dm_total / Bo * m'/m for sd_id in range(multiplicity.shape[0]): - moles_light = signed_water_mass[sd_id] / molar_mass_light - mass_ratio_heavy_to_light = (moles_heavy[sd_id] * molar_mass_heavy) / ( - moles_light * molar_mass_light - ) - dm_heavy_approx = ( - dm_total[sd_id] / bolin_number * mass_ratio_heavy_to_light + mass_ratio_heavy_to_total = ( + moles_heavy[sd_id] * molar_mass_heavy + ) / signed_water_mass[sd_id] + dm_heavy = ( + dm_total[sd_id] / bolin_number[sd_id] * mass_ratio_heavy_to_total ) - moles_heavy[sd_id] += dm_heavy_approx / molar_mass_heavy - mass_of_dry_air = dry_air_density[cell_id[sd_id]] * cell_volume + moles_heavy[sd_id] += dm_heavy / molar_mass_heavy + mass_of_dry_air = ( + dry_air_density[cell_id[sd_id]] * cell_volume + ) # TODO: pass from outside ambient_isotope_mixing_ratio[cell_id[sd_id]] -= ( - dm_heavy_approx * multiplicity[sd_id] / mass_of_dry_air + dm_heavy * multiplicity[sd_id] / mass_of_dry_air ) return body diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 4faa2c3b34..b5badc0211 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -451,9 +451,7 @@ def isotopic_fractionation(self, heavy_isotopes: tuple): "18O": "M_1H2_18O", }[isotope], ), - dm_total=self.attributes[ - f"diffusional growth mass change" - ], # TODO: total vs. heavy + dm_total=self.attributes["diffusional growth mass change"], bolin_number=self.attributes[f"Bolin number for {isotope}"], signed_water_mass=self.attributes["signed water mass"], multiplicity=self.attributes["multiplicity"], diff --git a/tests/unit_tests/attributes/test_diffusional_growth_mass_change.py b/tests/unit_tests/attributes/test_diffusional_growth_mass_change.py new file mode 100644 index 0000000000..38f2cabb25 --- /dev/null +++ b/tests/unit_tests/attributes/test_diffusional_growth_mass_change.py @@ -0,0 +1,41 @@ +import numpy as np +from PySDM import Builder +from PySDM.dynamics import Condensation, AmbientThermodynamics +from PySDM.physics import si +from ..dynamics.test_vapour_deposition_on_ice import MoistBox + + +def test_diffusional_growth_mass_change(backend_instance): + # arrange + n_sd = 1 + dt = 1 * si.s + + builder = Builder( + backend=backend_instance, n_sd=n_sd, environment=MoistBox(dt=dt, dv=np.nan) + ) + builder.add_dynamic(AmbientThermodynamics()) + builder.add_dynamic(Condensation()) + builder.request_attribute("diffusional growth mass change") + particulator = builder.build( + attributes={ + "water mass": np.ones(n_sd) * si.ng, + "multiplicity": np.ones(n_sd), + "dry volume": (dry_volume := np.ones(n_sd) * si.nm**3), + "kappa times dry volume": np.ones(n_sd) * dry_volume * (kappa := 0.6), + } + ) + particulator.environment["rhod"] = 1 * si.kg / si.m**3 + particulator.environment["thd"] = 300 * si.K + particulator.environment["water_vapour_mixing_ratio"] = 10 * si.g / si.kg + + # act + particulator.run(steps=1) + + # assert + diffusional_growth_mass_change = particulator.attributes[ + "diffusional growth mass change" + ].get() + assert abs(diffusional_growth_mass_change) > 0 + + # test happens if coagulation added -< check if throws assertion + # diff --git a/tests/unit_tests/dynamics/test_isotopic_fractionation.py b/tests/unit_tests/dynamics/test_isotopic_fractionation.py index 914ce24a09..736e116391 100644 --- a/tests/unit_tests/dynamics/test_isotopic_fractionation.py +++ b/tests/unit_tests/dynamics/test_isotopic_fractionation.py @@ -125,7 +125,9 @@ def test_no_isotope_fractionation_if_droplet_size_unchanged(backend_class): # act particulator.attributes["diffusional growth mass change"].data[:] = 0 - particulator.dynamics["IsotopicFractionation"]() + particulator.dynamics[ + "IsotopicFractionation" + ]() # TODO: call condensation as well! # assert assert ( diff --git a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py index 5533f9427f..6b1e337c8b 100644 --- a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py +++ b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py @@ -24,6 +24,7 @@ class MoistBox(Box, Moist): def __init__(self, dt: float, dv: float, mixed_phase: bool = False): Box.__init__(self, dt, dv) Moist.__init__(self, dt, self.mesh, variables=["rhod"], mixed_phase=mixed_phase) + self.dv = self.mesh.dv def register(self, builder): Moist.register(self, builder) From 3ec2fbcb863d4d68ac5bd1afdf946a9b1561c84b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Wed, 13 Aug 2025 13:02:30 +0200 Subject: [PATCH 10/23] fix bugs; add todos for next tests --- .../physics/diffusional_growth_mass_change.py | 10 ++++++++-- PySDM/particulator.py | 10 +++++++++- .../dynamics/test_isotopic_fractionation.py | 12 +++++++++++- 3 files changed, 28 insertions(+), 4 deletions(-) diff --git a/PySDM/attributes/physics/diffusional_growth_mass_change.py b/PySDM/attributes/physics/diffusional_growth_mass_change.py index 63100a4b26..dce3c856a5 100644 --- a/PySDM/attributes/physics/diffusional_growth_mass_change.py +++ b/PySDM/attributes/physics/diffusional_growth_mass_change.py @@ -4,14 +4,20 @@ @register_attribute() class DiffusionalGrowthMassChange(DerivedAttribute): def __init__(self, builder): - self.water_mass = builder.get_attribute("water mass") + self.water_mass = builder.get_attribute("signed water mass") super().__init__( builder, name="diffusional growth mass change", dependencies=(self.water_mass,), ) - builder.particulator.observers.append(self) assert "Collision" not in builder.particulator.dynamics + for triggers in ( + builder.particulator.observers, + builder.particulator.initialisers, + ): + triggers.append(self) + + def setup(self): self.notify() def notify(self): diff --git a/PySDM/particulator.py b/PySDM/particulator.py index b5badc0211..bca1adc170 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -29,8 +29,9 @@ def __init__(self, n_sd, backend): self.dynamics = {} self.products = {} self.observers = [] + self.initialisers = [] - self.n_steps = 0 + self.n_steps = -1 self.sorting_scheme = "default" self.condensation_solver = None @@ -48,6 +49,9 @@ def __init__(self, n_sd, backend): self.null = self.Storage.empty(0, dtype=float) def run(self, steps): + if self.n_steps == -1: + self._notify_initialisers() + self.n_steps = 0 for _ in range(steps): for key, dynamic in self.dynamics.items(): with self.timers[key]: @@ -60,6 +64,10 @@ def _notify_observers(self): for observer in reversed_order_so_that_environment_is_last: observer.notify() + def _notify_initialisers(self): + for initialiser in self.initialisers: + initialiser.setup() + @property def Storage(self): return self.backend.Storage diff --git a/tests/unit_tests/dynamics/test_isotopic_fractionation.py b/tests/unit_tests/dynamics/test_isotopic_fractionation.py index 736e116391..e070879cce 100644 --- a/tests/unit_tests/dynamics/test_isotopic_fractionation.py +++ b/tests/unit_tests/dynamics/test_isotopic_fractionation.py @@ -101,7 +101,7 @@ def test_no_isotope_fractionation_if_droplet_size_unchanged(backend_class): attributes = DUMMY_ATTRIBUTES.copy() attributes["moles_2H"] = particle_initial_isotope_content - attributes["water mass"] = 1 * si.ng + attributes["signed water mass"] = 1 * si.ng attributes["multiplicity"] = np.ones(1) builder = Builder( @@ -137,3 +137,13 @@ def test_no_isotope_fractionation_if_droplet_size_unchanged(backend_class): assert ( particulator.attributes["moles_2H"][0] == particle_initial_isotope_content ) + + # TODO + @staticmethod + def test_scenario_from_gedzelman_fig_X(): + pass + + # TODO + @staticmethod + def test_heavy_isotope_changes_match_bolin_number(): + pass From 4f2e0e8445ec111c5cfcfba5a810b7d2236b1d23 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Mon, 18 Aug 2025 13:50:54 +0200 Subject: [PATCH 11/23] add test_scenario_from_gedzelman_fig_2 draft --- .../dynamics/test_isotopic_fractionation.py | 58 ++++++++++++++++++- 1 file changed, 55 insertions(+), 3 deletions(-) diff --git a/tests/unit_tests/dynamics/test_isotopic_fractionation.py b/tests/unit_tests/dynamics/test_isotopic_fractionation.py index e070879cce..b97563a149 100644 --- a/tests/unit_tests/dynamics/test_isotopic_fractionation.py +++ b/tests/unit_tests/dynamics/test_isotopic_fractionation.py @@ -12,7 +12,8 @@ from PySDM.dynamics.isotopic_fractionation import HEAVY_ISOTOPES from PySDM.environments import Box from PySDM.physics import si - +from PySDM.physics.constants_defaults import R_str +from examples.PySDM_examples.Fisher_1991.fig_2 import temperature DUMMY_ATTRIBUTES = { attr: np.asarray([np.nan if attr != "multiplicity" else 0]) @@ -139,9 +140,60 @@ def test_no_isotope_fractionation_if_droplet_size_unchanged(backend_class): ) # TODO + R_SMOW = Formulae().constants.VSMOW_R_2H + @staticmethod - def test_scenario_from_gedzelman_fig_X(): - pass + @pytest.mark.parametrize( + "RH, R_rain, sign_of_dR_vap, sign_of_dR_rain", + ( + (0.5, 0.86 * R_SMOW, -1, 1), + (0.5, 0.9 * R_SMOW, 1, 1), + (0.5, 0.98 * R_SMOW, 1, -1), + ), + ) + def test_scenario_from_gedzelman_fig_2( + RH, R_rain, sign_of_dR_vap, sign_of_dR_rain, backend_class + ): + # arrange + particle_initial_isotope_content = 666.0 * si.moles + cell_volume = 1 * si.m**3 + mass_dry_air = 1 * si.kg + temperature = trivia.C2K(10) * si.K + vap_isotopic_ratio = trivia.isotopic_delta_2_ratio(-200 * PER_MILLE, R_SMOW) + vap_moles_light = pressure * cell_volume * M_1H2_16O / R_str / temperature + vap_mass_heavy_isotope = vap_isotopic_ratio * vap_moles_light * M_2H + + ambient_initial_isotope_mixing_ratio = vap_mass_heavy_isotope / mass_dry_air + + attributes = DUMMY_ATTRIBUTES.copy() + attributes["moles_2H"] = particle_initial_isotope_content + attributes["signed water mass"] = 1 * si.ng + attributes["multiplicity"] = np.ones(1) + + builder = Builder( + n_sd=1, + backend=backend_class( + formulae=Formulae(isotope_relaxation_timescale="MiyakeEtAl1968"), + ), + environment=Box(dv=cell_volume, dt=-1 * si.s), + ) + builder.add_dynamic(Condensation()) + builder.add_dynamic(IsotopicFractionation(isotopes=("2H",))) + particulator = builder.build(attributes=attributes, products=()) + particulator.environment["RH"] = RH + particulator.environment["dry_air_density"] = mass_dry_air / cell_volume + particulator.environment["mixing_ratio_2H"] = ( + ambient_initial_isotope_mixing_ratio + ) + + # act + particulator.dynamics["IsotopicFractionation"]() + dR_vap = new_vap_isotope_ratio - vap_isotope_ratio + dR_rain = new_R_rain - R_rain + + # assert + assert np.sign(dR_vap) == sign_of_dR_vap + assert np.sign(dR_rain) == sign_of_dR_rain # TODO @staticmethod From 19f4ebc866f406c867758bb71dd70813e90d32af Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Mon, 18 Aug 2025 19:22:33 +0200 Subject: [PATCH 12/23] clean imports --- .../dynamics/test_isotopic_fractionation.py | 21 ++++++++++--------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/tests/unit_tests/dynamics/test_isotopic_fractionation.py b/tests/unit_tests/dynamics/test_isotopic_fractionation.py index b97563a149..1b6f926168 100644 --- a/tests/unit_tests/dynamics/test_isotopic_fractionation.py +++ b/tests/unit_tests/dynamics/test_isotopic_fractionation.py @@ -12,8 +12,7 @@ from PySDM.dynamics.isotopic_fractionation import HEAVY_ISOTOPES from PySDM.environments import Box from PySDM.physics import si -from PySDM.physics.constants_defaults import R_str -from examples.PySDM_examples.Fisher_1991.fig_2 import temperature +from PySDM.physics.constants_defaults import R_str, PER_MILLE DUMMY_ATTRIBUTES = { attr: np.asarray([np.nan if attr != "multiplicity" else 0]) @@ -26,6 +25,9 @@ ) } +# TODO +R_SMOW = Formulae().constants.VSMOW_R_2H + class TestIsotopicFractionation: @staticmethod @@ -139,9 +141,6 @@ def test_no_isotope_fractionation_if_droplet_size_unchanged(backend_class): particulator.attributes["moles_2H"][0] == particle_initial_isotope_content ) - # TODO - R_SMOW = Formulae().constants.VSMOW_R_2H - @staticmethod @pytest.mark.parametrize( "RH, R_rain, sign_of_dR_vap, sign_of_dR_rain", @@ -155,11 +154,15 @@ def test_scenario_from_gedzelman_fig_2( RH, R_rain, sign_of_dR_vap, sign_of_dR_rain, backend_class ): # arrange + formulae = Formulae(isotope_relaxation_timescale="MiyakeEtAl1968") + particle_initial_isotope_content = 666.0 * si.moles cell_volume = 1 * si.m**3 mass_dry_air = 1 * si.kg - temperature = trivia.C2K(10) * si.K - vap_isotopic_ratio = trivia.isotopic_delta_2_ratio(-200 * PER_MILLE, R_SMOW) + temperature = formulae.trivia.C2K(10) * si.K + vap_isotopic_ratio = formulae.trivia.isotopic_delta_2_ratio( + -200 * PER_MILLE, R_SMOW + ) vap_moles_light = pressure * cell_volume * M_1H2_16O / R_str / temperature vap_mass_heavy_isotope = vap_isotopic_ratio * vap_moles_light * M_2H @@ -172,9 +175,7 @@ def test_scenario_from_gedzelman_fig_2( builder = Builder( n_sd=1, - backend=backend_class( - formulae=Formulae(isotope_relaxation_timescale="MiyakeEtAl1968"), - ), + backend=backend_class(formulae=formulae), environment=Box(dv=cell_volume, dt=-1 * si.s), ) builder.add_dynamic(Condensation()) From fb675d9d2ad2e217e26850dfdd790158092d0770 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Tue, 19 Aug 2025 22:49:31 +0200 Subject: [PATCH 13/23] populate test --- .../dynamics/test_isotopic_fractionation.py | 41 +++++++++++-------- 1 file changed, 25 insertions(+), 16 deletions(-) diff --git a/tests/unit_tests/dynamics/test_isotopic_fractionation.py b/tests/unit_tests/dynamics/test_isotopic_fractionation.py index 1b6f926168..8ac43b93a6 100644 --- a/tests/unit_tests/dynamics/test_isotopic_fractionation.py +++ b/tests/unit_tests/dynamics/test_isotopic_fractionation.py @@ -12,13 +12,13 @@ from PySDM.dynamics.isotopic_fractionation import HEAVY_ISOTOPES from PySDM.environments import Box from PySDM.physics import si -from PySDM.physics.constants_defaults import R_str, PER_MILLE +from PySDM.physics.constants_defaults import PER_MILLE DUMMY_ATTRIBUTES = { attr: np.asarray([np.nan if attr != "multiplicity" else 0]) for attr in ( "multiplicity", - "water mass", + "signed water mass", "dry volume", "kappa times dry volume", *[f"moles_{isotope}" for isotope in HEAVY_ISOTOPES], @@ -155,21 +155,21 @@ def test_scenario_from_gedzelman_fig_2( ): # arrange formulae = Formulae(isotope_relaxation_timescale="MiyakeEtAl1968") + const = formulae.constants - particle_initial_isotope_content = 666.0 * si.moles + n_2H_liq = 666.0 * si.moles + n_liq = n_2H_liq / R_rain cell_volume = 1 * si.m**3 mass_dry_air = 1 * si.kg temperature = formulae.trivia.C2K(10) * si.K - vap_isotopic_ratio = formulae.trivia.isotopic_delta_2_ratio( - -200 * PER_MILLE, R_SMOW - ) - vap_moles_light = pressure * cell_volume * M_1H2_16O / R_str / temperature - vap_mass_heavy_isotope = vap_isotopic_ratio * vap_moles_light * M_2H - - ambient_initial_isotope_mixing_ratio = vap_mass_heavy_isotope / mass_dry_air + R_vap = formulae.trivia.isotopic_delta_2_ratio(-200 * PER_MILLE, R_SMOW) + e = RH * formulae.saturation_vapour_pressure.pvs_water(temperature) + n_vap = (e + const.p_STP) * cell_volume / const.R_str / temperature + mass_2H_vap = R_vap * n_vap * const.M_2H_1H_16O + r_v = mass_2H_vap / mass_dry_air attributes = DUMMY_ATTRIBUTES.copy() - attributes["moles_2H"] = particle_initial_isotope_content + attributes["moles_2H"] = n_2H_liq attributes["signed water mass"] = 1 * si.ng attributes["multiplicity"] = np.ones(1) @@ -183,14 +183,23 @@ def test_scenario_from_gedzelman_fig_2( particulator = builder.build(attributes=attributes, products=()) particulator.environment["RH"] = RH particulator.environment["dry_air_density"] = mass_dry_air / cell_volume - particulator.environment["mixing_ratio_2H"] = ( - ambient_initial_isotope_mixing_ratio - ) + particulator.environment["mixing_ratio_2H"] = r_v # act + dn_liq_dt = -1 + dm_dt = dn_liq_dt * const.M_1H2_16O + particulator.attributes["diffusional growth mass change"].data[:] = dm_dt particulator.dynamics["IsotopicFractionation"]() - dR_vap = new_vap_isotope_ratio - vap_isotope_ratio - dR_rain = new_R_rain - R_rain + + new_n_liq = dn_liq_dt + n_liq + new_r_v = particulator.environment["mixing_ratio_2H"][0] + new_n_2H_liq = particulator.attributes["moles_2H"][0] + + new_R_vap = new_r_v * mass_dry_air / n_vap / const.M_2H_1H_16O # n_vap?? + new_R_liq = new_n_2H_liq / new_n_liq + + dR_vap = new_R_vap - R_vap + dR_rain = new_R_liq - R_rain # assert assert np.sign(dR_vap) == sign_of_dR_vap From e8bf6052b3f4b2e60e0ce4a7c041f72a5107d84b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Wed, 20 Aug 2025 12:47:29 +0200 Subject: [PATCH 14/23] update test gedzelman --- .../dynamics/test_isotopic_fractionation.py | 74 ++++++++++++++----- 1 file changed, 56 insertions(+), 18 deletions(-) diff --git a/tests/unit_tests/dynamics/test_isotopic_fractionation.py b/tests/unit_tests/dynamics/test_isotopic_fractionation.py index 8ac43b93a6..291e1deb15 100644 --- a/tests/unit_tests/dynamics/test_isotopic_fractionation.py +++ b/tests/unit_tests/dynamics/test_isotopic_fractionation.py @@ -150,28 +150,42 @@ def test_no_isotope_fractionation_if_droplet_size_unchanged(backend_class): (0.5, 0.98 * R_SMOW, 1, -1), ), ) - def test_scenario_from_gedzelman_fig_2( - RH, R_rain, sign_of_dR_vap, sign_of_dR_rain, backend_class + def test_scenario_from_gedzelman_fig_2_for_single_superdroplet_and_single_isotope( + RH, + R_rain, + sign_of_dR_vap, + sign_of_dR_rain, + backend_class, + # TODO: multiplicity, dt, ... ): # arrange formulae = Formulae(isotope_relaxation_timescale="MiyakeEtAl1968") const = formulae.constants - n_2H_liq = 666.0 * si.moles - n_liq = n_2H_liq / R_rain + multiplicity = np.ones(1) + m_t = formulae.trivia.volume(radius=1 * si.mm) * const.rho_w + cell_volume = 1 * si.m**3 mass_dry_air = 1 * si.kg + d_mixing_ratio_env = 0.01 * si.g / si.kg temperature = formulae.trivia.C2K(10) * si.K R_vap = formulae.trivia.isotopic_delta_2_ratio(-200 * PER_MILLE, R_SMOW) + e = RH * formulae.saturation_vapour_pressure.pvs_water(temperature) - n_vap = (e + const.p_STP) * cell_volume / const.R_str / temperature - mass_2H_vap = R_vap * n_vap * const.M_2H_1H_16O + n_vap_total = e * cell_volume / const.R_str / temperature + n_vap_heavy = n_vap_total * R_vap / (1 + R_vap) + mass_2H_vap = n_vap_heavy * const.M_2H r_v = mass_2H_vap / mass_dry_air attributes = DUMMY_ATTRIBUTES.copy() - attributes["moles_2H"] = n_2H_liq + attributes["moles_2H"] = ( + m_t / (1 + const.M_1H / (R_rain * const.M_2H)) / const.M_2H + ) + for isotope in HEAVY_ISOTOPES: + if isotope != "2H": + attributes[f"moles_{isotope}"] = 0 attributes["signed water mass"] = 1 * si.ng - attributes["multiplicity"] = np.ones(1) + attributes["multiplicity"] = multiplicity builder = Builder( n_sd=1, @@ -180,27 +194,51 @@ def test_scenario_from_gedzelman_fig_2( ) builder.add_dynamic(Condensation()) builder.add_dynamic(IsotopicFractionation(isotopes=("2H",))) + builder.request_attribute("delta_2H") particulator = builder.build(attributes=attributes, products=()) particulator.environment["RH"] = RH particulator.environment["dry_air_density"] = mass_dry_air / cell_volume particulator.environment["mixing_ratio_2H"] = r_v + # sanity check + assert ( + formulae.trivia.isotopic_delta_2_ratio( + particulator.attributes["delta_2H"][0], R_SMOW + ) + == R_rain + ) + # act - dn_liq_dt = -1 - dm_dt = dn_liq_dt * const.M_1H2_16O - particulator.attributes["diffusional growth mass change"].data[:] = dm_dt - particulator.dynamics["IsotopicFractionation"]() - new_n_liq = dn_liq_dt + n_liq - new_r_v = particulator.environment["mixing_ratio_2H"][0] - new_n_2H_liq = particulator.attributes["moles_2H"][0] + droplet_dm = -d_mixing_ratio_env * mass_dry_air / multiplicity + particulator.attributes["diffusional growth mass change"].data[:] = droplet_dm - new_R_vap = new_r_v * mass_dry_air / n_vap / const.M_2H_1H_16O # n_vap?? - new_R_liq = new_n_2H_liq / new_n_liq + print(particulator.environment["mixing_ratio_2H"][0]) + print(particulator.attributes["delta_2H"][0]) + particulator.dynamics["IsotopicFractionation"]() + + print(particulator.environment["mixing_ratio_2H"][0]) + print(particulator.attributes["delta_2H"][0]) + + # new_n_liq = dn_liq_dt + n_liq + new_r2H_v = particulator.environment["mixing_ratio_2H"][0] + new_delta_2H = particulator.attributes["delta_2H"][0] + # new_n_2H_liq = particulator.attributes["moles_2H"][0] + # + new_n_heavy_vap = new_r2H_v * mass_dry_air / const.M_2H + molar_mass_vap_new = np.nan # f(mixing_ratio_total, mixing_ratio_2H) + new_n_total_vap = ( + n_vap_total + d_mixing_ratio_env * mass_dry_air / molar_mass_vap_new + ) # or old? + new_n_light_vap = new_n_total_vap - new_n_heavy_vap + new_R_vap = new_n_heavy_vap / new_n_light_vap + new_R_liq = formulae.trivia.isotopic_delta_2_ratio(new_delta_2H, R_SMOW) + # dR_vap = new_R_vap - R_vap dR_rain = new_R_liq - R_rain - + print(dR_rain, dR_vap) + assert False # assert assert np.sign(dR_vap) == sign_of_dR_vap assert np.sign(dR_rain) == sign_of_dR_rain From 7ee99f3635097c47a5ecfea1582d3c328c7ad1ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Fri, 19 Sep 2025 11:58:13 +0200 Subject: [PATCH 15/23] update Gedzelman test; add moles_heavy_atom to trivia --- PySDM/physics/trivia.py | 23 ++++++++ .../dynamics/test_isotopic_fractionation.py | 53 ++++++++++--------- tests/unit_tests/physics/test_trivia.py | 47 +++++++++++++++- 3 files changed, 96 insertions(+), 27 deletions(-) diff --git a/PySDM/physics/trivia.py b/PySDM/physics/trivia.py index 93f53ce6d8..0d1e6bde1a 100644 --- a/PySDM/physics/trivia.py +++ b/PySDM/physics/trivia.py @@ -130,6 +130,29 @@ def isotopic_enrichment_to_delta_SMOW(E, delta_0_SMOW): """ return (E + 1) * (delta_0_SMOW + 1) - 1 + @staticmethod + def moles_heavy_atom( + const, + delta, + mass_total, + molar_mass_heavy_molecule, + R_STD, + light_atoms_per_light_molecule, + ): + return mass_total / ( + ( + 1 + + const.M_1H2_16O + / ( + light_atoms_per_light_molecule + * (delta + 1) + * R_STD + * molar_mass_heavy_molecule + ) + ) + * molar_mass_heavy_molecule + ) + @staticmethod def mixing_ratio_to_specific_content(mixing_ratio): return mixing_ratio / (1 + mixing_ratio) diff --git a/tests/unit_tests/dynamics/test_isotopic_fractionation.py b/tests/unit_tests/dynamics/test_isotopic_fractionation.py index 291e1deb15..84e141997b 100644 --- a/tests/unit_tests/dynamics/test_isotopic_fractionation.py +++ b/tests/unit_tests/dynamics/test_isotopic_fractionation.py @@ -25,9 +25,6 @@ ) } -# TODO -R_SMOW = Formulae().constants.VSMOW_R_2H - class TestIsotopicFractionation: @staticmethod @@ -143,16 +140,16 @@ def test_no_isotope_fractionation_if_droplet_size_unchanged(backend_class): @staticmethod @pytest.mark.parametrize( - "RH, R_rain, sign_of_dR_vap, sign_of_dR_rain", + "RH, delta_rain, sign_of_dR_vap, sign_of_dR_rain", ( - (0.5, 0.86 * R_SMOW, -1, 1), - (0.5, 0.9 * R_SMOW, 1, 1), - (0.5, 0.98 * R_SMOW, 1, -1), + (0.5, 0.86 - 1, -1, 1), + (0.5, 0.9 - 1, 1, 1), + (0.5, 0.98 - 1, 1, -1), ), ) def test_scenario_from_gedzelman_fig_2_for_single_superdroplet_and_single_isotope( RH, - R_rain, + delta_rain, sign_of_dR_vap, sign_of_dR_rain, backend_class, @@ -163,28 +160,32 @@ def test_scenario_from_gedzelman_fig_2_for_single_superdroplet_and_single_isotop const = formulae.constants multiplicity = np.ones(1) - m_t = formulae.trivia.volume(radius=1 * si.mm) * const.rho_w - + m_t = 1 * si.ng cell_volume = 1 * si.m**3 mass_dry_air = 1 * si.kg d_mixing_ratio_env = 0.01 * si.g / si.kg temperature = formulae.trivia.C2K(10) * si.K - R_vap = formulae.trivia.isotopic_delta_2_ratio(-200 * PER_MILLE, R_SMOW) + R_vap = formulae.trivia.isotopic_delta_2_ratio( + -200 * PER_MILLE, const.VSMOW_R_2H + ) e = RH * formulae.saturation_vapour_pressure.pvs_water(temperature) n_vap_total = e * cell_volume / const.R_str / temperature n_vap_heavy = n_vap_total * R_vap / (1 + R_vap) mass_2H_vap = n_vap_heavy * const.M_2H - r_v = mass_2H_vap / mass_dry_air attributes = DUMMY_ATTRIBUTES.copy() - attributes["moles_2H"] = ( - m_t / (1 + const.M_1H / (R_rain * const.M_2H)) / const.M_2H + attributes["moles_2H"] = formulae.trivia.moles_heavy_atom( + delta=delta_rain, + mass_total=m_t, + molar_mass_heavy_molecule=const.M_2H_1H_16O, + R_STD=const.VSMOW_R_2H, + light_atoms_per_light_molecule=2, ) for isotope in HEAVY_ISOTOPES: if isotope != "2H": attributes[f"moles_{isotope}"] = 0 - attributes["signed water mass"] = 1 * si.ng + attributes["signed water mass"] = m_t attributes["multiplicity"] = multiplicity builder = Builder( @@ -198,14 +199,11 @@ def test_scenario_from_gedzelman_fig_2_for_single_superdroplet_and_single_isotop particulator = builder.build(attributes=attributes, products=()) particulator.environment["RH"] = RH particulator.environment["dry_air_density"] = mass_dry_air / cell_volume - particulator.environment["mixing_ratio_2H"] = r_v + particulator.environment["mixing_ratio_2H"] = mass_2H_vap / mass_dry_air - # sanity check - assert ( - formulae.trivia.isotopic_delta_2_ratio( - particulator.attributes["delta_2H"][0], R_SMOW - ) - == R_rain + # sanity check for initial condition + np.testing.assert_approx_equal( + particulator.attributes["delta_2H"][0], delta_rain, significant=3 ) # act @@ -233,12 +231,15 @@ def test_scenario_from_gedzelman_fig_2_for_single_superdroplet_and_single_isotop ) # or old? new_n_light_vap = new_n_total_vap - new_n_heavy_vap new_R_vap = new_n_heavy_vap / new_n_light_vap - new_R_liq = formulae.trivia.isotopic_delta_2_ratio(new_delta_2H, R_SMOW) + new_R_liq = formulae.trivia.isotopic_delta_2_ratio( + new_delta_2H, const.VSMOW_R_2H + ) # dR_vap = new_R_vap - R_vap - dR_rain = new_R_liq - R_rain - print(dR_rain, dR_vap) - assert False + dR_rain = new_R_liq - formulae.trivia.isotopic_delta_2_ratio( + delta_rain, const.VSMOW_R_2H + ) + # assert assert np.sign(dR_vap) == sign_of_dR_vap assert np.sign(dR_rain) == sign_of_dR_rain diff --git a/tests/unit_tests/physics/test_trivia.py b/tests/unit_tests/physics/test_trivia.py index b3157e00f0..239df1de25 100644 --- a/tests/unit_tests/physics/test_trivia.py +++ b/tests/unit_tests/physics/test_trivia.py @@ -5,7 +5,7 @@ from PySDM import Formulae from PySDM.physics.dimensional_analysis import DimensionalAnalysis -from PySDM.physics import constants_defaults +from PySDM.physics import constants_defaults, si class TestTrivia: @@ -113,3 +113,48 @@ def test_celsius_to_kelvin(): # assert assert temperature_in_kelvin == temperature_in_celsius + 273.15 + + @staticmethod + @pytest.mark.parametrize("delta", (0.86 - 1, 0.9 - 1, 0.98 - 1)) + def test_moles_heavy_atom( + backend_class, + delta, + m_t=1 * si.ng, + ): + # arrange + formulae = Formulae() + const = formulae.constants + + from PySDM.dynamics.isotopic_fractionation import HEAVY_ISOTOPES + from PySDM import Builder + from PySDM.environments import Box + + attributes = { + "multiplicity": np.asarray([0]), + "signed water mass": m_t, + "moles_2H": formulae.trivia.moles_heavy_atom( + delta=delta, + mass_total=m_t, + molar_mass_heavy_molecule=const.M_2H_1H_16O, + R_STD=const.VSMOW_R_2H, + light_atoms_per_light_molecule=2, + ), + } + for isotope in HEAVY_ISOTOPES: + if isotope != "2H": + attributes[f"moles_{isotope}"] = 0 + + builder = Builder( + n_sd=1, + backend=backend_class(formulae=formulae), + environment=Box(dv=np.nan, dt=-1 * si.s), + ) + builder.request_attribute("delta_2H") + particulator = builder.build(attributes=attributes, products=()) + # ratio = moles_2H / + # formulae.trivia.isotopic_delta_2_ratio(ratio, const.VSMOW_R_2H) + + # sanity check for initial condition + np.testing.assert_approx_equal( + particulator.attributes["delta_2H"][0], delta, significant=10 + ) From 9d7540752edc50f89bd04d0335c04674268d530b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Fri, 19 Sep 2025 12:34:00 +0200 Subject: [PATCH 16/23] add test for atribute moles light water --- tests/unit_tests/attributes/test_isotopes.py | 36 ++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/tests/unit_tests/attributes/test_isotopes.py b/tests/unit_tests/attributes/test_isotopes.py index aaafe3c73d..a592515099 100644 --- a/tests/unit_tests/attributes/test_isotopes.py +++ b/tests/unit_tests/attributes/test_isotopes.py @@ -156,3 +156,39 @@ def test_bolin_number_attribute( # assert assert value == expected_tau + + @staticmethod + def test_moles( + backend_class, + m_t=1 * si.ng, + ): + # arrange + formulae = Formulae() + + from PySDM.dynamics.isotopic_fractionation import HEAVY_ISOTOPES + from PySDM import Builder + from PySDM.environments import Box + + attributes = { + "multiplicity": np.asarray([0]), + "signed water mass": m_t, + } + for isotope in HEAVY_ISOTOPES: + # if isotope != "2H": + attributes[f"moles_{isotope}"] = 44 + + builder = Builder( + n_sd=1, + backend=backend_class(formulae=formulae), + environment=Box(dv=np.nan, dt=-1 * si.s), + ) + builder.request_attribute("moles light water") + builder.request_attribute("moles_16O") + particulator = builder.build(attributes=attributes, products=()) + + # sanity check for initial condition + np.testing.assert_approx_equal( + particulator.attributes["moles light water"][0], + particulator.attributes["moles_16O"][0], + significant=10, + ) From 5a18a53336afa9b77e8bb6bbd4a615aa1f1eed7e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Fri, 19 Sep 2025 12:48:26 +0200 Subject: [PATCH 17/23] trivia test does not depend on backend class --- PySDM/attributes/isotopes/moles.py | 4 ++-- tests/unit_tests/attributes/test_isotopes.py | 5 ++-- tests/unit_tests/physics/test_trivia.py | 24 ++++++-------------- 3 files changed, 12 insertions(+), 21 deletions(-) diff --git a/PySDM/attributes/isotopes/moles.py b/PySDM/attributes/isotopes/moles.py index 4bf71cb938..f510a337de 100644 --- a/PySDM/attributes/isotopes/moles.py +++ b/PySDM/attributes/isotopes/moles.py @@ -81,8 +81,8 @@ def __init__(self, builder): builder=builder, name="moles_16O", attrs_to_multiplier={ - builder.get_attribute("moles_2H"): 0.5, - builder.get_attribute("moles_3H"): 0.5, + builder.get_attribute("moles_2H"): 1.0, + builder.get_attribute("moles_3H"): 1.0, builder.get_attribute("moles light water"): 1.0, }, ) diff --git a/tests/unit_tests/attributes/test_isotopes.py b/tests/unit_tests/attributes/test_isotopes.py index a592515099..246ee6d310 100644 --- a/tests/unit_tests/attributes/test_isotopes.py +++ b/tests/unit_tests/attributes/test_isotopes.py @@ -174,7 +174,6 @@ def test_moles( "signed water mass": m_t, } for isotope in HEAVY_ISOTOPES: - # if isotope != "2H": attributes[f"moles_{isotope}"] = 44 builder = Builder( @@ -189,6 +188,8 @@ def test_moles( # sanity check for initial condition np.testing.assert_approx_equal( particulator.attributes["moles light water"][0], - particulator.attributes["moles_16O"][0], + particulator.attributes["moles_16O"][0] + - particulator.attributes["moles_2H"][0] + - particulator.attributes["moles_3H"][0], significant=10, ) diff --git a/tests/unit_tests/physics/test_trivia.py b/tests/unit_tests/physics/test_trivia.py index 239df1de25..d85e92c119 100644 --- a/tests/unit_tests/physics/test_trivia.py +++ b/tests/unit_tests/physics/test_trivia.py @@ -117,7 +117,6 @@ def test_celsius_to_kelvin(): @staticmethod @pytest.mark.parametrize("delta", (0.86 - 1, 0.9 - 1, 0.98 - 1)) def test_moles_heavy_atom( - backend_class, delta, m_t=1 * si.ng, ): @@ -126,8 +125,6 @@ def test_moles_heavy_atom( const = formulae.constants from PySDM.dynamics.isotopic_fractionation import HEAVY_ISOTOPES - from PySDM import Builder - from PySDM.environments import Box attributes = { "multiplicity": np.asarray([0]), @@ -143,18 +140,11 @@ def test_moles_heavy_atom( for isotope in HEAVY_ISOTOPES: if isotope != "2H": attributes[f"moles_{isotope}"] = 0 - - builder = Builder( - n_sd=1, - backend=backend_class(formulae=formulae), - environment=Box(dv=np.nan, dt=-1 * si.s), - ) - builder.request_attribute("delta_2H") - particulator = builder.build(attributes=attributes, products=()) - # ratio = moles_2H / - # formulae.trivia.isotopic_delta_2_ratio(ratio, const.VSMOW_R_2H) - - # sanity check for initial condition - np.testing.assert_approx_equal( - particulator.attributes["delta_2H"][0], delta, significant=10 + attributes["moles light water"] = ( + attributes["signed water mass"] / const.M_1H2_16O ) + ratio = attributes["moles_2H"] / attributes["moles light water"] + sut = formulae.trivia.isotopic_ratio_2_delta(ratio, const.VSMOW_R_2H) + + # assert + np.testing.assert_approx_equal(actual=sut, desired=delta, significant=10) From d46ed4c64d271a9a91406546a94f9152279474d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Fri, 19 Sep 2025 14:10:22 +0200 Subject: [PATCH 18/23] update delta function --- PySDM/attributes/isotopes/delta.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/PySDM/attributes/isotopes/delta.py b/PySDM/attributes/isotopes/delta.py index 2baacc0256..b96544a2a2 100644 --- a/PySDM/attributes/isotopes/delta.py +++ b/PySDM/attributes/isotopes/delta.py @@ -10,15 +10,9 @@ class DeltaImpl(DerivedAttribute): def __init__(self, builder, *, heavy_isotope): assert heavy_isotope[:-1].isnumeric() - if heavy_isotope[-1] == "H": - light_isotope = "1H" - elif heavy_isotope[-1] == "O": - light_isotope = "16O" - else: - raise NotImplementedError() - - self.heavy_isotope = builder.get_attribute("moles_" + heavy_isotope) - self.light_isotope = builder.get_attribute(f"moles_{light_isotope}") + + self.heavy_isotope = builder.get_attribute(f"moles_{heavy_isotope}") + self.light_isotope = builder.get_attribute(f"moles light water") super().__init__( builder, name="delta_" + heavy_isotope, From 75c37cf5a622ffc4084cde642ea00ff720f7ac6e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Fri, 19 Sep 2025 14:27:51 +0200 Subject: [PATCH 19/23] add test for new trivia function; remove unused parameter --- PySDM/physics/trivia.py | 12 +-------- tests/unit_tests/physics/test_trivia.py | 36 +++++++++++-------------- 2 files changed, 17 insertions(+), 31 deletions(-) diff --git a/PySDM/physics/trivia.py b/PySDM/physics/trivia.py index 0d1e6bde1a..9a163a3aa3 100644 --- a/PySDM/physics/trivia.py +++ b/PySDM/physics/trivia.py @@ -137,19 +137,9 @@ def moles_heavy_atom( mass_total, molar_mass_heavy_molecule, R_STD, - light_atoms_per_light_molecule, ): return mass_total / ( - ( - 1 - + const.M_1H2_16O - / ( - light_atoms_per_light_molecule - * (delta + 1) - * R_STD - * molar_mass_heavy_molecule - ) - ) + (1 + const.M_1H2_16O / ((delta + 1) * R_STD * molar_mass_heavy_molecule)) * molar_mass_heavy_molecule ) diff --git a/tests/unit_tests/physics/test_trivia.py b/tests/unit_tests/physics/test_trivia.py index d85e92c119..ca0465e25a 100644 --- a/tests/unit_tests/physics/test_trivia.py +++ b/tests/unit_tests/physics/test_trivia.py @@ -116,7 +116,7 @@ def test_celsius_to_kelvin(): @staticmethod @pytest.mark.parametrize("delta", (0.86 - 1, 0.9 - 1, 0.98 - 1)) - def test_moles_heavy_atom( + def test_moles_2H_atom( delta, m_t=1 * si.ng, ): @@ -126,25 +126,21 @@ def test_moles_heavy_atom( from PySDM.dynamics.isotopic_fractionation import HEAVY_ISOTOPES - attributes = { - "multiplicity": np.asarray([0]), - "signed water mass": m_t, - "moles_2H": formulae.trivia.moles_heavy_atom( - delta=delta, - mass_total=m_t, - molar_mass_heavy_molecule=const.M_2H_1H_16O, - R_STD=const.VSMOW_R_2H, - light_atoms_per_light_molecule=2, - ), - } - for isotope in HEAVY_ISOTOPES: - if isotope != "2H": - attributes[f"moles_{isotope}"] = 0 - attributes["moles light water"] = ( - attributes["signed water mass"] / const.M_1H2_16O + signed_water_mass = m_t + moles_2H = formulae.trivia.moles_heavy_atom( + delta=delta, + mass_total=signed_water_mass, + molar_mass_heavy_molecule=const.M_2H_1H_16O, + R_STD=const.VSMOW_R_2H, + ) + moles_light_water = ( + signed_water_mass - const.M_2H_1H_16O * moles_2H + ) / const.M_1H2_16O + + # act + sut = formulae.trivia.isotopic_ratio_2_delta( + ratio=moles_2H / moles_light_water, reference_ratio=const.VSMOW_R_2H ) - ratio = attributes["moles_2H"] / attributes["moles light water"] - sut = formulae.trivia.isotopic_ratio_2_delta(ratio, const.VSMOW_R_2H) # assert - np.testing.assert_approx_equal(actual=sut, desired=delta, significant=10) + np.testing.assert_approx_equal(actual=sut, desired=delta, significant=5) From 7d4fbfb70ade31db4541330ec4afc1a4db6584ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Fri, 19 Sep 2025 14:28:16 +0200 Subject: [PATCH 20/23] reuse trivia fn --- tests/unit_tests/attributes/test_isotopes.py | 4 ++-- tests/unit_tests/dynamics/test_isotopic_fractionation.py | 7 ++++--- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/tests/unit_tests/attributes/test_isotopes.py b/tests/unit_tests/attributes/test_isotopes.py index 246ee6d310..bf7949b3e2 100644 --- a/tests/unit_tests/attributes/test_isotopes.py +++ b/tests/unit_tests/attributes/test_isotopes.py @@ -185,11 +185,11 @@ def test_moles( builder.request_attribute("moles_16O") particulator = builder.build(attributes=attributes, products=()) - # sanity check for initial condition + # assert np.testing.assert_approx_equal( particulator.attributes["moles light water"][0], particulator.attributes["moles_16O"][0] - particulator.attributes["moles_2H"][0] - particulator.attributes["moles_3H"][0], - significant=10, + significant=5, ) diff --git a/tests/unit_tests/dynamics/test_isotopic_fractionation.py b/tests/unit_tests/dynamics/test_isotopic_fractionation.py index 84e141997b..53fac14d1a 100644 --- a/tests/unit_tests/dynamics/test_isotopic_fractionation.py +++ b/tests/unit_tests/dynamics/test_isotopic_fractionation.py @@ -171,7 +171,9 @@ def test_scenario_from_gedzelman_fig_2_for_single_superdroplet_and_single_isotop e = RH * formulae.saturation_vapour_pressure.pvs_water(temperature) n_vap_total = e * cell_volume / const.R_str / temperature - n_vap_heavy = n_vap_total * R_vap / (1 + R_vap) + n_vap_heavy = n_vap_total * formulae.trivia.mixing_ratio_to_specific_content( + R_vap + ) mass_2H_vap = n_vap_heavy * const.M_2H attributes = DUMMY_ATTRIBUTES.copy() @@ -180,7 +182,6 @@ def test_scenario_from_gedzelman_fig_2_for_single_superdroplet_and_single_isotop mass_total=m_t, molar_mass_heavy_molecule=const.M_2H_1H_16O, R_STD=const.VSMOW_R_2H, - light_atoms_per_light_molecule=2, ) for isotope in HEAVY_ISOTOPES: if isotope != "2H": @@ -203,7 +204,7 @@ def test_scenario_from_gedzelman_fig_2_for_single_superdroplet_and_single_isotop # sanity check for initial condition np.testing.assert_approx_equal( - particulator.attributes["delta_2H"][0], delta_rain, significant=3 + particulator.attributes["delta_2H"][0], delta_rain, significant=5 ) # act From 80204cb06f7ad51fa23c9ce2c538787b3d1087f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Fri, 19 Sep 2025 14:51:50 +0200 Subject: [PATCH 21/23] clean test - remove unused calculations --- .../dynamics/test_isotopic_fractionation.py | 37 ++++--------------- 1 file changed, 7 insertions(+), 30 deletions(-) diff --git a/tests/unit_tests/dynamics/test_isotopic_fractionation.py b/tests/unit_tests/dynamics/test_isotopic_fractionation.py index 53fac14d1a..9d8fc7f7f2 100644 --- a/tests/unit_tests/dynamics/test_isotopic_fractionation.py +++ b/tests/unit_tests/dynamics/test_isotopic_fractionation.py @@ -208,43 +208,20 @@ def test_scenario_from_gedzelman_fig_2_for_single_superdroplet_and_single_isotop ) # act - droplet_dm = -d_mixing_ratio_env * mass_dry_air / multiplicity particulator.attributes["diffusional growth mass change"].data[:] = droplet_dm - - print(particulator.environment["mixing_ratio_2H"][0]) - print(particulator.attributes["delta_2H"][0]) - particulator.dynamics["IsotopicFractionation"]() - print(particulator.environment["mixing_ratio_2H"][0]) - print(particulator.attributes["delta_2H"][0]) - - # new_n_liq = dn_liq_dt + n_liq - new_r2H_v = particulator.environment["mixing_ratio_2H"][0] - new_delta_2H = particulator.attributes["delta_2H"][0] - # new_n_2H_liq = particulator.attributes["moles_2H"][0] - # - new_n_heavy_vap = new_r2H_v * mass_dry_air / const.M_2H - molar_mass_vap_new = np.nan # f(mixing_ratio_total, mixing_ratio_2H) - new_n_total_vap = ( - n_vap_total + d_mixing_ratio_env * mass_dry_air / molar_mass_vap_new - ) # or old? - new_n_light_vap = new_n_total_vap - new_n_heavy_vap - new_R_vap = new_n_heavy_vap / new_n_light_vap - new_R_liq = formulae.trivia.isotopic_delta_2_ratio( - new_delta_2H, const.VSMOW_R_2H + # assert + assert ( + np.sign(particulator.environment["mixing_ratio_2H"][0] - R_vap) + == sign_of_dR_vap ) - # - dR_vap = new_R_vap - R_vap - dR_rain = new_R_liq - formulae.trivia.isotopic_delta_2_ratio( - delta_rain, const.VSMOW_R_2H + assert ( + np.sign(particulator.attributes["delta_2H"][0] - delta_rain) + == sign_of_dR_rain ) - # assert - assert np.sign(dR_vap) == sign_of_dR_vap - assert np.sign(dR_rain) == sign_of_dR_rain - # TODO @staticmethod def test_heavy_isotope_changes_match_bolin_number(): From 86e8a21c0105ed1c893bd3d05563a87b241d3175 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Fri, 19 Sep 2025 17:59:02 +0200 Subject: [PATCH 22/23] revert testing changes; --- PySDM/attributes/isotopes/delta.py | 8 +++++++- PySDM/attributes/isotopes/moles.py | 10 +++++----- PySDM/physics/trivia.py | 12 +++++++++++- .../dynamics/test_isotopic_fractionation.py | 18 +++++++----------- tests/unit_tests/physics/test_trivia.py | 1 + 5 files changed, 31 insertions(+), 18 deletions(-) diff --git a/PySDM/attributes/isotopes/delta.py b/PySDM/attributes/isotopes/delta.py index b96544a2a2..df8c19ca79 100644 --- a/PySDM/attributes/isotopes/delta.py +++ b/PySDM/attributes/isotopes/delta.py @@ -10,9 +10,15 @@ class DeltaImpl(DerivedAttribute): def __init__(self, builder, *, heavy_isotope): assert heavy_isotope[:-1].isnumeric() + if heavy_isotope[-1] == "H": + light_isotope = "1H" + elif heavy_isotope[-1] == "O": + light_isotope = "16O" + else: + raise NotImplementedError() self.heavy_isotope = builder.get_attribute(f"moles_{heavy_isotope}") - self.light_isotope = builder.get_attribute(f"moles light water") + self.light_isotope = builder.get_attribute(f"moles_{light_isotope}") super().__init__( builder, name="delta_" + heavy_isotope, diff --git a/PySDM/attributes/isotopes/moles.py b/PySDM/attributes/isotopes/moles.py index f510a337de..1b4a812886 100644 --- a/PySDM/attributes/isotopes/moles.py +++ b/PySDM/attributes/isotopes/moles.py @@ -64,11 +64,11 @@ def __init__(self, builder): super().__init__( builder=builder, name="moles_1H", - attrs_to_multiplier={ - builder.get_attribute("moles_17O"): 2.0, - builder.get_attribute("moles_18O"): 2.0, - builder.get_attribute("moles_2H"): 1.0, - builder.get_attribute("moles_3H"): 1.0, + attrs_to_multiplier={ # TODO + # builder.get_attribute("moles_17O"): 2.0, + # builder.get_attribute("moles_18O"): 2.0, + # builder.get_attribute("moles_2H"): 1.0, + # builder.get_attribute("moles_3H"): 1.0, builder.get_attribute("moles light water"): 2.0, }, ) diff --git a/PySDM/physics/trivia.py b/PySDM/physics/trivia.py index 9a163a3aa3..0d1e6bde1a 100644 --- a/PySDM/physics/trivia.py +++ b/PySDM/physics/trivia.py @@ -137,9 +137,19 @@ def moles_heavy_atom( mass_total, molar_mass_heavy_molecule, R_STD, + light_atoms_per_light_molecule, ): return mass_total / ( - (1 + const.M_1H2_16O / ((delta + 1) * R_STD * molar_mass_heavy_molecule)) + ( + 1 + + const.M_1H2_16O + / ( + light_atoms_per_light_molecule + * (delta + 1) + * R_STD + * molar_mass_heavy_molecule + ) + ) * molar_mass_heavy_molecule ) diff --git a/tests/unit_tests/dynamics/test_isotopic_fractionation.py b/tests/unit_tests/dynamics/test_isotopic_fractionation.py index 9d8fc7f7f2..e2797c2a05 100644 --- a/tests/unit_tests/dynamics/test_isotopic_fractionation.py +++ b/tests/unit_tests/dynamics/test_isotopic_fractionation.py @@ -165,9 +165,7 @@ def test_scenario_from_gedzelman_fig_2_for_single_superdroplet_and_single_isotop mass_dry_air = 1 * si.kg d_mixing_ratio_env = 0.01 * si.g / si.kg temperature = formulae.trivia.C2K(10) * si.K - R_vap = formulae.trivia.isotopic_delta_2_ratio( - -200 * PER_MILLE, const.VSMOW_R_2H - ) + R_vap = formulae.trivia.isotopic_delta_2_ratio(-200, const.VSMOW_R_2H) e = RH * formulae.saturation_vapour_pressure.pvs_water(temperature) n_vap_total = e * cell_volume / const.R_str / temperature @@ -182,6 +180,7 @@ def test_scenario_from_gedzelman_fig_2_for_single_superdroplet_and_single_isotop mass_total=m_t, molar_mass_heavy_molecule=const.M_2H_1H_16O, R_STD=const.VSMOW_R_2H, + light_atoms_per_light_molecule=2, ) for isotope in HEAVY_ISOTOPES: if isotope != "2H": @@ -212,15 +211,12 @@ def test_scenario_from_gedzelman_fig_2_for_single_superdroplet_and_single_isotop particulator.attributes["diffusional growth mass change"].data[:] = droplet_dm particulator.dynamics["IsotopicFractionation"]() + new_R_vap = particulator.environment["mixing_ratio_2H"][0] + new_delta_rain = particulator.attributes["delta_2H"][0] + # assert - assert ( - np.sign(particulator.environment["mixing_ratio_2H"][0] - R_vap) - == sign_of_dR_vap - ) - assert ( - np.sign(particulator.attributes["delta_2H"][0] - delta_rain) - == sign_of_dR_rain - ) + assert np.sign(new_R_vap - R_vap) == sign_of_dR_vap + assert np.sign(new_delta_rain - delta_rain) == sign_of_dR_rain # TODO @staticmethod diff --git a/tests/unit_tests/physics/test_trivia.py b/tests/unit_tests/physics/test_trivia.py index ca0465e25a..00f4d19d94 100644 --- a/tests/unit_tests/physics/test_trivia.py +++ b/tests/unit_tests/physics/test_trivia.py @@ -132,6 +132,7 @@ def test_moles_2H_atom( mass_total=signed_water_mass, molar_mass_heavy_molecule=const.M_2H_1H_16O, R_STD=const.VSMOW_R_2H, + light_atoms_per_light_molecule=2, ) moles_light_water = ( signed_water_mass - const.M_2H_1H_16O * moles_2H From db20fe9e6c1da1ad028a6cdc98546da915939ec6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Sat, 20 Sep 2025 16:25:31 +0200 Subject: [PATCH 23/23] test for moles heavy in trivia --- tests/unit_tests/physics/test_trivia.py | 43 +++++++++++++++++-------- 1 file changed, 29 insertions(+), 14 deletions(-) diff --git a/tests/unit_tests/physics/test_trivia.py b/tests/unit_tests/physics/test_trivia.py index 00f4d19d94..181cbfe39c 100644 --- a/tests/unit_tests/physics/test_trivia.py +++ b/tests/unit_tests/physics/test_trivia.py @@ -116,32 +116,47 @@ def test_celsius_to_kelvin(): @staticmethod @pytest.mark.parametrize("delta", (0.86 - 1, 0.9 - 1, 0.98 - 1)) - def test_moles_2H_atom( + @pytest.mark.parametrize("water_mass", np.linspace(10**-2, 10**3, 9) * si.ng) + @pytest.mark.parametrize( + "heavy_isotope_name, heavy_isotope_molecule", + (("2H", "2H_1H_16O"), ("17O", "1H2_17O"), ("18O", "1H2_17O")), + ) + def test_moles_heavy_atom( delta, - m_t=1 * si.ng, + water_mass, + heavy_isotope_name, + heavy_isotope_molecule, ): # arrange formulae = Formulae() const = formulae.constants - from PySDM.dynamics.isotopic_fractionation import HEAVY_ISOTOPES + molar_mass_heavy_molecule = getattr(const, f"M_{heavy_isotope_molecule}") + R_STD = getattr(const, f"VSMOW_R_{heavy_isotope_name}") + if heavy_isotope_name[-1] == "O": + light_atoms_per_light_molecule = 1 + elif heavy_isotope_name[-1] == "H": + light_atoms_per_light_molecule = 2 + else: + light_atoms_per_light_molecule = 0 - signed_water_mass = m_t - moles_2H = formulae.trivia.moles_heavy_atom( + moles_atom = formulae.trivia.moles_heavy_atom( delta=delta, - mass_total=signed_water_mass, - molar_mass_heavy_molecule=const.M_2H_1H_16O, - R_STD=const.VSMOW_R_2H, - light_atoms_per_light_molecule=2, + mass_total=water_mass, + molar_mass_heavy_molecule=molar_mass_heavy_molecule, + R_STD=R_STD, + light_atoms_per_light_molecule=light_atoms_per_light_molecule, ) moles_light_water = ( - signed_water_mass - const.M_2H_1H_16O * moles_2H - ) / const.M_1H2_16O + moles_atom + / light_atoms_per_light_molecule + / formulae.trivia.isotopic_delta_2_ratio(delta=delta, reference_ratio=R_STD) + ) # act - sut = formulae.trivia.isotopic_ratio_2_delta( - ratio=moles_2H / moles_light_water, reference_ratio=const.VSMOW_R_2H + sut = ( + moles_atom * molar_mass_heavy_molecule + moles_light_water * const.M_1H2_16O ) # assert - np.testing.assert_approx_equal(actual=sut, desired=delta, significant=5) + np.testing.assert_approx_equal(actual=sut, desired=water_mass, significant=5)