From a29b033b426caa351e33537c4351e6d7eb731e06 Mon Sep 17 00:00:00 2001 From: sfonxu Date: Thu, 23 Jan 2025 12:40:15 +0100 Subject: [PATCH 1/2] Implementation of Twomey (1959) eq'n for CCN activation spectra --- .../ccn_activation_spectrum/__init__.py | 2 ++ .../ccn_activation_spectrum/twomey_1959.py | 15 ++++++++++ .../physics/test_ccn_activation_spectrum.py | 30 +++++++++++++++++++ 3 files changed, 47 insertions(+) create mode 100644 PySDM/physics/ccn_activation_spectrum/__init__.py create mode 100644 PySDM/physics/ccn_activation_spectrum/twomey_1959.py create mode 100644 tests/unit_tests/physics/test_ccn_activation_spectrum.py diff --git a/PySDM/physics/ccn_activation_spectrum/__init__.py b/PySDM/physics/ccn_activation_spectrum/__init__.py new file mode 100644 index 0000000000..237fc2a9b0 --- /dev/null +++ b/PySDM/physics/ccn_activation_spectrum/__init__.py @@ -0,0 +1,2 @@ +from PySDM.impl.null_physics_class import Null +from .twomey_1959 import Twomey1959 diff --git a/PySDM/physics/ccn_activation_spectrum/twomey_1959.py b/PySDM/physics/ccn_activation_spectrum/twomey_1959.py new file mode 100644 index 0000000000..5560fff7ac --- /dev/null +++ b/PySDM/physics/ccn_activation_spectrum/twomey_1959.py @@ -0,0 +1,15 @@ +import numpy as np +""" +[Twomey 1959](https://doi.org/10.1007/BF01993560). Note that supersaturation is expressed in temperature unit as the elevation of the dew point or as percent of supersaturation; and concentrations are reported for 10C and 800 mb. +""" + +class Twomey1959: + def __init__(self, const): + assert np.isfinite(const.TWOMEY_K) + assert np.isfinite(const.TWOMEY_N0) + + @staticmethod + def ccn_concentration(const, saturation_ratio): + return const.TWOMEY_N0*np.power(saturation_ratio-1, const.TWOMEY_K) + + diff --git a/tests/unit_tests/physics/test_ccn_activation_spectrum.py b/tests/unit_tests/physics/test_ccn_activation_spectrum.py new file mode 100644 index 0000000000..d68fb23671 --- /dev/null +++ b/tests/unit_tests/physics/test_ccn_activation_spectrum.py @@ -0,0 +1,30 @@ +import numpy as np +from matplotlib import pyplot + +from PySDM import Formulae +from PySDM.physics.constants import PER_CENT, si, in_unit + +def test_twomey_and_wojciechowski_1969_fig1(plot=True): + """[Twomey, Wojciechowski 1969](https://doi.org/10.1175/1520-0469(1969)26%3C648:OOTGVO%3E2.0.CO;2) + """ + #arange + for k, N in zip([0.5,0.7], [100,500]): + formulae = Formulae(ccn_activation_spectrum="Twomey1959",constants={"TWOMEY_K": k, "TWOMEY_N0": N/si.cm**3}) + supersaturation = np.logspace(np.log10(.2), np.log10(9))*PER_CENT + #act + activated_nuclei_concentration = formulae.ccn_activation_spectrum.ccn_concentration(saturation_ratio=supersaturation+1) + #plot + pyplot.plot(in_unit(supersaturation,PER_CENT), + in_unit(activated_nuclei_concentration,si.cm**-3), + label=f"{k=}" + ) + pyplot.xlim(0.1,10) + pyplot.ylim(1,1000) + pyplot.xscale("log") + pyplot.yscale("log") + pyplot.xlabel("Percent supersaturation") + pyplot.ylabel("Nuclei [cm$^{-3}$]") + pyplot.grid() + pyplot.legend() + pyplot.show() + #assert From 72ccc3f0ec4d6fd1e63cc2669a87108ca815e1df Mon Sep 17 00:00:00 2001 From: sfonxu Date: Thu, 23 Jan 2025 12:43:48 +0100 Subject: [PATCH 2/2] formulae and constants changes --- PySDM/formulae.py | 4 +++- PySDM/physics/__init__.py | 1 + PySDM/physics/constants_defaults.py | 2 ++ 3 files changed, 6 insertions(+), 1 deletion(-) diff --git a/PySDM/formulae.py b/PySDM/formulae.py index af8dc10c62..a3cd8af6df 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -56,6 +56,7 @@ def __init__( # pylint: disable=too-many-locals terminal_velocity: str = "GunnKinzer1949", air_dynamic_viscosity: str = "ZografosEtAl1987", bulk_phase_partitioning: str = "Null", + ccn_activation_spectrum: str = "Null", handle_all_breakups: bool = False, ): # initialisation of the fields below is just to silence pylint and to enable code hints @@ -88,7 +89,8 @@ def __init__( # pylint: disable=too-many-locals self.air_dynamic_viscosity = air_dynamic_viscosity self.terminal_velocity = terminal_velocity self.bulk_phase_partitioning = bulk_phase_partitioning - + self.ccn_activation_spectrum = ccn_activation_spectrum + self._components = tuple( i for i in dir(self) diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py index 93b1c4f3d0..f251794d28 100644 --- a/PySDM/physics/__init__.py +++ b/PySDM/physics/__init__.py @@ -46,5 +46,6 @@ air_dynamic_viscosity, terminal_velocity, bulk_phase_partitioning, + ccn_activation_spectrum, ) from .constants import convert_to, in_unit, si diff --git a/PySDM/physics/constants_defaults.py b/PySDM/physics/constants_defaults.py index aee675d08f..45e178e769 100644 --- a/PySDM/physics/constants_defaults.py +++ b/PySDM/physics/constants_defaults.py @@ -437,3 +437,5 @@ def compute_derived_values(c: dict): isotopologue; in the paper timescale is calculated for tritium with assumption of no tritium in the environment around the drop (Table 1). """ +TWOMEY_K = np.nan +TWOMEY_N0 = np.nan