Homogeneous freezing example + FIX for mark updated attributes in homogeneous freezing and thaw + expanded unit test of record freezing temperature (#1682)

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

Author: 3 people

.github/workflows/tests.yml

@@ -185,7 +185,7 @@ jobs:
             test-suite: "tests/smoke_tests/parcel_d"
       fail-fast: true
     runs-on: ${{ matrix.platform }}
-    timeout-minutes: ${{ startsWith(matrix.platform, 'windows-') && 35 || 22 }}
+    timeout-minutes: ${{ startsWith(matrix.platform, 'windows-') && 35 || 25 }}
     steps:
       - uses: actions/checkout@v4.1.6
         with:
@@ -212,7 +212,8 @@ jobs:
         python-version: ["3.12"]
         test-suite:
           - "isotopes_chemistry_extraterrestrial"
-          - "coagulation_freezing"
+          - "homogeneous_freezing"
+          - "coagulation_immersion_freezing"
           - "condensation_a"
           - "condensation_b"
           - "condensation_c"
@@ -237,7 +238,12 @@ jobs:
             test-suite: "multi-process_d"
           - platform: "macos-14"
             test-suite: "multi-process_e"
-      fail-fast: true
+          # TODO #1870
+          - platform: "ubuntu-24.04"
+            test-suite: "homogeneous_freezing"
+          - platform: "windows-latest"
+            test-suite: "homogeneous_freezing"
+      fail-fast: false
     runs-on: ${{ matrix.platform }}
     timeout-minutes: ${{ startsWith(matrix.platform, 'windows-') && 35 || 30 }}
     steps:
@@ -288,7 +294,7 @@ jobs:
           python -m pytest ${{ env.pytest_options }} tests/examples_tests/test_tests_completeness.py
       - run: |
           python -m pip install pytest-timeout
-          python -m pytest ${{ env.pytest_options }} --basetemp=/tmp/pytest --timeout_method=thread --timeout=${{ startsWith(matrix.platform, 'windows-') && 1000 || 900 }} tests/examples_tests/test_run* --suite ${{ matrix.test-suite }}
+          python -m pytest ${{ env.pytest_options }} --basetemp=/tmp/pytest --timeout_method=thread --timeout=${{ startsWith(matrix.platform, 'windows-') && 1200 || 900 }} tests/examples_tests/test_run* --suite ${{ matrix.test-suite }}
 
       - if: ( ! startsWith(matrix.platform, 'windows-') ) && contains(matrix.test-suite, env.anim_test-suite) && matrix.python-version == env.anim_python-version
         run: |

PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py

@@ -83,7 +83,7 @@ def _condensation(
         RH_max=RH_max.data,
         success=success.data,
     )
-    particulator.attributes.mark_updated("water mass")
+    particulator.attributes.mark_updated("signed water mass")
 
 
 @lru_cache()

PySDM/particulator.py

@@ -573,6 +573,7 @@ def homogeneous_freezing_time_dependent(self, *, rand: Storage):
             temperature=self.environment["T"],
             relative_humidity_ice=self.environment["RH_ice"],
         )
+        self.attributes.mark_updated("signed water mass")
 
     def homogeneous_freezing_threshold(self):
         self.backend.homogeneous_freezing_threshold(
@@ -583,6 +584,7 @@ def homogeneous_freezing_threshold(self):
             temperature=self.environment["T"],
             relative_humidity_ice=self.environment["RH_ice"],
         )
+        self.attributes.mark_updated("signed water mass")
 
     def thaw_instantaneous(self):
         self.backend.thaw_instantaneous(
@@ -592,3 +594,4 @@ def thaw_instantaneous(self):
             cell=self.attributes["cell id"],
             temperature=self.environment["T"],
         )
+        self.attributes.mark_updated("signed water mass")

PySDM/physics/constants_defaults.py

@@ -356,6 +356,16 @@
 KOOP_MURRAY_C6 = -7.46953e-7 / si.K**6
 """ 〃 """
 
