other isotopes in isotope methods; cache D_ratio and alpha in bolin_number; add diffusivity ratios using Grahams law for other isotopes; tests for Bolin number (#1834)

Co-authored-by: Sylwester Arabas <sylwester.arabas@agh.edu.pl>

Author: AgnieszkaZaba

PySDM/attributes/isotopes/bolin_numbers.py

@@ -1,6 +1,6 @@
 """
 Derived droplet attribute representing the Bolin number for heavy
-water isotopologues (D, 17O, 18O substitutions).
+water isotopologues (2H, 17O, 18O substitutions).
 
 The Bolin number is a dimensionless coefficient relating the
 heavy-isotope mole tendency to the bulk liquid-water mass tendency
@@ -24,6 +24,7 @@ def __init__(self, builder, *, heavy_isotope: str):
         heavy_isotope : str
             Heavy isotopologue identifier (entry of ``HEAVY_ISOTOPES``).
         """
+        self.isotope = heavy_isotope
         self.moles_heavy = builder.get_attribute(f"moles_{heavy_isotope}")
         self.moles_light = builder.get_attribute("moles light water")  # TODO #1787
         self.cell_id = builder.get_attribute("cell id")
@@ -41,9 +42,10 @@ def recalculate(self):
         self.particulator.backend.bolin_number(
             output=self.data,
             cell_id=self.cell_id.data,
+            isotope=self.isotope,
             relative_humidity=self.particulator.environment["RH"],
             temperature=self.particulator.environment["T"],
-            density_dry_air=self.particulator.environment["dry_air_density"],
+            density_dry_air=self.particulator.environment["rhod"],
             moles_light_molecule=self.moles_light.data,
             moles_heavy=self.moles_heavy.data,
             molality_in_dry_air=self.molality_in_dry_air,

PySDM/backends/impl_numba/methods/isotope_methods.py

@@ -2,7 +2,7 @@
 CPU implementation of isotope-relates backend methods
 """
 
-from functools import cached_property
+from functools import cached_property, lru_cache
 
 import numba
 
@@ -76,9 +76,7 @@ def body(
                     dm_heavy = dm_total[sd_id] / Bo * mass_ratio_heavy_to_total
                 dn_heavy_molecule = dm_heavy / molar_mass_heavy_molecule
                 moles_heavy_molecule[sd_id] += dn_heavy_molecule
-                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
                 molality_in_dry_air[cell_id[sd_id]] -= (
                     dn_heavy_molecule * multiplicity[sd_id] / mass_of_dry_air
                 )
@@ -102,7 +100,7 @@ def isotopic_fractionation(
         """Update heavy-isotope composition during droplet growth/evaporation."""
         self._isotopic_fractionation_body(
             cell_id=cell_id.data,
-            cell_volume=cell_volume.data,
+            cell_volume=cell_volume,
             multiplicity=multiplicity.data,
             dm_total=dm_total.data,
             signed_water_mass=signed_water_mass.data,
@@ -133,6 +131,8 @@ def body(
             *,
             output,
             cell_id,
+            alpha,
+            D_ratio,
             relative_humidity,
             temperature,
             density_dry_air,
@@ -148,38 +148,51 @@ def body(
                 conc_vap_total = (
                     pvs_water * relative_humidity[cell_id[i]] / ff.constants.R_str / T
                 )
-                rho_v = pvs_water / T / ff.constants.Rv
-
                 isotopic_fraction = ff.trivia__isotopic_fraction(
                     molality_in_dry_air=molality_in_dry_air[cell_id[i]],
                     density_dry_air=density_dry_air[cell_id[i]],
                     total_vap_concentration=conc_vap_total,
                 )
-                D_ratio_heavy_to_light = (
-                    ff.isotope_diffusivity_ratios__ratio_2H_heavy_to_light(T)
-                )
                 output[i] = ff.isotope_relaxation_timescale__bolin_number(
-                    D_ratio_heavy_to_light=D_ratio_heavy_to_light,
-                    alpha=ff.isotope_equilibrium_fractionation_factors__alpha_l_2H(T),
-                    D_light=ff.constants.D0,
+                    D_ratio_heavy_to_light=D_ratio(T),
+                    alpha=alpha(T),
                     Fk=ff.drop_growth__Fk(
                         T=T, K=ff.constants.K0, lv=ff.constants.l_tri
                     ),
+                    Fd=ff.drop_growth__Fd(
+                        T=T,
+                        D=ff.constants.D0,
+                        pvs=pvs_water,
+                    ),
                     R_vap=ff.trivia__isotopic_ratio_assuming_single_heavy_isotope(
                         isotopic_fraction
                     ),
                     R_liq=moles_heavy_atom / moles_light_isotope,
                     relative_humidity=relative_humidity[cell_id[i]],
-                    rho_v=rho_v,
                 )
 
         return body
 
+    @lru_cache
+    def alpha_l(self, isotope):
+        return getattr(
+            self.formulae_flattened,
+            f"isotope_equilibrium_fractionation_factors__alpha_l_{isotope}",
+        )
+
+    @lru_cache
+    def D_ratio(self, isotope):
+        return getattr(
+            self.formulae_flattened,
+            f"isotope_diffusivity_ratios__ratio_{isotope}_heavy_to_light",
+        )
+
     def bolin_number(
         self,
         *,
         output,
         cell_id,
+        isotope,
         relative_humidity,
         temperature,
         density_dry_air,
@@ -191,6 +204,8 @@ def bolin_number(
         self._bolin_number_body(
             output=output.data,
             cell_id=cell_id.data,
+            D_ratio=self.D_ratio(isotope),
+            alpha=self.alpha_l(isotope),
             relative_humidity=relative_humidity.data,
             temperature=temperature.data,
             density_dry_air=density_dry_air.data,

PySDM/backends/impl_thrust_rtc/methods/isotope_methods.py

@@ -63,7 +63,7 @@ def __isotopic_fractionation(self):
 
             moles_heavy_molecule[i] += dn_heavy;
 
-            real_type mass_of_dry_air = dry_air_density[cid] * cell_volume[cid];
+            real_type mass_of_dry_air = dry_air_density[cid] * cell_volume;
 
             atomicAdd((real_type*) &molality_in_dry_air[cid], -dn_heavy * multiplicity[i] / mass_of_dry_air);
             """.replace("real_type", self._get_c_type()),
@@ -87,7 +87,7 @@ def isotopic_fractionation(
             n=len(multiplicity),
             args=(
                 cell_id.data,
-                cell_volume.data,
+                self._get_floating_point(cell_volume),
                 multiplicity.data,
                 dm_total.data,
                 signed_water_mass.data,

PySDM/dynamics/isotopic_fractionation.py

@@ -41,8 +41,10 @@ def register(self, builder):
                 raise AssertionError(
                     f"Isotopic fractionation not implemented for {isotope}"
                 )
-            builder.request_attribute(f"moles_{isotope}")
             builder.request_attribute(f"Bolin number for {isotope}")
 
+        for isotope in HEAVY_ISOTOPES:
+            builder.request_attribute(f"moles_{isotope}")
+
     def __call__(self):
         self.particulator.isotopic_fractionation(self.isotopes)

PySDM/particulator.py

@@ -456,7 +456,7 @@ def isotopic_fractionation(self, heavy_isotopes: tuple):
                 multiplicity=self.attributes["multiplicity"],
                 dm_total=self.attributes["diffusional growth mass change"],
                 signed_water_mass=self.attributes["signed water mass"],
-                dry_air_density=self.environment["dry_air_density"],
+                dry_air_density=self.environment["rhod"],
                 molar_mass_heavy_molecule=getattr(
                     self.formulae.constants,
                     {

PySDM/physics/drop_growth/howell_1949.py

@@ -4,7 +4,7 @@
 same as in [Mason 1951](https://doi.org/10.1088/0370-1301/64/9/307)
 
 The notation for terms associated with heat conduction and diffusion are from eq. 7.17
- in [Rogers & Yau 1971](https://archive.org/details/shortcourseinclo0000roge_m3k2).
+ in [Rogers & Yau 1989](https://archive.org/details/shortcourseinclo0000roge_m3k2).
 """
 
 from .fick import Fick

PySDM/physics/isotope_diffusivity_ratios/grahams_law.py

@@ -18,3 +18,15 @@ def ratio_3H_heavy_to_light(const, temperature):  # pylint: disable=unused-argum
         return (
             (2 * const.M_1H + const.M_16O) / (const.M_3H + const.M_1H + const.M_16O)
         ) ** const.ONE_HALF
+
+    @staticmethod
+    def ratio_17O_heavy_to_light(const, temperature):  # pylint: disable=unused-argument
+        return (
+            (2 * const.M_1H + const.M_16O) / (2 * const.M_1H + const.M_17O)
+        ) ** const.ONE_HALF
+
+    @staticmethod
+    def ratio_18O_heavy_to_light(const, temperature):  # pylint: disable=unused-argument
+        return (
+            (2 * const.M_1H + const.M_16O) / (2 * const.M_1H + const.M_18O)
+        ) ** const.ONE_HALF

PySDM/physics/isotope_relaxation_timescale/zaba_et_al.py

@@ -22,15 +22,13 @@ def tau(
 
     @staticmethod
     def bolin_number(
-        const,
+        relative_humidity,
         D_ratio_heavy_to_light,
         alpha,
-        D_light,
-        Fk,
         R_vap,
         R_liq,
-        relative_humidity,
-        rho_v,
+        Fd,
+        Fk,
     ):  # pylint: disable=too-many-arguments,too-many-positional-arguments
         # TODO #1809 Numba can't compile when * in bolin_number
         """Heavy to total isotopic-timescales ratio (tau_heavy/tau_total).
@@ -43,9 +41,7 @@ def bolin_number(
             * (1 - relative_humidity)
             / D_ratio_heavy_to_light
             / (
-                (1 + rho_v * D_light * Fk / const.rho_w)
-                * relative_humidity
-                * (1 - alpha * R_vap / R_liq)
+                (1 + Fk / Fd) * relative_humidity * (1 - alpha * R_vap / R_liq)
                 - relative_humidity
                 + 1
             )

docs/bibliography.json

@@ -757,7 +757,8 @@
       "PySDM/physics/drop_growth/howell_1949.py",
       "PySDM/physics/drop_growth/mason_1971.py",
       "PySDM/dynamics/terminal_velocity/rogers_and_yau.py",
-      "PySDM/dynamics/terminal_velocity/power_series.py"
+      "PySDM/dynamics/terminal_velocity/power_series.py",
+      "examples/PySDM_examples/Gedzelman_and_Arnold_1994/fig_2.ipynb"
     ],
     "title": "A short course in clouds physics",
     "label": "Rogers & Yau 1989 (Pergamon Press)"