new attribute: freezing supersaturation (recorded for the last freezing event) (#1878)

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

Author: tluettm

PySDM/attributes/ice/__init__.py

@@ -4,4 +4,5 @@
 
 from .cooling_rate import CoolingRate
 from .freezing_temperature import FreezingTemperature
+from .freezing_supersaturation import SupersaturationOfLastFreezing
 from .immersed_surface_area import ImmersedSurfaceArea

PySDM/attributes/ice/freezing_supersaturation.py

@@ -0,0 +1,36 @@
+"""
+particle freezing supersaturation wrt ice
+"""
+
+from PySDM.attributes.impl import register_attribute
+from ..impl import DerivedAttribute
+
+
+@register_attribute()
+class SupersaturationOfLastFreezing(DerivedAttribute):
+    """time-dependent variant: assigned upon freezing"""
+
+    def __init__(self, builder):
+        assert "Freezing" in builder.particulator.dynamics
+        assert (
+            builder.particulator.dynamics["Freezing"].immersion_freezing != "singular"
+        )
+        self.signed_water_mass = builder.get_attribute("signed water mass")
+        self.cell_id = builder.get_attribute("cell id")
+        super().__init__(
+            builder,
+            name="supersaturation of last freezing",
+            dependencies=(self.signed_water_mass, self.cell_id),
+        )
+        builder.particulator.observers.append(self)
+
+    def notify(self):
+        self.update()
+
+    def recalculate(self):
+        self.particulator.backend.record_freezing_supersaturations(
+            data=self.data,
+            cell_id=self.cell_id.data,
+            relative_humidity_ice=self.particulator.environment["RH_ice"],
+            signed_water_mass=self.signed_water_mass.data,
+        )

PySDM/backends/impl_numba/methods/freezing_methods.py

@@ -299,3 +299,29 @@ def record_freezing_temperatures(
             temperature=temperature.data,
             signed_water_mass=signed_water_mass.data,
         )
+
+    @cached_property
+    def _record_freezing_supersaturations_body(self):
+        ff = self.formulae_flattened
+
+        @numba.njit(**{**self.default_jit_flags, "fastmath": False})
+        def body(data, cell_id, relative_humidity_ice, signed_water_mass):
+            for drop_id in numba.prange(len(data)):  # pylint: disable=not-an-iterable
+                if ff.trivia__unfrozen(signed_water_mass[drop_id]):
+                    if data[drop_id] > 0:
+                        data[drop_id] = np.nan
+                else:
+                    if np.isnan(data[drop_id]):
+                        data[drop_id] = relative_humidity_ice[cell_id[drop_id]]
+
+        return body
+
+    def record_freezing_supersaturations(
+        self, *, data, cell_id, relative_humidity_ice, signed_water_mass
+    ):
+        self._record_freezing_supersaturations_body(
+            data=data.data,
+            cell_id=cell_id.data,
+            relative_humidity_ice=relative_humidity_ice.data,
+            signed_water_mass=signed_water_mass.data,
+        )

examples/PySDM_examples/Luettmer_homogeneous_freezing/plot.py

@@ -16,18 +16,15 @@
 
 
 def cumulative_histogram(data, bins, reverse=False, density=True):
-    # Compute regular histogram using given bins
     hist, bin_edges = np.histogram(data, bins=bins, density=False)
 
-    # Cumulative sum
     if reverse:
         cum_hist = np.cumsum(hist[::-1])[::-1]
         cum_hist_0 = cum_hist[0]
     else:
         cum_hist = np.cumsum(hist)
         cum_hist_0 = cum_hist[-1]
 
-    # Normalize
     if density:
         cum_hist = cum_hist / cum_hist_0
 
@@ -56,8 +53,7 @@ def plot_thermodynamics_and_bulk(
     qv = np.asarray(output["qv"])
     T_frz = np.asarray(output["T_frz"])
     if show_conc:
-        nc = np.asarray(output["ns"])
-        ni = np.asarray(output["ni"])
+        nc, ni = np.asarray(output["ns"]), np.asarray(output["ni"])
 
     if t_lim is None:
         t_lim = np.amax(time)
@@ -67,8 +63,7 @@ def plot_thermodynamics_and_bulk(
     first_T_frz = T[first_ice_idx]
 
     if not show_tf:
-        rc = np.asarray(output["rs"])
-        ri = np.asarray(output["ri"])
+        rc, ri = np.asarray(output["rs"]), np.asarray(output["ri"])
 
     if show_jhom:
         svp = Formulae(
@@ -90,9 +85,7 @@ def plot_thermodynamics_and_bulk(
         radius = np.asarray(output["radius"])
         multiplicity = np.asarray(output["multiplicity"])
 
-    _, axs = pyplot.subplots(
-        1, 4, figsize=(20, 5), sharex=False, constrained_layout=True
-    )
+    _, axs = pyplot.subplots(1, 4, figsize=(20, 5), constrained_layout=True)
 
     # Temperture profile
     iax = 0
@@ -247,30 +240,40 @@ def plot_thermodynamics_and_bulk(
     ax.grid(visible=True)
 
 
-def plot_freezing_temperatures_histogram(ax, simulation):
+def plot_freezing_temperatures_histogram(ax, simulation, plot_rhi=False):
 
     number_of_ensemble_runs = simulation["settings"]["number_of_ensemble_runs"]
 
     for i in range(number_of_ensemble_runs):
         output = simulation["ensemble_member_outputs"][i]
-        T_frz = np.asarray(output["T_frz"])
+        if plot_rhi:
+            var = np.asarray(output["RHi_frz"])
+            bins = np.linspace(1, 1.6, num=60, endpoint=True)
+            cumulative = 1
+        else:
+            var = np.asarray(output["T_frz"])
+            bins = T_frz_bins_kelvin
+            cumulative = -1
         title = "Nucleation rate=" + simulation["settings"]["hom_freezing"]
 
-        # Freezing temperatures
         ax.hist(
-            T_frz,
-            bins=T_frz_bins_kelvin,
+            var,
+            bins=bins,
             density=True,
-            cumulative=-1,
+            cumulative=cumulative,
             alpha=1.0,
             histtype="step",
             linewidth=1.5,
         )
 
-        ax.set_xlim(left=234, right=239)
-        ax.axvline(x=235, color="k", linestyle="--")
+        if plot_rhi:
+            ax.set_xlim(left=1.0, right=1.6)
+            ax.set_xlabel("freezing supersaturation wrt ice", fontsize=ax_lab_fsize)
+        else:
+            ax.set_xlim(left=234, right=239)
+            ax.axvline(x=235, color="k", linestyle="--")
+            ax.set_xlabel("freezing temperature [K]", fontsize=ax_lab_fsize)
         ax.set_title(title, fontsize=ax_lab_fsize)
-        ax.set_xlabel("freezing temperature [K]", fontsize=ax_lab_fsize)
         ax.set_ylabel("frozen fraction", fontsize=ax_lab_fsize)
         ax.tick_params(labelsize=tick_fsize)
     return ax

examples/PySDM_examples/Luettmer_homogeneous_freezing/simple_homogenous_freezing_example.ipynb

@@ -79,6 +79,7 @@ def __init__(self, settings):
             "signed water mass": self.initial_mass,
         }
         builder.request_attribute("temperature of last freezing")
+        builder.request_attribute("supersaturation of last freezing")
         builder.request_attribute("radius")
         builder.request_attribute("wet to critical volume ratio")
 
@@ -140,6 +141,9 @@ def save(self, output):
         output["T_frz"] = self.particulator.attributes[
             "temperature of last freezing"
         ].data.tolist()
+        output["RHi_frz"] = self.particulator.attributes[
+            "supersaturation of last freezing"
+        ].data.tolist()
         if not output["radius"]:
             output["radius"] = self.particulator.attributes["radius"].data.tolist()
             output["multiplicity"] = self.particulator.attributes[
@@ -153,6 +157,7 @@ def run(self):
 
         output = {
             "T_frz": [],
+            "RHi_frz": [],
             "radius": [],
             "multiplicity": [],
         }

examples/PySDM_examples/Luettmer_homogeneous_freezing/simulation.py

@@ -38,6 +38,10 @@ def test_record_freezing_temperature_on_time_dependent_freeze(
             pytest.skip("TODO #1495")
 
         # arrange
+        attr_names = (
+            "temperature of last freezing",
+            "supersaturation of last freezing",
+        )
         formulae = Formulae(
             particle_shape_and_density="MixedPhaseSpheres",
             heterogeneous_ice_nucleation_rate="Constant",
@@ -52,15 +56,18 @@ def test_record_freezing_temperature_on_time_dependent_freeze(
             backend=backend_class(formulae=formulae),
             environment=Box(dt=1 * si.s, dv=1 * si.m**3),
             dynamics=(
-                Freezing(
-                    immersion_freezing=freezing_mode["immersion_freezing"],
-                    homogeneous_freezing=freezing_mode["homogeneous_freezing"],
-                    thaw="instantaneous",
-                ),
+                [
+                    Freezing(
+                        immersion_freezing=freezing_mode["immersion_freezing"],
+                        homogeneous_freezing=freezing_mode["homogeneous_freezing"],
+                        thaw="instantaneous",
+                    )
+                ]
             ),
         )
         if record_freezing_temperature:
-            builder.request_attribute("temperature of last freezing")
+            for attr_name in attr_names:
+                builder.request_attribute(attr_name)
         particulator = builder.build(
             attributes={
                 "multiplicity": np.asarray([1]),
@@ -69,46 +76,59 @@ def test_record_freezing_temperature_on_time_dependent_freeze(
             }
         )
 
-        temp_1 = 200 * si.K
-        temp_2 = 250 * si.K
-        temp_3 = 220 * si.K
+        temp_list = np.asarray([200 * si.K, 300 * si.K, 220 * si.K])
+        rh_list = np.asarray([0.9, 1.0, 1.1])
         particulator.environment["a_w_ice"] = np.nan
-        particulator.environment["T"] = temp_2
+        particulator.environment["T"] = temp_list[1]
+        particulator.environment["RH_ice"] = rh_list[0]
 
         # act & assert
-        attr_name = "temperature of last freezing"
         if not record_freezing_temperature:
-            assert attr_name not in particulator.attributes
+            for attr_name in attr_names:
+                assert attr_name not in particulator.attributes
         else:
             # never frozen yet
-            np.isnan(particulator.attributes[attr_name].to_ndarray()).all()
+            for attr_name in attr_names:
+                np.isnan(particulator.attributes[attr_name].to_ndarray()).all()
 
             # still not frozen since RH not greater than 100%
-            particulator.environment["RH"] = 1.0
-            particulator.environment["RH_ice"] = 1.0
+            particulator.environment["RH"] = rh_list[1]
+            particulator.environment["RH_ice"] = rh_list[1]
             particulator.run(steps=1)
-            np.isnan(particulator.attributes[attr_name].to_ndarray()).all()
+            for attr_name in attr_names:
+                np.isnan(particulator.attributes[attr_name].to_ndarray()).all()
             assert all(particulator.attributes["signed water mass"].to_ndarray() > 0)
 
             # should freeze and record T1
-            particulator.environment["RH"] += EPSILON_RH
-            particulator.environment["RH_ice"] += EPSILON_RH
-            particulator.environment["T"] = temp_1
+            particulator.environment["RH"] = rh_list[2]
+            particulator.environment["RH_ice"] = rh_list[2]
+            particulator.environment["T"] = temp_list[0]
             particulator.run(steps=1)
             assert all(particulator.attributes["signed water mass"].to_ndarray() < 0)
-            assert all(particulator.attributes[attr_name].to_ndarray() == temp_1)
+            assert all(
+                particulator.attributes[attr_names[0]].to_ndarray() == temp_list[0]
+            )
+            assert all(
+                particulator.attributes[attr_names[1]].to_ndarray() == rh_list[2]
+            )
 
             # should thaw
-            particulator.environment["T"] = 300 * si.K
+            particulator.environment["T"] = temp_list[1]
             particulator.run(steps=1)
             assert all(particulator.attributes["signed water mass"].to_ndarray() > 0)
-            np.isnan(particulator.attributes[attr_name].to_ndarray()).all()
+            for attr_name in attr_names:
+                np.isnan(particulator.attributes[attr_name].to_ndarray()).all()
 
             # should re-freeze and record T2
-            particulator.environment["T"] = temp_3
+            particulator.environment["T"] = temp_list[2]
             particulator.run(steps=1)
             assert all(particulator.attributes["signed water mass"].to_ndarray() < 0)
-            assert all(particulator.attributes[attr_name].to_ndarray() == temp_3)
+            assert all(
+                particulator.attributes[attr_names[0]].to_ndarray() == temp_list[2]
+            )
+            assert all(
+                particulator.attributes[attr_names[1]].to_ndarray() == rh_list[2]
+            )
 
     # TODO #599
     def test_no_subsaturated_freezing(self):