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

PySDM/particulator.py

@@ -546,6 +546,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(
@@ -556,6 +557,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(
@@ -565,3 +567,4 @@ def thaw_instantaneous(self):
             cell=self.attributes["cell id"],
             temperature=self.environment["T"],
         )
+        self.attributes.mark_updated("signed water mass")

examples/PySDM_examples/Luettmer_homogeneous_freezing/__init__.py

@@ -0,0 +1,5 @@
+"""
+Homogeneous freezing example
+"""
+
+from .settings import Settings

examples/PySDM_examples/Luettmer_homogeneous_freezing/plot.py

@@ -0,0 +1,254 @@
+import matplotlib.pyplot as plt
+from matplotlib import pyplot
+import numpy as np
+from scipy.ndimage import histogram
+import seaborn as sns
+
+from PySDM import Formulae
+
+formulae = Formulae(
+    particle_shape_and_density="MixedPhaseSpheres",
+)
+
+# general plot settings
+ax_lab_fsize = 15
+tick_fsize = 15
+# title_fsize = 15
+# line_width = 2.5
+T_frz_bins = np.linspace(-38.5, -33, num=7, endpoint=True)
+
+
+def plot_thermodynamics_and_bulk(simulation, title_add=""):
+
+    output = simulation["ensemble_member_outputs"][0]
+    time = output["t"]
+    T = np.asarray(output["T"])
+    RH = np.asarray(output["RH"])
+    RHi = np.asarray(output["RHi"])
+    qc = np.asarray(output["LWC"])
+    qi = np.asarray(output["IWC"])
+    qv = np.asarray(output["qv"])
+    qt = qc + qv + qi
+    rc = np.asarray(output["rs"])
+    ri = np.asarray(output["ri"])
+
+    svp = Formulae(
+        saturation_vapour_pressure="FlatauWalkoCotton"
+    ).saturation_vapour_pressure
+    a_w_ice = svp.pvs_ice(T) / svp.pvs_water(T)
+    d_a_w_ice = (RHi / 100 - 1) * a_w_ice
+
+    j_hom_rate = Formulae(
+        homogeneous_ice_nucleation_rate="Koop2000"
+    ).homogeneous_ice_nucleation_rate
+    koop_2000 = j_hom_rate.j_hom(T, d_a_w_ice)
+    j_hom_rate = Formulae(
+        homogeneous_ice_nucleation_rate="KoopMurray2016"
+    ).homogeneous_ice_nucleation_rate
+    koop_murray_2016 = j_hom_rate.j_hom(T, d_a_w_ice)
+    j_hom_rate = Formulae(
+        homogeneous_ice_nucleation_rate="Koop_Correction"
+    ).homogeneous_ice_nucleation_rate
+    spichtinger_2023 = j_hom_rate.j_hom(T, d_a_w_ice)
+
+    fig, axs = pyplot.subplots(
+        2, 2, figsize=(10, 10), sharex=False, constrained_layout=True
+    )
+
+    title = (
+        "Freezing method="
+        + simulation["settings"]["hom_freezing"]
+        + " n_sd="
+        + str(simulation["settings"]["n_sd"])
+        + " w="
+        + str(simulation["settings"]["w_updraft"])
+        + title_add
+    )
+    fig.suptitle(title, fontsize=16)
+
+    """ Temperture profile """
+    ax = axs[0, 0]
+    ax.plot(time, formulae.trivia.K2C(T), color="black", linestyle="-", label="T")
+    ax.set_xlabel("time [s]", fontsize=ax_lab_fsize)
+    ax.set_ylabel("temperature [°C]", fontsize=ax_lab_fsize)
+    ax.legend(loc="upper right", fontsize=ax_lab_fsize)
+    ax.tick_params(labelsize=tick_fsize)
+
+    twin = ax.twinx()
+    twin.plot(time, RH, color="red", linestyle="-", label="RH")
+    twin.plot(time, RHi, color="blue", linestyle="-", label="RHi")
+    twin.set_ylabel("relative humidity [%]", fontsize=ax_lab_fsize)
+    twin.legend(loc="upper left", fontsize=ax_lab_fsize)
+    twin.tick_params(labelsize=tick_fsize)
+
+    """ Water activity difference profile """
+    ax = axs[0, 1]
+    ax.plot(time, d_a_w_ice, color="gray", linestyle="-")
+    ax.set_xlabel("time [s]", fontsize=ax_lab_fsize)
+    ax.set_ylabel("water activity difference", fontsize=ax_lab_fsize)
+    ax.set_ylim(0.2, 0.35)
+    ax.tick_params(labelsize=tick_fsize)
+    twin = ax.twinx()
+    twin.plot(time, koop_2000, color="black", linestyle="-", label="Koop2000")
+    twin.plot(
+        time, koop_murray_2016, color="blue", linestyle="-", label="KoopMurray2016"
+    )
+    twin.plot(
+        time, spichtinger_2023, color="red", linestyle="-", label="Spichtinger2023"
+    )
+    twin.set_yscale("log")
+    twin.set_ylabel("nucleation rate", fontsize=ax_lab_fsize)
+    twin.tick_params(labelsize=tick_fsize)
+    twin.legend(fontsize=ax_lab_fsize)
+
+    """ Mass content """
+    ax = axs[1, 0]
+    ax.plot(time, qc, color="red", linestyle="-", label="water")
+    ax.plot(time, qi, color="blue", linestyle="-", label="ice")
+    ax.plot(time, qv, color="black", linestyle="-", label="vapor")
+    ax.plot(time, qt, color="black", linestyle="dotted", label="total")
+    ax.set_yscale("log")
+    ax.set_ylim(1e-5, 1e-2)
+    ax.set_xlabel("time [s]", fontsize=ax_lab_fsize)
+    ax.set_ylabel("mass content [kg/kg]", fontsize=ax_lab_fsize)
+    ax.legend(fontsize=ax_lab_fsize)
+    ax.tick_params(labelsize=tick_fsize)
+
+    """ Mean radius """
+    ax = axs[1, 1]
+    ax.plot(time, rc * 1e6, color="red", linestyle="-", label="water")
+    ax.plot(time, ri * 1e6, color="blue", linestyle="-", label="ice")
+    ax.set_yscale("log")
+    ax.set_ylim(1e-2, 1e2)
+    ax.set_xlabel("time [s]", fontsize=ax_lab_fsize)
+    ax.set_ylabel("mean radius [µm]", fontsize=ax_lab_fsize)
+    ax.legend(fontsize=ax_lab_fsize)
+    ax.tick_params(labelsize=tick_fsize)
+
+
+def plot_freezing_temperatures_histogram(ax, simulation):
+
+    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"][-1])
+
+        title = "Freezing method=" + simulation["settings"]["hom_freezing"]
+
+        """ Freezing temperatures """
+        hist = ax.hist(
+            formulae.trivia.K2C(T_frz),
+            bins=T_frz_bins,
+            density=True,
+            cumulative=-1,
+            alpha=1.0,
+            histtype="step",
+            linewidth=1.5,
+        )
+
+    ax.set_title(title, fontsize=ax_lab_fsize)
+    ax.set_xlabel("freezing temperature [°C]", fontsize=ax_lab_fsize)
+    ax.set_ylabel("frequency", fontsize=ax_lab_fsize)
+    ax.tick_params(labelsize=tick_fsize)
+
+    return ax
+
+
+def plot_freezing_temperatures_2d_histogram(histogram_data_dict):
+
+    vertical_updrafts_bins = np.geomspace(0.05, 15, num=6, endpoint=True)
+
+    hom_freezing_types = histogram_data_dict.keys()
+
+    fig, axs = pyplot.subplots(2, 2, figsize=(10, 10), constrained_layout=True)
+    axs = axs.ravel()
+    i = 0
+    for hom_freezing_type in hom_freezing_types:
+
+        title = "Freezing method=" + hom_freezing_type
+
+        ax = axs[i]
+        T_frz = formulae.trivia.K2C(
+            np.asarray(histogram_data_dict[hom_freezing_type]["T_frz_histogram_list"])
+        )
+
+        hist, x, y = np.histogram2d(
+            T_frz,
+            histogram_data_dict[hom_freezing_type]["w_updraft_histogram_list"],
+            bins=(T_frz_bins, vertical_updrafts_bins),
+        )
+        y = np.log10(y)
+        X, Y = np.meshgrid(x, y)
+
+        hist = hist.T / sum(hist.flatten())
+
+        c = ax.pcolor(X, Y, hist)
+        fig.colorbar(c, ax=ax)
+
+        ax.set_title(title, fontsize=ax_lab_fsize)
+        ax.set_xlabel("freezing temperature [°C]", fontsize=ax_lab_fsize)
+        ax.set_ylabel("vertical updraft [m/s]", fontsize=ax_lab_fsize)
+
+        i += 1
+
+
+def plot_freezing_temperatures_2d_histogram_seaborn(histogram_data_dict, title_add=""):
+
+    # hom_freezing_types = ["KoopMurray2016", "Koop_Correction", "Koop2000"]
+    sns.set_theme(style="ticks")
+
+    hom_freezing_type = histogram_data_dict["hom_freezing_type"]
+    # for i, hom_freezing_type in enumerate(hom_freezing_types):
+
+    T_frz = formulae.trivia.K2C(
+        (np.asarray(histogram_data_dict["T_frz_histogram_list"]))
+    )
+    if "w_updraft_histogram_list" in histogram_data_dict:
+        w = histogram_data_dict["w_updraft_histogram_list"]
+        y_label = "vertical updraft [m/s]"
+    elif "n_ccn_histogram_list" in histogram_data_dict:
+        w = histogram_data_dict["n_ccn_histogram_list"]
+        y_label = "ccn concentration [1/m^³]"
+    elif "rc_max_histogram_list" in histogram_data_dict:
+        w = (
+            np.asarray(
+                histogram_data_dict["rc_max_histogram_list"]
+            )
+            * 1e6
+        )
+        y_label = "(maximum) radius [µm]"
+
+    xlim = (-38.5, -32)
+    h = sns.JointGrid(
+        x=T_frz,
+        y=w,
+        xlim=xlim,
+    )
+    h.ax_joint.set(yscale="log")
+    if hom_freezing_type == "KoopMurray2016":
+        x_pos_cbar = 0.75
+    else:
+        x_pos_cbar = 0.15
+    cax = h.figure.add_axes([x_pos_cbar, 0.55, 0.02, 0.2])
+    h.plot_joint(
+        sns.histplot,
+        stat="density",
+        binwidth=0.25,
+        discrete=(False, False),
+        pmax=0.8,
+        cbar=True,
+        cbar_ax=cax,
+    )
+
+    h.plot_marginals(
+        sns.histplot,
+        element="step",
+    )
+    h.set_axis_labels("freezing temperature [°C]", y_label, fontsize=ax_lab_fsize)
+    h.ax_joint.set_title(
+        "Freezing method=" + hom_freezing_type + title_add, pad=70, fontsize=ax_lab_fsize
+    )
+    h.ax_marg_y.remove()
+
+    return h

examples/PySDM_examples/Luettmer_homogeneous_freezing/settings.py

@@ -0,0 +1,111 @@
+import time
+
+import numpy as np
+from pystrict import strict
+
+from PySDM import Formulae
+from PySDM.physics.constants import si
+from PySDM.initialisation.spectra import Lognormal
+from PySDM.initialisation.sampling import spectral_sampling
+from tests.unit_tests.backends.test_oxidation import formulae
+
+
+class Settings:
+    def __init__(
+        self,
+        *,
+        hom_freezing: str,
+        n_sd: int,
+        w_updraft: float,
+        T0: float,
+        # dt: float,
+        dz: float,
+        N_dv_droplet_distribution: float,
+        r_mean_droplet_distribution: float,
+        sigma_droplet_distribution: float = None,
+        type_droplet_distribution: str,
+        p0: float = 200 * si.hectopascals,
+        RH_0: float = 1.0,
+        kappa: float = 0.64,
+        condensation_enable: bool = True,
+        deposition_enable: bool = True,
+        deposition_adaptive: bool = False,
+        backend=None,
+        number_of_ensemble_runs: None,
+    ):
+        self.backend = backend
+        print(
+            "Setting up simulation for "
+            + hom_freezing
+            + " with wpdraft="
+            + str(w_updraft)
+            + " and n_sd="
+            + str(n_sd)
+            + " and n_dv="
+            + str(N_dv_droplet_distribution)
+        )
+        self.n_sd = n_sd
+        self.w_updraft = w_updraft
+        self.N_dv_droplet_distribution = N_dv_droplet_distribution
+        self.r_mean_droplet_distribution = r_mean_droplet_distribution
+        self.sigma_droplet_distribution = sigma_droplet_distribution
+        self.type_droplet_distribution = type_droplet_distribution
+
+        self.mass_of_dry_air = 1000 * si.kilogram
+        self.initial_pressure = p0
+        self.initial_water_supersaturation = RH_0
+        # self.initial_ice_supersaturation = RHi_0
+        self.kappa = kappa
+        self.initial_temperature = T0
+
+        self.condensation_enable = condensation_enable
+        self.deposition_enable = deposition_enable
+        self.deposition_adaptive = deposition_adaptive
+
+        if hom_freezing == "threshold":
+            self.hom_freezing_type = "threshold"
+        else:
+            self.hom_freezing_type = "time-dependent"
+
+        self.formulae = backend.formulae
+
+        const = self.formulae.constants
+        pvs_w = self.formulae.saturation_vapour_pressure.pvs_water(
+            self.initial_temperature
+        )
+        self.initial_water_vapour_mixing_ratio = const.eps / (
+            self.initial_pressure / self.initial_water_supersaturation / pvs_w - 1
+        )
+
+        dry_air_density = (
+            self.formulae.trivia.p_d(
+                self.initial_pressure, self.initial_water_vapour_mixing_ratio
+            )
+            / self.initial_temperature
+            / const.Rd
+        )
+
+        if self.type_droplet_distribution == ("monodisperse"):
+            self.r_dry = np.ones(self.n_sd) * r_mean_droplet_distribution
+            self.specific_concentration = (
+                np.ones(self.n_sd)
+                * N_dv_droplet_distribution
+                / self.n_sd
+                / dry_air_density
+            )
+
+        elif self.type_droplet_distribution == ("lognormal"):
+            spectrum = Lognormal(
+                norm_factor=N_dv_droplet_distribution / dry_air_density,
+                m_mode=r_mean_droplet_distribution,
+                s_geom=sigma_droplet_distribution,
+            )
+            self.r_dry, self.specific_concentration = spectral_sampling.Linear(
+                spectrum
+            ).sample(n_sd)
+
+        self.dz = dz
+        self.t_max_duration = 20000
+        self.dt = dz / self.w_updraft
+        print( self.dt, dz, self.w_updraft)
+        self.n_output = 1