From bdc1fd858e37bd33db0066cde287727e3eb617e5 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Tue, 18 Mar 2025 16:06:23 +0100 Subject: [PATCH 01/12] modified sedimentation unit test to test for ice particles --- .../dynamics/displacement/displacement_settings.py | 4 ++-- .../dynamics/displacement/test_sedimentation.py | 11 +++++++---- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/tests/unit_tests/dynamics/displacement/displacement_settings.py b/tests/unit_tests/dynamics/displacement/displacement_settings.py index e16d481ac..3d9ad159b 100644 --- a/tests/unit_tests/dynamics/displacement/displacement_settings.py +++ b/tests/unit_tests/dynamics/displacement/displacement_settings.py @@ -9,9 +9,9 @@ class DisplacementSettings: # pylint: disable=too-few-public-methods - def __init__(self, n_sd=1, grid=None, positions=None, courant_field_data=None): + def __init__(self, n_sd=1, volume=None, grid=None, positions=None, courant_field_data=None): self.n = np.ones(n_sd, dtype=np.int64) - self.volume = np.ones(n_sd, dtype=np.float64) + self.volume = volume or np.ones(n_sd, dtype=np.float64) self.grid = grid or (1, 1) self.courant_field_data = courant_field_data or ( np.array([[0, 0]]).T, diff --git a/tests/unit_tests/dynamics/displacement/test_sedimentation.py b/tests/unit_tests/dynamics/displacement/test_sedimentation.py index b1c96dd29..5cfaad7ab 100644 --- a/tests/unit_tests/dynamics/displacement/test_sedimentation.py +++ b/tests/unit_tests/dynamics/displacement/test_sedimentation.py @@ -1,6 +1,6 @@ # pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring,no-member import numpy as np - +import pytest from .displacement_settings import DisplacementSettings @@ -11,12 +11,15 @@ def __init__(self, backend, particles): def get(self): return self.values - +VOLUMES = (1., -1.) class TestSedimentation: # pylint: disable=too-few-public-methods @staticmethod - def test_boundary_condition(backend_class): + @pytest.mark.parametrize( + "volume", [np.asarray((v,)) for v in VOLUMES] + ) + def test_boundary_condition(backend_class, volume): # Arrange - settings = DisplacementSettings() + settings = DisplacementSettings(n_sd=len(volume), volume=volume) settings.dt = 1 settings.sedimentation = True sut, particulator = settings.get_displacement( From 3745b928ab5ea5250a817f9e01e227c4cb6173fe Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Mon, 24 Mar 2025 17:27:06 +0100 Subject: [PATCH 02/12] added empirical parametrisation for columnar ice crystals --- PySDM/dynamics/terminal_velocity/__init__.py | 1 + .../terminal_velocity/columnar_ice_crystal.py | 16 +++++++++++ PySDM/physics/constants_defaults.py | 26 +++++++++++++++++ PySDM/physics/terminal_velocity/__init__.py | 1 + .../terminal_velocity/columnar_ice_crystal.py | 28 +++++++++++++++++++ 5 files changed, 72 insertions(+) create mode 100644 PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py create mode 100644 PySDM/physics/terminal_velocity/columnar_ice_crystal.py diff --git a/PySDM/dynamics/terminal_velocity/__init__.py b/PySDM/dynamics/terminal_velocity/__init__.py index c94a06bca..16b5119a6 100644 --- a/PySDM/dynamics/terminal_velocity/__init__.py +++ b/PySDM/dynamics/terminal_velocity/__init__.py @@ -5,3 +5,4 @@ from PySDM.dynamics.terminal_velocity.gunn_and_kinzer import GunnKinzer1949 from PySDM.dynamics.terminal_velocity.power_series import PowerSeries from PySDM.dynamics.terminal_velocity.rogers_and_yau import RogersYau +from PySDM.dynamics.terminal_velocity.columnar_ice_crystal import columnar_ice_crystal \ No newline at end of file diff --git a/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py b/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py new file mode 100644 index 000000000..6c6fc4305 --- /dev/null +++ b/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py @@ -0,0 +1,16 @@ +""" +[Spichtinger & Gierens 2009](https://doi.org/10.5194/acp-9-685-2009) +Eq. (18) .Assumed shape is columnar based on empirical parameterizations of +[Heymsfield & Iaquinta (2000)](https://doi.org/10.1175/1520-0469(2000)057<0916:CCTV>2.0.CO;2) +[Barthazy & Schefold (2006)](https://doi.org/10.1016/j.atmosres.2005.12.009) +""" + +class columnar_ice_crystal: # pylint: disable=too-few-public-methods + def __init__(self, particulator): + self.particulator = particulator + + def __call__(self, output, mass): + self.particulator.backend.terminal_velocity( + values=output.data, + mass=mass.data, + ) \ No newline at end of file diff --git a/PySDM/physics/constants_defaults.py b/PySDM/physics/constants_defaults.py index 5247d27b2..46f8d10ad 100644 --- a/PySDM/physics/constants_defaults.py +++ b/PySDM/physics/constants_defaults.py @@ -554,6 +554,32 @@ ROGERS_YAU_TERM_VEL_MEDIUM_R_LIMIT = 600 * si.um """ 〃 """ +SPICHTINGER_GIERENS_TERM_VEL_LIMIT_0 = 2.146e-13 * si.kg +""" empirical terminal velocity formulation +for columnar ice crystals from Table 2. in +[Spichtinger & Gierens 2009](https://doi.org/10.5194/acp-9-685-2009) """ +SPICHTINGER_GIERENS_TERM_VEL_LIMIT_1 = 2.166e-9 * si.kg +""" 〃 """ +SPICHTINGER_GIERENS_TERM_VEL_LIMIT_2 = 4.264e-8 * si.kg +""" 〃 """ +SPICHTINGER_TERM_DELTA_COEFF0 = 0.42 +""" 〃 """ +SPICHTINGER_TERM_DELTA_COEFF1 = 0.57 +""" 〃 """ +SPICHTINGER_TERM_DELTA_COEFF2 = 0.31 +""" 〃 """ +SPICHTINGER_TERM_DELTA_COEFF3 = 0.096 +""" 〃 """ +SPICHTINGER_TERM_GAMMA_COEFF0 = 735.4 * si.m / si.s / si.kg**SPICHTINGER_TERM_DELTA_COEFF0 +""" 〃 """ +SPICHTINGER_TERM_GAMMA_COEFF1 = 63292.4 * si.m / si.s / si.kg**SPICHTINGER_TERM_DELTA_COEFF1 +""" 〃 """ +SPICHTINGER_TERM_GAMMA_COEFF2 = 329.8 * si.m / si.s / si.kg**SPICHTINGER_TERM_DELTA_COEFF2 +""" 〃 """ +SPICHTINGER_TERM_GAMMA_COEFF3 = 8.8 * si.m / si.s / si.kg**SPICHTINGER_TERM_DELTA_COEFF3 +""" 〃 """ + + W76W_G0 = -2.9912729e3 * si.K**2 """ [Wexler 1976](https://doi.org/10.6028/jres.080A.071) saturation vapour pressure """ W76W_G1 = -6.0170128e3 * si.K diff --git a/PySDM/physics/terminal_velocity/__init__.py b/PySDM/physics/terminal_velocity/__init__.py index 2e012dc47..5c4a59bc0 100644 --- a/PySDM/physics/terminal_velocity/__init__.py +++ b/PySDM/physics/terminal_velocity/__init__.py @@ -3,3 +3,4 @@ from .rogers_yau import RogersYau from .gunn_kinzer_1949 import GunnKinzer1949 from .power_series import PowerSeries +from .columnar_ice_crystal import columnar_ice_crystal \ No newline at end of file diff --git a/PySDM/physics/terminal_velocity/columnar_ice_crystal.py b/PySDM/physics/terminal_velocity/columnar_ice_crystal.py new file mode 100644 index 000000000..a36636edd --- /dev/null +++ b/PySDM/physics/terminal_velocity/columnar_ice_crystal.py @@ -0,0 +1,28 @@ +""" +[Spichtinger & Gierens 2009](https://doi.org/10.5194/acp-9-685-2009) +Eq. (18) .Assumed shape is columnar based on empirical parameterizations of +[Heymsfield & Iaquinta (2000)](https://doi.org/10.1175/1520-0469(2000)057<0916:CCTV>2.0.CO;2) +[Barthazy & Schefold (2006)](https://doi.org/10.1016/j.atmosres.2005.12.009) +""" + +import numpy as np + + +class columnar_ice_crystal: # pylint: disable=too-few-public-methods + def __init__(self, _): + pass + + @staticmethod + def v_term(const, mass): + return np.where( + mass < const.SPICHTINGER_GIERENS_TERM_VEL_LIMIT_0, + const.SPICHTINGER_TERM_GAMMA_COEFF0 * mass**const.SPICHTINGER_TERM_DELTA_COEFF0, + np.where( + mass < const.SPICHTINGER_GIERENS_TERM_VEL_LIMIT_1, + const.SPICHTINGER_TERM_GAMMA_COEFF1 * mass ** const.SPICHTINGER_TERM_DELTA_COEFF1, + np.where( mass < const.SPICHTINGER_GIERENS_TERM_VEL_LIMIT_2, + const.SPICHTINGER_TERM_GAMMA_COEFF2 * mass ** const.SPICHTINGER_TERM_DELTA_COEFF2, + const.SPICHTINGER_TERM_GAMMA_COEFF3 * mass ** const.SPICHTINGER_TERM_DELTA_COEFF3 + ) + ), + ) From e7f8feeeaf1961f2191b2ab57adb2a8c18ece54f Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Fri, 28 Mar 2025 15:41:39 +0100 Subject: [PATCH 03/12] added atmospheric correction factor for columnar ice crystals terminal velocity --- PySDM/physics/constants_defaults.py | 9 ++++++++- .../terminal_velocity/columnar_ice_crystal.py | 14 ++++++++++---- 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/PySDM/physics/constants_defaults.py b/PySDM/physics/constants_defaults.py index 46f8d10ad..05c74cdab 100644 --- a/PySDM/physics/constants_defaults.py +++ b/PySDM/physics/constants_defaults.py @@ -578,7 +578,14 @@ """ 〃 """ SPICHTINGER_TERM_GAMMA_COEFF3 = 8.8 * si.m / si.s / si.kg**SPICHTINGER_TERM_DELTA_COEFF3 """ 〃 """ - +SPICHTINGER_CORRECTION_P0 = 300 * si.hectopascal +""" 〃 """ +SPICHTINGER_CORRECTION_P_EXPO = -0.178 +""" 〃 """ +SPICHTINGER_CORRECTION_T0 = 233 * si.kelvin +""" 〃 """ +SPICHTINGER_CORRECTION_T_EXPO = -0.394 +""" 〃 """ W76W_G0 = -2.9912729e3 * si.K**2 """ [Wexler 1976](https://doi.org/10.6028/jres.080A.071) saturation vapour pressure """ diff --git a/PySDM/physics/terminal_velocity/columnar_ice_crystal.py b/PySDM/physics/terminal_velocity/columnar_ice_crystal.py index a36636edd..bbe3749e0 100644 --- a/PySDM/physics/terminal_velocity/columnar_ice_crystal.py +++ b/PySDM/physics/terminal_velocity/columnar_ice_crystal.py @@ -1,6 +1,6 @@ """ [Spichtinger & Gierens 2009](https://doi.org/10.5194/acp-9-685-2009) -Eq. (18) .Assumed shape is columnar based on empirical parameterizations of +Eq. (18) and (19) Assumed shape is columnar based on empirical parameterizations of [Heymsfield & Iaquinta (2000)](https://doi.org/10.1175/1520-0469(2000)057<0916:CCTV>2.0.CO;2) [Barthazy & Schefold (2006)](https://doi.org/10.1016/j.atmosres.2005.12.009) """ @@ -19,10 +19,16 @@ def v_term(const, mass): const.SPICHTINGER_TERM_GAMMA_COEFF0 * mass**const.SPICHTINGER_TERM_DELTA_COEFF0, np.where( mass < const.SPICHTINGER_GIERENS_TERM_VEL_LIMIT_1, - const.SPICHTINGER_TERM_GAMMA_COEFF1 * mass ** const.SPICHTINGER_TERM_DELTA_COEFF1, + const.SPICHTINGER_TERM_GAMMA_COEFF1 * mass**const.SPICHTINGER_TERM_DELTA_COEFF1, np.where( mass < const.SPICHTINGER_GIERENS_TERM_VEL_LIMIT_2, - const.SPICHTINGER_TERM_GAMMA_COEFF2 * mass ** const.SPICHTINGER_TERM_DELTA_COEFF2, - const.SPICHTINGER_TERM_GAMMA_COEFF3 * mass ** const.SPICHTINGER_TERM_DELTA_COEFF3 + const.SPICHTINGER_TERM_GAMMA_COEFF2 * mass**const.SPICHTINGER_TERM_DELTA_COEFF2, + const.SPICHTINGER_TERM_GAMMA_COEFF3 * mass**const.SPICHTINGER_TERM_DELTA_COEFF3 ) ), ) + + @staticmethod + def atmospheric_correction_factor(const, temperature, pressure): + return( (pressure / const.SPICHTINGER_CORRECTION_P0)**const.SPICHTINGER_CORRECTION_P_EXPO * + (const.SPICHTINGER_CORRECTION_T0 / temperature)**const.SPICHTINGER_CORRECTION_T_EXPO + ) \ No newline at end of file From 4901282bdafaa01a0f56c2de4f5d151c8052a9ae Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Fri, 28 Mar 2025 16:27:25 +0100 Subject: [PATCH 04/12] name change of columnar ice crystal sedimentation class --- PySDM/dynamics/terminal_velocity/__init__.py | 2 +- PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py | 2 +- PySDM/formulae.py | 3 ++- PySDM/physics/terminal_velocity/__init__.py | 2 +- PySDM/physics/terminal_velocity/columnar_ice_crystal.py | 2 +- 5 files changed, 6 insertions(+), 5 deletions(-) diff --git a/PySDM/dynamics/terminal_velocity/__init__.py b/PySDM/dynamics/terminal_velocity/__init__.py index 16b5119a6..17d0bcf65 100644 --- a/PySDM/dynamics/terminal_velocity/__init__.py +++ b/PySDM/dynamics/terminal_velocity/__init__.py @@ -5,4 +5,4 @@ from PySDM.dynamics.terminal_velocity.gunn_and_kinzer import GunnKinzer1949 from PySDM.dynamics.terminal_velocity.power_series import PowerSeries from PySDM.dynamics.terminal_velocity.rogers_and_yau import RogersYau -from PySDM.dynamics.terminal_velocity.columnar_ice_crystal import columnar_ice_crystal \ No newline at end of file +from PySDM.dynamics.terminal_velocity.columnar_ice_crystal import ColumnarIceCrystal \ No newline at end of file diff --git a/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py b/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py index 6c6fc4305..92320c94f 100644 --- a/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py +++ b/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py @@ -5,7 +5,7 @@ [Barthazy & Schefold (2006)](https://doi.org/10.1016/j.atmosres.2005.12.009) """ -class columnar_ice_crystal: # pylint: disable=too-few-public-methods +class ColumnarIceCrystal: # pylint: disable=too-few-public-methods def __init__(self, particulator): self.particulator = particulator diff --git a/PySDM/formulae.py b/PySDM/formulae.py index 853000eb7..7329ca86f 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -19,7 +19,7 @@ from PySDM import physics from PySDM.backends.impl_numba import conf -from PySDM.dynamics.terminal_velocity import GunnKinzer1949, PowerSeries, RogersYau +from PySDM.dynamics.terminal_velocity import GunnKinzer1949, PowerSeries, RogersYau, ColumnarIceCrystal from PySDM.dynamics.terminal_velocity.gunn_and_kinzer import TpDependent @@ -55,6 +55,7 @@ def __init__( # pylint: disable=too-many-locals optical_depth: str = "Null", particle_shape_and_density: str = "LiquidSpheres", terminal_velocity: str = "GunnKinzer1949", + terminal_velocity_ice: str = "ColumnarIceCrystal", air_dynamic_viscosity: str = "ZografosEtAl1987", bulk_phase_partitioning: str = "Null", handle_all_breakups: bool = False, diff --git a/PySDM/physics/terminal_velocity/__init__.py b/PySDM/physics/terminal_velocity/__init__.py index 5c4a59bc0..5dad6a8f7 100644 --- a/PySDM/physics/terminal_velocity/__init__.py +++ b/PySDM/physics/terminal_velocity/__init__.py @@ -3,4 +3,4 @@ from .rogers_yau import RogersYau from .gunn_kinzer_1949 import GunnKinzer1949 from .power_series import PowerSeries -from .columnar_ice_crystal import columnar_ice_crystal \ No newline at end of file +from .columnar_ice_crystal import ColumnarIceCrystal \ No newline at end of file diff --git a/PySDM/physics/terminal_velocity/columnar_ice_crystal.py b/PySDM/physics/terminal_velocity/columnar_ice_crystal.py index bbe3749e0..d8487cb42 100644 --- a/PySDM/physics/terminal_velocity/columnar_ice_crystal.py +++ b/PySDM/physics/terminal_velocity/columnar_ice_crystal.py @@ -8,7 +8,7 @@ import numpy as np -class columnar_ice_crystal: # pylint: disable=too-few-public-methods +class ColumnarIceCrystal: # pylint: disable=too-few-public-methods def __init__(self, _): pass From 599c40853320557fd39b7856b831c9059de571f7 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Fri, 28 Mar 2025 17:02:27 +0100 Subject: [PATCH 05/12] changed mass attribute to signed mass attribute for relative fall velocity --- PySDM/attributes/physics/relative_fall_velocity.py | 8 ++++---- PySDM/dynamics/relaxed_velocity.py | 2 +- tests/unit_tests/attributes/test_fall_velocity.py | 6 +++--- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/PySDM/attributes/physics/relative_fall_velocity.py b/PySDM/attributes/physics/relative_fall_velocity.py index a5e170c00..d047cad4b 100644 --- a/PySDM/attributes/physics/relative_fall_velocity.py +++ b/PySDM/attributes/physics/relative_fall_velocity.py @@ -1,5 +1,5 @@ """ -Attributes for tracking droplet velocity +Attributes for tracking particle velocity """ from PySDM.attributes.impl import ( @@ -29,13 +29,13 @@ def __init__(self, builder): class RelativeFallVelocity(DerivedAttribute): def __init__(self, builder): self.momentum = builder.get_attribute("relative fall momentum") - self.water_mass = builder.get_attribute("water mass") + self.signed_water_mass = builder.get_attribute("signed water mass") super().__init__( builder, name="relative fall velocity", - dependencies=(self.momentum, self.water_mass), + dependencies=(self.momentum, self.signed_water_mass), ) def recalculate(self): - self.data.ratio(self.momentum.get(), self.water_mass.get()) + self.data.ratio(self.momentum.get(), self.signed_water_mass.get()) diff --git a/PySDM/dynamics/relaxed_velocity.py b/PySDM/dynamics/relaxed_velocity.py index 0bbf5e699..cd6d32a06 100644 --- a/PySDM/dynamics/relaxed_velocity.py +++ b/PySDM/dynamics/relaxed_velocity.py @@ -63,7 +63,7 @@ def register(self, builder): "relative fall momentum" ) self.terminal_vel_attr: Attribute = builder.get_attribute("terminal velocity") - self.water_mass_attr: Attribute = builder.get_attribute("water mass") + self.water_mass_attr: Attribute = builder.get_attribute("signed water mass") self.sqrt_radius_attr: Attribute = builder.get_attribute( "square root of radius" ) diff --git a/tests/unit_tests/attributes/test_fall_velocity.py b/tests/unit_tests/attributes/test_fall_velocity.py index 5a4489c34..8719a094f 100644 --- a/tests/unit_tests/attributes/test_fall_velocity.py +++ b/tests/unit_tests/attributes/test_fall_velocity.py @@ -77,7 +77,7 @@ def test_fall_velocity_calculation(default_attributes, backend_instance): assert np.allclose( particulator.attributes["relative fall velocity"].to_ndarray(), particulator.attributes["relative fall momentum"].to_ndarray() - / (particulator.attributes["water mass"].to_ndarray()), + / (particulator.attributes["signed water mass"].to_ndarray()), ) @staticmethod @@ -151,7 +151,7 @@ def test_attribute_selection(backend_instance): _ = builder.build( attributes={ "multiplicity": np.ones(1), - "water mass": np.zeros(1), + "signed water mass": np.zeros(1), "relative fall momentum": np.zeros(1), }, products=(), @@ -170,7 +170,7 @@ def test_attribute_selection(backend_instance): _ = builder.build( attributes={ "multiplicity": np.ones(1), - "water mass": np.zeros(1), + "signed water mass": np.zeros(1), }, products=(), ) From 99c41006f1acb58d31b486580151f8c7de3f9689 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Fri, 28 Mar 2025 17:58:39 +0100 Subject: [PATCH 06/12] separated ice and liquid terminal velocity classes --- PySDM/attributes/physics/terminal_velocity.py | 10 +++++++--- PySDM/formulae.py | 4 ++++ PySDM/physics/__init__.py | 1 + PySDM/physics/terminal_velocity/__init__.py | 3 +-- PySDM/physics/terminal_velocity_ice/__init__.py | 3 +++ .../columnar_ice_crystal.py | 0 tests/unit_tests/dynamics/test_terminal_velocity.py | 2 +- 7 files changed, 17 insertions(+), 6 deletions(-) create mode 100644 PySDM/physics/terminal_velocity_ice/__init__.py rename PySDM/physics/{terminal_velocity => terminal_velocity_ice}/columnar_ice_crystal.py (100%) diff --git a/PySDM/attributes/physics/terminal_velocity.py b/PySDM/attributes/physics/terminal_velocity.py index 43f1aaa52..4cac173a0 100644 --- a/PySDM/attributes/physics/terminal_velocity.py +++ b/PySDM/attributes/physics/terminal_velocity.py @@ -13,12 +13,16 @@ class TerminalVelocity(DerivedAttribute): def __init__(self, builder): self.radius = builder.get_attribute("radius") - dependencies = [self.radius] + self.signed_water_mass = builder.get_attribute("signed water mass") + dependencies = [self.radius,self.signed_water_mass] super().__init__(builder, name="terminal velocity", dependencies=dependencies) - self.approximation = builder.formulae.terminal_velocity_class( + self.approximation_liquid = builder.formulae.terminal_velocity_class( + builder.particulator + ) + self.approximation_ice = builder.formulae.terminal_velocity_ice_class( builder.particulator ) def recalculate(self): - self.approximation(self.data, self.radius.get()) + self.approximation_liquid(self.data, self.radius.get()) diff --git a/PySDM/formulae.py b/PySDM/formulae.py index 7329ca86f..581a06eef 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -90,6 +90,7 @@ def __init__( # pylint: disable=too-many-locals self.particle_shape_and_density = particle_shape_and_density self.air_dynamic_viscosity = air_dynamic_viscosity self.terminal_velocity = terminal_velocity + self.terminal_velocity_ice = terminal_velocity_ice self.bulk_phase_partitioning = bulk_phase_partitioning self._components = tuple( @@ -152,6 +153,9 @@ def __init__( # pylint: disable=too-many-locals "TpDependent": TpDependent, "PowerSeries": PowerSeries, }[terminal_velocity] + self.terminal_velocity_ice_class = { + "ColumnarIceCrystal": ColumnarIceCrystal, + }[terminal_velocity_ice] def __str__(self): description = [] diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py index 65ad336cb..6efeaab33 100644 --- a/PySDM/physics/__init__.py +++ b/PySDM/physics/__init__.py @@ -46,6 +46,7 @@ ventilation, air_dynamic_viscosity, terminal_velocity, + terminal_velocity_ice, bulk_phase_partitioning, ) from .constants import convert_to, in_unit, si diff --git a/PySDM/physics/terminal_velocity/__init__.py b/PySDM/physics/terminal_velocity/__init__.py index 5dad6a8f7..74f9ec8c7 100644 --- a/PySDM/physics/terminal_velocity/__init__.py +++ b/PySDM/physics/terminal_velocity/__init__.py @@ -1,6 +1,5 @@ -"""terminal velocity formulae""" +"""terminal velocity formulae for liquid droplets""" from .rogers_yau import RogersYau from .gunn_kinzer_1949 import GunnKinzer1949 from .power_series import PowerSeries -from .columnar_ice_crystal import ColumnarIceCrystal \ No newline at end of file diff --git a/PySDM/physics/terminal_velocity_ice/__init__.py b/PySDM/physics/terminal_velocity_ice/__init__.py new file mode 100644 index 000000000..1d9032019 --- /dev/null +++ b/PySDM/physics/terminal_velocity_ice/__init__.py @@ -0,0 +1,3 @@ +"""terminal velocity formulae for ice particles""" + +from .columnar_ice_crystal import ColumnarIceCrystal \ No newline at end of file diff --git a/PySDM/physics/terminal_velocity/columnar_ice_crystal.py b/PySDM/physics/terminal_velocity_ice/columnar_ice_crystal.py similarity index 100% rename from PySDM/physics/terminal_velocity/columnar_ice_crystal.py rename to PySDM/physics/terminal_velocity_ice/columnar_ice_crystal.py diff --git a/tests/unit_tests/dynamics/test_terminal_velocity.py b/tests/unit_tests/dynamics/test_terminal_velocity.py index 0b148bc76..7be2d1eae 100644 --- a/tests/unit_tests/dynamics/test_terminal_velocity.py +++ b/tests/unit_tests/dynamics/test_terminal_velocity.py @@ -73,7 +73,7 @@ def test_terminal_velocity_boundary_values( builder.request_attribute("terminal velocity") particulator = builder.build( attributes={ - "water mass": np.asarray([water_mass]), + "signed water mass": np.asarray([water_mass]), "multiplicity": np.asarray([-1]), } ) From d111cce38ec1940dfac6b71535fdc5c3b9913784 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Sat, 29 Mar 2025 16:08:43 +0100 Subject: [PATCH 07/12] added check for negative radius (indicating ice mass) to RogersYau and GunnKinzer terminal velocity --- PySDM/attributes/physics/terminal_velocity.py | 1 + .../methods/terminal_velocity_methods.py | 15 ++++++++++----- .../unit_tests/dynamics/test_terminal_velocity.py | 3 +++ 3 files changed, 14 insertions(+), 5 deletions(-) diff --git a/PySDM/attributes/physics/terminal_velocity.py b/PySDM/attributes/physics/terminal_velocity.py index 4cac173a0..5c0db0028 100644 --- a/PySDM/attributes/physics/terminal_velocity.py +++ b/PySDM/attributes/physics/terminal_velocity.py @@ -25,4 +25,5 @@ def __init__(self, builder): ) def recalculate(self): + print( "recalculate terminal_velocity" ) self.approximation_liquid(self.data, self.radius.get()) diff --git a/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py b/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py index ee0e14e6b..4d9ec275a 100644 --- a/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py +++ b/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py @@ -10,25 +10,28 @@ class TerminalVelocityMethods(BackendMethods): + # TODO: give terminal velocity functions name of the parametrisation @cached_property def _interpolation_body(self): @numba.njit(**self.default_jit_flags) def body(output, radius, factor, b, c): for i in numba.prange(len(radius)): # pylint: disable=not-an-iterable - if radius[i] < 0: - output[i] = 0 - else: + if radius[i] > 0: r_id = int(factor * radius[i]) r_rest = ((factor * radius[i]) % 1) / factor output[i] = b[r_id] + r_rest * c[r_id] + # TODO: check if output 0 for radius 0 is necessary + elif radius == 0: + output[i] = 0 return body + # TODO: give terminal velocity functions name of the parametrisation def interpolation(self, *, output, radius, factor, b, c): return self._interpolation_body( output.data, radius.data, factor, b.data, c.data ) - + # TODO: give terminal velocity functions name of the parametrisation @cached_property def _terminal_velocity_body(self): v_term = self.formulae.terminal_velocity.v_term @@ -36,10 +39,12 @@ def _terminal_velocity_body(self): @numba.njit(**self.default_jit_flags) def body(*, values, radius): for i in numba.prange(len(values)): # pylint: disable=not-an-iterable - values[i] = v_term(radius[i]) + if radius[i] >= 0.: + values[i] = v_term(radius[i]) return body + # TODO: give terminal velocity functions name of the parametrisation def terminal_velocity(self, *, values, radius): self._terminal_velocity_body(values=values, radius=radius) diff --git a/tests/unit_tests/dynamics/test_terminal_velocity.py b/tests/unit_tests/dynamics/test_terminal_velocity.py index 7be2d1eae..a2b6a4139 100644 --- a/tests/unit_tests/dynamics/test_terminal_velocity.py +++ b/tests/unit_tests/dynamics/test_terminal_velocity.py @@ -51,7 +51,9 @@ def test_approximation(backend_class, plot=False): ( ("GunnKinzer1949", 0 * si.g, None, 0), ("GunnKinzer1949", 1e10 * si.kg, pytest.raises(ValueError, match="Radii"), -1), + ("GunnKinzer1949", -1 * si.g, None, np.nan), ("RogersYau", 0 * si.g, None, 0), + ("RogersYau", -1 * si.g, None, np.nan), ("TpDependent", 0 * si.g, None, 0), ), ) @@ -83,6 +85,7 @@ def test_terminal_velocity_boundary_values( (v_term,) = particulator.attributes["terminal velocity"].to_ndarray() # assert + print( v_term, expected_v_term) np.testing.assert_approx_equal(v_term, expected_v_term) From 13fbcd3885dc111db839002101ff85656b55c1d8 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Sat, 29 Mar 2025 20:04:00 +0100 Subject: [PATCH 08/12] further work on columnar ice crystal terminal velocity and added a unit test --- PySDM/attributes/physics/terminal_velocity.py | 2 +- .../methods/terminal_velocity_methods.py | 18 ++++++++ .../terminal_velocity/columnar_ice_crystal.py | 6 +-- .../columnar_ice_crystal.py | 2 +- .../dynamics/test_terminal_velocity.py | 44 ++++++++++++++++++- 5 files changed, 66 insertions(+), 6 deletions(-) diff --git a/PySDM/attributes/physics/terminal_velocity.py b/PySDM/attributes/physics/terminal_velocity.py index 5c0db0028..b27885e6f 100644 --- a/PySDM/attributes/physics/terminal_velocity.py +++ b/PySDM/attributes/physics/terminal_velocity.py @@ -25,5 +25,5 @@ def __init__(self, builder): ) def recalculate(self): - print( "recalculate terminal_velocity" ) self.approximation_liquid(self.data, self.radius.get()) + self.approximation_ice(self.data, self.signed_water_mass.get()) diff --git a/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py b/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py index 4d9ec275a..ec6406073 100644 --- a/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py +++ b/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py @@ -67,3 +67,21 @@ def power_series(self, *, values, radius, num_terms, prefactors, powers): prefactors=prefactors, powers=powers, ) + + + def terminal_velocity_columnar_ice_crystals(self, *, values, signed_water_mass): + self._terminal_velocity_columnar_ice_crystals_body( values=values, + signed_water_mass=signed_water_mass, + ) + + @cached_property + def _terminal_velocity_columnar_ice_crystals_body(self): + v_base_term = self.formulae.terminal_velocity_ice.v_base_term + + # @numba.njit(**self.default_jit_flags) + def body(*, values, signed_water_mass): + for i in numba.prange(len(values)): # pylint: disable=not-an-iterable + if signed_water_mass[i] < 0: + values[i] = v_base_term(-signed_water_mass[i]) + + return body \ No newline at end of file diff --git a/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py b/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py index 92320c94f..968970697 100644 --- a/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py +++ b/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py @@ -9,8 +9,8 @@ class ColumnarIceCrystal: # pylint: disable=too-few-public-methods def __init__(self, particulator): self.particulator = particulator - def __call__(self, output, mass): - self.particulator.backend.terminal_velocity( + def __call__(self, output, signed_water_mass): + self.particulator.backend.terminal_velocity_columnar_ice_crystals( values=output.data, - mass=mass.data, + signed_water_mass=signed_water_mass.data, ) \ No newline at end of file diff --git a/PySDM/physics/terminal_velocity_ice/columnar_ice_crystal.py b/PySDM/physics/terminal_velocity_ice/columnar_ice_crystal.py index d8487cb42..89638000f 100644 --- a/PySDM/physics/terminal_velocity_ice/columnar_ice_crystal.py +++ b/PySDM/physics/terminal_velocity_ice/columnar_ice_crystal.py @@ -13,7 +13,7 @@ def __init__(self, _): pass @staticmethod - def v_term(const, mass): + def v_base_term(const, mass): return np.where( mass < const.SPICHTINGER_GIERENS_TERM_VEL_LIMIT_0, const.SPICHTINGER_TERM_GAMMA_COEFF0 * mass**const.SPICHTINGER_TERM_DELTA_COEFF0, diff --git a/tests/unit_tests/dynamics/test_terminal_velocity.py b/tests/unit_tests/dynamics/test_terminal_velocity.py index a2b6a4139..538a8b4b0 100644 --- a/tests/unit_tests/dynamics/test_terminal_velocity.py +++ b/tests/unit_tests/dynamics/test_terminal_velocity.py @@ -6,7 +6,7 @@ import pytest from PySDM import Builder, Formulae -from PySDM.dynamics.terminal_velocity import GunnKinzer1949, PowerSeries, RogersYau +from PySDM.dynamics.terminal_velocity import GunnKinzer1949, PowerSeries, RogersYau, ColumnarIceCrystal from PySDM.environments import Box from PySDM.physics import constants as const from PySDM.physics import si @@ -113,3 +113,45 @@ def test_power_series(backend_class, prefactors, powers): u_term_true = particulator.backend.Storage.from_ndarray(u) np.testing.assert_array_almost_equal(u, u_term_true) + + + +def test_columnar_ice_crystal_terminal_velocity(backend_class, plot=False): + """Fig. 3 in [Spichtinger & Gierens 2009](https://doi.org/10.5194/acp-9-685-2009)""" + if backend_class.__name__ == "ThrustRTC": + pytest.skip() + # arrange + water_mass = ( + np.logspace(base=10, start=-16, stop=-7, num=10) * si.kg + ) + terminal_velocity_reference = [1.4e-04, 3.7e-04, 9.7e-04, 2.5e-03, + 9.1e-03, 3.4e-02, 1.3e-01, 4.7e-01, + 1.1e+00, 1.9e+00] + + + particulator = DummyParticulator( + backend_class, n_sd=len(water_mass), formulae=Formulae(terminal_velocity_ice="ColumnarIceCrystal") + ) + + terminal_velocity = particulator.backend.Storage.empty((len(water_mass),), float) + ColumnarIceCrystal(particulator=particulator)(terminal_velocity, -water_mass) + + plt.xlabel("mass (kg)") + plt.ylabel("terminal velocity (m/s)") + plt.xlim(water_mass[0], water_mass[-1]) + plt.xscale("log") + plt.ylim(1e-4, 1e1) + plt.yscale("log") + plt.grid() + + plt.plot(water_mass, terminal_velocity_reference, color="black") + plt.plot(water_mass, terminal_velocity, color="red") + + # plot + if plot: + plt.show() + else: + plt.clf() + + # assert + np.testing.assert_almost_equal(terminal_velocity, terminal_velocity_reference, decimal=1) \ No newline at end of file From 43b35db46899d16e478e138e8aac31d6f36efc36 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Mon, 31 Mar 2025 16:46:26 +0200 Subject: [PATCH 09/12] added temperature and pressure dependence for columnar ice crystal terminal velocity --- PySDM/attributes/physics/terminal_velocity.py | 13 +++++-- .../methods/terminal_velocity_methods.py | 15 +++++--- .../terminal_velocity/columnar_ice_crystal.py | 5 ++- .../dynamics/test_terminal_velocity.py | 35 ++++++++++++++----- 4 files changed, 52 insertions(+), 16 deletions(-) diff --git a/PySDM/attributes/physics/terminal_velocity.py b/PySDM/attributes/physics/terminal_velocity.py index b27885e6f..106bc29f7 100644 --- a/PySDM/attributes/physics/terminal_velocity.py +++ b/PySDM/attributes/physics/terminal_velocity.py @@ -14,7 +14,8 @@ class TerminalVelocity(DerivedAttribute): def __init__(self, builder): self.radius = builder.get_attribute("radius") self.signed_water_mass = builder.get_attribute("signed water mass") - dependencies = [self.radius,self.signed_water_mass] + self.cell_id = builder.get_attribute("cell id") + dependencies = [self.radius,self.signed_water_mass,self.cell_id] super().__init__(builder, name="terminal velocity", dependencies=dependencies) self.approximation_liquid = builder.formulae.terminal_velocity_class( @@ -26,4 +27,12 @@ def __init__(self, builder): def recalculate(self): self.approximation_liquid(self.data, self.radius.get()) - self.approximation_ice(self.data, self.signed_water_mass.get()) + if self.formulae.particle_shape_and_density.supports_mixed_phase(): + self.approximation_ice(self.data, + self.signed_water_mass.get(), + self.cell_id.get(), + self.particulator.environment["T"], + self.particulator.environment["p"], + ) + + diff --git a/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py b/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py index ec6406073..9be8126ef 100644 --- a/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py +++ b/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py @@ -69,19 +69,26 @@ def power_series(self, *, values, radius, num_terms, prefactors, powers): ) - def terminal_velocity_columnar_ice_crystals(self, *, values, signed_water_mass): + def terminal_velocity_columnar_ice_crystals(self, *, values, signed_water_mass, cell_id, temperature, pressure): self._terminal_velocity_columnar_ice_crystals_body( values=values, signed_water_mass=signed_water_mass, + cell_id=cell_id, + temperature=temperature, + pressure=pressure, ) @cached_property def _terminal_velocity_columnar_ice_crystals_body(self): v_base_term = self.formulae.terminal_velocity_ice.v_base_term + atmospheric_correction_factor = self.formulae.terminal_velocity_ice.atmospheric_correction_factor - # @numba.njit(**self.default_jit_flags) - def body(*, values, signed_water_mass): + @numba.njit(**self.default_jit_flags) + def body(*, values, signed_water_mass,cell_id,temperature,pressure): for i in numba.prange(len(values)): # pylint: disable=not-an-iterable if signed_water_mass[i] < 0: - values[i] = v_base_term(-signed_water_mass[i]) + cid = cell_id[i] + correction = atmospheric_correction_factor(temperature[cid], pressure[cid]) + values[i] = v_base_term(-signed_water_mass[i]) * correction + return body \ No newline at end of file diff --git a/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py b/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py index 968970697..34d8ad30a 100644 --- a/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py +++ b/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py @@ -9,8 +9,11 @@ class ColumnarIceCrystal: # pylint: disable=too-few-public-methods def __init__(self, particulator): self.particulator = particulator - def __call__(self, output, signed_water_mass): + def __call__(self, output, signed_water_mass,cell_id,temperature,pressure): self.particulator.backend.terminal_velocity_columnar_ice_crystals( values=output.data, signed_water_mass=signed_water_mass.data, + cell_id=cell_id.data, + temperature=temperature.data, + pressure=pressure.data, ) \ No newline at end of file diff --git a/tests/unit_tests/dynamics/test_terminal_velocity.py b/tests/unit_tests/dynamics/test_terminal_velocity.py index 538a8b4b0..5be165340 100644 --- a/tests/unit_tests/dynamics/test_terminal_velocity.py +++ b/tests/unit_tests/dynamics/test_terminal_velocity.py @@ -85,7 +85,6 @@ def test_terminal_velocity_boundary_values( (v_term,) = particulator.attributes["terminal velocity"].to_ndarray() # assert - print( v_term, expected_v_term) np.testing.assert_approx_equal(v_term, expected_v_term) @@ -124,18 +123,37 @@ def test_columnar_ice_crystal_terminal_velocity(backend_class, plot=False): water_mass = ( np.logspace(base=10, start=-16, stop=-7, num=10) * si.kg ) - terminal_velocity_reference = [1.4e-04, 3.7e-04, 9.7e-04, 2.5e-03, - 9.1e-03, 3.4e-02, 1.3e-01, 4.7e-01, - 1.1e+00, 1.9e+00] + terminal_velocity_reference = ( + np.array([1.4e-04, 3.7e-04, 9.7e-04, 2.5e-03, + 9.1e-03, 3.4e-02, 1.3e-01, 4.7e-01, + 1.1e+00, 1.9e+00]) * si.m / si.s + ) + ambient_temperature = 233 * si.K + ambient_pressure = 300 * si.hectopascal + env = Box(dt=None, dv=None) + formulae_enabling_terminal_velocity_ice_calculation = Formulae( + particle_shape_and_density="MixedPhaseSpheres", + terminal_velocity_ice="ColumnarIceCrystal", + ) + builder = Builder( + backend=backend_class(formulae_enabling_terminal_velocity_ice_calculation), + n_sd=len(water_mass), + environment=env, + ) - particulator = DummyParticulator( - backend_class, n_sd=len(water_mass), formulae=Formulae(terminal_velocity_ice="ColumnarIceCrystal") + builder.request_attribute("terminal velocity") + particulator = builder.build( + attributes={"signed water mass": -water_mass, "multiplicity": np.ones_like(water_mass)} ) - terminal_velocity = particulator.backend.Storage.empty((len(water_mass),), float) - ColumnarIceCrystal(particulator=particulator)(terminal_velocity, -water_mass) + particulator.environment["T"] = ambient_temperature + particulator.environment["p"] = ambient_pressure + + # act + terminal_velocity = particulator.attributes["terminal velocity"].to_ndarray() + # plot plt.xlabel("mass (kg)") plt.ylabel("terminal velocity (m/s)") plt.xlim(water_mass[0], water_mass[-1]) @@ -147,7 +165,6 @@ def test_columnar_ice_crystal_terminal_velocity(backend_class, plot=False): plt.plot(water_mass, terminal_velocity_reference, color="black") plt.plot(water_mass, terminal_velocity, color="red") - # plot if plot: plt.show() else: From c185563ee5af4549c82d454dfb5fc807a552508a Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Mon, 31 Mar 2025 19:12:07 +0200 Subject: [PATCH 10/12] added todo --- PySDM/attributes/physics/terminal_velocity.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/PySDM/attributes/physics/terminal_velocity.py b/PySDM/attributes/physics/terminal_velocity.py index 106bc29f7..0e84c0db4 100644 --- a/PySDM/attributes/physics/terminal_velocity.py +++ b/PySDM/attributes/physics/terminal_velocity.py @@ -27,6 +27,8 @@ def __init__(self, builder): def recalculate(self): self.approximation_liquid(self.data, self.radius.get()) + # TODO: fix issue that order of functions calls changes result. approximation_liquid will override + # approximation_ice since r < 0 is not a suitable test for ice particles with mixed-phase spheres shape active if self.formulae.particle_shape_and_density.supports_mixed_phase(): self.approximation_ice(self.data, self.signed_water_mass.get(), From 62e119201e3013ae4a6373992726416c4cf4139e Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Tue, 1 Apr 2025 13:34:11 +0200 Subject: [PATCH 11/12] added terminal velocity for smooth ice spheres --- .../methods/terminal_velocity_methods.py | 28 ++++++++++++++++++- PySDM/dynamics/terminal_velocity/__init__.py | 3 +- .../terminal_velocity/columnar_ice_crystal.py | 2 +- .../dynamics/terminal_velocity/spheres_ice.py | 15 ++++++++++ PySDM/formulae.py | 3 +- .../physics/terminal_velocity_ice/__init__.py | 3 +- .../terminal_velocity_ice/spheres_ice.py | 20 +++++++++++++ 7 files changed, 69 insertions(+), 5 deletions(-) create mode 100644 PySDM/dynamics/terminal_velocity/spheres_ice.py create mode 100644 PySDM/physics/terminal_velocity_ice/spheres_ice.py diff --git a/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py b/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py index 9be8126ef..76c69f6a1 100644 --- a/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py +++ b/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py @@ -83,12 +83,38 @@ def _terminal_velocity_columnar_ice_crystals_body(self): atmospheric_correction_factor = self.formulae.terminal_velocity_ice.atmospheric_correction_factor @numba.njit(**self.default_jit_flags) - def body(*, values, signed_water_mass,cell_id,temperature,pressure): + def body(*, values, signed_water_mass, cell_id, temperature, pressure): for i in numba.prange(len(values)): # pylint: disable=not-an-iterable if signed_water_mass[i] < 0: cid = cell_id[i] correction = atmospheric_correction_factor(temperature[cid], pressure[cid]) values[i] = v_base_term(-signed_water_mass[i]) * correction + return body + + + + def terminal_velocity_ice_spheres(self, *, values, signed_water_mass, cell_id, temperature, pressure): + self._terminal_velocity_ice_spheres_body(values=values, + signed_water_mass=signed_water_mass, + cell_id=cell_id, + temperature=temperature, + ) + + @cached_property + def _terminal_velocity_ice_spheres_body(self): + v_base_term = self.formulae.terminal_velocity_ice.v_base_term + stokes_prefactor = self.formulae.terminal_velocity_ice.stokes_regime + formulae = self.formulae_flattened + + def body(*, values, signed_water_mass, cell_id, temperature): + for i in numba.prange(len(values)): # pylint: disable=not-an-iterable + if signed_water_mass[i] < 0: + cid = cell_id[i] + radius = formulae.particle_shape_and_density__mass_to_radius(signed_water_mass[i]) + dynamic_viscosity = formulae.air_dynamic_viscosity__eta_air(temperature[cid]) + prefactor = stokes_prefactor( radius, dynamic_viscosity) + values[i] = v_base_term(prefactor, radius) + print( dynamic_viscosity, signed_water_mass[i], radius, prefactor, values[i] ) return body \ No newline at end of file diff --git a/PySDM/dynamics/terminal_velocity/__init__.py b/PySDM/dynamics/terminal_velocity/__init__.py index 17d0bcf65..981c822bb 100644 --- a/PySDM/dynamics/terminal_velocity/__init__.py +++ b/PySDM/dynamics/terminal_velocity/__init__.py @@ -5,4 +5,5 @@ from PySDM.dynamics.terminal_velocity.gunn_and_kinzer import GunnKinzer1949 from PySDM.dynamics.terminal_velocity.power_series import PowerSeries from PySDM.dynamics.terminal_velocity.rogers_and_yau import RogersYau -from PySDM.dynamics.terminal_velocity.columnar_ice_crystal import ColumnarIceCrystal \ No newline at end of file +from PySDM.dynamics.terminal_velocity.columnar_ice_crystal import ColumnarIceCrystal +from PySDM.dynamics.terminal_velocity.spheres_ice import IceSphere \ No newline at end of file diff --git a/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py b/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py index 34d8ad30a..ab3b75f8b 100644 --- a/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py +++ b/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py @@ -9,7 +9,7 @@ class ColumnarIceCrystal: # pylint: disable=too-few-public-methods def __init__(self, particulator): self.particulator = particulator - def __call__(self, output, signed_water_mass,cell_id,temperature,pressure): + def __call__(self, output, signed_water_mass, cell_id, temperature, pressure): self.particulator.backend.terminal_velocity_columnar_ice_crystals( values=output.data, signed_water_mass=signed_water_mass.data, diff --git a/PySDM/dynamics/terminal_velocity/spheres_ice.py b/PySDM/dynamics/terminal_velocity/spheres_ice.py new file mode 100644 index 000000000..c7f3a82eb --- /dev/null +++ b/PySDM/dynamics/terminal_velocity/spheres_ice.py @@ -0,0 +1,15 @@ + + + +class IceSphere: # pylint: disable=too-few-public-methods + def __init__(self, particulator): + self.particulator = particulator + + def __call__(self, output, signed_water_mass, cell_id, temperature, pressure): + self.particulator.backend.terminal_velocity_ice_spheres( + values=output.data, + signed_water_mass=signed_water_mass.data, + cell_id=cell_id.data, + temperature=temperature.data, + pressure=pressure.data, + ) \ No newline at end of file diff --git a/PySDM/formulae.py b/PySDM/formulae.py index 801a6848f..6a07dc542 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -19,7 +19,7 @@ from PySDM import physics from PySDM.backends.impl_numba import conf -from PySDM.dynamics.terminal_velocity import GunnKinzer1949, PowerSeries, RogersYau, ColumnarIceCrystal +from PySDM.dynamics.terminal_velocity import GunnKinzer1949, PowerSeries, RogersYau, ColumnarIceCrystal, IceSphere from PySDM.dynamics.terminal_velocity.gunn_and_kinzer import TpDependent @@ -162,6 +162,7 @@ def __init__( # pylint: disable=too-many-locals }[terminal_velocity] self.terminal_velocity_ice_class = { "ColumnarIceCrystal": ColumnarIceCrystal, + "IceSphere": IceSphere, }[terminal_velocity_ice] def __str__(self): diff --git a/PySDM/physics/terminal_velocity_ice/__init__.py b/PySDM/physics/terminal_velocity_ice/__init__.py index 1d9032019..ca51d0e52 100644 --- a/PySDM/physics/terminal_velocity_ice/__init__.py +++ b/PySDM/physics/terminal_velocity_ice/__init__.py @@ -1,3 +1,4 @@ """terminal velocity formulae for ice particles""" -from .columnar_ice_crystal import ColumnarIceCrystal \ No newline at end of file +from .columnar_ice_crystal import ColumnarIceCrystal +from .spheres_ice import IceSphere \ No newline at end of file diff --git a/PySDM/physics/terminal_velocity_ice/spheres_ice.py b/PySDM/physics/terminal_velocity_ice/spheres_ice.py new file mode 100644 index 000000000..7a6bfe9be --- /dev/null +++ b/PySDM/physics/terminal_velocity_ice/spheres_ice.py @@ -0,0 +1,20 @@ +""" +terminal velocity of smooth ice spheres +""" + +class IceSphere: # pylint: disable=too-few-public-methods + def __init__(self, _): + pass + + @staticmethod + def v_base_term(radius, prefactor): + return( prefactor * 2 * radius) + + @staticmethod + def stokes_regime(const, radius, dynamic_viscosity): + return( const.g_std * const.rho_i * radius / 9 / dynamic_viscosity ) + + # TODO: find parametrisation for general functional relationship between reynolds number and drag coefficient + @staticmethod + def general_flow_regime(const): + pass \ No newline at end of file From 8b7557f24ca3360f666ef02fb28509fbde1338d4 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Tue, 1 Apr 2025 15:04:01 +0200 Subject: [PATCH 12/12] added new unit test for terminal velocity of ice particles --- .../dynamics/test_terminal_velocity.py | 49 ++++++++++++++++++- 1 file changed, 48 insertions(+), 1 deletion(-) diff --git a/tests/unit_tests/dynamics/test_terminal_velocity.py b/tests/unit_tests/dynamics/test_terminal_velocity.py index 5be165340..9efee0bac 100644 --- a/tests/unit_tests/dynamics/test_terminal_velocity.py +++ b/tests/unit_tests/dynamics/test_terminal_velocity.py @@ -115,7 +115,54 @@ def test_power_series(backend_class, prefactors, powers): -def test_columnar_ice_crystal_terminal_velocity(backend_class, plot=False): + +@pytest.mark.parametrize("ice_variant", ("ColumnarIceCrystal","IceSphere")) +def test_ice_particle_terminal_velocities_basics(backend_class, ice_variant): + if backend_class.__name__ == "ThrustRTC": + pytest.skip() + + # arrange + water_mass = ( + np.logspace(base=10, start=-16, stop=-7, num=10) * si.kg + ) + env = Box(dt=None, dv=None) + formulae_enabling_terminal_velocity_ice_calculation = Formulae( + particle_shape_and_density="MixedPhaseSpheres", + terminal_velocity_ice=ice_variant, + ) + builder = Builder( + backend=backend_class(formulae_enabling_terminal_velocity_ice_calculation), + n_sd=len(water_mass), + environment=env, + ) + builder.request_attribute("terminal velocity") + particulator = builder.build( + attributes={"signed water mass": -water_mass, "multiplicity": np.ones_like(water_mass)} + ) + atmospheric_settings = [ + {"temperature": 233 * si.kelvin, "pressure": 300 * si.hectopascal}, + {"temperature": 270 * si.kelvin, "pressure": 1000 * si.hectopascal}, + ] + for setting in atmospheric_settings: + + particulator.environment["T"] = setting["temperature"] + particulator.environment["p"] = setting["pressure"] + + # act + particulator.run(steps=1) + terminal_velocity = particulator.attributes["terminal velocity"].to_ndarray() + # TODO: find a way to update terminal velocity attribute when environment variables change + print(terminal_velocity) + + # assert + assert all(~np.isnan(terminal_velocity)) + assert all(terminal_velocity > 0.0) + assert all(np.diff(terminal_velocity) > 0.0) + + + + +def test_columnar_ice_crystal_terminal_velocity_against_spichtinger_and_gierens_2009_fig_3(backend_class, plot=False): """Fig. 3 in [Spichtinger & Gierens 2009](https://doi.org/10.5194/acp-9-685-2009)""" if backend_class.__name__ == "ThrustRTC": pytest.skip()