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 1/9] 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 2/9] 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 3/9] 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 4/9] 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 5/9] 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 6/9] 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 7/9] 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 8/9] 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 9/9] 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)