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/attributes/physics/terminal_velocity.py b/PySDM/attributes/physics/terminal_velocity.py index 43f1aaa52..0e84c0db4 100644 --- a/PySDM/attributes/physics/terminal_velocity.py +++ b/PySDM/attributes/physics/terminal_velocity.py @@ -13,12 +13,28 @@ 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") + 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 = 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()) + # 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(), + 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 ee0e14e6b..76c69f6a1 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) @@ -62,3 +67,54 @@ 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, 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, 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/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/PySDM/dynamics/terminal_velocity/__init__.py b/PySDM/dynamics/terminal_velocity/__init__.py index c94a06bca..981c822bb 100644 --- a/PySDM/dynamics/terminal_velocity/__init__.py +++ b/PySDM/dynamics/terminal_velocity/__init__.py @@ -5,3 +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 +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 new file mode 100644 index 000000000..ab3b75f8b --- /dev/null +++ b/PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py @@ -0,0 +1,19 @@ +""" +[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 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): + 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/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 252fedd04..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 +from PySDM.dynamics.terminal_velocity import GunnKinzer1949, PowerSeries, RogersYau, ColumnarIceCrystal, IceSphere from PySDM.dynamics.terminal_velocity.gunn_and_kinzer import TpDependent @@ -58,6 +58,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, @@ -96,6 +97,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( @@ -158,6 +160,10 @@ def __init__( # pylint: disable=too-many-locals "TpDependent": TpDependent, "PowerSeries": PowerSeries, }[terminal_velocity] + self.terminal_velocity_ice_class = { + "ColumnarIceCrystal": ColumnarIceCrystal, + "IceSphere": IceSphere, + }[terminal_velocity_ice] def __str__(self): description = [] diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py index a39b3c8b1..4dbef2c76 100644 --- a/PySDM/physics/__init__.py +++ b/PySDM/physics/__init__.py @@ -49,6 +49,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/constants_defaults.py b/PySDM/physics/constants_defaults.py index 241d9f693..bfd2c30e4 100644 --- a/PySDM/physics/constants_defaults.py +++ b/PySDM/physics/constants_defaults.py @@ -580,6 +580,39 @@ 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 +""" 〃 """ +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 """ W76W_G1 = -6.0170128e3 * si.K diff --git a/PySDM/physics/terminal_velocity/__init__.py b/PySDM/physics/terminal_velocity/__init__.py index 2e012dc47..74f9ec8c7 100644 --- a/PySDM/physics/terminal_velocity/__init__.py +++ b/PySDM/physics/terminal_velocity/__init__.py @@ -1,4 +1,4 @@ -"""terminal velocity formulae""" +"""terminal velocity formulae for liquid droplets""" from .rogers_yau import RogersYau from .gunn_kinzer_1949 import GunnKinzer1949 diff --git a/PySDM/physics/terminal_velocity_ice/__init__.py b/PySDM/physics/terminal_velocity_ice/__init__.py new file mode 100644 index 000000000..ca51d0e52 --- /dev/null +++ b/PySDM/physics/terminal_velocity_ice/__init__.py @@ -0,0 +1,4 @@ +"""terminal velocity formulae for ice particles""" + +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/columnar_ice_crystal.py b/PySDM/physics/terminal_velocity_ice/columnar_ice_crystal.py new file mode 100644 index 000000000..89638000f --- /dev/null +++ b/PySDM/physics/terminal_velocity_ice/columnar_ice_crystal.py @@ -0,0 +1,34 @@ +""" +[Spichtinger & Gierens 2009](https://doi.org/10.5194/acp-9-685-2009) +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) +""" + +import numpy as np + + +class ColumnarIceCrystal: # pylint: disable=too-few-public-methods + def __init__(self, _): + pass + + @staticmethod + 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, + 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 + ) + ), + ) + + @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 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 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=(), ) 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( diff --git a/tests/unit_tests/dynamics/test_terminal_velocity.py b/tests/unit_tests/dynamics/test_terminal_velocity.py index 0b148bc76..9efee0bac 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 @@ -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), ), ) @@ -73,7 +75,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]), } ) @@ -110,3 +112,110 @@ 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) + + + + +@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() + # arrange + water_mass = ( + np.logspace(base=10, start=-16, stop=-7, num=10) * si.kg + ) + 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, + ) + + builder.request_attribute("terminal velocity") + particulator = builder.build( + attributes={"signed water mass": -water_mass, "multiplicity": np.ones_like(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]) + 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") + + 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