diff --git a/PySDM/attributes/isotopes/__init__.py b/PySDM/attributes/isotopes/__init__.py index 26875f6fe3..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 +from ..isotopes import delta, bolin_numbers diff --git a/PySDM/attributes/isotopes/bolin_numbers.py b/PySDM/attributes/isotopes/bolin_numbers.py new file mode 100644 index 0000000000..9246052751 --- /dev/null +++ b/PySDM/attributes/isotopes/bolin_numbers.py @@ -0,0 +1,41 @@ +""" +# TODO: define Bolin 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 BolinNumberImpl(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 number for " + heavy_isotope, + dependencies=(self.moles_heavy_isotope,), + ) + + def recalculate(self): + self.particulator.bolin_number( + output=self.data, + molar_mass=self.molar_mass, + moles_heavy_isotope=self.moles_heavy_isotope.data, + ) + + +def make_bolin_number_factory(heavy_isotope: str): + def _factory(builder): + return BolinNumberImpl(builder, heavy_isotope=heavy_isotope) + + return _factory + + +for heavy_isotope in HEAVY_ISOTOPES: + 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 new file mode 100644 index 0000000000..9246052751 --- /dev/null +++ b/PySDM/attributes/isotopes/bolins_numbers.py @@ -0,0 +1,41 @@ +""" +# TODO: define Bolin 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 BolinNumberImpl(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 number for " + heavy_isotope, + dependencies=(self.moles_heavy_isotope,), + ) + + def recalculate(self): + self.particulator.bolin_number( + output=self.data, + molar_mass=self.molar_mass, + moles_heavy_isotope=self.moles_heavy_isotope.data, + ) + + +def make_bolin_number_factory(heavy_isotope: str): + def _factory(builder): + return BolinNumberImpl(builder, heavy_isotope=heavy_isotope) + + return _factory + + +for heavy_isotope in HEAVY_ISOTOPES: + register_attribute(name=f"Bolin number for {heavy_isotope}")( + make_bolin_number_factory(heavy_isotope) + ) 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/__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..63100a4b26 --- /dev/null +++ b/PySDM/attributes/physics/diffusional_growth_mass_change.py @@ -0,0 +1,21 @@ +from PySDM.attributes.impl import register_attribute, DerivedAttribute + + +@register_attribute() +class DiffusionalGrowthMassChange(DerivedAttribute): + def __init__(self, builder): + self.water_mass = builder.get_attribute("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 + self.notify() + + def notify(self): + self.data[:] = -self.water_mass.data + + def recalculate(self): + 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 f33e6761dc..27981151be 100644 --- a/PySDM/backends/impl_numba/methods/isotope_methods.py +++ b/PySDM/backends/impl_numba/methods/isotope_methods.py @@ -24,5 +24,98 @@ 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, "parallel": False}) + def body( + *, + multiplicity, + dm_total, + bolin_number, + signed_water_mass, + moles_heavy, + molar_mass_heavy, + ambient_isotope_mixing_ratio, + cell_id, + cell_volume, + dry_air_density, + ): + # 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 / Bo * m'/m + + for sd_id in range(multiplicity.shape[0]): + 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 / 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 * multiplicity[sd_id] / mass_of_dry_air + ) + + return body + + def isotopic_fractionation( + self, + *, + multiplicity, + dm_total, + bolin_number, + signed_water_mass, + moles_heavy, + molar_mass_heavy, + ambient_isotope_mixing_ratio, + cell_id, + cell_volume, + dry_air_density, + ): + self._isotopic_fractionation_body( + multiplicity=multiplicity.data, + 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, + dry_air_density=dry_air_density.data, + ) + + @cached_property + 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__bolin_number( + moles_heavy_isotope=moles_heavy_isotope[i], + relative_humidity=relative_humidity[cell_id[i]], + molar_mass=molar_mass, + ) + + return body + + def bolin_number( + self, *, output, molar_mass, cell_id, moles_heavy_isotope, relative_humidity + ): + self._bolin_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/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 32ae0697d5..b5badc0211 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -440,10 +440,40 @@ def calculate_displacement( ) def isotopic_fractionation(self, heavy_isotopes: tuple): - self.backend.isotopic_fractionation() for isotope in heavy_isotopes: + self.backend.isotopic_fractionation( + 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["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"], + moles_heavy=self.attributes[f"moles_{isotope}"], + 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[ + f"mixing_ratio_{isotope}" + ], + ) self.attributes.mark_updated(f"moles_{isotope}") + 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"], + 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 2d0913897d..e164884184 100644 --- a/PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py +++ b/PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py @@ -22,3 +22,9 @@ 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 bolin_number( + const, moles_heavy_isotope, relative_humidity, molar_mass + ): # pylint: disable=unused-argument + return 44.0 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..345f4560b3 --- /dev/null +++ b/examples/PySDM_examples/Strzabala_2026/kinematic_2d_isotopes.ipynb @@ -0,0 +1,184 @@ +{ + "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", + "id": "initial_id", + "metadata": { + "collapsed": true, + "ExecuteTime": { + "end_time": "2025-05-27T14:36:02.233992Z", + "start_time": "2025-05-27T14:35:58.798955Z" + } + }, + "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, 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": { + "ExecuteTime": { + "end_time": "2025-05-27T14:36:10.341002Z", + "start_time": "2025-05-27T14:36:02.342762Z" + } + }, + "cell_type": "code", + "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", + "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": { + "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..23e0a92673 100644 --- a/examples/PySDM_examples/utils/kinematic_2d/simulation.py +++ b/examples/PySDM_examples/utils/kinematic_2d/simulation.py @@ -1,4 +1,5 @@ import numpy as np + from PySDM_examples.utils.kinematic_2d.make_default_product_collection import ( make_default_product_collection, ) @@ -32,7 +33,11 @@ def __init__(self, settings, storage, SpinUp, backend_class=CPU): def products(self): return self.particulator.products - def reinit(self, products=None): + def reinit( # pylint: disable=too-many-locals, too-many-branches, too-many-statements + self, + products=None, + additional_advectees_initial_profiles: dict[str, np.ndarray] = None, + ): formulae = self.settings.formulae backend = self.backend_class(formulae=formulae) environment = Kinematic2D( @@ -70,6 +75,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( ( 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/attributes/test_isotopes.py b/tests/unit_tests/attributes/test_isotopes.py index 6cb3c2dd66..aaafe3c73d 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_bolin_number_attribute( + backend_class, + heavy_isotope: str, + moles_heavy_isotope: float, + relative_humidity: float, + expected_tau: float, + variant: str, + ): # pylint: disable=too-many-arguments + 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 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..736e116391 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 @@ -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 ( @@ -91,3 +90,50 @@ 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" + ]() # TODO: call condensation as well! + + # assert + assert ( + particulator.environment["mixing_ratio_2H"][0] + == ambient_initial_isotope_mixing_ratio + ) + assert ( + particulator.attributes["moles_2H"][0] == particle_initial_isotope_content + ) 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)