fix Mv constant calc (D and T abundances count twice in water molecular abundances); share initial parameters across J&N and Yang examples (#1258)

Author: slayoo

PySDM/physics/constants_defaults.py

@@ -38,6 +38,7 @@
 )
 
 # https://web.archive.org/web/20200729203147/https://nucleus.iaea.org/rpst/documents/VSMOW_SLAP.pdf
+# heavy-to-light isotope abundance ratios
 VSMOW_R_2H = 155.76 * PPM
 VSMOW_R_3H = 1.85e-11 * PPM
 VSMOW_R_18O = 2005.20 * PPM
@@ -53,19 +54,6 @@
 M_17O = 16.99913175651 * si.g / si.mole
 M_18O = 17.99915961287 * si.g / si.mole
 
-Mv = (
-    1
-    - Trivia.mixing_ratio_to_specific_content(
-        VSMOW_R_2H / 2 + VSMOW_R_3H / 2 + VSMOW_R_17O + VSMOW_R_18O
-    )
-) * (
-    (M_1H * 2 + M_16O)
-    + (M_2H + M_1H + M_16O) * VSMOW_R_2H
-    + (M_3H + M_1H + M_16O) * VSMOW_R_3H
-    + (M_1H * 2 + M_17O) * VSMOW_R_17O
-    + (M_1H * 2 + M_18O) * VSMOW_R_18O
-)
-
 R_str = sci.R * si.joule / si.kelvin / si.mole
 N_A = sci.N_A / si.mole
 
@@ -356,6 +344,39 @@
 
 
 def compute_derived_values(c: dict):
+    """
+    computes derived quantities such as molar mass ratios, etc.
+
+    water molar mass is computed from from molecular masses and VSMOW isotope abundances
+    (and neglecting molecular binding energies)
+    for discussion, see:
+    - caption of Table 2.1 in [Gat 2010](https://doi.org/10.1142/p027)
+    - [IAPWS Guidelines](http://www.iapws.org/relguide/fundam.pdf)
+    """
+
+    c["Mv"] = (
+        (
+            1
+            - 2 * Trivia.mixing_ratio_to_specific_content(c["VSMOW_R_2H"])
+            - 2 * Trivia.mixing_ratio_to_specific_content(c["VSMOW_R_3H"])
+            - 1 * Trivia.mixing_ratio_to_specific_content(c["VSMOW_R_17O"])
+            - 1 * Trivia.mixing_ratio_to_specific_content(c["VSMOW_R_18O"])
+        )
+        * (c["M_1H"] * 2 + c["M_16O"])
+        + 2
+        * Trivia.mixing_ratio_to_specific_content(c["VSMOW_R_2H"])
+        * (c["M_2H"] + c["M_1H"] + c["M_16O"])
+        + 2
+        * Trivia.mixing_ratio_to_specific_content(c["VSMOW_R_3H"])
+        * (c["M_3H"] + c["M_1H"] + c["M_16O"])
+        + 1
+        * Trivia.mixing_ratio_to_specific_content(c["VSMOW_R_17O"])
+        * (c["M_1H"] * 2 + c["M_17O"])
+        + 1
+        * Trivia.mixing_ratio_to_specific_content(c["VSMOW_R_18O"])
+        * (c["M_1H"] * 2 + c["M_18O"])
+    )
+
     c["eps"] = c["Mv"] / c["Md"]
     c["Rd"] = c["R_str"] / c["Md"]
     c["Rv"] = c["R_str"] / c["Mv"]

examples/PySDM_examples/Jensen_and_Nugent_2017/settings.py

@@ -3,15 +3,19 @@
 from PySDM.physics import si
 from PySDM.initialisation.spectra import Lognormal, Sum
 
+INITIAL_RELATIVE_HUMIDITY = 0.8561
+INITIAL_TEMPERATURE = 284.3 * si.K
+INITIAL_PRESSURE = 938.5 * si.hPa
+INITIAL_ALTITUDE = 600 * si.metres
+
 
 @strict
 class Settings:
     def __init__(self, *, aerosol: str, cloud_type: str):