+KOOP_MURRAY_DAW_C0 = -273.111686
+""" homogeneous ice nucleation rate for pure water droplets adapted
+for DAW ([Spichtinger et al. 2023](https://doi.org/10.5194/acp-23-2035-2023)) """
+KOOP_MURRAY_DAW_C1 = 1689.25134
+""" 〃 """
+KOOP_MURRAY_DAW_C2 = -2449.65042
+""" 〃 """
+KOOP_MURRAY_DAW_UNIT = 1 / si.m**3 / si.s
+""" 〃 """
+
 J_HET = np.nan
 J_HOM = np.nan
 """ constant ice nucleation rates """

PySDM/physics/homogeneous_ice_nucleation_rate/__init__.py

@@ -7,3 +7,4 @@
 from .koop import Koop2000
 from .koop_corr import Koop_Correction
 from .koop_murray import KoopMurray2016
+from .koop_murray_dwa_adaption import KoopMurray2016_DWA

PySDM/physics/homogeneous_ice_nucleation_rate/koop_murray_dwa_adaption.py

@@ -0,0 +1,34 @@
+"""
+Koop and Murray homogeneous nucleation rate parameterization for pure water droplets
+adapted for water activity formulation
+ ([Spichtinger et al. 2023](https://doi.org/10.5194/acp-23-2035-2023))
+"""
+
+import numpy as np
+
+
+class KoopMurray2016_DWA:
+    def __init__(self, const):
+        pass
+
+    @staticmethod
+    def d_a_w_ice_within_range(const, da_w_ice):
+        return da_w_ice >= const.KOOP_MIN_DA_W_ICE
+
+    @staticmethod
+    def d_a_w_ice_maximum(const, da_w_ice):
+        return np.where(
+            da_w_ice > const.KOOP_MAX_DA_W_ICE, const.KOOP_MAX_DA_W_ICE, da_w_ice
+        )
+
+    @staticmethod
+    def j_hom(const, T, da_w_ice):  # pylint: disable=unused-argument
+        return (
+            10
+            ** (
+                const.KOOP_MURRAY_DAW_C0
+                + const.KOOP_MURRAY_DAW_C1 * da_w_ice
+                + const.KOOP_MURRAY_DAW_C2 * da_w_ice**2.0
+            )
+            * const.KOOP_MURRAY_DAW_UNIT
+        )

docs/bibliography.json

@@ -948,6 +948,7 @@
       "PySDM/physics/constants_defaults.py",
       "PySDM/physics/diffusion_ice_capacity/columnar.py",
       "PySDM/physics/homogeneous_ice_nucleation_rate/koop_corr.py",
+      "PySDM/physics/homogeneous_ice_nucleation_rate/koop_murray_dwa_adaption.py",
       "examples/PySDM_examples/Spichtinger_et_al_2023/__init__.py",
       "examples/PySDM_examples/Spichtinger_et_al_2023/data/reference_bulk.py",
       "examples/PySDM_examples/Spichtinger_et_al_2023/fig_B1.ipynb",
@@ -1014,14 +1015,22 @@
     "label": "Bartman 2020 (MSc thesis)",
     "title": "PySDM v1.0: Pythonic particle-based cloud microphysics package"
   },
