Add isotopes to simulation kinematic 2d

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

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)
+    )

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)
+    )

PySDM/attributes/isotopes/delta.py

@@ -17,7 +17,7 @@ def __init__(self, builder, *, heavy_isotope):
         else:
             raise NotImplementedError()
 
-        self.heavy_isotope = builder.get_attribute("moles_" + heavy_isotope)
+        self.heavy_isotope = builder.get_attribute(f"moles_{heavy_isotope}")
         self.light_isotope = builder.get_attribute(f"moles_{light_isotope}")
         super().__init__(
             builder,

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,
             },
         )
 
@@ -73,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,
             },
         )
@@ -90,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,
             },
         )

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

PySDM/attributes/physics/diffusional_growth_mass_change.py

@@ -0,0 +1,27 @@
+from PySDM.attributes.impl import register_attribute, DerivedAttribute
+
+
+@register_attribute()
+class DiffusionalGrowthMassChange(DerivedAttribute):
+    def __init__(self, builder):
+        self.water_mass = builder.get_attribute("signed water mass")
+        super().__init__(
+            builder,
+            name="diffusional growth mass change",
+            dependencies=(self.water_mass,),
+        )
+        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):
+        self.data[:] = -self.water_mass.data
+
+    def recalculate(self):
+        self.data[:] += self.water_mass.data

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,
+        )

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)

PySDM/particulator.py

@@ -30,8 +30,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
@@ -49,6 +50,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]:
@@ -61,6 +65,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
@@ -441,10 +449,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,
         *,

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