diff --git a/PySDM/backends/impl_common/freezing_attributes.py b/PySDM/backends/impl_common/freezing_attributes.py index d94f6d48d..4ff665a9b 100644 --- a/PySDM/backends/impl_common/freezing_attributes.py +++ b/PySDM/backends/impl_common/freezing_attributes.py @@ -19,3 +19,11 @@ class TimeDependentAttributes( """groups attributes required in time-dependent regime""" __slots__ = () + + +class TimeDependentHomogeneousAttributes( + namedtuple("TimeDependentHomogeneousAttributes", ("volume", "water_mass")) +): + """groups attributes required in time-dependent regime for homogeneous freezing""" + + __slots__ = () diff --git a/PySDM/backends/impl_numba/methods/freezing_methods.py b/PySDM/backends/impl_numba/methods/freezing_methods.py index 3a973d7bb..d2cc71ec0 100644 --- a/PySDM/backends/impl_numba/methods/freezing_methods.py +++ b/PySDM/backends/impl_numba/methods/freezing_methods.py @@ -1,5 +1,5 @@ """ -CPU implementation of backend methods for freezing (singular and time-dependent immersion freezing) +CPU implementation of backend methods for homogeneous freezing and heterogeneous freezing (singular and time-dependent immersion freezing) """ import numba @@ -10,6 +10,7 @@ from ...impl_common.freezing_attributes import ( SingularAttributes, TimeDependentAttributes, + TimeDependentHomogeneousAttributes, ) @@ -17,6 +18,7 @@ class FreezingMethods(BackendMethods): def __init__(self): BackendMethods.__init__(self) unfrozen_and_saturated = self.formulae.trivia.unfrozen_and_saturated + unfrozen_and_ice_saturated = self.formulae.trivia.unfrozen_and_ice_saturated frozen_and_above_freezing_point = ( self.formulae.trivia.frozen_and_above_freezing_point ) @@ -92,6 +94,48 @@ def freeze_time_dependent_body( # pylint: disable=unused-argument,too-many-argu self.freeze_time_dependent_body = freeze_time_dependent_body + + j_hom = self.formulae.homogeneous_ice_nucleation_rate.j_hom + + @numba.njit(**self.default_jit_flags) + def freeze_time_dependent_homogeneous_body( # pylint: disable=unused-argument,too-many-arguments + rand, + attributes, + timestep, + cell, + a_w_ice, + temperature, + relative_humidity_ice, + record_freezing_temperature, + freezing_temperature, + thaw, + ): + + n_sd = len(attributes.water_mass) + for i in numba.prange(n_sd): # pylint: disable=not-an-iterable + cell_id = cell[i] + #relative_humidity_ice = relative_humidity[cell_id] / a_w_ice[cell_id] + if thaw and frozen_and_above_freezing_point( + attributes.water_mass[i], temperature[cell_id] + ): + _thaw(attributes.water_mass, i) + elif unfrozen_and_ice_saturated( + attributes.water_mass[i], relative_humidity_ice[cell_id] + ): + d_a_w_ice = (relative_humidity_ice[cell_id] - 1) * a_w_ice[cell_id] + rate = j_hom(temperature, d_a_w_ice) + # TODO #594: this assumes constant T throughout timestep, can we do better? + prob = 1 - np.exp( # TODO #599: common code for Poissonian prob + -rate * attributes.volume[i] * timestep + ) + if rand[i] < prob: + _freeze(attributes.water_mass, i) + # if record_freezing_temperature: + # freezing_temperature[i] = temperature[cell_id] + + + self.freeze_time_dependent_homogeneous_body = freeze_time_dependent_homogeneous_body + def freeze_singular( self, *, attributes, temperature, relative_humidity, cell, thaw: bool ): @@ -106,6 +150,8 @@ def freeze_singular( thaw=thaw, ) + + def freeze_time_dependent( self, *, @@ -137,3 +183,38 @@ def freeze_time_dependent( ), thaw=thaw, ) + + + + def freeze_time_dependent_homogeneous( + self, + *, + rand, + attributes, + timestep, + cell, + a_w_ice, + temperature, + relative_humidity_ice, + record_freezing_temperature, + freezing_temperature, + thaw: bool + ): + self.freeze_time_dependent_homogeneous_body( + rand.data, + TimeDependentHomogeneousAttributes( + volume=attributes.volume.data, + water_mass=attributes.water_mass.data, + ), + timestep, + cell.data, + a_w_ice.data, + temperature.data, + relative_humidity_ice.data, + record_freezing_temperature=record_freezing_temperature, + freezing_temperature=( + freezing_temperature.data if record_freezing_temperature else None + ), + thaw=thaw, + ) + diff --git a/PySDM/dynamics/freezing.py b/PySDM/dynamics/freezing.py index a7d134244..2956e641d 100644 --- a/PySDM/dynamics/freezing.py +++ b/PySDM/dynamics/freezing.py @@ -1,10 +1,11 @@ """ -immersion freezing using either singular or time-dependent formulation +homogeneous freezing and heterogeneous freezing (singular and time-dependent immersion freezing) """ from PySDM.backends.impl_common.freezing_attributes import ( SingularAttributes, TimeDependentAttributes, + TimeDependentHomogeneousAttributes, ) from PySDM.physics.heterogeneous_ice_nucleation_rate import Null from PySDM.dynamics.impl import register_dynamic @@ -12,9 +13,11 @@ @register_dynamic() class Freezing: - def __init__(self, *, singular=True, record_freezing_temperature=False, thaw=False): + def __init__(self, *, singular=True, homogeneous_freezing=False, immersion_freezing=True, record_freezing_temperature=False, thaw=False): assert not (record_freezing_temperature and singular) self.singular = singular + self.homogeneous_freezing = homogeneous_freezing + self.immersion_freezing = immersion_freezing self.record_freezing_temperature = record_freezing_temperature self.thaw = thaw self.enable = True @@ -35,7 +38,8 @@ def register(self, builder): assert not isinstance( builder.formulae.heterogeneous_ice_nucleation_rate, Null ) - builder.request_attribute("immersed surface area") + if self.immersion_freezing: + builder.request_attribute("immersed surface area") self.rand = self.particulator.Storage.empty( self.particulator.n_sd, dtype=float ) @@ -43,6 +47,8 @@ def register(self, builder): self.particulator.n_sd, self.particulator.backend.formulae.seed ) + + def __call__(self): if "Coalescence" in self.particulator.dynamics: # TODO #594 @@ -53,34 +59,60 @@ def __call__(self): if not self.enable: return - if self.singular: - self.particulator.backend.freeze_singular( - attributes=SingularAttributes( - freezing_temperature=self.particulator.attributes[ - "freezing temperature" - ], - water_mass=self.particulator.attributes["water mass"], - ), - temperature=self.particulator.environment["T"], - relative_humidity=self.particulator.environment["RH"], - cell=self.particulator.attributes["cell id"], - thaw=self.thaw, - ) - else: - self.rand.urand(self.rng) - self.particulator.backend.freeze_time_dependent( - rand=self.rand, - attributes=TimeDependentAttributes( + if self.immersion_freezing: + + if self.singular: + self.particulator.backend.freeze_singular( + attributes=SingularAttributes( + freezing_temperature=self.particulator.attributes[ + "freezing temperature" + ], + water_mass=self.particulator.attributes["water mass"], + ), + temperature=self.particulator.environment["T"], + relative_humidity=self.particulator.environment["RH"], + cell=self.particulator.attributes["cell id"], + thaw=self.thaw, + ) + else: + self.rand.urand(self.rng) + self.particulator.backend.freeze_time_dependent( + rand=self.rand, + attributes=TimeDependentAttributes( immersed_surface_area=self.particulator.attributes[ "immersed surface area" ], water_mass=self.particulator.attributes["water mass"], + ), + timestep=self.particulator.dt, + cell=self.particulator.attributes["cell id"], + a_w_ice=self.particulator.environment["a_w_ice"], + temperature=self.particulator.environment["T"], + relative_humidity=self.particulator.environment["RH"], + record_freezing_temperature=self.record_freezing_temperature, + freezing_temperature=( + self.particulator.attributes["freezing temperature"] + if self.record_freezing_temperature + else None + ), + thaw=self.thaw, + ) + + + if self.homogeneous_freezing: + + self.rand.urand(self.rng) + self.particulator.backend.freeze_time_dependent_homogeneous( + rand=self.rand, + attributes=TimeDependentHomogeneousAttributes( + volume=self.particulator.attributes["volume"], + water_mass=self.particulator.attributes["water mass"], ), timestep=self.particulator.dt, cell=self.particulator.attributes["cell id"], a_w_ice=self.particulator.environment["a_w_ice"], temperature=self.particulator.environment["T"], - relative_humidity=self.particulator.environment["RH"], + relative_humidity_ice=self.particulator.environment["RH_ice"], record_freezing_temperature=self.record_freezing_temperature, freezing_temperature=( self.particulator.attributes["freezing temperature"] @@ -90,6 +122,8 @@ def __call__(self): thaw=self.thaw, ) + + self.particulator.attributes.mark_updated("water mass") if self.record_freezing_temperature: self.particulator.attributes.mark_updated("freezing temperature") diff --git a/PySDM/formulae.py b/PySDM/formulae.py index af8dc10c6..2e4f69588 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -44,6 +44,7 @@ def __init__( # pylint: disable=too-many-locals hydrostatics: str = "Default", freezing_temperature_spectrum: str = "Null", heterogeneous_ice_nucleation_rate: str = "Null", + homogeneous_ice_nucleation_rate: str = "Null", fragmentation_function: str = "AlwaysN", isotope_equilibrium_fractionation_factors: str = "Null", isotope_meteoric_water_line_excess: str = "Null", @@ -76,6 +77,7 @@ def __init__( # pylint: disable=too-many-locals self.hydrostatics = hydrostatics self.freezing_temperature_spectrum = freezing_temperature_spectrum self.heterogeneous_ice_nucleation_rate = heterogeneous_ice_nucleation_rate + self.homogeneous_ice_nucleation_rate = homogeneous_ice_nucleation_rate self.fragmentation_function = fragmentation_function self.isotope_equilibrium_fractionation_factors = ( isotope_equilibrium_fractionation_factors diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py index 93b1c4f3d..d907e0919 100644 --- a/PySDM/physics/__init__.py +++ b/PySDM/physics/__init__.py @@ -25,6 +25,7 @@ fragmentation_function, freezing_temperature_spectrum, heterogeneous_ice_nucleation_rate, + homogeneous_ice_nucleation_rate, hydrostatics, hygroscopicity, impl, diff --git a/PySDM/physics/constants_defaults.py b/PySDM/physics/constants_defaults.py index 439738ce8..6dfcd6f9d 100644 --- a/PySDM/physics/constants_defaults.py +++ b/PySDM/physics/constants_defaults.py @@ -185,7 +185,26 @@ ABIFM_C = np.inf ABIFM_UNIT = 1 / si.cm**2 / si.s +KOOP_2000_C1 = -906.7 +KOOP_2000_C2 = 8502 +KOOP_2000_C3 = 26924 +KOOP_2000_C4 = 29180 +KOOP_UNIT = 1 / si.cm**3 / si.s + +KOOP_CORR = -1.522 + +KOOP_MURRAY_C0 = -3020.684 +KOOP_MURRAY_C1 = -425.921 +KOOP_MURRAY_C2 = -25.9779 +KOOP_MURRAY_C3 = -0.868451 +KOOP_MURRAY_C4 = -1.66203e-2 +KOOP_MURRAY_C5 = -1.71736e-4 +KOOP_MURRAY_C6 = -7.46953e-7 + +T_hom_freeze = 235. * si.kelvin + J_HET = np.nan +J_HOM = np.nan STRAUB_E_D1 = 0.04 * si.cm STRAUB_MU2 = 0.095 * si.cm diff --git a/PySDM/physics/homogeneous_ice_nucleation_rate/__init__.py b/PySDM/physics/homogeneous_ice_nucleation_rate/__init__.py new file mode 100644 index 000000000..3cb1eabeb --- /dev/null +++ b/PySDM/physics/homogeneous_ice_nucleation_rate/__init__.py @@ -0,0 +1,7 @@ +""" +homogeneous-freezing rate (aka J_hom) formulations +""" + +from .constant import Constant +from .null import Null +from .koop import KOOP diff --git a/PySDM/physics/homogeneous_ice_nucleation_rate/constant.py b/PySDM/physics/homogeneous_ice_nucleation_rate/constant.py new file mode 100644 index 000000000..f5f32353d --- /dev/null +++ b/PySDM/physics/homogeneous_ice_nucleation_rate/constant.py @@ -0,0 +1,14 @@ +""" +constant rate formulation (for tests) +""" + +import numpy as np + + +class Constant: # pylint: disable=too-few-public-methods + def __init__(self, const): + assert np.isfinite(const.J_HOM) + + @staticmethod + def j_hom(const, T, a_w_ice): # pylint: disable=unused-argument + return const.J_HOM diff --git a/PySDM/physics/homogeneous_ice_nucleation_rate/koop.py b/PySDM/physics/homogeneous_ice_nucleation_rate/koop.py new file mode 100644 index 000000000..144b42edd --- /dev/null +++ b/PySDM/physics/homogeneous_ice_nucleation_rate/koop.py @@ -0,0 +1,14 @@ +""" +Koop homogeneous nucleation rate parameterization for solution droplets +valid for 0.26 < da_w_ice < 0.34 + ([Koop et al. 2000](https://doi.org/10.1038/35020537)) +""" + + +class KOOP: # pylint: disable=too-few-public-methods + def __init__(self, const): + pass + + @staticmethod + def j_hom(const, T, da_w_ice): + return 10**(const.KOOP_2000_C1 + const.KOOP_2000_C2 * da_w_ice - const.KOOP_2000_C3 * da_w_ice**2. + const.KOOP_2000_C4 * da_w_ice**3.) * const.KOOP_UNIT diff --git a/PySDM/physics/homogeneous_ice_nucleation_rate/koop_corr.py b/PySDM/physics/homogeneous_ice_nucleation_rate/koop_corr.py new file mode 100644 index 000000000..145370a56 --- /dev/null +++ b/PySDM/physics/homogeneous_ice_nucleation_rate/koop_corr.py @@ -0,0 +1,15 @@ +""" +Koop homogeneous nucleation rate parameterization for solution droplets [Koop et al. 2000] corrected +such that it coincides with homogeneous nucleation rate parameterization for pure water droplets +[Koop and Murray 2016] at water saturation between 235K < T < 240K + ([Spichtinger et al. 2023](https://doi.org/10.5194/acp-23-2035-2023)) +""" + + +class KOOP_CORR: # pylint: disable=too-few-public-methods + def __init__(self, const): + pass + + @staticmethod + def j_hom(const, T, da_w_ice): + return 10**(const.KOOP_2000_C1 + const.KOOP_2000_C2 * da_w_ice - const.KOOP_2000_C3 * da_w_ice**2. + const.KOOP_2000_C4 * da_w_ice**3. + const.KOOP_CORR) * const.KOOP_UNIT diff --git a/PySDM/physics/homogeneous_ice_nucleation_rate/koop_murray.py b/PySDM/physics/homogeneous_ice_nucleation_rate/koop_murray.py new file mode 100644 index 000000000..0b5b4ce33 --- /dev/null +++ b/PySDM/physics/homogeneous_ice_nucleation_rate/koop_murray.py @@ -0,0 +1,15 @@ +""" +Koop and Murray homogeneous nucleation rate parameterization for pure water droplets +at water saturation + ([Koop and Murray 2016](https://doi.org/10.1063/1.4962355)) +""" + + +class KOOP_MURRAY: # pylint: disable=too-few-public-methods + def __init__(self, const): + pass + + @staticmethod + def j_hom(const, T, da_w_ice): + T_diff = T - const.T_tri + return 10**(const.KOOP_MURRAY_C0 + const.KOOP_MURRAY_C1*T_diff + const.KOOP_MURRAY_C2*T_diff**2. + const.KOOP_MURRAY_C3*T_diff**3. + const.KOOP_MURRAY_C4*T_diff**4. + const.KOOP_MURRAY_C5*T_diff**5. + const.KOOP_MURRAY_C6*T_diff**6.) * const.KOOP_UNIT diff --git a/PySDM/physics/homogeneous_ice_nucleation_rate/null.py b/PySDM/physics/homogeneous_ice_nucleation_rate/null.py new file mode 100644 index 000000000..81ceff95b --- /dev/null +++ b/PySDM/physics/homogeneous_ice_nucleation_rate/null.py @@ -0,0 +1,13 @@ +""" +do-nothing null formulation (needed as other formulations require parameters + to be set before instantiation of Formulae) +""" + + +class Null: # pylint: disable=too-few-public-methods,unused-argument + def __init__(self, _): + pass + + @staticmethod + def j_hom(const, T, d_a_w_ice): # pylint: disable=unused-argument + return 0 diff --git a/PySDM/physics/trivia.py b/PySDM/physics/trivia.py index 3bac176aa..7331a1047 100644 --- a/PySDM/physics/trivia.py +++ b/PySDM/physics/trivia.py @@ -79,6 +79,10 @@ def th_std(const, p, T): def unfrozen_and_saturated(water_mass, relative_humidity): return water_mass > 0 and relative_humidity > 1 + @staticmethod + def unfrozen_and_ice_saturated(water_mass, relative_humidity_ice): + return water_mass > 0 and relative_humidity_ice > 1 + @staticmethod def frozen_and_above_freezing_point(const, water_mass, temperature): return water_mass < 0 and temperature > const.T0 diff --git a/tests/unit_tests/backends/test_freezing_methods.py b/tests/unit_tests/backends/test_freezing_methods.py index 719a738eb..5b8484d3d 100644 --- a/tests/unit_tests/backends/test_freezing_methods.py +++ b/tests/unit_tests/backends/test_freezing_methods.py @@ -208,6 +208,128 @@ def low(t): np.testing.assert_array_less(low(arg), data) + @staticmethod + @pytest.mark.parametrize("double_precision", (True,False)) + # pylint: disable=too-many-locals + def test_homogeneous_freeze_time_dependent(backend_class, double_precision, plot=False): + if backend_class.__name__ == "Numba" and not double_precision: + pytest.skip() + if backend_class.__name__ == "ThrustRTC": + pytest.skip() + + + + # Arrange + seed = 44 + cases = ( + {"dt": 5e5, "N": 1}, + {"dt": 1e6, "N": 1}, + {"dt": 5e5, "N": 8}, + {"dt": 1e6, "N": 8}, + {"dt": 5e5, "N": 16}, + {"dt": 1e6, "N": 16}, + ) + rate = 1e-9 + + number_of_real_droplets = 1024 + total_time = ( + 0.25e9 # effectively interpreted here as seconds, i.e. cycle = 1 * si.s + ) + + # dummy (but must-be-set) values + initial_water_mass = ( + 1000 # for sign flip (ice water has negative volumes) + ) + d_v = 666 # products use conc., dividing there, multiplying here, value does not matter + + droplet_volume = initial_water_mass / 1000. + + def hgh(t): + return np.exp(-0.75 * rate * (t - total_time / 4)) + + def low(t): + return np.exp(-1.25 * rate * (t + total_time / 4)) + + + RHi = 1.5 + T = 230 + + # Act + output = {} + + formulae = Formulae( + particle_shape_and_density="MixedPhaseSpheres", + homogeneous_ice_nucleation_rate="Constant", + constants={"J_HOM": rate / droplet_volume}, + seed=seed, + ) + products = (IceWaterContent(name="qi"),) + + for case in cases: + n_sd = int(number_of_real_droplets // case["N"]) + assert n_sd == number_of_real_droplets / case["N"] + assert total_time // case["dt"] == total_time / case["dt"] + + key = f"{case['dt']}:{case['N']}" + output[key] = {"unfrozen_fraction": [], "dt": case["dt"], "N": case["N"]} + + env = Box(dt=case["dt"], dv=d_v) + builder = Builder( + n_sd=n_sd, + backend=backend_class( + formulae=formulae, double_precision=double_precision + ), + environment=env, + ) + builder.add_dynamic(Freezing(singular=False,homogeneous_freezing=True,immersion_freezing=False)) + attributes = { + "multiplicity": np.full(n_sd, int(case["N"])), + "water mass": np.full(n_sd, initial_water_mass), + } + particulator = builder.build(attributes=attributes, products=products) + pvs_ice = particulator.formulae.saturation_vapour_pressure.pvs_ice(T) + pvs_water = particulator.formulae.saturation_vapour_pressure.pvs_water(T) + particulator.environment["RH_ice"] = RHi + particulator.environment["a_w_ice"] = pvs_ice / pvs_water + particulator.environment["T"] = T + + cell_id = 0 + for i in range(int(total_time / case["dt"]) + 1): + particulator.run(0 if i == 0 else 1) + + ice_mass_per_volume = particulator.products["qi"].get()[cell_id] + ice_mass = ice_mass_per_volume * d_v + ice_number = ice_mass / initial_water_mass + unfrozen_fraction = 1 - ice_number / number_of_real_droplets + output[key]["unfrozen_fraction"].append(unfrozen_fraction) + + + # Plot + fit_x = np.linspace(0, total_time, num=100) + fit_y = np.exp(-rate * fit_x) + + for out in output.values(): + pyplot.step( + out["dt"] * np.arange(len(out["unfrozen_fraction"])), + out["unfrozen_fraction"], + label=f"dt={out['dt']:.2g} / N={out['N']}", + marker=".", + linewidth=1 + out["N"] // 8, + ) + + _plot_fit(fit_x, fit_y, low, hgh, total_time) + if plot: + pyplot.show() + + # Assert + for out in output.values(): + data = np.asarray(out["unfrozen_fraction"]) + arg = out["dt"] * np.arange(len(data)) + np.testing.assert_array_less(data, hgh(arg)) + np.testing.assert_array_less(low(arg), data) + + + def _plot_fit(fit_x, fit_y, low, hgh, total_time): pyplot.plot( fit_x, fit_y, color="black", linestyle="--", label="theory", linewidth=5