-    "https://web.archive.org/web/20220629123450/https://web.gps.caltech.edu/~als/research-articles/other_stuff/hayes-2004-3.pdf" : {
+  "https://web.archive.org/web/20220629123450/https://web.gps.caltech.edu/~als/research-articles/other_stuff/hayes-2004-3.pdf" : {
     "usages": [
       "tests/unit_tests/physics/test_constants.py"
     ],
     "label": "Hayes 2004",
     "title": "An Introduction to Isotopic Calculations"
   },
-  "https://doi.org/10.1175/1520-0469%281974%29031%3C1040%3ASTTDCE%3E2.0.CO%3B2" : {
+  "https://doi.org/10.1039/c003297b": {
+    "usages": [
+      "examples/PySDM_examples/Murray_et_al_2010/fig_3.ipynb",
+      "examples/PySDM_examples/Murray_et_al_2010/__init__.py"
+    ],
+    "title": "Kinetics of the homogeneous freezing of water",
+    "label": "Murray et al. 2010 (Phys. Chem. Chem. Phys. 12)"
+  },
+  "https://doi.org/10.1175/1520-0469%281974%29031%3C1040%3ASTTDCE%3E2.0.CO%3B2": {
     "usages": [
       "PySDM/dynamics/collisions/collision_kernels/long1974.py",
       "examples/PySDM_examples/Long_1974/__init__.py",

examples/PySDM_examples/Luettmer_homogeneous_freezing/__init__.py

@@ -0,0 +1,17 @@
+"""
+Homogeneous freezing example
+
+.. include:: ./fig_1_2_3.ipynb
+.. include:: ./fig_4_5_6_S3_S4.ipynb
+.. include:: ./fig_S1_S2.ipynb
+.. include:: ./simple_homogenous_freezing_example.ipynb
+"""
+
+from .settings import Settings
+from .simulation import Simulation
+from .plot import (
+    plot_thermodynamics_and_bulk,
+    plot_freezing_temperatures_histogram,
+    plot_freezing_temperatures_histogram_allinone,
+    plot_freezing_temperatures_2d_histogram_seaborn,
+)

examples/PySDM_examples/Luettmer_homogeneous_freezing/commons.py

@@ -0,0 +1,91 @@
+"""commons for homogeneous freezing notebooks"""
+
+from PySDM_examples.Luettmer_homogeneous_freezing.settings import Settings
+from PySDM_examples.Luettmer_homogeneous_freezing.simulation import Simulation
+from PySDM import Formulae
+from PySDM.physics.constants import si
+from PySDM.backends import CPU
+
+formulae = Formulae(
+    particle_shape_and_density="MixedPhaseSpheres",
+)
+
+
+def run_simulations(setting):
+
+    simulation = {
+        "settings": setting,
+        "ensemble_member_outputs": [],
+    }
+    for _ in range(setting["number_of_ensemble_runs"]):
+        model_setup = Settings(**simulation["settings"])
+        model_setup.formulae.seed += 1
+        model = Simulation(model_setup)
+        simulation["ensemble_member_outputs"].append(model.run())
+
+    return simulation
+
+
+def hom_pure_droplet_freezing_backend():
+    cmn = {"override_jit_flags": {"parallel": False}}
+    backends = {
+        "threshold": CPU(
+            formulae=Formulae(
+                particle_shape_and_density="MixedPhaseSpheres",
+                homogeneous_ice_nucleation_rate="Null",
+                saturation_vapour_pressure="MurphyKoop2005",
+            ),
+            **cmn,
+        ),
+        "KoopMurray2016": CPU(
+            formulae=Formulae(
+                particle_shape_and_density="MixedPhaseSpheres",
+                homogeneous_ice_nucleation_rate="KoopMurray2016",
+                saturation_vapour_pressure="MurphyKoop2005",
+            ),
+            **cmn,
+        ),
+        "Spichtinger2023": CPU(
+            formulae=Formulae(
+                particle_shape_and_density="MixedPhaseSpheres",
+                homogeneous_ice_nucleation_rate="Koop_Correction",
+                saturation_vapour_pressure="MurphyKoop2005",
+            ),
+            **cmn,
+        ),
+        "Koop2000": CPU(
+            formulae=Formulae(
+                particle_shape_and_density="MixedPhaseSpheres",
+                homogeneous_ice_nucleation_rate="Koop2000",
+                saturation_vapour_pressure="MurphyKoop2005",
+            ),
+            **cmn,
+        ),
+        "KoopMurray2016_DWA": CPU(
+            formulae=Formulae(
+                particle_shape_and_density="MixedPhaseSpheres",
+                homogeneous_ice_nucleation_rate="KoopMurray2016_DWA",
+                saturation_vapour_pressure="MurphyKoop2005",
+            ),
+            **cmn,
+        ),
+    }
+    return backends
+
+
+def hom_pure_droplet_freezing_standard_setup():
+    standard = {
+        "n_sd": int(1e3),
+        "w_updraft": 1.0 * si.meter / si.second,
+        "T0": 245 * si.kelvin,
+        "dz": 1.0 * si.meter,
+        "n_ccn": 750 / si.cm**3,
+        "r_ccn": 15 * si.nanometer,
+        "type_droplet_distribution": "monodisperse",
+        "RH_0": 0.995,
+        "p0": 500 * si.hectopascals,
+        "deposition_enable": True,
+        "number_of_ensemble_runs": 1,
+        "silent": True,
+    }
+    return standard