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

PySDM/particulator.py

@@ -550,6 +550,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(
@@ -560,6 +561,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(
@@ -569,3 +571,4 @@ def thaw_instantaneous(self):
             cell=self.attributes["cell id"],
             temperature=self.environment["T"],
         )
+        self.attributes.mark_updated("signed water mass")

examples/PySDM_examples/homogeneous_freezing/__init__.py

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

examples/PySDM_examples/homogeneous_freezing/settings.py

@@ -0,0 +1,107 @@
+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,
+        *,
+        n_sd: int,
+        w_updraft: float,
+        T0: float,
+        dt: 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,
+        hom_freezing: str=None,
+    ):
+
+        print("Setting up simulation with " + hom_freezing)
+        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
+
+
+
+        if hom_freezing == "threshold":
+            hom_nucleation_rate = "Null"
+            self.hom_freezing_type = "threshold"
+        else:
+            hom_nucleation_rate = hom_freezing
+            self.hom_freezing_type = "time-dependent"
+
+        formulae = Formulae(
+            particle_shape_and_density="MixedPhaseSpheres",
+            homogeneous_ice_nucleation_rate=hom_nucleation_rate,
+            seed=time.time_ns()
+        )
+
+        self.formulae = formulae
+
+        const = self.formulae.constants
+        # pvs_i = self.formulae.saturation_vapour_pressure.pvs_ice(self.initial_temperature)
+        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.t_max_duration = 7200  # 3600 * 1.5 # total duration of simulation
+        self.dt = dt
+        self.n_output = 10  # int(self.t_duration / 100) #100 # number of output steps

examples/PySDM_examples/homogeneous_freezing/simulation.py

