From 19921943a4356a4b75f4089bd3e7d5ad92c30d32 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Tue, 17 Sep 2024 12:33:29 +0200 Subject: [PATCH 01/27] boilerplate code for starting with vapour deposition on ice --- PySDM/backends/impl_numba/methods/__init__.py | 1 + .../impl_numba/methods/deposition_methods.py | 26 +++++++++++ PySDM/backends/numba.py | 2 + PySDM/dynamics/__init__.py | 1 + PySDM/dynamics/vapour_deposition_on_ice.py | 16 +++++++ PySDM/particulator.py | 6 +++ .../dynamics/test_vapour_deposition_on_ice.py | 44 +++++++++++++++++++ 7 files changed, 96 insertions(+) create mode 100644 PySDM/backends/impl_numba/methods/deposition_methods.py create mode 100644 PySDM/dynamics/vapour_deposition_on_ice.py create mode 100644 tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py diff --git a/PySDM/backends/impl_numba/methods/__init__.py b/PySDM/backends/impl_numba/methods/__init__.py index 73d7ceed3..8ca6c05e9 100644 --- a/PySDM/backends/impl_numba/methods/__init__.py +++ b/PySDM/backends/impl_numba/methods/__init__.py @@ -13,3 +13,4 @@ from .physics_methods import PhysicsMethods from .terminal_velocity_methods import TerminalVelocityMethods from .seeding_methods import SeedingMethods +from .deposition_methods import DepositionMethods diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py new file mode 100644 index 000000000..737503411 --- /dev/null +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -0,0 +1,26 @@ +from functools import cached_property +from PySDM.backends.impl_common.backend_methods import BackendMethods +import numba + + +class DepositionMethods(BackendMethods): + @cached_property + def _deposition(self): + assert self.formulae.particle_shape_and_density.supports_mixed_phase() + + mass_to_volume = self.formulae.particle_shape_and_density.mass_to_volume + + @numba.jit(**self.default_jit_flags) + def body(water_mass): + volume = mass_to_volume(water_mass[0]) + print(volume) + + if water_mass[0] > 0: + water_mass[0] *= 1.1 + + return body + + def deposition(self, water_mass): + self._deposition( + water_mass=water_mass.data + ) diff --git a/PySDM/backends/numba.py b/PySDM/backends/numba.py index a3b558d11..3bf2197b8 100644 --- a/PySDM/backends/numba.py +++ b/PySDM/backends/numba.py @@ -29,6 +29,7 @@ class Numba( # pylint: disable=too-many-ancestors,duplicate-code methods.TerminalVelocityMethods, methods.IsotopeMethods, methods.SeedingMethods, + methods.DepositionMethods ): Storage = ImportedStorage Random = ImportedRandom @@ -77,3 +78,4 @@ def __init__(self, formulae=None, double_precision=True, override_jit_flags=None methods.TerminalVelocityMethods.__init__(self) methods.IsotopeMethods.__init__(self) methods.SeedingMethods.__init__(self) + methods.DepositionMethods.__init__(self) \ No newline at end of file diff --git a/PySDM/dynamics/__init__.py b/PySDM/dynamics/__init__.py index 66056cf0d..d296bf472 100644 --- a/PySDM/dynamics/__init__.py +++ b/PySDM/dynamics/__init__.py @@ -16,3 +16,4 @@ from PySDM.dynamics.freezing import Freezing from PySDM.dynamics.relaxed_velocity import RelaxedVelocity from PySDM.dynamics.seeding import Seeding +from PySDM.dynamics.vapour_deposition_on_ice import VapourDepositionOnIce diff --git a/PySDM/dynamics/vapour_deposition_on_ice.py b/PySDM/dynamics/vapour_deposition_on_ice.py new file mode 100644 index 000000000..7426309e2 --- /dev/null +++ b/PySDM/dynamics/vapour_deposition_on_ice.py @@ -0,0 +1,16 @@ +from PySDM.dynamics.impl import register_dynamic + + +@register_dynamic() +class VapourDepositionOnIce: + def __init__(self): + """ called by the user while building a particulator """ + self.particulator = None + + def register(self, *, builder): + """ called by the builder """ + self.particulator = builder.particulator + + def __call__(self): + """ called by the particulator during simulation""" + self.particulator.deposition() diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 7ab1098ce..5f4f89b5a 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -479,3 +479,9 @@ def seeding( self.attributes.mark_updated("multiplicity") for key in self.attributes.get_extensive_attribute_keys(): self.attributes.mark_updated(key) + + def deposition(self): + self.backend.deposition( + water_mass=self.attributes["water mass"] + ) + self.attributes.mark_updated("water mass") \ No newline at end of file diff --git a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py new file mode 100644 index 000000000..341272e6d --- /dev/null +++ b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py @@ -0,0 +1,44 @@ +import numpy as np + +import pytest + +from PySDM.physics import si +from PySDM.backends import CPU +from PySDM import Builder +from PySDM import Formulae +from PySDM.environments import Box +from PySDM.dynamics import VapourDepositionOnIce +from PySDM.products import IceWaterContent + + +@pytest.mark.parametrize('dt', (1 * si.s, .1 * si.s)) +@pytest.mark.parametrize('water_mass', (-si.ng, -si.ug, -si.mg, si.mg)) +@pytest.mark.parametrize('fastmath', (True, False)) +def test_TODO(dt, water_mass, fastmath, dv=1*si.m**3): + # arrange + n_sd = 1 + builder = Builder( + n_sd=n_sd, + environment=Box(dt=dt, dv=dv), + backend=CPU(formulae=Formulae( + fastmath=fastmath, + particle_shape_and_density="MixedPhaseSpheres" + )) + ) + deposition = VapourDepositionOnIce() + builder.add_dynamic(deposition) + particulator = builder.build( + attributes={ + 'multiplicity': np.full(shape=(n_sd,), fill_value=1), + 'water mass': np.full(shape=(n_sd,), fill_value=water_mass), + }, + products=(IceWaterContent(),) + ) + + # act + iwc_old = particulator.products['ice water content'].get().copy() + particulator.run(steps=1) + iwc_new = particulator.products['ice water content'].get().copy() + + # assert + assert (iwc_new > iwc_old).all() \ No newline at end of file From b800e6493d0e3d7d0a1db90808d93818c3183870 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Tue, 17 Sep 2024 17:24:49 +0200 Subject: [PATCH 02/27] Added looping over super particles and implemented ambient parameters --- .../impl_numba/methods/deposition_methods.py | 44 ++++++++++++++++--- PySDM/backends/numba.py | 4 +- PySDM/dynamics/vapour_deposition_on_ice.py | 6 +-- PySDM/particulator.py | 8 +++- .../dynamics/test_vapour_deposition_on_ice.py | 32 +++++++------- 5 files changed, 65 insertions(+), 29 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index 737503411..ba1390913 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -9,18 +9,48 @@ def _deposition(self): assert self.formulae.particle_shape_and_density.supports_mixed_phase() mass_to_volume = self.formulae.particle_shape_and_density.mass_to_volume + diffusion_coefficient_function = self.formulae.diffusion_thermics.D @numba.jit(**self.default_jit_flags) - def body(water_mass): - volume = mass_to_volume(water_mass[0]) - print(volume) + def body( + water_mass, ambient_temperature, ambient_total_pressure, time_step, cell_id + ): - if water_mass[0] > 0: - water_mass[0] *= 1.1 + n_sd = len(water_mass) + for i in numba.prange(n_sd): # pylint: disable=not-an-iterable + + if water_mass[i] < 0: + + cid = cell_id[i] + + volume = mass_to_volume(water_mass[i]) + + temperature = ambient_temperature[cid] + pressure = ambient_total_pressure[cid] + + diffusion_coefficient = diffusion_coefficient_function( + temperature, pressure + ) + + print(volume, temperature, pressure, diffusion_coefficient) + + water_mass[i] *= 1.1 return body - def deposition(self, water_mass): + def deposition( + self, + *, + water_mass, + ambient_temperature, + ambient_total_pressure, + time_step, + cell_id + ): self._deposition( - water_mass=water_mass.data + water_mass=water_mass.data, + ambient_temperature=ambient_temperature.data, + ambient_total_pressure=ambient_total_pressure.data, + time_step=time_step, + cell_id=cell_id.data, ) diff --git a/PySDM/backends/numba.py b/PySDM/backends/numba.py index 3bf2197b8..d138d905f 100644 --- a/PySDM/backends/numba.py +++ b/PySDM/backends/numba.py @@ -29,7 +29,7 @@ class Numba( # pylint: disable=too-many-ancestors,duplicate-code methods.TerminalVelocityMethods, methods.IsotopeMethods, methods.SeedingMethods, - methods.DepositionMethods + methods.DepositionMethods, ): Storage = ImportedStorage Random = ImportedRandom @@ -78,4 +78,4 @@ def __init__(self, formulae=None, double_precision=True, override_jit_flags=None methods.TerminalVelocityMethods.__init__(self) methods.IsotopeMethods.__init__(self) methods.SeedingMethods.__init__(self) - methods.DepositionMethods.__init__(self) \ No newline at end of file + methods.DepositionMethods.__init__(self) diff --git a/PySDM/dynamics/vapour_deposition_on_ice.py b/PySDM/dynamics/vapour_deposition_on_ice.py index 7426309e2..f0b188686 100644 --- a/PySDM/dynamics/vapour_deposition_on_ice.py +++ b/PySDM/dynamics/vapour_deposition_on_ice.py @@ -4,13 +4,13 @@ @register_dynamic() class VapourDepositionOnIce: def __init__(self): - """ called by the user while building a particulator """ + """called by the user while building a particulator""" self.particulator = None def register(self, *, builder): - """ called by the builder """ + """called by the builder""" self.particulator = builder.particulator def __call__(self): - """ called by the particulator during simulation""" + """called by the particulator during simulation""" self.particulator.deposition() diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 5f4f89b5a..8648d1e52 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -482,6 +482,10 @@ def seeding( def deposition(self): self.backend.deposition( - water_mass=self.attributes["water mass"] + water_mass=self.attributes["water mass"], + ambient_temperature=self.environment["T"], + ambient_total_pressure=self.environment["P"], + time_step=self.dt, + cell_id=self.attributes["cell id"], ) - self.attributes.mark_updated("water mass") \ No newline at end of file + self.attributes.mark_updated("water mass") diff --git a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py index 341272e6d..d976c1f1c 100644 --- a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py +++ b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py @@ -11,34 +11,36 @@ from PySDM.products import IceWaterContent -@pytest.mark.parametrize('dt', (1 * si.s, .1 * si.s)) -@pytest.mark.parametrize('water_mass', (-si.ng, -si.ug, -si.mg, si.mg)) -@pytest.mark.parametrize('fastmath', (True, False)) -def test_TODO(dt, water_mass, fastmath, dv=1*si.m**3): +@pytest.mark.parametrize("dt", (1 * si.s, 0.1 * si.s)) +@pytest.mark.parametrize("water_mass", (-si.ng, -si.ug, -si.mg, si.mg)) +@pytest.mark.parametrize("fastmath", (True, False)) +def test_TODO(dt, water_mass, fastmath, dv=1 * si.m**3): # arrange n_sd = 1 builder = Builder( n_sd=n_sd, environment=Box(dt=dt, dv=dv), - backend=CPU(formulae=Formulae( - fastmath=fastmath, - particle_shape_and_density="MixedPhaseSpheres" - )) + backend=CPU( + formulae=Formulae( + fastmath=fastmath, particle_shape_and_density="MixedPhaseSpheres" + ) + ), ) deposition = VapourDepositionOnIce() builder.add_dynamic(deposition) particulator = builder.build( attributes={ - 'multiplicity': np.full(shape=(n_sd,), fill_value=1), - 'water mass': np.full(shape=(n_sd,), fill_value=water_mass), + "multiplicity": np.full(shape=(n_sd,), fill_value=1), + "water mass": np.full(shape=(n_sd,), fill_value=water_mass), }, - products=(IceWaterContent(),) + products=(IceWaterContent(),), ) - + particulator.environment["T"] = 250 * si.K + particulator.environment["P"] = 500 * si.hPa # act - iwc_old = particulator.products['ice water content'].get().copy() + iwc_old = particulator.products["ice water content"].get().copy() particulator.run(steps=1) - iwc_new = particulator.products['ice water content'].get().copy() + iwc_new = particulator.products["ice water content"].get().copy() # assert - assert (iwc_new > iwc_old).all() \ No newline at end of file + assert (iwc_new > iwc_old).all() if water_mass < 0 else (iwc_new == iwc_old).all() From 5be8563b4f386c07475a16036d2b698985951a66 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Tue, 17 Sep 2024 17:35:41 +0200 Subject: [PATCH 03/27] Added commentary for pylint --- .../backends/impl_numba/methods/deposition_methods.py | 10 +++++++--- PySDM/dynamics/vapour_deposition_on_ice.py | 2 ++ .../dynamics/test_vapour_deposition_on_ice.py | 4 +++- 3 files changed, 12 insertions(+), 4 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index ba1390913..52150e24a 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -1,9 +1,11 @@ +""" basic water vapor deposition on ice for cpu backend """ + from functools import cached_property -from PySDM.backends.impl_common.backend_methods import BackendMethods import numba +from PySDM.backends.impl_common.backend_methods import BackendMethods -class DepositionMethods(BackendMethods): +class DepositionMethods(BackendMethods): # pylint:disable=too-few-public-methods @cached_property def _deposition(self): assert self.formulae.particle_shape_and_density.supports_mixed_phase() @@ -32,7 +34,9 @@ def body( temperature, pressure ) - print(volume, temperature, pressure, diffusion_coefficient) + print( + volume, temperature, pressure, diffusion_coefficient, time_step + ) water_mass[i] *= 1.1 diff --git a/PySDM/dynamics/vapour_deposition_on_ice.py b/PySDM/dynamics/vapour_deposition_on_ice.py index f0b188686..5eba7d0c1 100644 --- a/PySDM/dynamics/vapour_deposition_on_ice.py +++ b/PySDM/dynamics/vapour_deposition_on_ice.py @@ -1,3 +1,5 @@ +""" basic water vapor deposition on ice """ + from PySDM.dynamics.impl import register_dynamic diff --git a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py index d976c1f1c..7ae15e80d 100644 --- a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py +++ b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py @@ -1,3 +1,5 @@ +"""basic water vapor deposition on ice test""" + import numpy as np import pytest @@ -14,7 +16,7 @@ @pytest.mark.parametrize("dt", (1 * si.s, 0.1 * si.s)) @pytest.mark.parametrize("water_mass", (-si.ng, -si.ug, -si.mg, si.mg)) @pytest.mark.parametrize("fastmath", (True, False)) -def test_TODO(dt, water_mass, fastmath, dv=1 * si.m**3): +def test_iwc_lower_after_timestep(dt, water_mass, fastmath, dv=1 * si.m**3): # arrange n_sd = 1 builder = Builder( From 894a39eb0f220216b9a5958f351e2eec95cbfaeb Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 18 Sep 2024 10:28:16 +0200 Subject: [PATCH 04/27] improve readability by using unfrozen trivia function --- PySDM/backends/impl_numba/methods/deposition_methods.py | 4 ++-- PySDM/physics/trivia.py | 4 ++++ 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index 52150e24a..06bd50b23 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -12,6 +12,7 @@ def _deposition(self): mass_to_volume = self.formulae.particle_shape_and_density.mass_to_volume diffusion_coefficient_function = self.formulae.diffusion_thermics.D + liquid = self.formulae.trivia.unzfroze @numba.jit(**self.default_jit_flags) def body( @@ -21,8 +22,7 @@ def body( n_sd = len(water_mass) for i in numba.prange(n_sd): # pylint: disable=not-an-iterable - if water_mass[i] < 0: - + if not liquid(water_mass[i]): cid = cell_id[i] volume = mass_to_volume(water_mass[i]) diff --git a/PySDM/physics/trivia.py b/PySDM/physics/trivia.py index 3bac176aa..97cad6129 100644 --- a/PySDM/physics/trivia.py +++ b/PySDM/physics/trivia.py @@ -75,6 +75,10 @@ def p_d(const, p, water_vapour_mixing_ratio): def th_std(const, p, T): return T * np.power(const.p1000 / p, const.Rd_over_c_pd) + @staticmethod + def unfrozen(water_mass): + return water_mass > 0 + @staticmethod def unfrozen_and_saturated(water_mass, relative_humidity): return water_mass > 0 and relative_humidity > 1 From 4e0f21289f791f7b84594e605d8b876377dbf7bd Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 18 Sep 2024 10:40:39 +0200 Subject: [PATCH 05/27] employ formula_flattened trick in deposition backend method --- PySDM/backends/impl_numba/methods/deposition_methods.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index 06bd50b23..ad01cedc6 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -10,9 +10,10 @@ class DepositionMethods(BackendMethods): # pylint:disable=too-few-public-method def _deposition(self): assert self.formulae.particle_shape_and_density.supports_mixed_phase() - mass_to_volume = self.formulae.particle_shape_and_density.mass_to_volume + formulae = self.formulae_flattened + liquid = formulae.trivia__unzfroze + diffusion_coefficient_function = self.formulae.diffusion_thermics.D - liquid = self.formulae.trivia.unzfroze @numba.jit(**self.default_jit_flags) def body( @@ -25,7 +26,9 @@ def body( if not liquid(water_mass[i]): cid = cell_id[i] - volume = mass_to_volume(water_mass[i]) + volume = formulae.particle_shape_and_density__mass_to_volume( + water_mass[i] + ) temperature = ambient_temperature[cid] pressure = ambient_total_pressure[cid] From 839f78ccd4359f64eae77d194e798ccf402e8bb3 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 18 Sep 2024 10:42:07 +0200 Subject: [PATCH 06/27] renaming condensation coordinate to diffusion coordinate --- .../impl_numba/methods/condensation_methods.py | 14 ++++++-------- PySDM/formulae.py | 4 ++-- PySDM/physics/__init__.py | 2 +- .../__init__.py | 2 +- .../volume.py | 2 +- .../volume_logarithm.py | 3 ++- 6 files changed, 13 insertions(+), 14 deletions(-) rename PySDM/physics/{condensation_coordinate => diffusion_coordinate}/__init__.py (52%) rename PySDM/physics/{condensation_coordinate => diffusion_coordinate}/volume.py (83%) rename PySDM/physics/{condensation_coordinate => diffusion_coordinate}/volume_logarithm.py (84%) diff --git a/PySDM/backends/impl_numba/methods/condensation_methods.py b/PySDM/backends/impl_numba/methods/condensation_methods.py index 1b8d33ce2..525f9bddc 100644 --- a/PySDM/backends/impl_numba/methods/condensation_methods.py +++ b/PySDM/backends/impl_numba/methods/condensation_methods.py @@ -390,7 +390,7 @@ def minfun( # pylint: disable=too-many-arguments,too-many-locals K, ventilation_factor, ): - volume = formulae.condensation_coordinate__volume(x_new) + volume = formulae.diffusion_coordinate__volume(x_new) RH_eq = formulae.hygroscopicity__RH_eq( formulae.trivia__radius(volume), temperature, @@ -413,7 +413,7 @@ def minfun( # pylint: disable=too-many-arguments,too-many-locals return ( x_old - x_new - + timestep * formulae.condensation_coordinate__dx_dt(x_new, r_dr_dt) + + timestep * formulae.diffusion_coordinate__dx_dt(x_new, r_dr_dt) ) @numba.njit(**jit_flags) @@ -445,11 +445,9 @@ def calculate_ml_new( # pylint: disable=too-many-branches,too-many-arguments,to ) if v_drop <= 0: continue - x_old = formulae.condensation_coordinate__x(v_drop) + x_old = formulae.diffusion_coordinate__x(v_drop) r_old = formulae.trivia__radius(v_drop) - x_insane = formulae.condensation_coordinate__x( - attributes.vdry[drop] / 100 - ) + x_insane = formulae.diffusion_coordinate__x(attributes.vdry[drop] / 100) rd3 = attributes.vdry[drop] / formulae.constants.PI_4_3 sgm = formulae.surface_tension__sigma( T, v_drop, attributes.vdry[drop], attributes.f_org[drop] @@ -492,7 +490,7 @@ def calculate_ml_new( # pylint: disable=too-many-branches,too-many-arguments,to Kr, ventilation_factor, ) - dx_old = timestep * formulae.condensation_coordinate__dx_dt( + dx_old = timestep * formulae.diffusion_coordinate__dx_dt( x_old, r_dr_dt_old ) else: @@ -562,7 +560,7 @@ def calculate_ml_new( # pylint: disable=too-many-branches,too-many-arguments,to x_new = x_old mass_new = formulae.particle_shape_and_density__volume_to_mass( - formulae.condensation_coordinate__volume(x_new) + formulae.diffusion_coordinate__volume(x_new) ) mass_cr = formulae.particle_shape_and_density__volume_to_mass( attributes.v_cr[drop] diff --git a/PySDM/formulae.py b/PySDM/formulae.py index af8dc10c6..2c470214e 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -30,7 +30,7 @@ def __init__( # pylint: disable=too-many-locals constants: Optional[dict] = None, seed: int = None, fastmath: bool = True, - condensation_coordinate: str = "VolumeLogarithm", + diffusion_coordinate: str = "VolumeLogarithm", saturation_vapour_pressure: str = "FlatauWalkoCotton", latent_heat: str = "Kirchhoff", hygroscopicity: str = "KappaKoehlerLeadingTerms", @@ -62,7 +62,7 @@ def __init__( # pylint: disable=too-many-locals # in PyCharm and alike, all these fields are later overwritten within this ctor self.optical_albedo = optical_albedo self.optical_depth = optical_depth - self.condensation_coordinate = condensation_coordinate + self.diffusion_coordinate = diffusion_coordinate self.saturation_vapour_pressure = saturation_vapour_pressure self.hygroscopicity = hygroscopicity self.drop_growth = drop_growth diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py index 93b1c4f3d..e443ad5c1 100644 --- a/PySDM/physics/__init__.py +++ b/PySDM/physics/__init__.py @@ -17,7 +17,7 @@ """ from . import ( - condensation_coordinate, + diffusion_coordinate, constants_defaults, diffusion_kinetics, diffusion_thermics, diff --git a/PySDM/physics/condensation_coordinate/__init__.py b/PySDM/physics/diffusion_coordinate/__init__.py similarity index 52% rename from PySDM/physics/condensation_coordinate/__init__.py rename to PySDM/physics/diffusion_coordinate/__init__.py index 593ff9050..292a3e132 100644 --- a/PySDM/physics/condensation_coordinate/__init__.py +++ b/PySDM/physics/diffusion_coordinate/__init__.py @@ -1,5 +1,5 @@ """ -definitions of particle-size coordinates for the condensation solver +definitions of particle-size coordinates for the vapour diffusion solvers """ from .volume import Volume diff --git a/PySDM/physics/condensation_coordinate/volume.py b/PySDM/physics/diffusion_coordinate/volume.py similarity index 83% rename from PySDM/physics/condensation_coordinate/volume.py rename to PySDM/physics/diffusion_coordinate/volume.py index 5e16da18f..7b8ea4be1 100644 --- a/PySDM/physics/condensation_coordinate/volume.py +++ b/PySDM/physics/diffusion_coordinate/volume.py @@ -1,5 +1,5 @@ """ -particle volume as condensation coordinate (i.e. no transformation) +particle volume as diffusion coordinate (i.e. no transformation) """ import numpy as np diff --git a/PySDM/physics/condensation_coordinate/volume_logarithm.py b/PySDM/physics/diffusion_coordinate/volume_logarithm.py similarity index 84% rename from PySDM/physics/condensation_coordinate/volume_logarithm.py rename to PySDM/physics/diffusion_coordinate/volume_logarithm.py index e5f3bf041..e4ea3d722 100644 --- a/PySDM/physics/condensation_coordinate/volume_logarithm.py +++ b/PySDM/physics/diffusion_coordinate/volume_logarithm.py @@ -1,5 +1,6 @@ """ -logarithm of particle volume as coordinate (ensures non-negative values) +logarithm of particle volume as diffusion coordinate +(ensures non-negative values) """ import numpy as np From 7777b0c5fbe7582db3c076e216f1a48e49dae4d9 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 18 Sep 2024 11:51:53 +0200 Subject: [PATCH 07/27] work in progress on mass-based diffusion coordinate for deposition --- .../impl_numba/methods/deposition_methods.py | 34 ++++++++++++++++--- PySDM/dynamics/vapour_deposition_on_ice.py | 1 + PySDM/particulator.py | 2 ++ .../physics/diffusion_coordinate/__init__.py | 2 ++ PySDM/physics/diffusion_coordinate/mass.py | 22 ++++++++++++ .../diffusion_coordinate/mass_logarithm.py | 23 +++++++++++++ .../mixed_phase_spheres.py | 4 +++ .../dynamics/test_vapour_deposition_on_ice.py | 13 +++++-- 8 files changed, 94 insertions(+), 7 deletions(-) create mode 100644 PySDM/physics/diffusion_coordinate/mass.py create mode 100644 PySDM/physics/diffusion_coordinate/mass_logarithm.py diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index ad01cedc6..38bbfeb3d 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -11,13 +11,19 @@ def _deposition(self): assert self.formulae.particle_shape_and_density.supports_mixed_phase() formulae = self.formulae_flattened - liquid = formulae.trivia__unzfroze + liquid = formulae.trivia__unfrozen diffusion_coefficient_function = self.formulae.diffusion_thermics.D @numba.jit(**self.default_jit_flags) def body( - water_mass, ambient_temperature, ambient_total_pressure, time_step, cell_id + water_mass, + ambient_temperature, + ambient_total_pressure, + time_step, + cell_id, + reynolds_number, + schmidt_number, ): n_sd = len(water_mass) @@ -32,6 +38,16 @@ def body( temperature = ambient_temperature[cid] pressure = ambient_total_pressure[cid] + capacity = formulae.particle_shape_and_density__ice_mass_to_radius( + water_mass[i] + ) + + ventilation_factor = formulae.ventilation__ventilation_coefficient( + sqrt_re_times_cbrt_sc=formulae.trivia__sqrt_re_times_cbrt_sc( + Re=reynolds_number[i], + Sc=schmidt_number[cid], + ) + ) diffusion_coefficient = diffusion_coefficient_function( temperature, pressure @@ -41,7 +57,13 @@ def body( volume, temperature, pressure, diffusion_coefficient, time_step ) - water_mass[i] *= 1.1 + dm_dt = 1.1 + + x_old = formulae.diffusion_coordinate__x(water_mass[i]) + dx_dt_old = formulae.diffusion_coordinate__dx_dt(x_old, dm_dt) + x_new = formulae.trivia__explicit_euler(x_old, time_step, dx_dt_old) + + water_mass[i] = formulae.diffusion_coordinate__mass(x_new) return body @@ -52,7 +74,9 @@ def deposition( ambient_temperature, ambient_total_pressure, time_step, - cell_id + cell_id, + reynolds_number, + schmidt_number, ): self._deposition( water_mass=water_mass.data, @@ -60,4 +84,6 @@ def deposition( ambient_total_pressure=ambient_total_pressure.data, time_step=time_step, cell_id=cell_id.data, + reynolds_number=reynolds_number.data, + schmidt_number=schmidt_number.data, ) diff --git a/PySDM/dynamics/vapour_deposition_on_ice.py b/PySDM/dynamics/vapour_deposition_on_ice.py index 5eba7d0c1..802099f80 100644 --- a/PySDM/dynamics/vapour_deposition_on_ice.py +++ b/PySDM/dynamics/vapour_deposition_on_ice.py @@ -12,6 +12,7 @@ def __init__(self): def register(self, *, builder): """called by the builder""" self.particulator = builder.particulator + builder.request_attribute("Reynolds number") def __call__(self): """called by the particulator during simulation""" diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 8648d1e52..6a1a98310 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -487,5 +487,7 @@ def deposition(self): ambient_total_pressure=self.environment["P"], time_step=self.dt, cell_id=self.attributes["cell id"], + reynolds_number=self.attributes["Reynolds number"], + schmidt_number=self.environment["Schmidt number"], ) self.attributes.mark_updated("water mass") diff --git a/PySDM/physics/diffusion_coordinate/__init__.py b/PySDM/physics/diffusion_coordinate/__init__.py index 292a3e132..c397584f9 100644 --- a/PySDM/physics/diffusion_coordinate/__init__.py +++ b/PySDM/physics/diffusion_coordinate/__init__.py @@ -4,3 +4,5 @@ from .volume import Volume from .volume_logarithm import VolumeLogarithm +from .mass import Mass +from .mass_logarithm import MassLogarithm diff --git a/PySDM/physics/diffusion_coordinate/mass.py b/PySDM/physics/diffusion_coordinate/mass.py new file mode 100644 index 000000000..5cb6f5176 --- /dev/null +++ b/PySDM/physics/diffusion_coordinate/mass.py @@ -0,0 +1,22 @@ +""" +particle mass as diffusion coordinate (i.e. no transformation) +""" + +import numpy as np + + +class Mass: + def __init__(self, _): + pass + + @staticmethod + def dx_dt(const, x, dm_dt): + return dm_dt + + @staticmethod + def mass(_, x): + return x + + @staticmethod + def x(_, mass): + return mass diff --git a/PySDM/physics/diffusion_coordinate/mass_logarithm.py b/PySDM/physics/diffusion_coordinate/mass_logarithm.py new file mode 100644 index 000000000..f2544ca3d --- /dev/null +++ b/PySDM/physics/diffusion_coordinate/mass_logarithm.py @@ -0,0 +1,23 @@ +""" +logarithm of particle mass as diffusion coordinate +(ensures non-negative values) +""" + +import numpy as np + + +class MassLogarithm: + def __init__(self, _): + pass + + @staticmethod + def dx_dt(const, x, dm_dt): + return dm_dt / np.exp(x) + + @staticmethod + def mass(x): + return np.exp(x) + + @staticmethod + def x(mass): + return np.log(mass) diff --git a/PySDM/physics/particle_shape_and_density/mixed_phase_spheres.py b/PySDM/physics/particle_shape_and_density/mixed_phase_spheres.py index a9c886d1f..204938185 100644 --- a/PySDM/physics/particle_shape_and_density/mixed_phase_spheres.py +++ b/PySDM/physics/particle_shape_and_density/mixed_phase_spheres.py @@ -30,3 +30,7 @@ def volume_to_mass(const, volume): @staticmethod def radius_to_mass(const, radius): raise NotImplementedError() + + @staticmethod + def ice_mass_to_radius(const, ice_mass): + return np.power(-ice_mass / const.PI_4_3 / const.rho_i, const.ONE_THIRD) diff --git a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py index 7ae15e80d..3bdc7df3b 100644 --- a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py +++ b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py @@ -4,7 +4,7 @@ import pytest -from PySDM.physics import si +from PySDM.physics import si, diffusion_coordinate from PySDM.backends import CPU from PySDM import Builder from PySDM import Formulae @@ -16,7 +16,10 @@ @pytest.mark.parametrize("dt", (1 * si.s, 0.1 * si.s)) @pytest.mark.parametrize("water_mass", (-si.ng, -si.ug, -si.mg, si.mg)) @pytest.mark.parametrize("fastmath", (True, False)) -def test_iwc_lower_after_timestep(dt, water_mass, fastmath, dv=1 * si.m**3): +@pytest.mark.parametrize("diffusion_coordinate", ("Mass", "MassLogarithm")) +def test_iwc_lower_after_timestep( + dt, water_mass, fastmath, diffusion_coordinate, dv=1 * si.m**3 +): # arrange n_sd = 1 builder = Builder( @@ -24,7 +27,9 @@ def test_iwc_lower_after_timestep(dt, water_mass, fastmath, dv=1 * si.m**3): environment=Box(dt=dt, dv=dv), backend=CPU( formulae=Formulae( - fastmath=fastmath, particle_shape_and_density="MixedPhaseSpheres" + fastmath=fastmath, + particle_shape_and_density="MixedPhaseSpheres", + diffusion_coordinate=diffusion_coordinate, ) ), ) @@ -39,6 +44,8 @@ def test_iwc_lower_after_timestep(dt, water_mass, fastmath, dv=1 * si.m**3): ) particulator.environment["T"] = 250 * si.K particulator.environment["P"] = 500 * si.hPa + particulator.environment["Schmidt number"] = 1 + # act iwc_old = particulator.products["ice water content"].get().copy() particulator.run(steps=1) From 0c879cb566c5fedf8798b46374995ed02eb41d5d Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Wed, 18 Sep 2024 15:43:00 +0200 Subject: [PATCH 08/27] Added ambient moisture and dm/dt --- .../impl_numba/methods/deposition_methods.py | 23 +++++++++++++++---- PySDM/particulator.py | 2 ++ .../dynamics/test_vapour_deposition_on_ice.py | 7 +++++- 3 files changed, 26 insertions(+), 6 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index 38bbfeb3d..34c5b105d 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -2,6 +2,7 @@ from functools import cached_property import numba +import numpy as np from PySDM.backends.impl_common.backend_methods import BackendMethods @@ -20,6 +21,8 @@ def body( water_mass, ambient_temperature, ambient_total_pressure, + ambient_humidity, + ambient_water_activity, time_step, cell_id, reynolds_number, @@ -30,6 +33,7 @@ def body( for i in numba.prange(n_sd): # pylint: disable=not-an-iterable if not liquid(water_mass[i]): + ice_mass = -water_mass[i] cid = cell_id[i] volume = formulae.particle_shape_and_density__mass_to_volume( @@ -53,17 +57,26 @@ def body( temperature, pressure ) - print( - volume, temperature, pressure, diffusion_coefficient, time_step + saturation_ratio_ice = 1.1 # formulae.saturation_vapour_pressure.pvs_ice( temperature ) + + dm_dt = ( + 4 + * np.pi + * capacity + * diffusion_coefficient + * (saturation_ratio_ice - 1) ) - dm_dt = 1.1 + print( + f" {volume=}, {temperature=}, {pressure=}, {diffusion_coefficient=}," + ) + print(f" {time_step=}, {saturation_ratio_ice=}, {dm_dt=}") - x_old = formulae.diffusion_coordinate__x(water_mass[i]) + x_old = formulae.diffusion_coordinate__x(ice_mass) dx_dt_old = formulae.diffusion_coordinate__dx_dt(x_old, dm_dt) x_new = formulae.trivia__explicit_euler(x_old, time_step, dx_dt_old) - water_mass[i] = formulae.diffusion_coordinate__mass(x_new) + water_mass[i] = -formulae.diffusion_coordinate__mass(x_new) return body diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 6a1a98310..8e9dc12ac 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -485,6 +485,8 @@ def deposition(self): water_mass=self.attributes["water mass"], ambient_temperature=self.environment["T"], ambient_total_pressure=self.environment["P"], + ambient_humidity=self.environment["RH"], + ambient_water_activity=self.environment["a_w_ice"], time_step=self.dt, cell_id=self.attributes["cell id"], reynolds_number=self.attributes["Reynolds number"], diff --git a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py index 3bdc7df3b..3e2d6a82b 100644 --- a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py +++ b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py @@ -42,8 +42,13 @@ def test_iwc_lower_after_timestep( }, products=(IceWaterContent(),), ) - particulator.environment["T"] = 250 * si.K + temperature = 250 * si.K + particulator.environment["T"] = temperature particulator.environment["P"] = 500 * si.hPa + particulator.environment["RH"] = 1.1 * si.dimensionless + pvs_ice = particulator.formulae.saturation_vapour_pressure.pvs_ice(temperature) + pvs_water = particulator.formulae.saturation_vapour_pressure.pvs_water(temperature) + particulator.environment["a_w_ice"] = pvs_ice / pvs_water particulator.environment["Schmidt number"] = 1 # act From 4d60bd2bd6265e41f159c4c7c184a6e6d33ac541 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Wed, 18 Sep 2024 15:52:25 +0200 Subject: [PATCH 09/27] Added saturation ratio over ice to deposition --- PySDM/backends/impl_numba/methods/deposition_methods.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index 34c5b105d..5d1ae9c2f 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -57,7 +57,9 @@ def body( temperature, pressure ) - saturation_ratio_ice = 1.1 # formulae.saturation_vapour_pressure.pvs_ice( temperature ) + saturation_ratio_ice = ( + ambient_humidity[cid] / ambient_water_activity[cid] + ) dm_dt = ( 4 @@ -86,6 +88,8 @@ def deposition( water_mass, ambient_temperature, ambient_total_pressure, + ambient_humidity, + ambient_water_activity, time_step, cell_id, reynolds_number, @@ -95,6 +99,8 @@ def deposition( water_mass=water_mass.data, ambient_temperature=ambient_temperature.data, ambient_total_pressure=ambient_total_pressure.data, + ambient_humidity=ambient_humidity.data, + ambient_water_activity=ambient_water_activity.data, time_step=time_step, cell_id=cell_id.data, reynolds_number=reynolds_number.data, From 857ca7d4cdddbdcdd59b78cf4bd3eaf3f737d660 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Thu, 19 Sep 2024 15:19:54 +0200 Subject: [PATCH 10/27] Added subsaturation and saturation to deposition on ice test --- .../impl_numba/methods/deposition_methods.py | 3 ++- .../dynamics/test_vapour_deposition_on_ice.py | 17 ++++++++++++++--- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index 5d1ae9c2f..e993a7e10 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -68,7 +68,8 @@ def body( * diffusion_coefficient * (saturation_ratio_ice - 1) ) - + if dm_dt == 0: + continue print( f" {volume=}, {temperature=}, {pressure=}, {diffusion_coefficient=}," ) diff --git a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py index 3e2d6a82b..cf40ef124 100644 --- a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py +++ b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py @@ -15,10 +15,11 @@ @pytest.mark.parametrize("dt", (1 * si.s, 0.1 * si.s)) @pytest.mark.parametrize("water_mass", (-si.ng, -si.ug, -si.mg, si.mg)) +@pytest.mark.parametrize("RHi", (1.1, 1.0, 0.9)) @pytest.mark.parametrize("fastmath", (True, False)) @pytest.mark.parametrize("diffusion_coordinate", ("Mass", "MassLogarithm")) def test_iwc_lower_after_timestep( - dt, water_mass, fastmath, diffusion_coordinate, dv=1 * si.m**3 + dt, water_mass, RHi, fastmath, diffusion_coordinate, dv=1 * si.m**3 ): # arrange n_sd = 1 @@ -45,16 +46,26 @@ def test_iwc_lower_after_timestep( temperature = 250 * si.K particulator.environment["T"] = temperature particulator.environment["P"] = 500 * si.hPa - particulator.environment["RH"] = 1.1 * si.dimensionless pvs_ice = particulator.formulae.saturation_vapour_pressure.pvs_ice(temperature) pvs_water = particulator.formulae.saturation_vapour_pressure.pvs_water(temperature) + vapour_pressure = RHi * pvs_ice + RH = vapour_pressure / pvs_water + particulator.environment["RH"] = RH # 1.1 * si.dimensionless particulator.environment["a_w_ice"] = pvs_ice / pvs_water particulator.environment["Schmidt number"] = 1 + print(f" {RHi=}, {RH=}, {vapour_pressure=}, {pvs_ice=}, {pvs_water=}") + # act iwc_old = particulator.products["ice water content"].get().copy() particulator.run(steps=1) iwc_new = particulator.products["ice water content"].get().copy() # assert - assert (iwc_new > iwc_old).all() if water_mass < 0 else (iwc_new == iwc_old).all() + if water_mass < 0 and RHi != 1: + if RHi > 1: + assert (iwc_new > iwc_old).all() + elif RHi < 1: + assert (iwc_new < iwc_old).all() + else: + assert (iwc_new == iwc_old).all() From d7783f466c4e553b434e0831f09a80897521933b Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Thu, 19 Sep 2024 16:10:40 +0200 Subject: [PATCH 11/27] Added change of water vapour mixing ratio due to deposition on ice --- .../impl_numba/methods/deposition_methods.py | 17 +++++++++++++++-- PySDM/particulator.py | 3 +++ .../dynamics/test_vapour_deposition_on_ice.py | 16 +++++++++++++++- 3 files changed, 33 insertions(+), 3 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index e993a7e10..6adc0bdb1 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -23,14 +23,18 @@ def body( ambient_total_pressure, ambient_humidity, ambient_water_activity, + ambient_vapour_mixing_ratio, + ambient_dry_air_density, + cell_volume, time_step, cell_id, reynolds_number, schmidt_number, ): - + delta_vapour_mass = 0.0 n_sd = len(water_mass) - for i in numba.prange(n_sd): # pylint: disable=not-an-iterable + # for i in numba.prange(n_sd): # pylint: disable=not-an-iterable + for i in range(n_sd): if not liquid(water_mass[i]): ice_mass = -water_mass[i] @@ -70,6 +74,9 @@ def body( ) if dm_dt == 0: continue + ambient_vapour_mixing_ratio[cell_id] -= ( + dm_dt * time_step / (cell_volume * ambient_dry_air_density) + ) print( f" {volume=}, {temperature=}, {pressure=}, {diffusion_coefficient=}," ) @@ -91,6 +98,9 @@ def deposition( ambient_total_pressure, ambient_humidity, ambient_water_activity, + ambient_vapour_mixing_ratio, + ambient_dry_air_density, + cell_volume, time_step, cell_id, reynolds_number, @@ -102,6 +112,9 @@ def deposition( ambient_total_pressure=ambient_total_pressure.data, ambient_humidity=ambient_humidity.data, ambient_water_activity=ambient_water_activity.data, + ambient_vapour_mixing_ratio=ambient_vapour_mixing_ratio.data, + ambient_dry_air_density=ambient_dry_air_density.data, + cell_volume=cell_volume, time_step=time_step, cell_id=cell_id.data, reynolds_number=reynolds_number.data, diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 8e9dc12ac..15e67d4b8 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -487,6 +487,9 @@ def deposition(self): ambient_total_pressure=self.environment["P"], ambient_humidity=self.environment["RH"], ambient_water_activity=self.environment["a_w_ice"], + ambient_vapour_mixing_ratio=self.environment["water_vapour_mixing_ratio"], + ambient_dry_air_density=self.environment["rhod"], + cell_volume=self.environment.mesh.dv, time_step=self.dt, cell_id=self.attributes["cell id"], reynolds_number=self.attributes["Reynolds number"], diff --git a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py index cf40ef124..a3dac4439 100644 --- a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py +++ b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py @@ -45,7 +45,8 @@ def test_iwc_lower_after_timestep( ) temperature = 250 * si.K particulator.environment["T"] = temperature - particulator.environment["P"] = 500 * si.hPa + pressure = 500 * si.hPa + particulator.environment["P"] = pressure pvs_ice = particulator.formulae.saturation_vapour_pressure.pvs_ice(temperature) pvs_water = particulator.formulae.saturation_vapour_pressure.pvs_water(temperature) vapour_pressure = RHi * pvs_ice @@ -53,6 +54,15 @@ def test_iwc_lower_after_timestep( particulator.environment["RH"] = RH # 1.1 * si.dimensionless particulator.environment["a_w_ice"] = pvs_ice / pvs_water particulator.environment["Schmidt number"] = 1 + rv0 = ( + particulator.formulae.constants.eps + * vapour_pressure + / (pressure - vapour_pressure) + ) + particulator.environment["water_vapour_mixing_ratio"] = rv0 + particulator.environment["rhod"] = (pressure - vapour_pressure) / ( + temperature * particulator.formulae.constants.Rd + ) print(f" {RHi=}, {RH=}, {vapour_pressure=}, {pvs_ice=}, {pvs_water=}") @@ -61,11 +71,15 @@ def test_iwc_lower_after_timestep( particulator.run(steps=1) iwc_new = particulator.products["ice water content"].get().copy() + rv_new = particulator.environment["water_vapour_mixing_ratio"][0] # assert if water_mass < 0 and RHi != 1: if RHi > 1: assert (iwc_new > iwc_old).all() + assert rv_new < rv0 elif RHi < 1: assert (iwc_new < iwc_old).all() + assert rv_new > rv0 else: assert (iwc_new == iwc_old).all() + assert rv_new == rv0 From 6e16036953e8b3ce24f3563ba6e5fe44dee39515 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 19 Sep 2024 18:32:46 +0200 Subject: [PATCH 12/27] temporarily make condensation and diffusion coordinates separate --- .../impl_numba/methods/condensation_methods.py | 14 ++++++++------ PySDM/physics/__init__.py | 1 + PySDM/physics/condensation_coordinate/__init__.py | 2 ++ PySDM/physics/diffusion_coordinate/__init__.py | 2 -- 4 files changed, 11 insertions(+), 8 deletions(-) create mode 100644 PySDM/physics/condensation_coordinate/__init__.py diff --git a/PySDM/backends/impl_numba/methods/condensation_methods.py b/PySDM/backends/impl_numba/methods/condensation_methods.py index 525f9bddc..1b8d33ce2 100644 --- a/PySDM/backends/impl_numba/methods/condensation_methods.py +++ b/PySDM/backends/impl_numba/methods/condensation_methods.py @@ -390,7 +390,7 @@ def minfun( # pylint: disable=too-many-arguments,too-many-locals K, ventilation_factor, ): - volume = formulae.diffusion_coordinate__volume(x_new) + volume = formulae.condensation_coordinate__volume(x_new) RH_eq = formulae.hygroscopicity__RH_eq( formulae.trivia__radius(volume), temperature, @@ -413,7 +413,7 @@ def minfun( # pylint: disable=too-many-arguments,too-many-locals return ( x_old - x_new - + timestep * formulae.diffusion_coordinate__dx_dt(x_new, r_dr_dt) + + timestep * formulae.condensation_coordinate__dx_dt(x_new, r_dr_dt) ) @numba.njit(**jit_flags) @@ -445,9 +445,11 @@ def calculate_ml_new( # pylint: disable=too-many-branches,too-many-arguments,to ) if v_drop <= 0: continue - x_old = formulae.diffusion_coordinate__x(v_drop) + x_old = formulae.condensation_coordinate__x(v_drop) r_old = formulae.trivia__radius(v_drop) - x_insane = formulae.diffusion_coordinate__x(attributes.vdry[drop] / 100) + x_insane = formulae.condensation_coordinate__x( + attributes.vdry[drop] / 100 + ) rd3 = attributes.vdry[drop] / formulae.constants.PI_4_3 sgm = formulae.surface_tension__sigma( T, v_drop, attributes.vdry[drop], attributes.f_org[drop] @@ -490,7 +492,7 @@ def calculate_ml_new( # pylint: disable=too-many-branches,too-many-arguments,to Kr, ventilation_factor, ) - dx_old = timestep * formulae.diffusion_coordinate__dx_dt( + dx_old = timestep * formulae.condensation_coordinate__dx_dt( x_old, r_dr_dt_old ) else: @@ -560,7 +562,7 @@ def calculate_ml_new( # pylint: disable=too-many-branches,too-many-arguments,to x_new = x_old mass_new = formulae.particle_shape_and_density__volume_to_mass( - formulae.diffusion_coordinate__volume(x_new) + formulae.condensation_coordinate__volume(x_new) ) mass_cr = formulae.particle_shape_and_density__volume_to_mass( attributes.v_cr[drop] diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py index e443ad5c1..e2b3d9e48 100644 --- a/PySDM/physics/__init__.py +++ b/PySDM/physics/__init__.py @@ -18,6 +18,7 @@ from . import ( diffusion_coordinate, + condensation_coordinate, constants_defaults, diffusion_kinetics, diffusion_thermics, diff --git a/PySDM/physics/condensation_coordinate/__init__.py b/PySDM/physics/condensation_coordinate/__init__.py new file mode 100644 index 000000000..d7292934d --- /dev/null +++ b/PySDM/physics/condensation_coordinate/__init__.py @@ -0,0 +1,2 @@ +from .volume import Volume +from .volume_logarithm import VolumeLogarithm diff --git a/PySDM/physics/diffusion_coordinate/__init__.py b/PySDM/physics/diffusion_coordinate/__init__.py index c397584f9..25eeb5fb0 100644 --- a/PySDM/physics/diffusion_coordinate/__init__.py +++ b/PySDM/physics/diffusion_coordinate/__init__.py @@ -2,7 +2,5 @@ definitions of particle-size coordinates for the vapour diffusion solvers """ -from .volume import Volume -from .volume_logarithm import VolumeLogarithm from .mass import Mass from .mass_logarithm import MassLogarithm From 7c2527b3afc0877cfe59d8afdb93d2e288c3b8dc Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 19 Sep 2024 18:33:06 +0200 Subject: [PATCH 13/27] fix for last commit --- .../{diffusion_coordinate => condensation_coordinate}/volume.py | 0 .../volume_logarithm.py | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename PySDM/physics/{diffusion_coordinate => condensation_coordinate}/volume.py (100%) rename PySDM/physics/{diffusion_coordinate => condensation_coordinate}/volume_logarithm.py (100%) diff --git a/PySDM/physics/diffusion_coordinate/volume.py b/PySDM/physics/condensation_coordinate/volume.py similarity index 100% rename from PySDM/physics/diffusion_coordinate/volume.py rename to PySDM/physics/condensation_coordinate/volume.py diff --git a/PySDM/physics/diffusion_coordinate/volume_logarithm.py b/PySDM/physics/condensation_coordinate/volume_logarithm.py similarity index 100% rename from PySDM/physics/diffusion_coordinate/volume_logarithm.py rename to PySDM/physics/condensation_coordinate/volume_logarithm.py From 928495907ff7d68c1454e11b7042384edcacf580 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 19 Sep 2024 18:33:39 +0200 Subject: [PATCH 14/27] check if Formulae.particle_shape_and_density supports mixed phase --- PySDM/dynamics/vapour_deposition_on_ice.py | 1 + 1 file changed, 1 insertion(+) diff --git a/PySDM/dynamics/vapour_deposition_on_ice.py b/PySDM/dynamics/vapour_deposition_on_ice.py index 802099f80..ac38946df 100644 --- a/PySDM/dynamics/vapour_deposition_on_ice.py +++ b/PySDM/dynamics/vapour_deposition_on_ice.py @@ -12,6 +12,7 @@ def __init__(self): def register(self, *, builder): """called by the builder""" self.particulator = builder.particulator + assert builder.formulae.particle_shape_and_density.supports_mixed_phase() builder.request_attribute("Reynolds number") def __call__(self): From abef7711de803ae116fe7cb508b99f4fe4184971 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 19 Sep 2024 18:33:53 +0200 Subject: [PATCH 15/27] add one TODO --- PySDM/physics/particle_shape_and_density/mixed_phase_spheres.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PySDM/physics/particle_shape_and_density/mixed_phase_spheres.py b/PySDM/physics/particle_shape_and_density/mixed_phase_spheres.py index 204938185..767c37cee 100644 --- a/PySDM/physics/particle_shape_and_density/mixed_phase_spheres.py +++ b/PySDM/physics/particle_shape_and_density/mixed_phase_spheres.py @@ -10,7 +10,7 @@ def __init__(self, _): pass @staticmethod - def supports_mixed_phase(_=None): + def supports_mixed_phase(_=None): # TODO: rename to "ice_phase?" return True @staticmethod From 055b7991b87845312f8504a7e585a690d7049beb Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 19 Sep 2024 18:34:18 +0200 Subject: [PATCH 16/27] one more commit for coordinates --- PySDM/formulae.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/PySDM/formulae.py b/PySDM/formulae.py index 2c470214e..47b26bdb4 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -30,7 +30,8 @@ def __init__( # pylint: disable=too-many-locals constants: Optional[dict] = None, seed: int = None, fastmath: bool = True, - diffusion_coordinate: str = "VolumeLogarithm", + diffusion_coordinate: str = "MassLogarithm", + condensation_coordinate: str = "VolumeLogarithm", saturation_vapour_pressure: str = "FlatauWalkoCotton", latent_heat: str = "Kirchhoff", hygroscopicity: str = "KappaKoehlerLeadingTerms", @@ -63,6 +64,7 @@ def __init__( # pylint: disable=too-many-locals self.optical_albedo = optical_albedo self.optical_depth = optical_depth self.diffusion_coordinate = diffusion_coordinate + self.condensation_coordinate = condensation_coordinate self.saturation_vapour_pressure = saturation_vapour_pressure self.hygroscopicity = hygroscopicity self.drop_growth = drop_growth From 86ebc89d4f40176ccd9493c5398050bf15db163a Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 19 Sep 2024 18:34:29 +0200 Subject: [PATCH 17/27] one more TODO --- PySDM/products/size_spectral/cloud_water_content.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/PySDM/products/size_spectral/cloud_water_content.py b/PySDM/products/size_spectral/cloud_water_content.py index 64b82b9e8..a9e2be629 100644 --- a/PySDM/products/size_spectral/cloud_water_content.py +++ b/PySDM/products/size_spectral/cloud_water_content.py @@ -5,6 +5,8 @@ LiquidWaterContent is just liquid IceWaterContent is just ice +TODO: why not to call ice a mixing ration then?! + """ import numpy as np From 0697f212657a55b1a2f8aea3ff4a81c363d4ddfd Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 19 Sep 2024 18:35:13 +0200 Subject: [PATCH 18/27] work in progress on A&A example: support for freezing and vapour deposition --- .../impl_numba/methods/deposition_methods.py | 29 +- .../Abade_and_Albuquerque_2024/fig_2.ipynb | 1826 ++++++++++++----- .../Abade_and_Albuquerque_2024/settings.py | 13 +- .../Abade_and_Albuquerque_2024/simulation.py | 32 +- 4 files changed, 1391 insertions(+), 509 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index 6adc0bdb1..c06716874 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -31,19 +31,13 @@ def body( reynolds_number, schmidt_number, ): - delta_vapour_mass = 0.0 n_sd = len(water_mass) # for i in numba.prange(n_sd): # pylint: disable=not-an-iterable for i in range(n_sd): - if not liquid(water_mass[i]): ice_mass = -water_mass[i] cid = cell_id[i] - volume = formulae.particle_shape_and_density__mass_to_volume( - water_mass[i] - ) - temperature = ambient_temperature[cid] pressure = ambient_total_pressure[cid] capacity = formulae.particle_shape_and_density__ice_mass_to_radius( @@ -72,20 +66,31 @@ def body( * diffusion_coefficient * (saturation_ratio_ice - 1) ) + # print("dm_dt", dm_dt) if dm_dt == 0: continue - ambient_vapour_mixing_ratio[cell_id] -= ( - dm_dt * time_step / (cell_volume * ambient_dry_air_density) - ) - print( - f" {volume=}, {temperature=}, {pressure=}, {diffusion_coefficient=}," + delta_rv_i = ( + -dm_dt + * time_step + / (cell_volume * ambient_dry_air_density[cid]) ) - print(f" {time_step=}, {saturation_ratio_ice=}, {dm_dt=}") + if -delta_rv_i > ambient_vapour_mixing_ratio[cid]: + assert False + # print('delta_rv_i', delta_rv_i) + # print('A', ambient_vapour_mixing_ratio[cid]) + ambient_vapour_mixing_ratio[cid] += delta_rv_i + # print('B', ambient_vapour_mixing_ratio[cid]) + # print( + # f" {volume=}, {temperature=}, {pressure=}, {diffusion_coefficient=}," + # ) + # print(f" {time_step=}, {saturation_ratio_ice=}, {dm_dt=}") x_old = formulae.diffusion_coordinate__x(ice_mass) dx_dt_old = formulae.diffusion_coordinate__dx_dt(x_old, dm_dt) x_new = formulae.trivia__explicit_euler(x_old, time_step, dx_dt_old) + # print(x_old, x_new, dm_dt) + water_mass[i] = -formulae.diffusion_coordinate__mass(x_new) return body diff --git a/examples/PySDM_examples/Abade_and_Albuquerque_2024/fig_2.ipynb b/examples/PySDM_examples/Abade_and_Albuquerque_2024/fig_2.ipynb index 88e11f4d2..b0bfbd745 100644 --- a/examples/PySDM_examples/Abade_and_Albuquerque_2024/fig_2.ipynb +++ b/examples/PySDM_examples/Abade_and_Albuquerque_2024/fig_2.ipynb @@ -44,13 +44,30 @@ "metadata": {}, "outputs": [], "source": [ - "simulation = Simulation(Settings(timestep=100 * si.s, n_sd=100))\n", - "output = simulation.run(nt=60, steps_per_output_interval=5)" + "simulations = {\n", + " 'Bulk': Simulation(Settings(\n", + " timestep=100 * si.s,\n", + " n_sd=100,\n", + " enable_immersion_freezing=False,\n", + " enable_vapour_deposition_on_ice=False,\n", + " )),\n", + " 'Homogeneous': Simulation(Settings(\n", + " timestep=100 * si.s,\n", + " n_sd=100,\n", + " enable_immersion_freezing=True,\n", + " enable_vapour_deposition_on_ice=True,\n", + " )),\n", + "}\n", + "\n", + "output = {\n", + " key: simulations[key].run(nt=60, steps_per_output_interval=5)\n", + " for key in simulations\n", + "}" ] }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "id": "b7abadbe-aa26-433d-9a9a-1984e3360db5", "metadata": {}, "outputs": [ @@ -65,7 +82,7 @@ " \n", " \n", " \n", - " 2024-07-24T18:43:17.015303\n", + " 2024-09-19T18:28:58.371636\n", " image/svg+xml\n", " \n", " \n", @@ -90,8 +107,8 @@ " \n", " \n", " \n", @@ -99,23 +116,23 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -201,43 +242,59 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", - " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", + " \n", " \n", " \n", " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", + "L 174.547576 385.821875 \n", + "\" clip-path=\"url(#p39676d5648)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -386,17 +531,17 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "L 174.547576 325.237917 \n", + "\" clip-path=\"url(#p39676d5648)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -406,17 +551,17 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "L 174.547576 264.653958 \n", + "\" clip-path=\"url(#p39676d5648)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -426,17 +571,17 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "L 174.547576 204.07 \n", + "\" clip-path=\"url(#p39676d5648)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -446,45 +591,19 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "L 174.547576 143.486042 \n", + "\" clip-path=\"url(#p39676d5648)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", " \n", " \n", " \n", @@ -492,17 +611,17 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "L 174.547576 82.902083 \n", + "\" clip-path=\"url(#p39676d5648)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -512,17 +631,17 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "L 174.547576 22.318125 \n", + "\" clip-path=\"url(#p39676d5648)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -565,7 +684,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -584,31 +703,6 @@ "L 628 4666 \n", "z\n", "\" transform=\"scale(0.015625)\"/>\n", - " \n", " \n", " \n", " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", + " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", " \n", - " \n", + " \n", + " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -927,113 +1204,113 @@ " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", " \n", - " \n", " \n", " \n", - " \n", " \n", " \n", - " \n", " \n", " \n", - " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1222,231 +1499,192 @@ " \n", " \n", " \n", - " \n", - " \n", + " \n", + "\" clip-path=\"url(#p27bc25918a)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + "\" clip-path=\"url(#p27bc25918a)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + "\" clip-path=\"url(#p27bc25918a)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + "\" clip-path=\"url(#p27bc25918a)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + "\" clip-path=\"url(#p27bc25918a)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + "\" clip-path=\"url(#p27bc25918a)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + "\" clip-path=\"url(#p27bc25918a)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", - " \n", " \n", " \n", @@ -1455,18 +1693,18 @@ "\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n", " \n", " \n", - " \n", " \n", " \n", - " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1605,14 +1843,14 @@ " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", " \n", "\n" @@ -1627,12 +1865,12 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "81ea5d5837394cd6960646bb7b1fd972", + "model_id": "00632a293ef541bb91beabbaea976d87", "version_major": 2, "version_minor": 0 }, "text/plain": [ - "HTML(value=\"./fig_2.pdf
\")" + "HBox(children=(HTML(value=\"./fig_2.pdf
\"), HTML(value=\"\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-09-19T18:28:58.643361\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.8.1, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pyplot.plot(output['Bulk']['T'], output['Bulk']['height'])\n", + "pyplot.grid()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "721e5606-9e19-44b9-87f3-cba86b54574e", + "metadata": {}, "outputs": [], "source": [] } diff --git a/examples/PySDM_examples/Abade_and_Albuquerque_2024/settings.py b/examples/PySDM_examples/Abade_and_Albuquerque_2024/settings.py index e909f2156..1386c9c66 100644 --- a/examples/PySDM_examples/Abade_and_Albuquerque_2024/settings.py +++ b/examples/PySDM_examples/Abade_and_Albuquerque_2024/settings.py @@ -7,9 +7,18 @@ @strict class Settings: - def __init__(self, *, n_sd: int, timestep: float): + def __init__( + self, + *, + n_sd: int, + timestep: float, + enable_immersion_freezing: bool = True, + enable_vapour_deposition_on_ice: bool = True, + ): self.n_sd = n_sd self.timestep = timestep + self.enable_immersion_freezing = enable_immersion_freezing + self.enable_vapour_deposition_on_ice = enable_vapour_deposition_on_ice self.initial_total_pressure = 1000 * si.hPa # note: not given in the paper @@ -17,6 +26,8 @@ def __init__(self, *, n_sd: int, timestep: float): self.formulae = Formulae( constants={"bulk_phase_partitioning_exponent": 0.1}, bulk_phase_partitioning="KaulEtAl2015", + particle_shape_and_density="MixedPhaseSpheres", + diffusion_coordinate="Mass", ) self.initial_water_vapour_mixing_ratio = 1.5 * si.g / si.kg self.parcel_linear_extent = 100 * si.m diff --git a/examples/PySDM_examples/Abade_and_Albuquerque_2024/simulation.py b/examples/PySDM_examples/Abade_and_Albuquerque_2024/simulation.py index ac1010b74..b1a329b2c 100644 --- a/examples/PySDM_examples/Abade_and_Albuquerque_2024/simulation.py +++ b/examples/PySDM_examples/Abade_and_Albuquerque_2024/simulation.py @@ -4,9 +4,19 @@ from PySDM import Builder from PySDM.backends import CPU -from PySDM.dynamics import Condensation, AmbientThermodynamics +from PySDM.dynamics import ( + Condensation, + AmbientThermodynamics, + VapourDepositionOnIce, + Freezing, +) from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity -from PySDM.products import AmbientTemperature, ParcelDisplacement, WaterMixingRatio +from PySDM.products import ( + AmbientTemperature, + ParcelDisplacement, + WaterMixingRatio, + SpecificIceWaterContent, +) from PySDM.environments import Parcel @@ -22,9 +32,16 @@ def __init__(self, settings): initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio, T0=settings.initial_temperature, w=settings.updraft, + mixed_phase=True, ), ) builder.add_dynamic(AmbientThermodynamics()) + + if settings.enable_immersion_freezing: + builder.add_dynamic(Freezing()) + if settings.enable_vapour_deposition_on_ice: + builder.add_dynamic(VapourDepositionOnIce()) + builder.add_dynamic(Condensation()) r_dry, n_in_dv = ConstantMultiplicity(settings.soluble_aerosol).sample( @@ -33,10 +50,17 @@ def __init__(self, settings): attributes = builder.particulator.environment.init_attributes( n_in_dv=n_in_dv, kappa=settings.kappa, r_dry=r_dry ) + + # + if settings.enable_immersion_freezing: + attributes["freezing temperature"] = np.full( + shape=(settings.n_sd,), fill_value=250 + ) + # + self.products = ( WaterMixingRatio(name="water", radius_range=(0, np.inf)), - WaterMixingRatio(name="ice", radius_range=(-np.inf, 0)), - WaterMixingRatio(name="total", radius_range=(-np.inf, np.inf)), + SpecificIceWaterContent(name="ice"), ParcelDisplacement(name="height"), AmbientTemperature(name="T"), ) From 70ed3a1fd4ece5a3dbfe3810ef9602845263845f Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 20 Sep 2024 00:26:35 +0200 Subject: [PATCH 19/27] missing multiplication by multiplicity; missing multiplication by density --- .../impl_numba/methods/deposition_methods.py | 23 +- PySDM/dynamics/vapour_deposition_on_ice.py | 1 + PySDM/particulator.py | 5 +- .../Abade_and_Albuquerque_2024/fig_2.ipynb | 2120 ++++++++++------- .../Abade_and_Albuquerque_2024/simulation.py | 7 +- 5 files changed, 1233 insertions(+), 923 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index c06716874..34df43a54 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -18,6 +18,7 @@ def _deposition(self): @numba.jit(**self.default_jit_flags) def body( + multiplicity, water_mass, ambient_temperature, ambient_total_pressure, @@ -59,38 +60,36 @@ def body( ambient_humidity[cid] / ambient_water_activity[cid] ) + rho_vs_ice = ( + formulae.saturation_vapour_pressure__pvs_ice(temperature) + / formulae.constants.Rv + / temperature + ) dm_dt = ( 4 * np.pi * capacity * diffusion_coefficient * (saturation_ratio_ice - 1) - ) - # print("dm_dt", dm_dt) + ) * rho_vs_ice + if dm_dt == 0: continue + delta_rv_i = ( -dm_dt + * multiplicity[i] * time_step / (cell_volume * ambient_dry_air_density[cid]) ) if -delta_rv_i > ambient_vapour_mixing_ratio[cid]: assert False - # print('delta_rv_i', delta_rv_i) - # print('A', ambient_vapour_mixing_ratio[cid]) ambient_vapour_mixing_ratio[cid] += delta_rv_i - # print('B', ambient_vapour_mixing_ratio[cid]) - # print( - # f" {volume=}, {temperature=}, {pressure=}, {diffusion_coefficient=}," - # ) - # print(f" {time_step=}, {saturation_ratio_ice=}, {dm_dt=}") x_old = formulae.diffusion_coordinate__x(ice_mass) dx_dt_old = formulae.diffusion_coordinate__dx_dt(x_old, dm_dt) x_new = formulae.trivia__explicit_euler(x_old, time_step, dx_dt_old) - # print(x_old, x_new, dm_dt) - water_mass[i] = -formulae.diffusion_coordinate__mass(x_new) return body @@ -98,6 +97,7 @@ def body( def deposition( self, *, + multiplicity, water_mass, ambient_temperature, ambient_total_pressure, @@ -112,6 +112,7 @@ def deposition( schmidt_number, ): self._deposition( + multiplicity=multiplicity.data, water_mass=water_mass.data, ambient_temperature=ambient_temperature.data, ambient_total_pressure=ambient_total_pressure.data, diff --git a/PySDM/dynamics/vapour_deposition_on_ice.py b/PySDM/dynamics/vapour_deposition_on_ice.py index ac38946df..a34eb0820 100644 --- a/PySDM/dynamics/vapour_deposition_on_ice.py +++ b/PySDM/dynamics/vapour_deposition_on_ice.py @@ -18,3 +18,4 @@ def register(self, *, builder): def __call__(self): """called by the particulator during simulation""" self.particulator.deposition() + self.particulator.update_TpRH() diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 15e67d4b8..0384650ca 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -482,12 +482,15 @@ def seeding( def deposition(self): self.backend.deposition( + multiplicity=self.attributes["multiplicity"], water_mass=self.attributes["water mass"], ambient_temperature=self.environment["T"], ambient_total_pressure=self.environment["P"], ambient_humidity=self.environment["RH"], ambient_water_activity=self.environment["a_w_ice"], - ambient_vapour_mixing_ratio=self.environment["water_vapour_mixing_ratio"], + ambient_vapour_mixing_ratio=self.environment.get_predicted( + "water_vapour_mixing_ratio" + ), ambient_dry_air_density=self.environment["rhod"], cell_volume=self.environment.mesh.dv, time_step=self.dt, diff --git a/examples/PySDM_examples/Abade_and_Albuquerque_2024/fig_2.ipynb b/examples/PySDM_examples/Abade_and_Albuquerque_2024/fig_2.ipynb index b0bfbd745..2a80b6d45 100644 --- a/examples/PySDM_examples/Abade_and_Albuquerque_2024/fig_2.ipynb +++ b/examples/PySDM_examples/Abade_and_Albuquerque_2024/fig_2.ipynb @@ -39,20 +39,20 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 6, "id": "478730ad-0c93-4adf-82a1-c606fde3c0b9", "metadata": {}, "outputs": [], "source": [ "simulations = {\n", " 'Bulk': Simulation(Settings(\n", - " timestep=100 * si.s,\n", + " timestep=1 * si.s,\n", " n_sd=100,\n", " enable_immersion_freezing=False,\n", " enable_vapour_deposition_on_ice=False,\n", " )),\n", " 'Homogeneous': Simulation(Settings(\n", - " timestep=100 * si.s,\n", + " timestep=1 * si.s,\n", " n_sd=100,\n", " enable_immersion_freezing=True,\n", " enable_vapour_deposition_on_ice=True,\n", @@ -60,14 +60,14 @@ "}\n", "\n", "output = {\n", - " key: simulations[key].run(nt=60, steps_per_output_interval=5)\n", + " key: simulations[key].run(nt=6000, steps_per_output_interval=500)\n", " for key in simulations\n", "}" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 7, "id": "b7abadbe-aa26-433d-9a9a-1984e3360db5", "metadata": {}, "outputs": [ @@ -82,7 +82,7 @@ " \n", " \n", " \n", - " 2024-09-19T18:28:58.371636\n", + " 2024-09-20T00:24:17.774939\n", " image/svg+xml\n", " \n", " \n", @@ -107,8 +107,8 @@ " \n", " \n", " \n", @@ -116,23 +116,23 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -242,59 +218,43 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", + "L 182.854167 385.821875 \n", + "\" clip-path=\"url(#pc02b7673f4)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -531,17 +403,17 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "L 182.854167 325.237917 \n", + "\" clip-path=\"url(#pc02b7673f4)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -551,17 +423,17 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "L 182.854167 264.653958 \n", + "\" clip-path=\"url(#pc02b7673f4)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -571,17 +443,17 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "L 182.854167 204.07 \n", + "\" clip-path=\"url(#pc02b7673f4)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -591,19 +463,45 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "L 182.854167 143.486042 \n", + "\" clip-path=\"url(#pc02b7673f4)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", @@ -611,17 +509,17 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "L 182.854167 82.902083 \n", + "\" clip-path=\"url(#pc02b7673f4)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -631,17 +529,17 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "L 182.854167 22.318125 \n", + "\" clip-path=\"url(#pc02b7673f4)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -684,7 +582,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -703,6 +601,31 @@ "L 628 4666 \n", "z\n", "\" transform=\"scale(0.015625)\"/>\n", + " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", @@ -945,23 +868,23 @@ "\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n", " \n", " \n", - " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", - " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", - " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1204,113 +1082,113 @@ " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", " \n", - " \n", " \n", " \n", - " \n", " \n", " \n", - " \n", " \n", " \n", - " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1499,192 +1377,192 @@ " \n", " \n", " \n", - " \n", - " \n", + " \n", + "\" clip-path=\"url(#pe710c15827)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + "\" clip-path=\"url(#pe710c15827)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + "\" clip-path=\"url(#pe710c15827)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + "\" clip-path=\"url(#pe710c15827)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + "\" clip-path=\"url(#pe710c15827)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + "\" clip-path=\"url(#pe710c15827)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + "\" clip-path=\"url(#pe710c15827)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", - " \n", " \n", " \n", @@ -1693,18 +1571,18 @@ "\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n", " \n", " \n", - " \n", " \n", " \n", - " \n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1783,16 +1661,16 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1820,16 +1698,16 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1843,14 +1721,14 @@ " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", " \n", "\n" @@ -1865,7 +1743,7 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "00632a293ef541bb91beabbaea976d87", + "model_id": "4bb849fb742d425c9c0b5fbd5690a8ea", "version_major": 2, "version_minor": 0 }, @@ -1881,7 +1759,7 @@ "fig, axs = pyplot.subplot_mosaic(\n", " (('Homogeneous', 'Stochastic', 'Bulk'),),\n", " figsize=(7, 6),\n", - " #sharex=True,\n", + " sharex=True,\n", " sharey=True,\n", " tight_layout=True,\n", ")\n", @@ -1931,7 +1809,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 8, "id": "892a1d98-c690-446e-9507-d40dc6269878", "metadata": {}, "outputs": [ @@ -1941,12 +1819,12 @@ "\n", "\n", - "\n", + "\n", " \n", " \n", " \n", " \n", - " 2024-09-19T18:28:58.643361\n", + " 2024-09-20T00:24:18.029570\n", " image/svg+xml\n", " \n", " \n", @@ -1961,9 +1839,9 @@ " \n", " \n", " \n", - " \n", @@ -1980,67 +1858,24 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", - " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", " \n", - " \n", - " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", - " \n", - " \n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", + " \n", " \n", " \n", " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", + "\" clip-path=\"url(#p007a01d97b)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -2286,17 +2333,17 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "\" clip-path=\"url(#p007a01d97b)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -2306,35 +2353,19 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "\" clip-path=\"url(#p007a01d97b)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", " \n", " \n", " \n", @@ -2343,17 +2374,17 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "\" clip-path=\"url(#p007a01d97b)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -2364,17 +2395,17 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "\" clip-path=\"url(#p007a01d97b)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -2385,17 +2416,17 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "\" clip-path=\"url(#p007a01d97b)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -2406,53 +2437,19 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", + "\" clip-path=\"url(#p007a01d97b)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n", " \n", - " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", " \n", " \n", " \n", @@ -2461,22 +2458,39 @@ " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -2516,7 +2811,14 @@ } ], "source": [ - "pyplot.plot(output['Bulk']['T'], output['Bulk']['height'])\n", + "for key in ('Homogeneous', 'Bulk'):\n", + " pyplot.plot(\n", + " np.asarray(output[key]['vapour']) + np.asarray(output[key]['water']) + np.asarray(output[key]['ice']),\n", + " output[key]['height'],\n", + " label=key,\n", + " linestyle=':' if key[0] == 'B' else '--'\n", + " )\n", + "pyplot.legend()\n", "pyplot.grid()" ] }, diff --git a/examples/PySDM_examples/Abade_and_Albuquerque_2024/simulation.py b/examples/PySDM_examples/Abade_and_Albuquerque_2024/simulation.py index b1a329b2c..1497d5dd9 100644 --- a/examples/PySDM_examples/Abade_and_Albuquerque_2024/simulation.py +++ b/examples/PySDM_examples/Abade_and_Albuquerque_2024/simulation.py @@ -13,6 +13,7 @@ from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity from PySDM.products import ( AmbientTemperature, + AmbientWaterVapourMixingRatio, ParcelDisplacement, WaterMixingRatio, SpecificIceWaterContent, @@ -36,14 +37,13 @@ def __init__(self, settings): ), ) builder.add_dynamic(AmbientThermodynamics()) + builder.add_dynamic(Condensation()) if settings.enable_immersion_freezing: builder.add_dynamic(Freezing()) if settings.enable_vapour_deposition_on_ice: builder.add_dynamic(VapourDepositionOnIce()) - builder.add_dynamic(Condensation()) - r_dry, n_in_dv = ConstantMultiplicity(settings.soluble_aerosol).sample( n_sd=settings.n_sd ) @@ -63,6 +63,9 @@ def __init__(self, settings): SpecificIceWaterContent(name="ice"), ParcelDisplacement(name="height"), AmbientTemperature(name="T"), + AmbientWaterVapourMixingRatio( + name="vapour", var="water_vapour_mixing_ratio" + ), ) super().__init__( particulator=builder.build(attributes=attributes, products=self.products) From da030a5c2a8f5c0fda0e034dfe31992031d5380b Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Wed, 25 Sep 2024 17:06:25 +0200 Subject: [PATCH 20/27] Removed updates to predicted water vapour mixing ration to make vapour deposition test work with box enviroment again --- PySDM/dynamics/vapour_deposition_on_ice.py | 1 - PySDM/particulator.py | 4 +--- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/PySDM/dynamics/vapour_deposition_on_ice.py b/PySDM/dynamics/vapour_deposition_on_ice.py index a34eb0820..ac38946df 100644 --- a/PySDM/dynamics/vapour_deposition_on_ice.py +++ b/PySDM/dynamics/vapour_deposition_on_ice.py @@ -18,4 +18,3 @@ def register(self, *, builder): def __call__(self): """called by the particulator during simulation""" self.particulator.deposition() - self.particulator.update_TpRH() diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 0384650ca..7cb69c512 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -488,9 +488,7 @@ def deposition(self): ambient_total_pressure=self.environment["P"], ambient_humidity=self.environment["RH"], ambient_water_activity=self.environment["a_w_ice"], - ambient_vapour_mixing_ratio=self.environment.get_predicted( - "water_vapour_mixing_ratio" - ), + ambient_vapour_mixing_ratio=self.environment["water_vapour_mixing_ratio"], ambient_dry_air_density=self.environment["rhod"], cell_volume=self.environment.mesh.dv, time_step=self.dt, From 6d3b3262a07de369bbb67c57143b76b58a818ec2 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Thu, 26 Sep 2024 17:20:43 +0200 Subject: [PATCH 21/27] Added setup for kinetic correction routines in ice deposition methods --- .../impl_numba/methods/deposition_methods.py | 9 +++---- PySDM/formulae.py | 2 ++ PySDM/physics/__init__.py | 1 + .../diffusion_ice_kinetics/__init__.py | 6 +++++ .../physics/diffusion_ice_kinetics/neglect.py | 24 +++++++++++++++++++ .../dynamics/test_vapour_deposition_on_ice.py | 7 +++++- 6 files changed, 44 insertions(+), 5 deletions(-) create mode 100644 PySDM/physics/diffusion_ice_kinetics/__init__.py create mode 100644 PySDM/physics/diffusion_ice_kinetics/neglect.py diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index 34df43a54..231e2bf87 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -14,7 +14,8 @@ def _deposition(self): formulae = self.formulae_flattened liquid = formulae.trivia__unfrozen - diffusion_coefficient_function = self.formulae.diffusion_thermics.D + # diffusion_coefficient_function = self.formulae.diffusion_thermics.D + # diffusion_kinetic_correction_function = self.formulae.diffusion_thermics.D @numba.jit(**self.default_jit_flags) def body( @@ -52,9 +53,9 @@ def body( ) ) - diffusion_coefficient = diffusion_coefficient_function( - temperature, pressure - ) + DTp = formulae.diffusion_thermics__D(temperature, pressure) + + diffusion_coefficient = formulae.diffusion_ice_kinetics__D(DTp, 0., temperature) saturation_ratio_ice = ( ambient_humidity[cid] / ambient_water_activity[cid] diff --git a/PySDM/formulae.py b/PySDM/formulae.py index 47b26bdb4..949054e89 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -38,6 +38,7 @@ def __init__( # pylint: disable=too-many-locals drop_growth: str = "Mason1971", surface_tension: str = "Constant", diffusion_kinetics: str = "FuchsSutugin", + diffusion_ice_kinetics: str = "Neglect", diffusion_thermics: str = "Neglect", ventilation: str = "Neglect", state_variable_triplet: str = "LibcloudphPlusPlus", @@ -70,6 +71,7 @@ def __init__( # pylint: disable=too-many-locals self.drop_growth = drop_growth self.surface_tension = surface_tension self.diffusion_kinetics = diffusion_kinetics + self.diffusion_ice_kinetics = diffusion_ice_kinetics self.latent_heat = latent_heat self.diffusion_thermics = diffusion_thermics self.ventilation = ventilation diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py index e2b3d9e48..9a4e126bb 100644 --- a/PySDM/physics/__init__.py +++ b/PySDM/physics/__init__.py @@ -21,6 +21,7 @@ condensation_coordinate, constants_defaults, diffusion_kinetics, + diffusion_ice_kinetics, diffusion_thermics, drop_growth, fragmentation_function, diff --git a/PySDM/physics/diffusion_ice_kinetics/__init__.py b/PySDM/physics/diffusion_ice_kinetics/__init__.py new file mode 100644 index 000000000..256b351c1 --- /dev/null +++ b/PySDM/physics/diffusion_ice_kinetics/__init__.py @@ -0,0 +1,6 @@ +""" +Formulae for handling transition-régime corrections in diffusional growh/evaporation of ice +""" + +from .neglect import Neglect + diff --git a/PySDM/physics/diffusion_ice_kinetics/neglect.py b/PySDM/physics/diffusion_ice_kinetics/neglect.py new file mode 100644 index 000000000..a357b594a --- /dev/null +++ b/PySDM/physics/diffusion_ice_kinetics/neglect.py @@ -0,0 +1,24 @@ +""" +no transition-regime corrections formulation +""" + + +class Neglect: + def __init__(self, _): + pass + + @staticmethod + def lambdaD(_, D, T): # pylint: disable=unused-argument + return -1 + + @staticmethod + def lambdaK(_, T, p): # pylint: disable=unused-argument + return -1 + + @staticmethod + def D(_, D, r, lmbd): # pylint: disable=unused-argument + return D + + @staticmethod + def K(_, K, r, lmbd): # pylint: disable=unused-argument + return K diff --git a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py index a3dac4439..330f45cdb 100644 --- a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py +++ b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py @@ -18,6 +18,11 @@ @pytest.mark.parametrize("RHi", (1.1, 1.0, 0.9)) @pytest.mark.parametrize("fastmath", (True, False)) @pytest.mark.parametrize("diffusion_coordinate", ("Mass", "MassLogarithm")) +# @pytest.mark.parametrize("dt", (1 * si.s,)) +# @pytest.mark.parametrize("water_mass", (-si.mg,)) +# @pytest.mark.parametrize("RHi", (1.1,)) +# @pytest.mark.parametrize("fastmath", (False,)) +# @pytest.mark.parametrize("diffusion_coordinate", ("Mass",)) def test_iwc_lower_after_timestep( dt, water_mass, RHi, fastmath, diffusion_coordinate, dv=1 * si.m**3 ): @@ -51,7 +56,7 @@ def test_iwc_lower_after_timestep( pvs_water = particulator.formulae.saturation_vapour_pressure.pvs_water(temperature) vapour_pressure = RHi * pvs_ice RH = vapour_pressure / pvs_water - particulator.environment["RH"] = RH # 1.1 * si.dimensionless + particulator.environment["RH"] = RH particulator.environment["a_w_ice"] = pvs_ice / pvs_water particulator.environment["Schmidt number"] = 1 rv0 = ( From 8e5f6dcbee2021e2f0a016244e3b4768b566e4f7 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Fri, 27 Sep 2024 18:10:48 +0200 Subject: [PATCH 22/27] Added kinetic correction for diffusivity and conductivity in ice vapour deposition --- .../impl_numba/methods/deposition_methods.py | 21 ++++++++++----- PySDM/formulae.py | 2 +- PySDM/physics/constants_defaults.py | 10 ++++++- .../diffusion_ice_kinetics/__init__.py | 2 +- .../diffusion_ice_kinetics/lamb_verlinde.py | 27 +++++++++++++++++++ .../physics/diffusion_ice_kinetics/neglect.py | 8 +++--- .../dynamics/test_vapour_deposition_on_ice.py | 2 +- 7 files changed, 57 insertions(+), 15 deletions(-) create mode 100644 PySDM/physics/diffusion_ice_kinetics/lamb_verlinde.py diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index 231e2bf87..b8b0bf87c 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -13,9 +13,7 @@ def _deposition(self): formulae = self.formulae_flattened liquid = formulae.trivia__unfrozen - - # diffusion_coefficient_function = self.formulae.diffusion_thermics.D - # diffusion_kinetic_correction_function = self.formulae.diffusion_thermics.D + @numba.jit(**self.default_jit_flags) def body( @@ -40,8 +38,13 @@ def body( ice_mass = -water_mass[i] cid = cell_id[i] + radius = formulae.particle_shape_and_density__ice_mass_to_radius( + water_mass[i] + ) + temperature = ambient_temperature[cid] pressure = ambient_total_pressure[cid] + rho = ambient_dry_air_density[cid] capacity = formulae.particle_shape_and_density__ice_mass_to_radius( water_mass[i] ) @@ -53,9 +56,13 @@ def body( ) ) - DTp = formulae.diffusion_thermics__D(temperature, pressure) - - diffusion_coefficient = formulae.diffusion_ice_kinetics__D(DTp, 0., temperature) + Dv_const = formulae.diffusion_thermics__D(temperature, pressure) + lambdaD = formulae.diffusion_ice_kinetics__lambdaD(temperature, pressure) + diffusion_coefficient = formulae.diffusion_ice_kinetics__D(Dv_const, radius, lambdaD, temperature) + + Ka_const = formulae.diffusion_thermics__K(temperature, pressure) + lambdaK = formulae.diffusion_ice_kinetics__lambdaK(temperature, pressure) + thermal_conductivity = formulae.diffusion_ice_kinetics__K(Ka_const, radius, lambdaK, temperature, rho) saturation_ratio_ice = ( ambient_humidity[cid] / ambient_water_activity[cid] @@ -81,7 +88,7 @@ def body( -dm_dt * multiplicity[i] * time_step - / (cell_volume * ambient_dry_air_density[cid]) + / (cell_volume * rho) ) if -delta_rv_i > ambient_vapour_mixing_ratio[cid]: assert False diff --git a/PySDM/formulae.py b/PySDM/formulae.py index 949054e89..ff6311212 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -38,7 +38,7 @@ def __init__( # pylint: disable=too-many-locals drop_growth: str = "Mason1971", surface_tension: str = "Constant", diffusion_kinetics: str = "FuchsSutugin", - diffusion_ice_kinetics: str = "Neglect", + diffusion_ice_kinetics: str = "LambVerlinde", diffusion_thermics: str = "Neglect", ventilation: str = "Neglect", state_variable_triplet: str = "LibcloudphPlusPlus", diff --git a/PySDM/physics/constants_defaults.py b/PySDM/physics/constants_defaults.py index 6a2feaf53..54ef3c6a4 100644 --- a/PySDM/physics/constants_defaults.py +++ b/PySDM/physics/constants_defaults.py @@ -62,10 +62,14 @@ K0 = 2.4e-2 * si.joules / si.metres / si.seconds / si.kelvins -# mass and heat accommodation coefficients +# mass and heat accommodation coefficients for condensation MAC = 1.0 HAC = 1.0 +# mass and heat accommodation coefficients for vapour deposition on ice +MAC_ice = 0.5 +HAC_ice = 0.7 + p1000 = 1000 * si.hectopascals c_pd = 1005 * si.joule / si.kilogram / si.kelvin c_pv = 1850 * si.joule / si.kilogram / si.kelvin @@ -135,6 +139,10 @@ # Delta v for diffusivity in Pruppacher & Klett eq. 13-14 dv_pk05 = 0.0 * si.metres +# mean free path of lambda of water molecules +# in pruppacher und klett +lmbd_w_0 = 6.6e-8 * si.metre + # Seinfeld & Pandis eq. 15.65 (Hall & Pruppacher 1976) d_l19_a = 0.211e-4 * si.metre**2 / si.second d_l19_b = 1.94 diff --git a/PySDM/physics/diffusion_ice_kinetics/__init__.py b/PySDM/physics/diffusion_ice_kinetics/__init__.py index 256b351c1..32f6bd5c4 100644 --- a/PySDM/physics/diffusion_ice_kinetics/__init__.py +++ b/PySDM/physics/diffusion_ice_kinetics/__init__.py @@ -3,4 +3,4 @@ """ from .neglect import Neglect - +from .lamb_verlinde import LambVerlinde diff --git a/PySDM/physics/diffusion_ice_kinetics/lamb_verlinde.py b/PySDM/physics/diffusion_ice_kinetics/lamb_verlinde.py new file mode 100644 index 000000000..303bc355c --- /dev/null +++ b/PySDM/physics/diffusion_ice_kinetics/lamb_verlinde.py @@ -0,0 +1,27 @@ +""" +transition-regime correction as in 'Physics and Chemistry of Clouds' by Lamb and Verlinde (2011), Chapter 8.2 +with free pathway of air/vapour (lambdaD) after Pruppacher and Klett (2010) +""" + +import numpy as np + + +class LambVerlinde: + def __init__(self, _): + pass + + @staticmethod + def lambdaD(_, T, p): # pylint: disable=unused-argument + return const.lmbd_w_0 * T / const.T_STP * const.p_STP / p + + @staticmethod + def lambdaK(_, T, p): # pylint: disable=unused-argument + return const.lmbd_w_0 * T / const.T_STP * const.p_STP / p + + @staticmethod + def D(_, D, r, lmbd, T): # pylint: disable=unused-argument + return D / ( r / (r + lmbd) + 4.0 * D / const.MAC_ice / np.sqrt(8. * const.Rv * T / const.PI) / r ) + + @staticmethod + def K(_, K, r, lmbd, T, rho): # pylint: disable=unused-argument + return K / ( r / (r + lmbd) + 4.0 * K / const.HAC_ice / np.sqrt(8. * const.Rd * T / const.PI) / const.c_pd / rho / r ) diff --git a/PySDM/physics/diffusion_ice_kinetics/neglect.py b/PySDM/physics/diffusion_ice_kinetics/neglect.py index a357b594a..5550bf036 100644 --- a/PySDM/physics/diffusion_ice_kinetics/neglect.py +++ b/PySDM/physics/diffusion_ice_kinetics/neglect.py @@ -6,9 +6,9 @@ class Neglect: def __init__(self, _): pass - + @staticmethod - def lambdaD(_, D, T): # pylint: disable=unused-argument + def lambdaD(_, T, p): # pylint: disable=unused-argument return -1 @staticmethod @@ -16,9 +16,9 @@ def lambdaK(_, T, p): # pylint: disable=unused-argument return -1 @staticmethod - def D(_, D, r, lmbd): # pylint: disable=unused-argument + def D(_, D, r, lmbd, T): # pylint: disable=unused-argument return D @staticmethod - def K(_, K, r, lmbd): # pylint: disable=unused-argument + def K(_, K, r, lmbd, T, rho): # pylint: disable=unused-argument return K diff --git a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py index 330f45cdb..90a8ea1ed 100644 --- a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py +++ b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py @@ -69,7 +69,7 @@ def test_iwc_lower_after_timestep( temperature * particulator.formulae.constants.Rd ) - print(f" {RHi=}, {RH=}, {vapour_pressure=}, {pvs_ice=}, {pvs_water=}") + # act iwc_old = particulator.products["ice water content"].get().copy() From f3a1bf8c3b89658370706fe2f34c72f209aded55 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Mon, 7 Oct 2024 13:26:45 +0200 Subject: [PATCH 23/27] Changed name of ice vapour deposition kinetic corection routine --- PySDM/formulae.py | 2 +- PySDM/physics/diffusion_ice_kinetics/__init__.py | 2 +- .../diffusion_ice_kinetics/{lamb_verlinde.py => standard.py} | 3 ++- 3 files changed, 4 insertions(+), 3 deletions(-) rename PySDM/physics/diffusion_ice_kinetics/{lamb_verlinde.py => standard.py} (94%) diff --git a/PySDM/formulae.py b/PySDM/formulae.py index ff6311212..61648f2c0 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -38,7 +38,7 @@ def __init__( # pylint: disable=too-many-locals drop_growth: str = "Mason1971", surface_tension: str = "Constant", diffusion_kinetics: str = "FuchsSutugin", - diffusion_ice_kinetics: str = "LambVerlinde", + diffusion_ice_kinetics: str = "Standard", diffusion_thermics: str = "Neglect", ventilation: str = "Neglect", state_variable_triplet: str = "LibcloudphPlusPlus", diff --git a/PySDM/physics/diffusion_ice_kinetics/__init__.py b/PySDM/physics/diffusion_ice_kinetics/__init__.py index 32f6bd5c4..b9aea5da1 100644 --- a/PySDM/physics/diffusion_ice_kinetics/__init__.py +++ b/PySDM/physics/diffusion_ice_kinetics/__init__.py @@ -3,4 +3,4 @@ """ from .neglect import Neglect -from .lamb_verlinde import LambVerlinde +from .standard import Standard diff --git a/PySDM/physics/diffusion_ice_kinetics/lamb_verlinde.py b/PySDM/physics/diffusion_ice_kinetics/standard.py similarity index 94% rename from PySDM/physics/diffusion_ice_kinetics/lamb_verlinde.py rename to PySDM/physics/diffusion_ice_kinetics/standard.py index 303bc355c..9dd7cb087 100644 --- a/PySDM/physics/diffusion_ice_kinetics/lamb_verlinde.py +++ b/PySDM/physics/diffusion_ice_kinetics/standard.py @@ -1,12 +1,13 @@ """ transition-regime correction as in 'Physics and Chemistry of Clouds' by Lamb and Verlinde (2011), Chapter 8.2 +or 13.1 in Pruppbacher and Klett (2010) with free pathway of air/vapour (lambdaD) after Pruppacher and Klett (2010) """ import numpy as np -class LambVerlinde: +class Standard: def __init__(self, _): pass From c229529f151a3041d405fbca5a7802154868c513 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Mon, 7 Oct 2024 16:01:19 +0200 Subject: [PATCH 24/27] Added capacity of ice crystals during deposition of vapour as its own physics routine --- .../impl_numba/methods/deposition_methods.py | 10 +++++----- PySDM/formulae.py | 2 ++ PySDM/physics/__init__.py | 1 + PySDM/physics/diffusion_ice_capacity/__init__.py | 5 +++++ PySDM/physics/diffusion_ice_capacity/spherical.py | 13 +++++++++++++ 5 files changed, 26 insertions(+), 5 deletions(-) create mode 100644 PySDM/physics/diffusion_ice_capacity/__init__.py create mode 100644 PySDM/physics/diffusion_ice_capacity/spherical.py diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index b8b0bf87c..131ac74da 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -32,7 +32,7 @@ def body( schmidt_number, ): n_sd = len(water_mass) - # for i in numba.prange(n_sd): # pylint: disable=not-an-iterable + ## for i in numba.prange(n_sd): # pylint: disable=not-an-iterable for i in range(n_sd): if not liquid(water_mass[i]): ice_mass = -water_mass[i] @@ -41,14 +41,14 @@ def body( radius = formulae.particle_shape_and_density__ice_mass_to_radius( water_mass[i] ) + diameter = radius * 2. temperature = ambient_temperature[cid] pressure = ambient_total_pressure[cid] rho = ambient_dry_air_density[cid] - capacity = formulae.particle_shape_and_density__ice_mass_to_radius( - water_mass[i] - ) - + + capacity = formulae.diffusion_ice_capacity__capacity( diameter ) + ventilation_factor = formulae.ventilation__ventilation_coefficient( sqrt_re_times_cbrt_sc=formulae.trivia__sqrt_re_times_cbrt_sc( Re=reynolds_number[i], diff --git a/PySDM/formulae.py b/PySDM/formulae.py index 61648f2c0..b03e43653 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -39,6 +39,7 @@ def __init__( # pylint: disable=too-many-locals surface_tension: str = "Constant", diffusion_kinetics: str = "FuchsSutugin", diffusion_ice_kinetics: str = "Standard", + diffusion_ice_capacity: str = "Spherical", diffusion_thermics: str = "Neglect", ventilation: str = "Neglect", state_variable_triplet: str = "LibcloudphPlusPlus", @@ -72,6 +73,7 @@ def __init__( # pylint: disable=too-many-locals self.surface_tension = surface_tension self.diffusion_kinetics = diffusion_kinetics self.diffusion_ice_kinetics = diffusion_ice_kinetics + self.diffusion_ice_capacity = diffusion_ice_capacity self.latent_heat = latent_heat self.diffusion_thermics = diffusion_thermics self.ventilation = ventilation diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py index 9a4e126bb..c172b3fad 100644 --- a/PySDM/physics/__init__.py +++ b/PySDM/physics/__init__.py @@ -22,6 +22,7 @@ constants_defaults, diffusion_kinetics, diffusion_ice_kinetics, + diffusion_ice_capacity, diffusion_thermics, drop_growth, fragmentation_function, diff --git a/PySDM/physics/diffusion_ice_capacity/__init__.py b/PySDM/physics/diffusion_ice_capacity/__init__.py new file mode 100644 index 000000000..d0082d718 --- /dev/null +++ b/PySDM/physics/diffusion_ice_capacity/__init__.py @@ -0,0 +1,5 @@ +""" +Formulae for capacity in diffusional growh/evaporation of ice +""" + +from .spherical import Spherical diff --git a/PySDM/physics/diffusion_ice_capacity/spherical.py b/PySDM/physics/diffusion_ice_capacity/spherical.py new file mode 100644 index 000000000..b8553f3c8 --- /dev/null +++ b/PySDM/physics/diffusion_ice_capacity/spherical.py @@ -0,0 +1,13 @@ +""" +capacity for approximation of ice crystals as spheres +""" + + +class Spherical: + + def __init__(self, _): + pass + + @staticmethod + def capacity(_, diameter, length=None ): # pylint: disable=unused-argument + return diameter / 2. From 75a61e53cc1fcc492a9e6cf08f68c26cbd0847db Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Mon, 7 Oct 2024 17:55:08 +0200 Subject: [PATCH 25/27] Added latent heat of sublimation routines --- .../impl_numba/methods/deposition_methods.py | 3 +++ PySDM/formulae.py | 2 ++ PySDM/physics/__init__.py | 1 + PySDM/physics/constants_defaults.py | 8 ++++++++ PySDM/physics/latent_heat_sublimation/__init__.py | 6 ++++++ PySDM/physics/latent_heat_sublimation/constant.py | 12 ++++++++++++ .../latent_heat_sublimation/murphy_koop_2005.py | 12 ++++++++++++ 7 files changed, 44 insertions(+) create mode 100644 PySDM/physics/latent_heat_sublimation/__init__.py create mode 100644 PySDM/physics/latent_heat_sublimation/constant.py create mode 100644 PySDM/physics/latent_heat_sublimation/murphy_koop_2005.py diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index 131ac74da..14047bb37 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -64,6 +64,9 @@ def body( lambdaK = formulae.diffusion_ice_kinetics__lambdaK(temperature, pressure) thermal_conductivity = formulae.diffusion_ice_kinetics__K(Ka_const, radius, lambdaK, temperature, rho) + latent_heat_sub = formulae.latent_heat_sublimation__ls(temperature) + + saturation_ratio_ice = ( ambient_humidity[cid] / ambient_water_activity[cid] ) diff --git a/PySDM/formulae.py b/PySDM/formulae.py index b03e43653..586553b67 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -34,6 +34,7 @@ def __init__( # pylint: disable=too-many-locals condensation_coordinate: str = "VolumeLogarithm", saturation_vapour_pressure: str = "FlatauWalkoCotton", latent_heat: str = "Kirchhoff", + latent_heat_sublimation: str = "MurphyKoop", hygroscopicity: str = "KappaKoehlerLeadingTerms", drop_growth: str = "Mason1971", surface_tension: str = "Constant", @@ -75,6 +76,7 @@ def __init__( # pylint: disable=too-many-locals self.diffusion_ice_kinetics = diffusion_ice_kinetics self.diffusion_ice_capacity = diffusion_ice_capacity self.latent_heat = latent_heat + self.latent_heat_sublimation = latent_heat_sublimation self.diffusion_thermics = diffusion_thermics self.ventilation = ventilation self.state_variable_triplet = state_variable_triplet diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py index c172b3fad..6a7c7f9ee 100644 --- a/PySDM/physics/__init__.py +++ b/PySDM/physics/__init__.py @@ -37,6 +37,7 @@ isotope_diffusivity_ratios, isotope_relaxation_timescale, latent_heat, + latent_heat_sublimation, optical_albedo, optical_depth, particle_advection, diff --git a/PySDM/physics/constants_defaults.py b/PySDM/physics/constants_defaults.py index 54ef3c6a4..34bb957e2 100644 --- a/PySDM/physics/constants_defaults.py +++ b/PySDM/physics/constants_defaults.py @@ -131,6 +131,7 @@ l_l19_a = 0.167 * si.dimensionless l_l19_b = 3.65e-4 / si.kelvin + # Thermal diffusivity constants from Lowe et al. (2019) k_l19_a = 4.2e-3 * si.joules / si.metres / si.seconds / si.kelvins k_l19_b = 1.0456 * si.dimensionless @@ -168,6 +169,13 @@ MK05_LIQ_C12 = 1 * si.K MK05_LIQ_C13 = 0.014025 / si.K +# Latent heat of sublimation from Murphy and Koop (2005) +MK05_SUB_C1 = 46782.5 * si.joule / si.mole +MK05_SUB_C2 = 35.8925 * si.joule / si.mole / si.kelvin +MK05_SUB_C3 = 0.07414 * si.joule / si.mole / si.kelvin**2 +MK05_SUB_C4 = 541.5 * si.joule / si.mole +MK05_SUB_C5 = 123.75 * si.kelvin + # standard pressure and temperature (ICAO) T_STP = (sci.zero_Celsius + 15) * si.kelvin p_STP = 101325 * si.pascal diff --git a/PySDM/physics/latent_heat_sublimation/__init__.py b/PySDM/physics/latent_heat_sublimation/__init__.py new file mode 100644 index 000000000..e42dda007 --- /dev/null +++ b/PySDM/physics/latent_heat_sublimation/__init__.py @@ -0,0 +1,6 @@ +""" +formulations of latent heat of sublimation temperature dependence (or lack thereof) +""" + +from .constant import Constant +from .murphy_koop_2005 import MurphyKoop diff --git a/PySDM/physics/latent_heat_sublimation/constant.py b/PySDM/physics/latent_heat_sublimation/constant.py new file mode 100644 index 000000000..cd1fae92b --- /dev/null +++ b/PySDM/physics/latent_heat_sublimation/constant.py @@ -0,0 +1,12 @@ +""" +temperature-independent latent heat of sublimation +""" + + +class Constant: # pylint: disable=too-few-public-methods + def __init__(self, _): + pass + + @staticmethod + def ls(const, T): # pylint: disable=unused-argument + return const.l_tri diff --git a/PySDM/physics/latent_heat_sublimation/murphy_koop_2005.py b/PySDM/physics/latent_heat_sublimation/murphy_koop_2005.py new file mode 100644 index 000000000..d50204525 --- /dev/null +++ b/PySDM/physics/latent_heat_sublimation/murphy_koop_2005.py @@ -0,0 +1,12 @@ +""" +temperature-dependent latent heat of sublimation from Murphy and Koop (2005) Eq. (5) +""" + + +class MurphyKoop: # pylint: disable=too-few-public-methods + def __init__(self, _): + pass + + @staticmethod + def ls(const, T): # pylint: disable=unused-argument + return (const.MK05_SUB_C1 + const.MK05_SUB_C2 * T - const.MK05_SUB_C3 * T**2. + const.MK05_SUB_C4 * np.exp(-(T/const.MK05_SUB_C5)**2.)) / const.Mv From 3a8cabcd1d505c179e08bfef82aa1a1d8ec35129 Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Tue, 8 Oct 2024 15:43:20 +0200 Subject: [PATCH 26/27] Added ice crystal temperature correction factor for vapour deposition --- .../impl_numba/methods/deposition_methods.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index 14047bb37..d14628e35 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -46,6 +46,9 @@ def body( temperature = ambient_temperature[cid] pressure = ambient_total_pressure[cid] rho = ambient_dry_air_density[cid] + Rv = formulae.constants.Rv + pvs_ice = formulae.saturation_vapour_pressure__pvs_ice(temperature) + latent_heat_sub = formulae.latent_heat_sublimation__ls(temperature) capacity = formulae.diffusion_ice_capacity__capacity( diameter ) @@ -64,18 +67,20 @@ def body( lambdaK = formulae.diffusion_ice_kinetics__lambdaK(temperature, pressure) thermal_conductivity = formulae.diffusion_ice_kinetics__K(Ka_const, radius, lambdaK, temperature, rho) - latent_heat_sub = formulae.latent_heat_sublimation__ls(temperature) - + + howell_factor = 1. / ((latent_heat_sub / Rv / temperature - 1.) * latent_heat_sub * diffusion_coefficient / temperature / thermal_conductivity + Rv * temperature / pvs_ice) + # print( f" {howell_factor=}, {latent_heat_sub= }, {Rv=}, {Dv_const=}, {diffusion_coefficient=}, {Ka_const=}, {thermal_conductivity=}, {pvs_ice=}") saturation_ratio_ice = ( ambient_humidity[cid] / ambient_water_activity[cid] ) rho_vs_ice = ( - formulae.saturation_vapour_pressure__pvs_ice(temperature) - / formulae.constants.Rv + pvs_ice + / Rv / temperature ) + dm_dt = ( 4 * np.pi From 95cf76e1dc819bbf70a965a3a0bfce80fb3adabf Mon Sep 17 00:00:00 2001 From: Tim Luettmer Date: Wed, 9 Oct 2024 13:19:32 +0200 Subject: [PATCH 27/27] Added latent heating for ice vapour deposition and included it in unit test --- PySDM/backends/impl_numba/methods/deposition_methods.py | 3 +++ tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py | 7 ++++++- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/PySDM/backends/impl_numba/methods/deposition_methods.py b/PySDM/backends/impl_numba/methods/deposition_methods.py index d14628e35..0ac5e4a6e 100644 --- a/PySDM/backends/impl_numba/methods/deposition_methods.py +++ b/PySDM/backends/impl_numba/methods/deposition_methods.py @@ -102,6 +102,9 @@ def body( assert False ambient_vapour_mixing_ratio[cid] += delta_rv_i + delta_T = -delta_rv_i * latent_heat_sub / formulae.constants.c_pd + ambient_temperature[cid] += delta_T + x_old = formulae.diffusion_coordinate__x(ice_mass) dx_dt_old = formulae.diffusion_coordinate__dx_dt(x_old, dm_dt) x_new = formulae.trivia__explicit_euler(x_old, time_step, dx_dt_old) diff --git a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py index 90a8ea1ed..bac8cf971 100644 --- a/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py +++ b/tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py @@ -72,19 +72,24 @@ def test_iwc_lower_after_timestep( # act + T0 = temperature iwc_old = particulator.products["ice water content"].get().copy() particulator.run(steps=1) + T_new = particulator.environment["T"][0] iwc_new = particulator.products["ice water content"].get().copy() - + rv_new = particulator.environment["water_vapour_mixing_ratio"][0] # assert if water_mass < 0 and RHi != 1: if RHi > 1: assert (iwc_new > iwc_old).all() assert rv_new < rv0 + assert T_new > T0 elif RHi < 1: assert (iwc_new < iwc_old).all() assert rv_new > rv0 + assert T_new < T0 else: assert (iwc_new == iwc_old).all() assert rv_new == rv0 + assert T_new == T0