-        # TODO #1266 reuse these values in the Yang et al. '18 example which is based on J&N'16
-        self.p0 = 938.5 * si.hPa
-        self.RH0 = 0.8561
-        self.T0 = 284.3 * si.K
-        self.z0 = 600 * si.m
+        self.p0 = INITIAL_PRESSURE
+        self.RH0 = INITIAL_RELATIVE_HUMIDITY
+        self.T0 = INITIAL_TEMPERATURE
+        self.z0 = INITIAL_ALTITUDE
         self.t_end_of_ascent = 1500 * si.s if cloud_type == "Sc" else None
         self.dt = 1 * si.s  # TODO #1266: not found in the paper yet
 

examples/PySDM_examples/Yang_et_al_2018/fig_2.ipynb

@@ -1,6 +1,8 @@
 import numpy as np
 from pystrict import strict
 
+from PySDM_examples import Jensen_and_Nugent_2017
+
 from PySDM.backends import CPU
 from PySDM.dynamics import condensation
 from PySDM.initialisation import spectra
@@ -47,10 +49,10 @@ def __init__(
         self.rtol_thd = condensation.DEFAULTS.rtol_thd
         self.dt_cond_range = condensation.DEFAULTS.cond_range
 
-        self.T0 = 284.3 * si.kelvin
-        self.initial_water_vapour_mixing_ratio = 7.6 * si.grams / si.kilogram
-        self.p0 = 938.5 * si.hectopascals
-        self.z0 = 600 * si.metres
+        self.T0 = Jensen_and_Nugent_2017.settings.INITIAL_TEMPERATURE
+        self.RH0 = Jensen_and_Nugent_2017.settings.INITIAL_RELATIVE_HUMIDITY
+        self.p0 = Jensen_and_Nugent_2017.settings.INITIAL_PRESSURE
+        self.z0 = Jensen_and_Nugent_2017.settings.INITIAL_ALTITUDE
         self.kappa = 0.53  # Petters and S. M. Kreidenweis mean growth-factor derived
 
         self.t0 = 1200 * si.second

examples/PySDM_examples/Yang_et_al_2018/settings.py

@@ -26,7 +26,15 @@ def __init__(self, settings, backend=CPU):
             dt=dt_output / self.n_substeps,
             mass_of_dry_air=settings.mass_of_dry_air,
             p0=settings.p0,
-            initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio,
+            initial_water_vapour_mixing_ratio=self.formulae.constants.eps
+            / (
+                settings.p0
+                / settings.RH0
+                / self.formulae.saturation_vapour_pressure.pvs_Celsius(
+                    settings.T0 - self.formulae.constants.T0
+                )
+                - 1
+            ),
             T0=settings.T0,
             w=settings.w,
             z0=settings.z0,

examples/PySDM_examples/Yang_et_al_2018/simulation.py

@@ -80,4 +80,4 @@ def test_relative_humidity():
     simulation = Simulation(settings)
 
     # Assert
-    assert round(simulation.particulator.environment["RH"][0], 3) == 0.859
+    assert round(simulation.particulator.environment["RH"][0], 3) == 0.856

tests/smoke_tests/parcel/yang_et_al_2018/test_initialisation.py

@@ -4,8 +4,8 @@
 from chempy import Substance
 
 # from PySDM.initialisation import spectra
+from PySDM import Formulae
 from PySDM.initialisation.aerosol_composition import DryAerosolMixture
-from PySDM.physics import constants_defaults as const
 from PySDM.physics import si
 
 
@@ -19,6 +19,7 @@
 )
 def test_volume_weighted_kappa_with_insoluble_compound(mass_fractions):
     # Arrange
+    const = Formulae().constants
     water_molar_volume = const.Mv / const.rho_w
     compounds = ("(NH4)2SO4", "insoluble")
     molar_masses = {

tests/unit_tests/initialisation/test_aerosol_init.py

@@ -118,9 +118,9 @@ def test_isotope_molar_masses_vs_wikipedia_or_scipy(item, value):
     @staticmethod
     def test_isotope_molar_masses_vsmow_vs_mean_water_molar_mass():
         np.testing.assert_approx_equal(
-            desired=constants_defaults.Mv,
-            actual=Substance.from_formula("H2O").mass * si.gram / si.mole,
-            significant=4.5,
+            actual=Formulae().constants.Mv,
+            desired=Substance.from_formula("H2O").mass * si.gram / si.mole,
+            significant=5.5,
         )
 
     @staticmethod