@@ -0,0 +1,174 @@
+import numpy as np
+
+import PySDM.products as PySDM_products
+from PySDM.backends import CPU
+from PySDM.builder import Builder
+from PySDM.dynamics import (
+    AmbientThermodynamics,
+    Condensation,
+    Freezing,
+    VapourDepositionOnIce,
+)
+from PySDM.environments import Parcel
+from PySDM.physics import constants as const
+from PySDM.initialisation import discretise_multiplicities, equilibrate_wet_radii
+
+
+class Simulation:
+    def __init__(self, settings, backend=CPU):
+
+        dt = settings.dt
+
+        formulae = settings.formulae
+
+        env = Parcel(
+            mixed_phase=True,
+            dt=dt,
+            mass_of_dry_air=settings.mass_of_dry_air,
+            p0=settings.initial_pressure,
+            initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio,
+            T0=settings.initial_temperature,
+            w=settings.w_updraft,
+        )
+
+        builder = Builder(
+            backend=backend(
+                formulae=settings.formulae,
+                **(
+                    {"override_jit_flags": {"parallel": False}}
+                    if backend == CPU
+                    else {}
+                ),
+            ),
+            n_sd=settings.n_sd,
+            environment=env,
+        )
+
+        builder.add_dynamic(AmbientThermodynamics())
+        if settings.condensation_enable:
+            builder.add_dynamic(Condensation())
+        if settings.deposition_enable:
+            builder.add_dynamic(VapourDepositionOnIce())
+        builder.add_dynamic(
+            Freezing(
+                homogeneous_freezing=settings.hom_freezing_type, immersion_freezing=None
+            )
+        )
+
+        self.n_sd = settings.n_sd
+        self.multiplicities = discretise_multiplicities(
+            settings.specific_concentration * env.mass_of_dry_air
+        )
+        self.r_dry = settings.r_dry
+        v_dry = settings.formulae.trivia.volume(radius=self.r_dry)
+        kappa = settings.kappa
+
+        self.r_wet = equilibrate_wet_radii(
+            r_dry=self.r_dry,
+            environment=builder.particulator.environment,
+            kappa_times_dry_volume=kappa * v_dry,
+        )
+        attributes = {
+            "multiplicity": self.multiplicities,
+            "dry volume": v_dry,
+            "kappa times dry volume": kappa * v_dry,
+            "signed water mass": formulae.particle_shape_and_density.radius_to_mass(
+                self.r_wet
+            ),
+        }
+        builder.request_attribute("temperature of last freezing")
+
+        products = [
+            PySDM_products.ParcelDisplacement(name="z"),
+            PySDM_products.Time(name="t"),
+            PySDM_products.AmbientRelativeHumidity(name="RH", unit="%"),
+            PySDM_products.AmbientRelativeHumidity(name="RH_ice", unit="%"),
+            PySDM_products.AmbientTemperature(name="T"),
+            PySDM_products.AmbientPressure(name="p", unit="hPa"),
+            PySDM_products.WaterMixingRatio(name="water", radius_range=(0, np.inf)),
+            PySDM_products.WaterMixingRatio(name="ice", radius_range=(-np.inf, 0)),
+            PySDM_products.WaterMixingRatio(
+                name="total", radius_range=(-np.inf, np.inf)
+            ),
+            PySDM_products.AmbientWaterVapourMixingRatio(
+                name="vapour", var="water_vapour_mixing_ratio"
+            ),
+            PySDM_products.ParticleConcentration(
+                name="n_s", unit="1/cm**3", radius_range=(0, np.inf)
+            ),
+            PySDM_products.ParticleConcentration(
+                name="n_i", unit="1/cm**3", radius_range=(-np.inf, 0)
+            ),
+            PySDM_products.MeanRadius(name="r_s", unit="µm", radius_range=(0, np.inf)),
+            PySDM_products.MeanRadius(name="r_i", unit="µm", radius_range=(-np.inf, 0)),
+        ]
+
+        self.particulator = builder.build(attributes, products)
+
+        self.n_output = settings.n_output
+        self.n_substeps = int(self.n_output / dt)
+        self.t_max_duration = settings.t_max_duration
+
+    def save(self, output):
+        cell_id = 0
+
+        output["z"].append(self.particulator.products["z"].get()[cell_id])
+        output["t"].append(self.particulator.products["t"].get())
+        output["RH"].append(self.particulator.products["RH"].get()[cell_id])
+        output["RHi"].append(self.particulator.products["RH_ice"].get()[cell_id])
+        output["T"].append(self.particulator.products["T"].get()[cell_id])
+        output["P"].append(self.particulator.products["p"].get()[cell_id])
+        output["LWC"].append(self.particulator.products["water"].get()[cell_id])
+        output["IWC"].append(self.particulator.products["ice"].get()[cell_id])
+        #   output["TWC"].append(self.particulator.products["total"].get()[cell_id])
+        output["qv"].append(self.particulator.products["vapour"].get()[cell_id])
+        output["ns"].append(self.particulator.products["n_s"].get()[cell_id])
+        output["ni"].append(self.particulator.products["n_i"].get()[cell_id])
+        output["rs"].append(self.particulator.products["r_s"].get()[cell_id])
+        output["ri"].append(self.particulator.products["r_i"].get()[cell_id])
+        output["water_mass"].append(
+            self.particulator.attributes["signed water mass"].data.tolist()
+        )
+        output["T_frz"].append(
+            self.particulator.attributes["temperature of last freezing"].data.tolist()
+        )
+
+    def run(self):
+
+        print( "Starting simulation..." )
+
+        output = {
+            "t": [],
+            "z": [],
+            "RH": [],
+            "RHi": [],
+            "T": [],
+            "P": [],
+            "LWC": [],
+            "IWC": [],
+            "qv": [],
+            "ns": [],
+            "ni": [],
+            "rs": [],
+            "ri": [],
+            "water_mass": [],
+            "T_frz": [],
+        }
+
+        self.save(output)
+
+        while True:
+
+            self.particulator.run(self.n_substeps)
+            self.save(output)
+
+            # print( output["t"][-1], output["T"][-1], output["LWC"][-1],  output["IWC"][-1] )
+
+            if output["LWC"][-1] == 0:
+                print("break due to LWC")
+                break
+            if output["t"][-1] >= self.t_max_duration:
+                print("time exceeded")
+                break
+
+        return output