From a07b9575055666b5e5fc1299b7bf4e27e254b4a8 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sun, 19 Dec 2021 11:40:08 -0600 Subject: [PATCH 01/33] adding test_freezing for ICMW (with both singular & time-dependent runs and assert for non-zero ice-water content) --- PySDM/initialisation/__init__.py | 2 +- .../arabas_et_al_2015/test_freezing.py | 48 +++++++++++++++++++ 2 files changed, 49 insertions(+), 1 deletion(-) create mode 100644 tests/smoke_tests/arabas_et_al_2015/test_freezing.py diff --git a/PySDM/initialisation/__init__.py b/PySDM/initialisation/__init__.py index ec15d47e5..2b2ad9abc 100644 --- a/PySDM/initialisation/__init__.py +++ b/PySDM/initialisation/__init__.py @@ -1,5 +1,5 @@ """ -Initialisation logic particle size spectra, sampling methods and +initialisation logic, particle size spectra, sampling methods and wet radii equilibration """ from .discretise_multiplicities import discretise_multiplicities diff --git a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py new file mode 100644 index 000000000..330d1215a --- /dev/null +++ b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py @@ -0,0 +1,48 @@ +# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring +import numpy as np +from matplotlib import pyplot +import pytest +from PySDM_examples.Szumowski_et_al_1998 import Simulation +from PySDM_examples.Arabas_et_al_2015 import Settings, SpinUp +from PySDM.physics import si +from PySDM.backends import CPU +from .dummy_storage import DummyStorage + + + +@pytest.mark.parametrize("singular", ( + pytest.param(False, id="singular: False"), + pytest.param(True, id="singular: True") +)) +# pylint: disable=redefined-outer-name +def test_freezing(singular, plot=False): + # Arrange + settings = Settings() + settings.dt = .5 * si.second + settings.grid = (3, 25) + + settings.simulation_time = 100 * settings.dt + settings.spin_up_time = 10 * settings.dt + + settings.output_interval = settings.dt # settings.simulation_time + + settings.processes['freezing'] = True + settings.processes['coalescence'] = False + + settings.freezing_singular = singular + settings.th_std0 -= 35 * si.K + settings.qv0 -= 7.15 * si.g/si.kg + + storage = DummyStorage() + simulation = Simulation(settings, storage, SpinUp=SpinUp, backend_class=CPU) + simulation.reinit() + + # Act + simulation.run() + + # Plot + if plot: + pass + + # Assert + assert (simulation.products['ice water content'].get() > 0).any() From 7726138fa470030590b32d0060e3509320e9d6b7 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sun, 19 Dec 2021 11:48:53 -0600 Subject: [PATCH 02/33] adding dummy_storage (refactored to make it common across tests) --- .../smoke_tests/arabas_et_al_2015/dummy_storage.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) create mode 100644 tests/smoke_tests/arabas_et_al_2015/dummy_storage.py diff --git a/tests/smoke_tests/arabas_et_al_2015/dummy_storage.py b/tests/smoke_tests/arabas_et_al_2015/dummy_storage.py new file mode 100644 index 000000000..2c0ccfee1 --- /dev/null +++ b/tests/smoke_tests/arabas_et_al_2015/dummy_storage.py @@ -0,0 +1,13 @@ +import numpy as np + + +class DummyStorage: + def __init__(self): + self.profiles = [] + + def init(*_): # pylint: disable=no-method-argument + pass + + def save(self, data: np.ndarray, step: int, name: str): # pylint: disable=unused-argument + if name == "qv_env": + self.profiles.append({"qv_env": np.mean(data, axis=0)}) From 82af378b88d9b02f6e69cd41abeb2cb3d8552e09 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sun, 19 Dec 2021 12:13:42 -0600 Subject: [PATCH 03/33] removing unused code, marking test for non-docstring treatment under pylint --- tests/smoke_tests/arabas_et_al_2015/dummy_storage.py | 1 + tests/smoke_tests/arabas_et_al_2015/test_freezing.py | 8 +------- 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/tests/smoke_tests/arabas_et_al_2015/dummy_storage.py b/tests/smoke_tests/arabas_et_al_2015/dummy_storage.py index 2c0ccfee1..df9665ffc 100644 --- a/tests/smoke_tests/arabas_et_al_2015/dummy_storage.py +++ b/tests/smoke_tests/arabas_et_al_2015/dummy_storage.py @@ -1,3 +1,4 @@ +# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring import numpy as np diff --git a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py index 330d1215a..e1448bcbd 100644 --- a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py +++ b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py @@ -1,6 +1,4 @@ # pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring -import numpy as np -from matplotlib import pyplot import pytest from PySDM_examples.Szumowski_et_al_1998 import Simulation from PySDM_examples.Arabas_et_al_2015 import Settings, SpinUp @@ -15,7 +13,7 @@ pytest.param(True, id="singular: True") )) # pylint: disable=redefined-outer-name -def test_freezing(singular, plot=False): +def test_freezing(singular): # Arrange settings = Settings() settings.dt = .5 * si.second @@ -40,9 +38,5 @@ def test_freezing(singular, plot=False): # Act simulation.run() - # Plot - if plot: - pass - # Assert assert (simulation.products['ice water content'].get() > 0).any() From d38a29fba3c57c092bb9dcb9df4111e85d3e296c Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Mon, 20 Dec 2021 13:50:33 -0600 Subject: [PATCH 04/33] bump examples req --- test-time-requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test-time-requirements.txt b/test-time-requirements.txt index 8af10bf0b..77493fde4 100644 --- a/test-time-requirements.txt +++ b/test-time-requirements.txt @@ -4,4 +4,4 @@ ghapi pytest # note: if cloning both PySDM and PySDM examples, consider "pip install -e" -PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@8267691#egg=PySDM-examples +PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@0806ab0#egg=PySDM-examples From bc238893c04d7a4fd816a3a168cb80278ebb968e Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 22 Dec 2021 11:49:56 -0600 Subject: [PATCH 05/33] bump examples req --- test-time-requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test-time-requirements.txt b/test-time-requirements.txt index 77493fde4..d1f04f86c 100644 --- a/test-time-requirements.txt +++ b/test-time-requirements.txt @@ -4,4 +4,4 @@ ghapi pytest # note: if cloning both PySDM and PySDM examples, consider "pip install -e" -PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@0806ab0#egg=PySDM-examples +PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@14ce005#egg=PySDM-examples From d874a69e6d2a1ffd3a580008eadd57d5f17d0d6c Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 23 Dec 2021 00:33:24 -0600 Subject: [PATCH 06/33] mixed-phase logic in moist environments --- PySDM/environments/impl/moist.py | 4 +++- PySDM/environments/kinematic_2d.py | 22 ++++++++-------------- PySDM/environments/parcel.py | 2 +- 3 files changed, 12 insertions(+), 16 deletions(-) diff --git a/PySDM/environments/impl/moist.py b/PySDM/environments/impl/moist.py index 00f0df0a4..58ea3a0dc 100644 --- a/PySDM/environments/impl/moist.py +++ b/PySDM/environments/impl/moist.py @@ -7,8 +7,10 @@ class Moist: - def __init__(self, dt, mesh, variables): + def __init__(self, dt, mesh, variables, mixed_phase=False): variables += ['qv', 'thd', 'T', 'p', 'RH'] + if mixed_phase: + variables += ['a_w_ice'] self.particulator = None self.dt = dt self.mesh = mesh diff --git a/PySDM/environments/kinematic_2d.py b/PySDM/environments/kinematic_2d.py index 837a679b6..2af7a67ac 100644 --- a/PySDM/environments/kinematic_2d.py +++ b/PySDM/environments/kinematic_2d.py @@ -7,13 +7,14 @@ from PySDM.impl.mesh import Mesh from PySDM.initialisation.equilibrate_wet_radii import equilibrate_wet_radii, default_rtol from PySDM.initialisation.discretise_multiplicities import discretise_multiplicities +from PySDM.initialisation.sampling import spectral_sampling from PySDM.environments.impl.moist import Moist from ..impl import arakawa_c class Kinematic2D(Moist): - def __init__(self, dt, grid, size, rhod_of): - super().__init__(dt, Mesh(grid, size), []) + def __init__(self, dt, grid, size, rhod_of, mixed_phase=False): + super().__init__(dt, Mesh(grid, size), [], mixed_phase=mixed_phase) self.rhod_of = rhod_of self.formulae = None @@ -32,28 +33,21 @@ def dv(self): def init_attributes(self, *, spatial_discretisation, kappa, - spectral_discretisation = None, - spectro_glacial_discretisation = None, + dry_radius_spectrum, rtol=default_rtol ): super().sync() self.notify() - assert spectro_glacial_discretisation is None or spectral_discretisation is None - attributes = {} with np.errstate(all='raise'): positions = spatial_discretisation.sample(self.mesh.grid, self.particulator.n_sd) attributes['cell id'], attributes['cell origin'], attributes['position in cell'] = \ self.mesh.cellular_attributes(positions) - if spectral_discretisation: - r_dry, n_per_kg = spectral_discretisation.sample(self.particulator.n_sd) - elif spectro_glacial_discretisation: - r_dry, T_fz, n_per_kg = spectro_glacial_discretisation.sample( - self.particulator.n_sd) - attributes['freezing temperature'] = T_fz - else: - raise NotImplementedError() + + r_dry, n_per_kg = spectral_sampling.ConstantMultiplicity( + spectrum=dry_radius_spectrum + ).sample(n_sd=self.particulator.n_sd) attributes['dry volume'] = self.formulae.trivia.volume(radius=r_dry) attributes['kappa times dry volume'] = kappa * attributes['dry volume'] diff --git a/PySDM/environments/parcel.py b/PySDM/environments/parcel.py index d31e8ec87..5726c6560 100644 --- a/PySDM/environments/parcel.py +++ b/PySDM/environments/parcel.py @@ -23,7 +23,7 @@ def __init__( super().__init__( dt, Mesh.mesh_0d(), - ['rhod', 'z', 't'] + (['a_w_ice'] if mixed_phase else []) + mixed_phase ) self.p0 = p0 From 5d89dc6e85ded3ffd062d636c53916b09062dec1 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 23 Dec 2021 00:49:04 -0600 Subject: [PATCH 07/33] removing MG08 leftovers --- .../morrison_and_grabowski_2008/__init__.py | 0 .../test_just_do_it.py | 37 ------------------- 2 files changed, 37 deletions(-) delete mode 100644 tests/smoke_tests/morrison_and_grabowski_2008/__init__.py delete mode 100644 tests/smoke_tests/morrison_and_grabowski_2008/test_just_do_it.py diff --git a/tests/smoke_tests/morrison_and_grabowski_2008/__init__.py b/tests/smoke_tests/morrison_and_grabowski_2008/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/tests/smoke_tests/morrison_and_grabowski_2008/test_just_do_it.py b/tests/smoke_tests/morrison_and_grabowski_2008/test_just_do_it.py deleted file mode 100644 index 7b2827097..000000000 --- a/tests/smoke_tests/morrison_and_grabowski_2008/test_just_do_it.py +++ /dev/null @@ -1,37 +0,0 @@ -# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring -from contextlib import AbstractContextManager -import numpy as np -from PySDM_examples.Szumowski_et_al_1998 import Simulation, Storage -from PySDM_examples.Morrison_and_Grabowski_2008 import ColdCumulus - - -class FlowFieldAsserts(AbstractContextManager): - def __init__(self, simulation): - self.particulator = simulation.particulator - self.panic = None - - def set_percent(self, percent): - if percent == 0: - return - advector = None - solvers = self.particulator.dynamics['EulerianAdvection'].solvers.mpdatas.values() - for solver in solvers: - assert advector is None or advector is solver.advector - advector = solver.advector - np.testing.assert_allclose(advector.get_component(1)[:, 0], 0) - np.testing.assert_allclose(advector.get_component(1)[:, -1], 0, atol=1e-15) - - def __exit__(self, exc_type, exc_val, exc_tb): - pass - - -def test_just_do_it(): - settings = ColdCumulus() - settings.kappa = 0 - for process in settings.processes: - settings.processes[process] = False - settings.processes['particle advection'] += 1 - settings.processes['fluid advection'] += 1 - simulation = Simulation(settings, Storage(), None) - simulation.reinit() - simulation.run(controller=FlowFieldAsserts(simulation)) From fa149598f48c24df559191ec315d80eba42b729c Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 23 Dec 2021 00:55:19 -0600 Subject: [PATCH 08/33] missing parcel vars --- PySDM/environments/parcel.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/PySDM/environments/parcel.py b/PySDM/environments/parcel.py index 5726c6560..a1cb70772 100644 --- a/PySDM/environments/parcel.py +++ b/PySDM/environments/parcel.py @@ -23,7 +23,8 @@ def __init__( super().__init__( dt, Mesh.mesh_0d(), - mixed_phase + ['rhod', 'z', 't'], + mixed_phase=mixed_phase ) self.p0 = p0 From 55727569bbb10e9fb0c98e77c6b1a677a153710a Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 23 Dec 2021 02:13:38 -0600 Subject: [PATCH 09/33] bump examples req --- test-time-requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test-time-requirements.txt b/test-time-requirements.txt index d1f04f86c..eb634731c 100644 --- a/test-time-requirements.txt +++ b/test-time-requirements.txt @@ -4,4 +4,4 @@ ghapi pytest # note: if cloning both PySDM and PySDM examples, consider "pip install -e" -PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@14ce005#egg=PySDM-examples +PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@17d9855#egg=PySDM-examples From 7a22e6fe7598c7a72d7bdefc3bdb58c2cc82fd88 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 23 Dec 2021 04:43:52 -0600 Subject: [PATCH 10/33] bump examples req --- test-time-requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test-time-requirements.txt b/test-time-requirements.txt index eb634731c..96022b757 100644 --- a/test-time-requirements.txt +++ b/test-time-requirements.txt @@ -4,4 +4,4 @@ ghapi pytest # note: if cloning both PySDM and PySDM examples, consider "pip install -e" -PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@17d9855#egg=PySDM-examples +PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@4900e90#egg=PySDM-examples From 48f837abd0967fcecc34b3440080792cf8cff4eb Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 23 Dec 2021 15:52:43 +0100 Subject: [PATCH 11/33] adapt test_freezing to new ICMW settings init + bump examples req --- test-time-requirements.txt | 2 +- .../arabas_et_al_2015/test_freezing.py | 15 ++++++++++++++- 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/test-time-requirements.txt b/test-time-requirements.txt index 96022b757..2b0ad5849 100644 --- a/test-time-requirements.txt +++ b/test-time-requirements.txt @@ -4,4 +4,4 @@ ghapi pytest # note: if cloning both PySDM and PySDM examples, consider "pip install -e" -PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@4900e90#egg=PySDM-examples +PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@02b5216#egg=PySDM-examples diff --git a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py index e1448bcbd..ef6e28d7d 100644 --- a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py +++ b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py @@ -2,10 +2,18 @@ import pytest from PySDM_examples.Szumowski_et_al_1998 import Simulation from PySDM_examples.Arabas_et_al_2015 import Settings, SpinUp +from PySDM import Formulae from PySDM.physics import si from PySDM.backends import CPU from .dummy_storage import DummyStorage +from PySDM.physics.freezing_temperature_spectrum import niemand_et_al_2012 +from PySDM.physics.heterogeneous_ice_nucleation_rate import abifm +# TODO #599 +niemand_et_al_2012.a = -0.517 +niemand_et_al_2012.b = 8.934 +abifm.m = 28.13797 +abifm.c = -2.92414 @pytest.mark.parametrize("singular", ( @@ -15,7 +23,12 @@ # pylint: disable=redefined-outer-name def test_freezing(singular): # Arrange - settings = Settings() + settings = Settings(Formulae( + condensation_coordinate='VolumeLogarithm', + fastmath=True, + freezing_temperature_spectrum='Niemand_et_al_2012', + heterogeneous_ice_nucleation_rate='ABIFM' + )) settings.dt = .5 * si.second settings.grid = (3, 25) From 1fa98b69be3476a5f059ed644e080be7cf8d61d2 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 24 Dec 2021 01:13:31 +0100 Subject: [PATCH 12/33] bump examples req; fixing imports in test_freezing --- test-time-requirements.txt | 2 +- tests/smoke_tests/arabas_et_al_2015/test_freezing.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test-time-requirements.txt b/test-time-requirements.txt index 2b0ad5849..cf49cf5c5 100644 --- a/test-time-requirements.txt +++ b/test-time-requirements.txt @@ -4,4 +4,4 @@ ghapi pytest # note: if cloning both PySDM and PySDM examples, consider "pip install -e" -PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@02b5216#egg=PySDM-examples +PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@ed59c6e#egg=PySDM-examples diff --git a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py index ef6e28d7d..123e495b5 100644 --- a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py +++ b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py @@ -5,9 +5,9 @@ from PySDM import Formulae from PySDM.physics import si from PySDM.backends import CPU -from .dummy_storage import DummyStorage from PySDM.physics.freezing_temperature_spectrum import niemand_et_al_2012 from PySDM.physics.heterogeneous_ice_nucleation_rate import abifm +from .dummy_storage import DummyStorage # TODO #599 niemand_et_al_2012.a = -0.517 From bf5aca697b9a5c61e4fadfa63f21d53164e4e3fd Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 24 Dec 2021 01:22:36 +0100 Subject: [PATCH 13/33] sync test_spin_up --- .../arabas_et_al_2015/test_spin_up.py | 17 +++-------------- 1 file changed, 3 insertions(+), 14 deletions(-) diff --git a/tests/smoke_tests/arabas_et_al_2015/test_spin_up.py b/tests/smoke_tests/arabas_et_al_2015/test_spin_up.py index fec14d217..10c905991 100644 --- a/tests/smoke_tests/arabas_et_al_2015/test_spin_up.py +++ b/tests/smoke_tests/arabas_et_al_2015/test_spin_up.py @@ -5,23 +5,12 @@ from PySDM_examples.Szumowski_et_al_1998 import Simulation from PySDM_examples.Arabas_et_al_2015 import Settings, SpinUp from PySDM.physics import si - +from PySDM.formulae import Formulae +from .dummy_storage import DummyStorage from ...backends_fixture import backend_class assert hasattr(backend_class, '_pytestfixturefunction') -class DummyStorage: - def __init__(self): - self.profiles = [] - - def init(*_): # pylint: disable=no-method-argument - pass - - def save(self, data: np.ndarray, step: int, name: str): # pylint: disable=unused-argument - if name == "qv_env": - self.profiles.append({"qv_env": np.mean(data, axis=0)}) - - @pytest.mark.parametrize("fastmath", ( pytest.param(False, id="fastmath: False"), pytest.param(True, id="fastmath: True") @@ -29,7 +18,7 @@ def save(self, data: np.ndarray, step: int, name: str): # pylint: disable=unuse # pylint: disable=redefined-outer-name def test_spin_up(backend_class, fastmath, plot=False): # Arrange - settings = Settings(fastmath=fastmath) + settings = Settings(Formulae(fastmath=fastmath)) settings.dt = .5 * si.second settings.grid = (3, 25) settings.simulation_time = 20 * settings.dt From e6bf65289fe1794d15f7d3156944497b5a4b6cfe Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 24 Dec 2021 09:59:12 +0100 Subject: [PATCH 14/33] move Rogers & Yau terminal velocity to a separate file --- .../terminal_velocity/gunn_and_kinzer.py | 26 ++----------------- .../terminal_velocity/rogers_and_yau.py | 24 +++++++++++++++++ .../attributes/test_terminal_velocity.py | 3 ++- 3 files changed, 28 insertions(+), 25 deletions(-) create mode 100644 PySDM/physics/terminal_velocity/rogers_and_yau.py diff --git a/PySDM/physics/terminal_velocity/gunn_and_kinzer.py b/PySDM/physics/terminal_velocity/gunn_and_kinzer.py index 09a2d17e8..21617067b 100644 --- a/PySDM/physics/terminal_velocity/gunn_and_kinzer.py +++ b/PySDM/physics/terminal_velocity/gunn_and_kinzer.py @@ -1,6 +1,7 @@ """ [Gunn & Kinzer 1949](https://doi.org/10.1175/1520-0469(1949)006%3C0243:TTVOFF%3E2.0.CO;2) - terminal velocities + terminal velocities used for both coalescence kernel evaluation as well as for particle + displacement """ import numba import numpy as np @@ -43,29 +44,6 @@ def __call__(self, output, radius): self.particulator.backend.interpolation(output, radius, self.factor, self.a, self.b) -class RogersYau: - """ - Rogers & Yau, equations: (8.5), (8.6), (8.8) - """ - def __init__(self, particles, - small_k=None, medium_k=None, large_k=None, - small_r_limit=None, medium_r_limit=None): - si = const.si - self.particles = particles - self.small_k = small_k or 1.19e6 / si.cm / si.s - self.medium_k = medium_k or 8e3 / si.s - self.large_k = large_k or 2.01e3 * si.cm ** (1 / 2) / si.s - self.small_r_limit = small_r_limit or 35 * si.um - self.medium_r_limit = medium_r_limit or 600 * si.um - - def __call__(self, output, radius): - self.particles.backend.terminal_velocity( - output.data, radius.data, - self.small_k, self.medium_k, self.large_k, - self.small_r_limit, self.medium_r_limit - ) - - # TODO #348 implement in backend logic class TpDependent: def __init__(self, _, small_r_limit): diff --git a/PySDM/physics/terminal_velocity/rogers_and_yau.py b/PySDM/physics/terminal_velocity/rogers_and_yau.py new file mode 100644 index 000000000..90909fb74 --- /dev/null +++ b/PySDM/physics/terminal_velocity/rogers_and_yau.py @@ -0,0 +1,24 @@ +""" +Rogers & Yau, equations: (8.5), (8.6), (8.8) +""" +from PySDM.physics import constants as const + + +class RogersYau: + def __init__(self, particles, + small_k=None, medium_k=None, large_k=None, + small_r_limit=None, medium_r_limit=None): + si = const.si + self.particles = particles + self.small_k = small_k or 1.19e6 / si.cm / si.s + self.medium_k = medium_k or 8e3 / si.s + self.large_k = large_k or 2.01e3 * si.cm ** (1 / 2) / si.s + self.small_r_limit = small_r_limit or 35 * si.um + self.medium_r_limit = medium_r_limit or 600 * si.um + + def __call__(self, output, radius): + self.particles.backend.terminal_velocity( + output.data, radius.data, + self.small_k, self.medium_k, self.large_k, + self.small_r_limit, self.medium_r_limit + ) diff --git a/tests/unit_tests/attributes/test_terminal_velocity.py b/tests/unit_tests/attributes/test_terminal_velocity.py index 8e056e75b..9446b3e16 100644 --- a/tests/unit_tests/attributes/test_terminal_velocity.py +++ b/tests/unit_tests/attributes/test_terminal_velocity.py @@ -2,7 +2,8 @@ import matplotlib.pyplot as plt import numpy as np from PySDM.physics import constants as const -from PySDM.physics.terminal_velocity.gunn_and_kinzer import RogersYau, Interpolation +from PySDM.physics.terminal_velocity.gunn_and_kinzer import Interpolation +from PySDM.physics.terminal_velocity.rogers_and_yau import RogersYau from tests.backends_fixture import backend_class from tests.unit_tests.dummy_particulator import DummyParticulator From d2623ff2a4b9031b2a2e8429c6d6e3332b74e2e2 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 24 Dec 2021 10:01:05 +0100 Subject: [PATCH 15/33] negative volumes handling --- PySDM/attributes/physics/volume.py | 2 ++ PySDM/backends/impl_numba/methods/collisions_methods.py | 9 ++++++--- PySDM/backends/impl_numba/storage_impl.py | 3 ++- .../impl_thrust_rtc/methods/collisions_methods.py | 1 + 4 files changed, 11 insertions(+), 4 deletions(-) diff --git a/PySDM/attributes/physics/volume.py b/PySDM/attributes/physics/volume.py index 00b95ae9b..ae8192278 100644 --- a/PySDM/attributes/physics/volume.py +++ b/PySDM/attributes/physics/volume.py @@ -1,5 +1,7 @@ """ particle (wet) volume, key attribute for coalescence +in simmulation involving mixed-phase clouds, positive values correspond to +liquid water and nagative values to ice """ from PySDM.attributes.impl.extensive_attribute import ExtensiveAttribute diff --git a/PySDM/backends/impl_numba/methods/collisions_methods.py b/PySDM/backends/impl_numba/methods/collisions_methods.py index fb522ca32..cdcd8e21a 100644 --- a/PySDM/backends/impl_numba/methods/collisions_methods.py +++ b/PySDM/backends/impl_numba/methods/collisions_methods.py @@ -177,9 +177,12 @@ def linear_collection_efficiency(self, params, output, radii, is_first_in_pair, @numba.njit(**conf.JIT_FLAGS) def interpolation_body(output, radius, factor, b, c): for i in numba.prange(len(radius)): # pylint: disable=not-an-iterable - r_id = int(factor * radius[i]) - r_rest = ((factor * radius[i]) % 1) / factor - output[i] = b[r_id] + r_rest * c[r_id] + if radius[i] < 0: + output[i] = 0 + else: + r_id = int(factor * radius[i]) + r_rest = ((factor * radius[i]) % 1) / factor + output[i] = b[r_id] + r_rest * c[r_id] def interpolation(self, output, radius, factor, b, c): return self.interpolation_body( diff --git a/PySDM/backends/impl_numba/storage_impl.py b/PySDM/backends/impl_numba/storage_impl.py index be6c439e4..77a7f494f 100644 --- a/PySDM/backends/impl_numba/storage_impl.py +++ b/PySDM/backends/impl_numba/storage_impl.py @@ -55,7 +55,8 @@ def sum_out_of_place(output, a, b): @numba.njit(**conf.JIT_FLAGS) def power(output, exponent): - output[:] = np.power(output, exponent) + # TODO #599 (was: output[:] = np.power(output, exponent)) + output[:] = np.sign(output) * np.power(np.abs(output), exponent) @numba.njit(**{**conf.JIT_FLAGS, **{'parallel': False}}) diff --git a/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py b/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py index 1b45149b0..9953dcb9c 100644 --- a/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py +++ b/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py @@ -164,6 +164,7 @@ def __init__(self): ''' ) + # TODO #599 r<0 self.__interpolation_body = trtc.For( ('output', 'radius', 'factor', 'a', 'b'), 'i', From a9bfd75179dc8202b21485b828514419f67526ad Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 24 Dec 2021 10:02:04 +0100 Subject: [PATCH 16/33] TODO labels, one array copy less in IWC product --- PySDM/backends/impl_numba/methods/freezing_methods.py | 5 +++-- PySDM/products/freezing/ice_water_content.py | 7 ++++--- .../freezing/total_unfrozen_immersed_surface_area.py | 1 + 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/freezing_methods.py b/PySDM/backends/impl_numba/methods/freezing_methods.py index 7def58270..677e38113 100644 --- a/PySDM/backends/impl_numba/methods/freezing_methods.py +++ b/PySDM/backends/impl_numba/methods/freezing_methods.py @@ -25,10 +25,11 @@ def _freeze(volume, i): @numba.njit(**{**conf.JIT_FLAGS, 'fastmath': self.formulae.fastmath}) def freeze_singular_body(attributes, temperature, relative_humidity, cell): - for i in numba.prange(len(attributes.freezing_temperature)): # pylint: disable=not-an-iterable + n_sd = len(attributes.freezing_temperature) + for i in numba.prange(n_sd): # pylint: disable=not-an-iterable if ( _unfrozen(attributes.wet_volume, i) and - relative_humidity[cell[i]] > 1 and + relative_humidity[cell[i]] > 1 and # TODO #599 - it is in Shima's formulation, is it needed? temperature[cell[i]] <= attributes.freezing_temperature[i] ): _freeze(attributes.wet_volume, i) diff --git a/PySDM/products/freezing/ice_water_content.py b/PySDM/products/freezing/ice_water_content.py index e5528454c..b9266e94a 100644 --- a/PySDM/products/freezing/ice_water_content.py +++ b/PySDM/products/freezing/ice_water_content.py @@ -12,11 +12,12 @@ def __init__(self, unit='kg/kg', name=None, __specific=False): self.specific = __specific def _impl(self, **kwargs): - self._download_moment_to_buffer('volume', rank=0, filter_range=(-np.inf, 0)) - conc = self.buffer.copy() - self._download_moment_to_buffer('volume', rank=1, filter_range=(-np.inf, 0)) result = self.buffer.copy() + + self._download_moment_to_buffer('volume', rank=0, filter_range=(-np.inf, 0)) + conc = self.buffer + result[:] *= -const.rho_i * conc / self.particulator.mesh.dv if self.specific: diff --git a/PySDM/products/freezing/total_unfrozen_immersed_surface_area.py b/PySDM/products/freezing/total_unfrozen_immersed_surface_area.py index 008372f7f..e9b55e1df 100644 --- a/PySDM/products/freezing/total_unfrozen_immersed_surface_area.py +++ b/PySDM/products/freezing/total_unfrozen_immersed_surface_area.py @@ -19,4 +19,5 @@ def _impl(self, **kwargs): result = np.copy(self.buffer) self._download_moment_to_buffer(**params, rank=0) result[:] *= self.buffer + # TODO #599 per volume / per gridbox ? return result From 73ce0be22f309f42f8cd5319578f6b3089728d5d Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 24 Dec 2021 10:08:04 +0100 Subject: [PATCH 17/33] shorter line --- PySDM/backends/impl_numba/methods/freezing_methods.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PySDM/backends/impl_numba/methods/freezing_methods.py b/PySDM/backends/impl_numba/methods/freezing_methods.py index 677e38113..19a718a18 100644 --- a/PySDM/backends/impl_numba/methods/freezing_methods.py +++ b/PySDM/backends/impl_numba/methods/freezing_methods.py @@ -29,7 +29,7 @@ def freeze_singular_body(attributes, temperature, relative_humidity, cell): for i in numba.prange(n_sd): # pylint: disable=not-an-iterable if ( _unfrozen(attributes.wet_volume, i) and - relative_humidity[cell[i]] > 1 and # TODO #599 - it is in Shima's formulation, is it needed? + relative_humidity[cell[i]] > 1 and # TODO #599 as in Shima, but is it needed? temperature[cell[i]] <= attributes.freezing_temperature[i] ): _freeze(attributes.wet_volume, i) From 3aa566fa5549260dfafe096b4d4513cd020dacab Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 24 Dec 2021 10:14:32 +0100 Subject: [PATCH 18/33] moving terminal velocity logic out of collision methods --- .../impl_numba/methods/collisions_methods.py | 54 +----------- .../methods/terminal_velocity_methods.py | 57 ++++++++++++ .../methods/collisions_methods.py | 75 ---------------- .../methods/terminal_velocity_methods.py | 87 +++++++++++++++++++ PySDM/backends/numba.py | 10 ++- PySDM/backends/thrust_rtc.py | 5 +- 6 files changed, 157 insertions(+), 131 deletions(-) create mode 100644 PySDM/backends/impl_numba/methods/terminal_velocity_methods.py create mode 100644 PySDM/backends/impl_thrust_rtc/methods/terminal_velocity_methods.py diff --git a/PySDM/backends/impl_numba/methods/collisions_methods.py b/PySDM/backends/impl_numba/methods/collisions_methods.py index cdcd8e21a..b42e6b43d 100644 --- a/PySDM/backends/impl_numba/methods/collisions_methods.py +++ b/PySDM/backends/impl_numba/methods/collisions_methods.py @@ -16,7 +16,7 @@ def pair_indices(i, idx, is_first_in_pair): return j, k -class AlgorithmicMethods(BackendMethods): +class CollisionsMethods(BackendMethods): @staticmethod @numba.njit(**{**conf.JIT_FLAGS, **{'parallel': False}}) def adaptive_sdm_end_body(dt_left, n_cell, cell_start): @@ -141,54 +141,6 @@ def compute_gamma(self, gamma, rand, multiplicity, cell_id, cell_id.data, collision_rate_deficit.data, collision_rate.data, is_first_in_pair.indicator.data) - @staticmethod - @numba.njit(**conf.JIT_FLAGS) - # pylint: disable=too-many-arguments,too-many-locals - def linear_collection_efficiency_body( - params, output, radii, is_first_in_pair, idx, length, unit - ): - A, B, D1, D2, E1, E2, F1, F2, G1, G2, G3, Mf, Mg = params - output[:] = 0 - for i in numba.prange(length - 1): # pylint: disable=not-an-iterable - if is_first_in_pair[i]: - if radii[idx[i]] > radii[idx[i + 1]]: - r = radii[idx[i]] / unit - r_s = radii[idx[i + 1]] / unit - else: - r = radii[idx[i + 1]] / unit - r_s = radii[idx[i]] / unit - p = r_s / r - if p not in (0, 1): - G = (G1 / r) ** Mg + G2 + G3 * r - Gp = (1 - p) ** G - if Gp != 0: - D = D1 / r ** D2 - E = E1 / r ** E2 - F = (F1 / r) ** Mf + F2 - output[i // 2] = A + B * p + D / p ** F + E / Gp - output[i // 2] = max(0, output[i // 2]) - - def linear_collection_efficiency(self, params, output, radii, is_first_in_pair, unit): - return self.linear_collection_efficiency_body( - params, output.data, radii.data, is_first_in_pair.indicator.data, - radii.idx.data, len(is_first_in_pair), unit) - - @staticmethod - @numba.njit(**conf.JIT_FLAGS) - def interpolation_body(output, radius, factor, b, c): - for i in numba.prange(len(radius)): # pylint: disable=not-an-iterable - if radius[i] < 0: - output[i] = 0 - else: - r_id = int(factor * radius[i]) - r_rest = ((factor * radius[i]) % 1) / factor - output[i] = b[r_id] + r_rest * c[r_id] - - def interpolation(self, output, radius, factor, b, c): - return self.interpolation_body( - output.data, radius.data, factor, b.data, c.data - ) - @staticmethod def make_cell_caretaker(idx, cell_start, scheme="default"): class CellCaretaker: # pylint: disable=too-few-public-methods @@ -210,11 +162,11 @@ def __init__(self, idx, cell_start, scheme): def __call__(self, cell_id, cell_idx, cell_start, idx): length = len(idx) if self.scheme == "counting_sort": - AlgorithmicMethods._counting_sort_by_cell_id_and_update_cell_start( + CollisionsMethods._counting_sort_by_cell_id_and_update_cell_start( self.tmp_idx.data, idx.data, cell_id.data, cell_idx.data, length, cell_start.data) elif self.scheme == "counting_sort_parallel": - AlgorithmicMethods._parallel_counting_sort_by_cell_id_and_update_cell_start( + CollisionsMethods._parallel_counting_sort_by_cell_id_and_update_cell_start( self.tmp_idx.data, idx.data, cell_id.data, cell_idx.data, length, cell_start.data, self.cell_starts.data) idx.data, self.tmp_idx.data = self.tmp_idx.data, idx.data diff --git a/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py b/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py new file mode 100644 index 000000000..9c3380257 --- /dev/null +++ b/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py @@ -0,0 +1,57 @@ +""" +CPU implementation of backend methods for terminal velocities +""" +import numba +import numpy as np +from PySDM.backends.impl_numba import conf +from PySDM.backends.impl_common.backend_methods import BackendMethods + + +class TerminalVelocityMethods(BackendMethods): + @staticmethod + @numba.njit(**conf.JIT_FLAGS) + # pylint: disable=too-many-arguments,too-many-locals + def linear_collection_efficiency_body( + params, output, radii, is_first_in_pair, idx, length, unit + ): + A, B, D1, D2, E1, E2, F1, F2, G1, G2, G3, Mf, Mg = params + output[:] = 0 + for i in numba.prange(length - 1): # pylint: disable=not-an-iterable + if is_first_in_pair[i]: + if radii[idx[i]] > radii[idx[i + 1]]: + r = radii[idx[i]] / unit + r_s = radii[idx[i + 1]] / unit + else: + r = radii[idx[i + 1]] / unit + r_s = radii[idx[i]] / unit + p = r_s / r + if p not in (0, 1): + G = (G1 / r) ** Mg + G2 + G3 * r + Gp = (1 - p) ** G + if Gp != 0: + D = D1 / r ** D2 + E = E1 / r ** E2 + F = (F1 / r) ** Mf + F2 + output[i // 2] = A + B * p + D / p ** F + E / Gp + output[i // 2] = max(0, output[i // 2]) + + def linear_collection_efficiency(self, params, output, radii, is_first_in_pair, unit): + return self.linear_collection_efficiency_body( + params, output.data, radii.data, is_first_in_pair.indicator.data, + radii.idx.data, len(is_first_in_pair), unit) + + @staticmethod + @numba.njit(**conf.JIT_FLAGS) + def interpolation_body(output, radius, factor, b, c): + for i in numba.prange(len(radius)): # pylint: disable=not-an-iterable + if radius[i] < 0: + output[i] = 0 + else: + r_id = int(factor * radius[i]) + r_rest = ((factor * radius[i]) % 1) / factor + output[i] = b[r_id] + r_rest * c[r_id] + + def interpolation(self, output, radius, factor, b, c): + return self.interpolation_body( + output.data, radius.data, factor, b.data, c.data + ) diff --git a/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py b/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py index 9953dcb9c..4e4bcecfa 100644 --- a/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py +++ b/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py @@ -10,40 +10,6 @@ class CollisionsMethods(ThrustRTCBackendMethods): def __init__(self): super().__init__() - self.__linear_collection_efficiency_body = trtc.For( - ( - 'A', 'B', 'D1', 'D2', 'E1', 'E2', 'F1', 'F2', 'G1', 'G2', 'G3', 'Mf', 'Mg', - 'output', 'radii', 'is_first_in_pair', 'idx', 'unit' - ), - "i", - ''' - if (is_first_in_pair[i]) { - real_type r = 0; - real_type r_s = 0; - if (radii[idx[i]] > radii[idx[i + 1]]) { - r = radii[idx[i]] / unit; - r_s = radii[idx[i + 1]] / unit; - } - else { - r = radii[idx[i + 1]] / unit; - r_s = radii[idx[i]] / unit; - } - real_type p = r_s / r; - if (p != 0 && p != 1) { - real_type G = pow((G1 / r), Mg) + G2 + G3 * r; - real_type Gp = pow((1 - p), G); - if (Gp != 0) { - real_type D = D1 / pow(r, D2); - real_type E = E1 / pow(r, E2); - real_type F = pow((F1 / r), Mf) + F2; - output[int(i / 2)] = A + B * p + D / pow(p, F) + E / Gp; - if (output[int(i / 2)] < 0) { - output[int(i / 2)] = 0; - } - } - } - } - '''.replace("real_type", self._get_c_type())) self.__adaptive_sdm_gamma_body_1 = trtc.For( ('dt_todo', 'dt_left', 'dt_max'), @@ -164,17 +130,6 @@ def __init__(self): ''' ) - # TODO #599 r<0 - self.__interpolation_body = trtc.For( - ('output', 'radius', 'factor', 'a', 'b'), - 'i', - ''' - auto r_id = (int64_t)(factor * radius[i]); - auto r_rest = (factor * radius[i] - r_id) / factor; - output[i] = a[r_id] + r_rest * b[r_id]; - ''' - ) - self.__normalize_body_0 = trtc.For( ('cell_start', 'norm_factor', 'dt_div_dv'), "i", @@ -288,36 +243,6 @@ def compute_gamma(self, gamma, rand, multiplicity, cell_id, (gamma.data, rand.data, multiplicity.idx.data, multiplicity.data, cell_id.data, collision_rate_deficit.data, collision_rate.data)) - def linear_collection_efficiency(self, params, output, radii, is_first_in_pair, unit): - A, B, D1, D2, E1, E2, F1, F2, G1, G2, G3, Mf, Mg = params - dA = self._get_floating_point(A) - dB = self._get_floating_point(B) - dD1 = self._get_floating_point(D1) - dD2 = self._get_floating_point(D2) - dE1 = self._get_floating_point(E1) - dE2 = self._get_floating_point(E2) - dF1 = self._get_floating_point(F1) - dF2 = self._get_floating_point(F2) - dG1 = self._get_floating_point(G1) - dG2 = self._get_floating_point(G2) - dG3 = self._get_floating_point(G3) - dMf = self._get_floating_point(Mf) - dMg = self._get_floating_point(Mg) - dunit = self._get_floating_point(unit) - trtc.Fill(output.data, trtc.DVDouble(0)) - self.__linear_collection_efficiency_body.launch_n( - len(is_first_in_pair) - 1, - [dA, dB, dD1, dD2, dE1, dE2, dF1, dF2, dG1, dG2, dG3, dMf, dMg, - output.data, radii.data, is_first_in_pair.indicator.data, radii.idx.data, dunit]) - - @nice_thrust(**NICE_THRUST_FLAGS) - def interpolation(self, output, radius, factor, b, c): - factor_device = trtc.DVInt64(factor) - self.__interpolation_body.launch_n( - len(radius), - (output.data, radius.data, factor_device, b.data, c.data) - ) - # pylint: disable=unused-argument def make_cell_caretaker(self, idx, cell_start, scheme=None): return self._sort_by_cell_id_and_update_cell_start diff --git a/PySDM/backends/impl_thrust_rtc/methods/terminal_velocity_methods.py b/PySDM/backends/impl_thrust_rtc/methods/terminal_velocity_methods.py new file mode 100644 index 000000000..b23520618 --- /dev/null +++ b/PySDM/backends/impl_thrust_rtc/methods/terminal_velocity_methods.py @@ -0,0 +1,87 @@ +""" +GPU implementation of backend methods for terminal velocities +""" +from PySDM.backends.impl_thrust_rtc.conf import NICE_THRUST_FLAGS +from PySDM.backends.impl_thrust_rtc.nice_thrust import nice_thrust +from ..conf import trtc +from ..methods.thrust_rtc_backend_methods import ThrustRTCBackendMethods + + +class TerminalVelocityMethods(ThrustRTCBackendMethods): + def __init__(self): + super().__init__() + self.__linear_collection_efficiency_body = trtc.For( + ( + 'A', 'B', 'D1', 'D2', 'E1', 'E2', 'F1', 'F2', 'G1', 'G2', 'G3', 'Mf', 'Mg', + 'output', 'radii', 'is_first_in_pair', 'idx', 'unit' + ), + "i", + ''' + if (is_first_in_pair[i]) { + real_type r = 0; + real_type r_s = 0; + if (radii[idx[i]] > radii[idx[i + 1]]) { + r = radii[idx[i]] / unit; + r_s = radii[idx[i + 1]] / unit; + } + else { + r = radii[idx[i + 1]] / unit; + r_s = radii[idx[i]] / unit; + } + real_type p = r_s / r; + if (p != 0 && p != 1) { + real_type G = pow((G1 / r), Mg) + G2 + G3 * r; + real_type Gp = pow((1 - p), G); + if (Gp != 0) { + real_type D = D1 / pow(r, D2); + real_type E = E1 / pow(r, E2); + real_type F = pow((F1 / r), Mf) + F2; + output[int(i / 2)] = A + B * p + D / pow(p, F) + E / Gp; + if (output[int(i / 2)] < 0) { + output[int(i / 2)] = 0; + } + } + } + } + '''.replace("real_type", self._get_c_type())) + + # TODO #599 r<0 + self.__interpolation_body = trtc.For( + ('output', 'radius', 'factor', 'a', 'b'), + 'i', + ''' + auto r_id = (int64_t)(factor * radius[i]); + auto r_rest = (factor * radius[i] - r_id) / factor; + output[i] = a[r_id] + r_rest * b[r_id]; + ''' + ) + + def linear_collection_efficiency(self, params, output, radii, is_first_in_pair, unit): + A, B, D1, D2, E1, E2, F1, F2, G1, G2, G3, Mf, Mg = params + dA = self._get_floating_point(A) + dB = self._get_floating_point(B) + dD1 = self._get_floating_point(D1) + dD2 = self._get_floating_point(D2) + dE1 = self._get_floating_point(E1) + dE2 = self._get_floating_point(E2) + dF1 = self._get_floating_point(F1) + dF2 = self._get_floating_point(F2) + dG1 = self._get_floating_point(G1) + dG2 = self._get_floating_point(G2) + dG3 = self._get_floating_point(G3) + dMf = self._get_floating_point(Mf) + dMg = self._get_floating_point(Mg) + dunit = self._get_floating_point(unit) + trtc.Fill(output.data, trtc.DVDouble(0)) + self.__linear_collection_efficiency_body.launch_n( + len(is_first_in_pair) - 1, + [dA, dB, dD1, dD2, dE1, dE2, dF1, dF2, dG1, dG2, dG3, dMf, dMg, + output.data, radii.data, is_first_in_pair.indicator.data, radii.idx.data, dunit]) + + @nice_thrust(**NICE_THRUST_FLAGS) + def interpolation(self, output, radius, factor, b, c): + factor_device = trtc.DVInt64(factor) + self.__interpolation_body.launch_n( + len(radius), + (output.data, radius.data, factor_device, b.data, c.data) + ) diff --git a/PySDM/backends/numba.py b/PySDM/backends/numba.py index 8015e4d37..e91d0ce96 100644 --- a/PySDM/backends/numba.py +++ b/PySDM/backends/numba.py @@ -2,7 +2,7 @@ Multi-threaded CPU backend using LLVM-powered just-in-time compilation """ -from PySDM.backends.impl_numba.methods.collisions_methods import AlgorithmicMethods +from PySDM.backends.impl_numba.methods.collisions_methods import CollisionsMethods from PySDM.backends.impl_numba.methods.pair_methods import PairMethods from PySDM.backends.impl_numba.methods.physics_methods import PhysicsMethods from PySDM.backends.impl_numba.methods.index_methods import IndexMethods @@ -11,13 +11,14 @@ from PySDM.backends.impl_numba.methods.freezing_methods import FreezingMethods from PySDM.backends.impl_numba.methods.chemistry_methods import ChemistryMethods from PySDM.backends.impl_numba.methods.displacement_methods import DisplacementMethods +from PySDM.backends.impl_numba.methods.terminal_velocity_methods import TerminalVelocityMethods from PySDM.backends.impl_numba.random import Random as ImportedRandom from PySDM.backends.impl_numba.storage import Storage as ImportedStorage from PySDM.formulae import Formulae class Numba( # pylint: disable=too-many-ancestors,duplicate-code - AlgorithmicMethods, + CollisionsMethods, PairMethods, IndexMethods, PhysicsMethods, @@ -25,7 +26,8 @@ class Numba( # pylint: disable=too-many-ancestors,duplicate-code ChemistryMethods, MomentsMethods, FreezingMethods, - DisplacementMethods + DisplacementMethods, + TerminalVelocityMethods ): Storage = ImportedStorage Random = ImportedRandom @@ -34,7 +36,7 @@ class Numba( # pylint: disable=too-many-ancestors,duplicate-code def __init__(self, formulae=None): self.formulae = formulae or Formulae() - AlgorithmicMethods.__init__(self) + CollisionsMethods.__init__(self) PairMethods.__init__(self) IndexMethods.__init__(self) PhysicsMethods.__init__(self) diff --git a/PySDM/backends/thrust_rtc.py b/PySDM/backends/thrust_rtc.py index d2c0b8fb1..7c097dcc6 100644 --- a/PySDM/backends/thrust_rtc.py +++ b/PySDM/backends/thrust_rtc.py @@ -12,6 +12,8 @@ from PySDM.backends.impl_thrust_rtc.methods.moments_methods import MomentsMethods from PySDM.backends.impl_thrust_rtc.methods.condensation_methods import CondensationMethods from PySDM.backends.impl_thrust_rtc.methods.displacement_methods import DisplacementMethods +from PySDM.backends.impl_thrust_rtc.methods.terminal_velocity_methods import \ + TerminalVelocityMethods from PySDM.backends.impl_thrust_rtc.storage import make_storage_class from PySDM.backends.impl_thrust_rtc.random import Random as ImportedRandom from PySDM.formulae import Formulae @@ -25,7 +27,8 @@ class ThrustRTC( # pylint: disable=duplicate-code,too-many-ancestors PhysicsMethods, CondensationMethods, MomentsMethods, - DisplacementMethods + DisplacementMethods, + TerminalVelocityMethods ): ENABLE = True Random = ImportedRandom From 40348fda278396c1f99747423176afda36d776fc Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 24 Dec 2021 13:17:40 +0100 Subject: [PATCH 19/33] ctor calls for terminal vel methods --- PySDM/backends/impl_numba/methods/terminal_velocity_methods.py | 1 - PySDM/backends/numba.py | 1 + PySDM/backends/thrust_rtc.py | 1 + 3 files changed, 2 insertions(+), 1 deletion(-) diff --git a/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py b/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py index 9c3380257..aee4fb463 100644 --- a/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py +++ b/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py @@ -2,7 +2,6 @@ CPU implementation of backend methods for terminal velocities """ import numba -import numpy as np from PySDM.backends.impl_numba import conf from PySDM.backends.impl_common.backend_methods import BackendMethods diff --git a/PySDM/backends/numba.py b/PySDM/backends/numba.py index e91d0ce96..b6ff67151 100644 --- a/PySDM/backends/numba.py +++ b/PySDM/backends/numba.py @@ -45,3 +45,4 @@ def __init__(self, formulae=None): MomentsMethods.__init__(self) FreezingMethods.__init__(self) DisplacementMethods.__init__(self) + TerminalVelocityMethods.__init__(self) diff --git a/PySDM/backends/thrust_rtc.py b/PySDM/backends/thrust_rtc.py index 7c097dcc6..9883d9dd7 100644 --- a/PySDM/backends/thrust_rtc.py +++ b/PySDM/backends/thrust_rtc.py @@ -51,6 +51,7 @@ def __init__(self, formulae=None, double_precision=False, debug=False, verbose=F CondensationMethods.__init__(self) MomentsMethods.__init__(self) DisplacementMethods.__init__(self) + TerminalVelocityMethods.__init__(self) trtc.Set_Kernel_Debug(debug) trtc.Set_Verbose(verbose) From cd07ec7fdad0a18733525fecb3d6b036199e1776 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sun, 26 Dec 2021 08:52:00 +0100 Subject: [PATCH 20/33] specify seed in test_freezing --- tests/smoke_tests/arabas_et_al_2015/test_freezing.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py index 123e495b5..cc65c4cb3 100644 --- a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py +++ b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py @@ -24,6 +24,7 @@ def test_freezing(singular): # Arrange settings = Settings(Formulae( + seed=44, condensation_coordinate='VolumeLogarithm', fastmath=True, freezing_temperature_spectrum='Niemand_et_al_2012', From 7076f4180e9764b13a91e1d3c6d831ab94ef1e9f Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sun, 26 Dec 2021 09:51:39 +0100 Subject: [PATCH 21/33] adding the RadiusBinnedNumberAveragedTerminalVelocity product --- ...inned_number_averaged_terminal_velocity.py | 35 +++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 PySDM/products/size_spectral/radius_binned_number_averaged_terminal_velocity.py diff --git a/PySDM/products/size_spectral/radius_binned_number_averaged_terminal_velocity.py b/PySDM/products/size_spectral/radius_binned_number_averaged_terminal_velocity.py new file mode 100644 index 000000000..672a7026c --- /dev/null +++ b/PySDM/products/size_spectral/radius_binned_number_averaged_terminal_velocity.py @@ -0,0 +1,35 @@ +import numpy as np +from PySDM.products.impl.spectrum_moment_product import SpectrumMomentProduct + +ATTR = 'terminal velocity' +RANK = 1 + +class RadiusBinnedNumberAveragedTerminalVelocity(SpectrumMomentProduct): + def __init__(self, radius_bin_edges, name=None, unit='m/s'): + super().__init__(name=name, unit=unit) + self.radius_bin_edges = radius_bin_edges + + def register(self, builder): + builder.request_attribute(ATTR) + + volume_bin_edges = builder.particulator.formulae.trivia.volume(self.radius_bin_edges) + self.attr_bins_edges = builder.particulator.backend.Storage.from_ndarray(volume_bin_edges) + + super().register(builder) + + self.shape = (*builder.particulator.mesh.grid, len(self.attr_bins_edges) - 1) + + + def _impl(self, **kwargs): + vals = np.empty([self.particulator.mesh.n_cell, len(self.attr_bins_edges) - 1]) + + self._recalculate_spectrum_moment( + attr=ATTR, + rank=RANK, + ) + + for i in range(vals.shape[1]): + self._download_spectrum_moment_to_buffer(rank=RANK, bin_number=i) + vals[:, i] = self.buffer.ravel() + + return np.squeeze(vals.reshape(self.shape)) From f8e096ba7decd91ad753a9beeed31a5973975c78 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sun, 26 Dec 2021 09:51:58 +0100 Subject: [PATCH 22/33] adding test_gui_settings for ICMW --- .../arabas_et_al_2015/test_gui_settings.py | 25 +++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 tests/smoke_tests/arabas_et_al_2015/test_gui_settings.py diff --git a/tests/smoke_tests/arabas_et_al_2015/test_gui_settings.py b/tests/smoke_tests/arabas_et_al_2015/test_gui_settings.py new file mode 100644 index 000000000..aa5e01789 --- /dev/null +++ b/tests/smoke_tests/arabas_et_al_2015/test_gui_settings.py @@ -0,0 +1,25 @@ +# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring +from PySDM_examples.Arabas_et_al_2015 import Settings +from PySDM_examples.Szumowski_et_al_1998.gui_settings import GUISettings + + +class TestGUISettings: + @staticmethod + def test_instantiate(): + gui_settings = GUISettings(Settings()) + + @staticmethod + def test_stream_function(): + # arrange + gui_settings = GUISettings(Settings()) + gui_settings.ui_rhod_w_max = None + failed = False + + # act + try: + _ = gui_settings.stream_function(0, 0, 0) + except AttributeError: + failed = True + + # assert + assert failed From 6a96a4e9221429d921ba2d6b5531fad1d258577c Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sun, 26 Dec 2021 10:01:25 +0100 Subject: [PATCH 23/33] pylint fixes --- .../radius_binned_number_averaged_terminal_velocity.py | 3 +++ tests/smoke_tests/arabas_et_al_2015/test_gui_settings.py | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/PySDM/products/size_spectral/radius_binned_number_averaged_terminal_velocity.py b/PySDM/products/size_spectral/radius_binned_number_averaged_terminal_velocity.py index 672a7026c..18402568f 100644 --- a/PySDM/products/size_spectral/radius_binned_number_averaged_terminal_velocity.py +++ b/PySDM/products/size_spectral/radius_binned_number_averaged_terminal_velocity.py @@ -1,3 +1,6 @@ +""" +Provides radius bin-resolved average terminal velocity (average is particle-number weighted) +""" import numpy as np from PySDM.products.impl.spectrum_moment_product import SpectrumMomentProduct diff --git a/tests/smoke_tests/arabas_et_al_2015/test_gui_settings.py b/tests/smoke_tests/arabas_et_al_2015/test_gui_settings.py index aa5e01789..a474a9497 100644 --- a/tests/smoke_tests/arabas_et_al_2015/test_gui_settings.py +++ b/tests/smoke_tests/arabas_et_al_2015/test_gui_settings.py @@ -6,7 +6,7 @@ class TestGUISettings: @staticmethod def test_instantiate(): - gui_settings = GUISettings(Settings()) + _ = GUISettings(Settings()) @staticmethod def test_stream_function(): From ebe9da19a025fd91d27b2eee549a24c0aa200d35 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sun, 26 Dec 2021 10:02:24 +0100 Subject: [PATCH 24/33] adding missing import in products' __init__ --- PySDM/products/size_spectral/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/PySDM/products/size_spectral/__init__.py b/PySDM/products/size_spectral/__init__.py index 208b35f23..f622f398c 100644 --- a/PySDM/products/size_spectral/__init__.py +++ b/PySDM/products/size_spectral/__init__.py @@ -7,3 +7,5 @@ from .total_particle_concentration import TotalParticleConcentration from .total_particle_specific_concentration import TotalParticleSpecificConcentration from .water_mixing_ratio import WaterMixingRatio +from .radius_binned_number_averaged_terminal_velocity \ + import RadiusBinnedNumberAveragedTerminalVelocity From 7c950b1612bf73320bfc167056c8eaac37ca354e Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sun, 26 Dec 2021 11:16:52 +0100 Subject: [PATCH 25/33] handling ctor arg for RadiusBinnedNumberAveragedTerminalVelocity in test_products --- tests/unit_tests/test_products.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_products.py b/tests/unit_tests/test_products.py index 62f82240b..dc03476e5 100644 --- a/tests/unit_tests/test_products.py +++ b/tests/unit_tests/test_products.py @@ -12,7 +12,8 @@ from PySDM.products import (AqueousMassSpectrum, AqueousMoleFraction, TotalDryMassMixingRatio, ParticleSizeSpectrumPerMass, GaseousMoleFraction, FreezableSpecificConcentration, DynamicWallTime, - ParticleSizeSpectrumPerVolume, ParticlesVolumeSpectrum) + ParticleSizeSpectrumPerVolume, ParticlesVolumeSpectrum, + RadiusBinnedNumberAveragedTerminalVelocity) _ARGUMENTS = { AqueousMassSpectrum: {'key': 'S_VI', 'dry_radius_bins_edges': (0, np.inf)}, @@ -23,7 +24,8 @@ FreezableSpecificConcentration: {'temperature_bins_edges': (0, 300)}, DynamicWallTime: {'dynamic': 'Condensation'}, ParticleSizeSpectrumPerVolume: {'radius_bins_edges': (0, np.inf)}, - ParticlesVolumeSpectrum: {'radius_bins_edges': (0, np.inf)} + ParticlesVolumeSpectrum: {'radius_bins_edges': (0, np.inf)}, + RadiusBinnedNumberAveragedTerminalVelocity: {'radius_bin_edges': (0, np.inf)} } From 02b826138b96c6eec0af283ab4346ebeb05097ad Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sun, 26 Dec 2021 20:12:21 +0100 Subject: [PATCH 26/33] test_freezing settings update --- tests/smoke_tests/arabas_et_al_2015/test_freezing.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py index cc65c4cb3..2180a4e85 100644 --- a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py +++ b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py @@ -31,7 +31,8 @@ def test_freezing(singular): heterogeneous_ice_nucleation_rate='ABIFM' )) settings.dt = .5 * si.second - settings.grid = (3, 25) + settings.grid = (5, 15) + settings.n_sd_per_gridbox = 16 settings.simulation_time = 100 * settings.dt settings.spin_up_time = 10 * settings.dt From 2b69bb986c6ea462310517e7cc34c83a2911c03c Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Mon, 27 Dec 2021 09:44:15 +0100 Subject: [PATCH 27/33] better handling of non-spatial dimensions in netCDF exporter (e.g., terminal velocity bins). closes #340 --- PySDM/exporters/netcdf_exporter.py | 24 ++++++++++++++----- .../aqueous_mass_spectrum.py | 2 +- .../freezable_specific_concentration.py | 2 +- PySDM/products/impl/product.py | 6 ++--- .../products/impl/spectrum_moment_product.py | 4 +++- .../size_spectral/particle_size_spectrum.py | 2 +- .../particles_volume_spectrum.py | 2 +- ...inned_number_averaged_terminal_velocity.py | 2 +- 8 files changed, 29 insertions(+), 15 deletions(-) diff --git a/PySDM/exporters/netcdf_exporter.py b/PySDM/exporters/netcdf_exporter.py index 532d54b49..b3ed5012b 100644 --- a/PySDM/exporters/netcdf_exporter.py +++ b/PySDM/exporters/netcdf_exporter.py @@ -3,8 +3,11 @@ """ import numpy as np from scipy.io.netcdf import netcdf_file +from PySDM.products.impl.spectrum_moment_product import SpectrumMomentProduct +DIM_SUFFIX = '_bin_left_edges' + class NetCDFExporter: def __init__(self, storage, settings, simulator, filename): self.storage = storage @@ -20,12 +23,16 @@ def _write_settings(self, ncdf): def _create_dimensions(self, ncdf): ncdf.createDimension("T", len(self.settings.output_steps)) + for index, label in enumerate(self.XZ): ncdf.createDimension(label, self.settings.grid[index]) - ncdf.createDimension( - "ParticleVolume", - len(self.settings.formulae.trivia.volume(self.settings.r_bins_edges)) - 1 - ) + + for name, instance in self.simulator.products.items(): + if isinstance(instance, SpectrumMomentProduct): + ncdf.createDimension( + f"{name}{DIM_SUFFIX}", + len(instance.attr_bins_edges) - 1 + ) def _create_variables(self, ncdf): self.vars = {} @@ -38,7 +45,12 @@ def _create_variables(self, ncdf): * (1 / 2 + np.arange(self.settings.grid[index]))) self.vars[label].units = "metres" - # TODO #340 ParticleVolume var + for name, instance in self.simulator.products.items(): + if isinstance(instance, SpectrumMomentProduct): + label = f"{name}{DIM_SUFFIX}" + self.vars[label] = ncdf.createVariable(label, 'f', (label,)) + self.vars[label][:] = instance.attr_bins_edges.to_ndarray()[:-1] + self.vars[label].units = instance.attr_unit for name, instance in self.simulator.products.items(): if name in self.vars: @@ -46,7 +58,7 @@ def _create_variables(self, ncdf): n_dimensions = len(instance.shape) if n_dimensions == 3: - dimensions = ("T", "X", "Z", "ParticleVolume") + dimensions = ("T", "X", "Z", f"{name}{DIM_SUFFIX}") elif n_dimensions == 2: dimensions = ("T", "X", "Z") elif n_dimensions == 0: diff --git a/PySDM/products/aqueous_chemistry/aqueous_mass_spectrum.py b/PySDM/products/aqueous_chemistry/aqueous_mass_spectrum.py index 84935aeb7..b91086cf8 100644 --- a/PySDM/products/aqueous_chemistry/aqueous_mass_spectrum.py +++ b/PySDM/products/aqueous_chemistry/aqueous_mass_spectrum.py @@ -11,7 +11,7 @@ class AqueousMassSpectrum(SpectrumMomentProduct): def __init__(self, key, dry_radius_bins_edges, specific=False, name=None, unit='kg/m^3'): - super().__init__(name=name, unit=unit) + super().__init__(name=name, unit=unit, attr_unit='m') self.key = key self.dry_radius_bins_edges = dry_radius_bins_edges self.molar_mass = Substance.from_formula(AQUEOUS_COMPOUNDS[key][0]).mass * si.g / si.mole diff --git a/PySDM/products/freezing/freezable_specific_concentration.py b/PySDM/products/freezing/freezable_specific_concentration.py index 98a56013a..3e7f57ea4 100644 --- a/PySDM/products/freezing/freezable_specific_concentration.py +++ b/PySDM/products/freezing/freezable_specific_concentration.py @@ -7,7 +7,7 @@ class FreezableSpecificConcentration(SpectrumMomentProduct): def __init__(self, temperature_bins_edges, name=None, unit='kg^-1 K^-1'): - super().__init__(name=name, unit=unit) + super().__init__(name=name, unit=unit, attr_unit='K') self.attr_bins_edges = temperature_bins_edges def register(self, builder): diff --git a/PySDM/products/impl/product.py b/PySDM/products/impl/product.py index 385a09c5b..8975728b5 100644 --- a/PySDM/products/impl/product.py +++ b/PySDM/products/impl/product.py @@ -17,7 +17,7 @@ class Product: def __init__(self, *, unit: str, name: Optional[str] = None): self.name = name or self._camel_case_to_words(self.__class__.__name__) - self._unit = self.__parse_unit(unit) + self._unit = self._parse_unit(unit) self.unit_magnitude_in_base_units = self._unit.to_base_units().magnitude self.__check_unit() @@ -36,7 +36,7 @@ def _download_to_buffer(self, storage): storage.download(self.buffer.ravel()) @staticmethod - def __parse_unit(unit: str): + def _parse_unit(unit: str): if unit in ('%', 'percent'): return .01 * _UNIT_REGISTRY.dimensionless if unit in ('PPB', 'ppb'): @@ -65,7 +65,7 @@ def __check_unit(self): raise AssertionError(f"unit parameter of {type(self).__name__}.__init__" f" is expected to have a non-empty default value") - default_unit = self.__parse_unit(default_unit_arg) + default_unit = self._parse_unit(default_unit_arg) if default_unit.to_base_units().magnitude != 1: raise AssertionError(f'default value "{default_unit_arg}"' diff --git a/PySDM/products/impl/spectrum_moment_product.py b/PySDM/products/impl/spectrum_moment_product.py index 1296a32f2..201d49214 100644 --- a/PySDM/products/impl/spectrum_moment_product.py +++ b/PySDM/products/impl/spectrum_moment_product.py @@ -7,9 +7,11 @@ class SpectrumMomentProduct(ABC, Product): - def __init__(self, name, unit): + def __init__(self, name, unit, attr_unit): super().__init__(name=name, unit=unit) + _ = self._parse_unit(unit) self.attr_bins_edges = None + self.attr_unit = attr_unit self.moment_0 = None self.moments = None diff --git a/PySDM/products/size_spectral/particle_size_spectrum.py b/PySDM/products/size_spectral/particle_size_spectrum.py index cf7685cd5..3d25f9e17 100644 --- a/PySDM/products/size_spectral/particle_size_spectrum.py +++ b/PySDM/products/size_spectral/particle_size_spectrum.py @@ -11,7 +11,7 @@ def __init__(self, radius_bins_edges, name, unit, dry=False, normalise_by_dv=Fal self.volume_attr = 'dry volume' if dry else 'volume' self.radius_bins_edges = radius_bins_edges self.normalise_by_dv = normalise_by_dv - super().__init__(name=name, unit=unit) + super().__init__(name=name, unit=unit, attr_unit='m') def register(self, builder): builder.request_attribute(self.volume_attr) diff --git a/PySDM/products/size_spectral/particles_volume_spectrum.py b/PySDM/products/size_spectral/particles_volume_spectrum.py index 13c7a56c0..362404b6b 100644 --- a/PySDM/products/size_spectral/particles_volume_spectrum.py +++ b/PySDM/products/size_spectral/particles_volume_spectrum.py @@ -7,7 +7,7 @@ class ParticlesVolumeSpectrum(SpectrumMomentProduct): def __init__(self, radius_bins_edges, name=None, unit='dimensionless'): - super().__init__(name=name, unit=unit) + super().__init__(name=name, unit=unit, attr_unit='m') self.radius_bins_edges = radius_bins_edges self.moment_0 = None self.moments = None diff --git a/PySDM/products/size_spectral/radius_binned_number_averaged_terminal_velocity.py b/PySDM/products/size_spectral/radius_binned_number_averaged_terminal_velocity.py index 18402568f..da47177f3 100644 --- a/PySDM/products/size_spectral/radius_binned_number_averaged_terminal_velocity.py +++ b/PySDM/products/size_spectral/radius_binned_number_averaged_terminal_velocity.py @@ -9,7 +9,7 @@ class RadiusBinnedNumberAveragedTerminalVelocity(SpectrumMomentProduct): def __init__(self, radius_bin_edges, name=None, unit='m/s'): - super().__init__(name=name, unit=unit) + super().__init__(name=name, unit=unit, attr_unit='m') self.radius_bin_edges = radius_bin_edges def register(self, builder): From 1631584bc31e597e9251f0d26c9e10710daffc62 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Mon, 27 Dec 2021 12:31:33 +0100 Subject: [PATCH 28/33] more super droplets for test_freezing --- tests/smoke_tests/arabas_et_al_2015/test_freezing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py index 2180a4e85..60f40b9fc 100644 --- a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py +++ b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py @@ -32,7 +32,7 @@ def test_freezing(singular): )) settings.dt = .5 * si.second settings.grid = (5, 15) - settings.n_sd_per_gridbox = 16 + settings.n_sd_per_gridbox = 64 settings.simulation_time = 100 * settings.dt settings.spin_up_time = 10 * settings.dt From b539cff1544614c098534560151de59381c3a006 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Mon, 27 Dec 2021 12:49:11 +0100 Subject: [PATCH 29/33] renaming ParticlesVolumeSpectrum -> ParticleVolumeVersusRadiusLogarithmSpectrum --- PySDM/products/size_spectral/__init__.py | 2 +- ... => particle_volume_versus_radius_logarithm_spectrum.py} | 6 ++++-- README.md | 6 +++--- tests/unit_tests/test_products.py | 4 ++-- 4 files changed, 10 insertions(+), 8 deletions(-) rename PySDM/products/size_spectral/{particles_volume_spectrum.py => particle_volume_versus_radius_logarithm_spectrum.py} (84%) diff --git a/PySDM/products/size_spectral/__init__.py b/PySDM/products/size_spectral/__init__.py index f622f398c..9753b58fa 100644 --- a/PySDM/products/size_spectral/__init__.py +++ b/PySDM/products/size_spectral/__init__.py @@ -3,7 +3,7 @@ from .mean_radius import MeanRadius from .particle_size_spectrum import ParticleSizeSpectrumPerMass, ParticleSizeSpectrumPerVolume from .particles_concentration import ParticleConcentration, ParticleSpecificConcentration -from .particles_volume_spectrum import ParticlesVolumeSpectrum +from .particle_volume_versus_radius_logarithm_spectrum import ParticleVolumeVersusRadiusLogarithmSpectrum from .total_particle_concentration import TotalParticleConcentration from .total_particle_specific_concentration import TotalParticleSpecificConcentration from .water_mixing_ratio import WaterMixingRatio diff --git a/PySDM/products/size_spectral/particles_volume_spectrum.py b/PySDM/products/size_spectral/particle_volume_versus_radius_logarithm_spectrum.py similarity index 84% rename from PySDM/products/size_spectral/particles_volume_spectrum.py rename to PySDM/products/size_spectral/particle_volume_versus_radius_logarithm_spectrum.py index 362404b6b..c9f681dc3 100644 --- a/PySDM/products/size_spectral/particles_volume_spectrum.py +++ b/PySDM/products/size_spectral/particle_volume_versus_radius_logarithm_spectrum.py @@ -1,11 +1,13 @@ """ -the dn-dlnr particle volume spectrum per volume of air (uses natural logarithm) +n_V^e(ln(r)) particle volume spectrum per volume of air (uses natural logarithm), +i.e. volume of particles per volume of air having in the size range ln(r) to +ln(r) + dln(r) """ import numpy as np from PySDM.products.impl.spectrum_moment_product import SpectrumMomentProduct -class ParticlesVolumeSpectrum(SpectrumMomentProduct): +class ParticleVolumeVersusRadiusLogarithmSpectrum(SpectrumMomentProduct): def __init__(self, radius_bins_edges, name=None, unit='dimensionless'): super().__init__(name=name, unit=unit, attr_unit='m') self.radius_bins_edges = radius_bins_edges diff --git a/README.md b/README.md index 3f2d930bb..a25bae6bd 100644 --- a/README.md +++ b/README.md @@ -213,14 +213,14 @@ from PySDM.environments import Box from PySDM.dynamics import Coalescence from PySDM.physics.coalescence_kernels import Golovin from PySDM.backends import CPU -from PySDM.products import ParticlesVolumeSpectrum +from PySDM.products import ParticleVolumeVersusRadiusLogarithmSpectrum radius_bins_edges = np.logspace(np.log10(10 * si.um), np.log10(5e3 * si.um), num=32) builder = Builder(n_sd=n_sd, backend=CPU()) -builder.set_environment(Box(dt=1 * si.s, dv=1e6 * si.m**3)) +builder.set_environment(Box(dt=1 * si.s, dv=1e6 * si.m ** 3)) builder.add_dynamic(Coalescence(kernel=Golovin(b=1.5e3 / si.s))) -products = [ParticlesVolumeSpectrum(radius_bins_edges=radius_bins_edges, name='dv/dlnr')] +products = [ParticleVolumeVersusRadiusLogarithmSpectrum(radius_bins_edges=radius_bins_edges, name='dv/dlnr')] particulator = builder.build(attributes, products) ``` diff --git a/tests/unit_tests/test_products.py b/tests/unit_tests/test_products.py index dc03476e5..4276c4b5b 100644 --- a/tests/unit_tests/test_products.py +++ b/tests/unit_tests/test_products.py @@ -12,7 +12,7 @@ from PySDM.products import (AqueousMassSpectrum, AqueousMoleFraction, TotalDryMassMixingRatio, ParticleSizeSpectrumPerMass, GaseousMoleFraction, FreezableSpecificConcentration, DynamicWallTime, - ParticleSizeSpectrumPerVolume, ParticlesVolumeSpectrum, + ParticleSizeSpectrumPerVolume, ParticleVolumeVersusRadiusLogarithmSpectrum, RadiusBinnedNumberAveragedTerminalVelocity) _ARGUMENTS = { @@ -24,7 +24,7 @@ FreezableSpecificConcentration: {'temperature_bins_edges': (0, 300)}, DynamicWallTime: {'dynamic': 'Condensation'}, ParticleSizeSpectrumPerVolume: {'radius_bins_edges': (0, np.inf)}, - ParticlesVolumeSpectrum: {'radius_bins_edges': (0, np.inf)}, + ParticleVolumeVersusRadiusLogarithmSpectrum: {'radius_bins_edges': (0, np.inf)}, RadiusBinnedNumberAveragedTerminalVelocity: {'radius_bin_edges': (0, np.inf)} } From 924c7d0961c2de144d35fe3c31e8b06a2ae4b1f4 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Mon, 27 Dec 2021 13:04:39 +0100 Subject: [PATCH 30/33] switching from exdown to pytest-codeblocks --- .github/workflows/julia.yml | 4 ++-- .github/workflows/matlab.yml | 4 ++-- .github/workflows/python.yml | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/.github/workflows/julia.yml b/.github/workflows/julia.yml index 88a0c8c0a..d508fea45 100644 --- a/.github/workflows/julia.yml +++ b/.github/workflows/julia.yml @@ -21,8 +21,8 @@ jobs: with: python-version: 3.9 - run: pip install -e . - - run: pip install exdown pytest - - run: python -c "import exdown; code=exdown.extract('README.md'); f=open('readme.jl', 'w'); f.writelines(block.code for block in code if block.syntax=='Julia'); f.close()" + - run: pip install pytest-codeblocks pytest + - run: python -c "import pytest-codeblocks; code=pytest-codeblocks.extract_from_file('README.md'); f=open('readme.jl', 'w'); f.writelines(block.code for block in code if block.syntax=='Julia'); f.close()" - uses: julia-actions/setup-julia@v1 - run: cat -n readme.jl - run: julia readme.jl diff --git a/.github/workflows/matlab.yml b/.github/workflows/matlab.yml index fa6aba883..cc0962eeb 100644 --- a/.github/workflows/matlab.yml +++ b/.github/workflows/matlab.yml @@ -21,8 +21,8 @@ jobs: with: python-version: 3.8 - run: pip install -e . - - run: pip install exdown pytest - - run: python -c "import exdown; code=exdown.extract('README.md'); f=open('readme.m', 'w'); f.writelines(block.code for block in code if block.syntax=='Matlab'); f.close()" + - run: pip install pytest-codeblocks pytest + - run: python -c "import pytest-codeblocks; code=pytest-codeblocks.extract_from_file('README.md'); f=open('readme.m', 'w'); f.writelines(block.code for block in code if block.syntax=='Matlab'); f.close()" - run: cat -n readme.m - uses: matlab-actions/setup-matlab@v0 with: diff --git a/.github/workflows/python.yml b/.github/workflows/python.yml index 38b1d969b..96f21a67c 100644 --- a/.github/workflows/python.yml +++ b/.github/workflows/python.yml @@ -21,8 +21,8 @@ jobs: with: python-version: 3.9 - run: pip install -e . - - run: pip install exdown pytest - - run: python -c "import exdown; code=exdown.extract('README.md'); f=open('readme.py', 'w'); f.writelines(block.code for block in code if block.syntax=='Python'); f.close()" + - run: pip install pytest-codeblocks pytest + - run: python -c "import pytest-codeblocks; code=pytest-codeblocks.extract_from_file('README.md'); f=open('readme.py', 'w'); f.writelines(block.code for block in code if block.syntax=='Python'); f.close()" - run: python -We readme.py - run: sed -i 's/CPU/GPU/g' readme.py - run: python -We readme.py From 022cd870bc687401448d05c5f586360cc573d908 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Mon, 27 Dec 2021 13:28:29 +0100 Subject: [PATCH 31/33] fix pytest-codeblocks import --- .github/workflows/julia.yml | 2 +- .github/workflows/matlab.yml | 2 +- .github/workflows/python.yml | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/julia.yml b/.github/workflows/julia.yml index d508fea45..c0315dc2a 100644 --- a/.github/workflows/julia.yml +++ b/.github/workflows/julia.yml @@ -22,7 +22,7 @@ jobs: python-version: 3.9 - run: pip install -e . - run: pip install pytest-codeblocks pytest - - run: python -c "import pytest-codeblocks; code=pytest-codeblocks.extract_from_file('README.md'); f=open('readme.jl', 'w'); f.writelines(block.code for block in code if block.syntax=='Julia'); f.close()" + - run: python -c "import pytest_codeblocks; code=pytest_codeblocks.extract_from_file('README.md'); f=open('readme.jl', 'w'); f.writelines(block.code for block in code if block.syntax=='Julia'); f.close()" - uses: julia-actions/setup-julia@v1 - run: cat -n readme.jl - run: julia readme.jl diff --git a/.github/workflows/matlab.yml b/.github/workflows/matlab.yml index cc0962eeb..ac68da5a8 100644 --- a/.github/workflows/matlab.yml +++ b/.github/workflows/matlab.yml @@ -22,7 +22,7 @@ jobs: python-version: 3.8 - run: pip install -e . - run: pip install pytest-codeblocks pytest - - run: python -c "import pytest-codeblocks; code=pytest-codeblocks.extract_from_file('README.md'); f=open('readme.m', 'w'); f.writelines(block.code for block in code if block.syntax=='Matlab'); f.close()" + - run: python -c "import pytest_codeblocks; code=pytest_codeblocks.extract_from_file('README.md'); f=open('readme.m', 'w'); f.writelines(block.code for block in code if block.syntax=='Matlab'); f.close()" - run: cat -n readme.m - uses: matlab-actions/setup-matlab@v0 with: diff --git a/.github/workflows/python.yml b/.github/workflows/python.yml index 96f21a67c..55ed28b41 100644 --- a/.github/workflows/python.yml +++ b/.github/workflows/python.yml @@ -22,7 +22,7 @@ jobs: python-version: 3.9 - run: pip install -e . - run: pip install pytest-codeblocks pytest - - run: python -c "import pytest-codeblocks; code=pytest-codeblocks.extract_from_file('README.md'); f=open('readme.py', 'w'); f.writelines(block.code for block in code if block.syntax=='Python'); f.close()" + - run: python -c "import pytest_codeblocks; code=pytest_codeblocks.extract_from_file('README.md'); f=open('readme.py', 'w'); f.writelines(block.code for block in code if block.syntax=='Python'); f.close()" - run: python -We readme.py - run: sed -i 's/CPU/GPU/g' readme.py - run: python -We readme.py From 770d1d8b2726c5cf0e031976b685d80eb7b83b58 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Mon, 27 Dec 2021 15:38:18 +0100 Subject: [PATCH 32/33] shorter lines --- PySDM/products/size_spectral/__init__.py | 3 ++- tests/unit_tests/test_products.py | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/PySDM/products/size_spectral/__init__.py b/PySDM/products/size_spectral/__init__.py index 9753b58fa..4b1a20110 100644 --- a/PySDM/products/size_spectral/__init__.py +++ b/PySDM/products/size_spectral/__init__.py @@ -3,7 +3,8 @@ from .mean_radius import MeanRadius from .particle_size_spectrum import ParticleSizeSpectrumPerMass, ParticleSizeSpectrumPerVolume from .particles_concentration import ParticleConcentration, ParticleSpecificConcentration -from .particle_volume_versus_radius_logarithm_spectrum import ParticleVolumeVersusRadiusLogarithmSpectrum +from .particle_volume_versus_radius_logarithm_spectrum import \ + ParticleVolumeVersusRadiusLogarithmSpectrum from .total_particle_concentration import TotalParticleConcentration from .total_particle_specific_concentration import TotalParticleSpecificConcentration from .water_mixing_ratio import WaterMixingRatio diff --git a/tests/unit_tests/test_products.py b/tests/unit_tests/test_products.py index 4276c4b5b..357837491 100644 --- a/tests/unit_tests/test_products.py +++ b/tests/unit_tests/test_products.py @@ -12,7 +12,8 @@ from PySDM.products import (AqueousMassSpectrum, AqueousMoleFraction, TotalDryMassMixingRatio, ParticleSizeSpectrumPerMass, GaseousMoleFraction, FreezableSpecificConcentration, DynamicWallTime, - ParticleSizeSpectrumPerVolume, ParticleVolumeVersusRadiusLogarithmSpectrum, + ParticleSizeSpectrumPerVolume, + ParticleVolumeVersusRadiusLogarithmSpectrum, RadiusBinnedNumberAveragedTerminalVelocity) _ARGUMENTS = { From 6901491a7469a8d75661b9a54d68ba39b668d323 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Mon, 27 Dec 2021 22:09:47 +0100 Subject: [PATCH 33/33] ParticlesVolumeSpectrum -> ParticleVolumeVersusRadiusLogarithmSpectrum in README --- README.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index a25bae6bd..1be872931 100644 --- a/README.md +++ b/README.md @@ -169,14 +169,14 @@ Box = pyimport("PySDM.environments").Box Coalescence = pyimport("PySDM.dynamics").Coalescence Golovin = pyimport("PySDM.physics.coalescence_kernels").Golovin CPU = pyimport("PySDM.backends").CPU -ParticlesVolumeSpectrum = pyimport("PySDM.products").ParticlesVolumeSpectrum +ParticleVolumeVersusRadiusLogarithmSpectrum = pyimport("PySDM.products").ParticleVolumeVersusRadiusLogarithmSpectrum radius_bins_edges = 10 .^ range(log10(10*si.um), log10(5e3*si.um), length=32) builder = Builder(n_sd=n_sd, backend=CPU()) builder.set_environment(Box(dt=1 * si.s, dv=1e6 * si.m^3)) builder.add_dynamic(Coalescence(kernel=Golovin(b=1.5e3 / si.s))) -products = [ParticlesVolumeSpectrum(radius_bins_edges=radius_bins_edges, name="dv/dlnr")] +products = [ParticleVolumeVersusRadiusLogarithmSpectrum(radius_bins_edges=radius_bins_edges, name="dv/dlnr")] particulator = builder.build(attributes, products) ``` @@ -189,14 +189,14 @@ Box = py.importlib.import_module('PySDM.environments').Box; Coalescence = py.importlib.import_module('PySDM.dynamics').Coalescence; Golovin = py.importlib.import_module('PySDM.physics.coalescence_kernels').Golovin; CPU = py.importlib.import_module('PySDM.backends').CPU; -ParticlesVolumeSpectrum = py.importlib.import_module('PySDM.products').ParticlesVolumeSpectrum; +ParticleVolumeVersusRadiusLogarithmSpectrum = py.importlib.import_module('PySDM.products').ParticleVolumeVersusRadiusLogarithmSpectrum; radius_bins_edges = logspace(log10(10 * si.um), log10(5e3 * si.um), 32); builder = Builder(pyargs('n_sd', int32(n_sd), 'backend', CPU())); builder.set_environment(Box(pyargs('dt', 1 * si.s, 'dv', 1e6 * si.m ^ 3))); builder.add_dynamic(Coalescence(pyargs('kernel', Golovin(1.5e3 / si.s)))); -products = py.list({ ParticlesVolumeSpectrum(pyargs( ... +products = py.list({ ParticleVolumeVersusRadiusLogarithmSpectrum(pyargs( ... 'radius_bins_edges', py.numpy.array(radius_bins_edges), ... 'name', 'dv/dlnr' ... )) });