From 09fad83d37e5b1225a74000b0d46f93eb89dd1c3 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Tue, 30 Sep 2025 04:37:05 +0200 Subject: [PATCH 1/3] early draft of deterministic deposition removal --- PySDM/dynamics/__init__.py | 1 + PySDM/dynamics/deposition_removal.py | 26 +++++++++ PySDM/particulator.py | 34 +++++++++++ docs/bibliography.json | 7 +++ .../dynamics/test_deposition_removal.py | 56 +++++++++++++++++++ 5 files changed, 124 insertions(+) create mode 100644 PySDM/dynamics/deposition_removal.py create mode 100644 tests/unit_tests/dynamics/test_deposition_removal.py diff --git a/PySDM/dynamics/__init__.py b/PySDM/dynamics/__init__.py index d296bf4720..7a1fb1c354 100644 --- a/PySDM/dynamics/__init__.py +++ b/PySDM/dynamics/__init__.py @@ -17,3 +17,4 @@ from PySDM.dynamics.relaxed_velocity import RelaxedVelocity from PySDM.dynamics.seeding import Seeding from PySDM.dynamics.vapour_deposition_on_ice import VapourDepositionOnIce +from PySDM.dynamics.deposition_removal import DepositionRemoval diff --git a/PySDM/dynamics/deposition_removal.py b/PySDM/dynamics/deposition_removal.py new file mode 100644 index 0000000000..3b1a582734 --- /dev/null +++ b/PySDM/dynamics/deposition_removal.py @@ -0,0 +1,26 @@ +"""deposition removal logic for zero-dimensional environments""" + +from PySDM.dynamics.impl import register_dynamic +import numpy as np + + +@register_dynamic() +class DepositionRemoval: + def __init__(self, *, all_or_nothing: bool): + """stochastic ("all or nothing") or deterministic (multiplicity altering) removal""" + self.all_or_nothing = all_or_nothing + self.length_scale = None + + def register(self, builder): + builder.request_attribute("relative fall velocity") + assert builder.particulator.environment.mesh.n_dims == 0 + self.particulator = builder.particulator + self.length_scale = np.cbrt(self.particulator.environment.mesh.dv) + + def __call__(self): + """see, e.g., the naive scheme in Algorithm 1 in + [Curtis et al. 2016](https://doi.org/10.1016/j.jcp.2016.06.029)""" + self.particulator.deposition_removal( + all_or_nothing=self.all_or_nothing, + length_scale=self.length_scale, + ) diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 5919feb805..710f889d66 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -565,3 +565,37 @@ def thaw_instantaneous(self): cell=self.attributes["cell id"], temperature=self.environment["T"], ) + + def deposition_removal(self, *, all_or_nothing, length_scale): + # TODO: to be moved into backend (and attributes?) + + prob_zero_events = self.formulae.trivia.poissonian_avoidance_function + + def backend_deposition_removal_all_or_nothing( + *, relative_fall_velocity, length_scale, timestep, u01 + ): + pass + + def backend_deposition_removal_deterministic( + *, relative_fall_velocity, multiplicity, length_scale, timestep + ): + for i, velocity in enumerate(relative_fall_velocity): + removal_rate = velocity / length_scale + survive_prob = prob_zero_events(r=removal_rate, dt=timestep) + assert 0 <= survive_prob <= 1 + multiplicity[i] *= survive_prob + + if all_or_nothing: + backend_deposition_removal_all_or_nothing( + relative_fall_velocity=self.attributes["relative fall velocity"].data, + length_scale=length_scale, + timestep=self.dt, + u01=None, + ) + else: + backend_deposition_removal_deterministic( + relative_fall_velocity=self.attributes["relative fall velocity"].data, + multiplicity=self.attributes["multiplicity"].data, + length_scale=length_scale, + timestep=self.dt, + ) diff --git a/docs/bibliography.json b/docs/bibliography.json index deb122ac9e..0aad561a4c 100644 --- a/docs/bibliography.json +++ b/docs/bibliography.json @@ -973,5 +973,12 @@ ], "label": "Barthazy and Schefold 2006 (Atmos. Res. 82)", "title": "Fall velocity of snowflakes of different riming degree and crystal types" + }, + "https://doi.org/10.1016/j.jcp.2016.06.029": { + "usages": [ + "PySDM/dynamics/deposition_removal.py" + ], + "label": "Curtis et al. 2016 (J. Comp. Phys. 322)", + "title": "Accelerated simulation of stochastic particle removal\nprocesses in particle-resolved aerosol models" } } diff --git a/tests/unit_tests/dynamics/test_deposition_removal.py b/tests/unit_tests/dynamics/test_deposition_removal.py new file mode 100644 index 0000000000..58d127c027 --- /dev/null +++ b/tests/unit_tests/dynamics/test_deposition_removal.py @@ -0,0 +1,56 @@ +from PySDM.dynamics import DepositionRemoval +from PySDM.physics import si +from PySDM.environments import Box +from PySDM.builder import Builder +from PySDM.backends import CPU +from PySDM.products import ParticleConcentration, SuperDropletCountPerGridbox, Time +from matplotlib import pyplot +import pytest +import numpy as np + + +class TestDepositionRemoval: + @staticmethod + @pytest.mark.parametrize("all_or_nothing", (True, False)) + def test_convergence_wrt_dt(all_or_nothing, plot=True): + # arrange + dt = 1 * si.s + dv = 666 * si.m**3 + n_steps = 100 + multiplicity = [1, 10, 100, 1000] + water_mass = [1 * si.ug, 2 * si.ug, 3 * si.ug, 4 * si.ug] + backend_instance = CPU() + + builder = Builder( + n_sd=len(multiplicity), + environment=Box(dv=dv, dt=dt), + backend=backend_instance, + ) + builder.add_dynamic(DepositionRemoval(all_or_nothing=all_or_nothing)) + particulator = builder.build( + attributes={ + "multiplicity": np.asarray(multiplicity), + "signed water mass": np.asarray(water_mass), + }, + products=(ParticleConcentration(), SuperDropletCountPerGridbox(), Time()), + ) + + # act + output = {name: [] for name in particulator.products} + for step in range(n_steps): + if step != 0: + particulator.run(steps=1) + for name, product in particulator.products.items(): + output[name].append(product.get() + 0) + + # plot + pyplot.title(f"{all_or_nothing=}") + pyplot.xlabel("time [s]") + pyplot.ylabel("particle concentration [m$^{-3}$]") + pyplot.plot(output["time"], output["particle concentration"]) + if plot: + pyplot.show() + else: + pyplot.clf() + + # assert From 6f9a3045750a2672f15fc11c03c58a42574e68a8 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Tue, 30 Sep 2025 07:08:02 +0200 Subject: [PATCH 2/3] renaming deposition -> sedimentation; plot=False not to hang CI --- PySDM/dynamics/__init__.py | 2 +- ...eposition_removal.py => sedimentation_removal.py} | 4 ++-- PySDM/particulator.py | 10 +++++----- ...tion_removal.py => test_sedimentation_removal.py} | 12 +++++++----- 4 files changed, 15 insertions(+), 13 deletions(-) rename PySDM/dynamics/{deposition_removal.py => sedimentation_removal.py} (92%) rename tests/unit_tests/dynamics/{test_deposition_removal.py => test_sedimentation_removal.py} (83%) diff --git a/PySDM/dynamics/__init__.py b/PySDM/dynamics/__init__.py index 7a1fb1c354..45725a2dac 100644 --- a/PySDM/dynamics/__init__.py +++ b/PySDM/dynamics/__init__.py @@ -17,4 +17,4 @@ from PySDM.dynamics.relaxed_velocity import RelaxedVelocity from PySDM.dynamics.seeding import Seeding from PySDM.dynamics.vapour_deposition_on_ice import VapourDepositionOnIce -from PySDM.dynamics.deposition_removal import DepositionRemoval +from PySDM.dynamics.sedimentation_removal import SedimentationRemoval diff --git a/PySDM/dynamics/deposition_removal.py b/PySDM/dynamics/sedimentation_removal.py similarity index 92% rename from PySDM/dynamics/deposition_removal.py rename to PySDM/dynamics/sedimentation_removal.py index 3b1a582734..bdb3414f7e 100644 --- a/PySDM/dynamics/deposition_removal.py +++ b/PySDM/dynamics/sedimentation_removal.py @@ -5,7 +5,7 @@ @register_dynamic() -class DepositionRemoval: +class SedimentationRemoval: def __init__(self, *, all_or_nothing: bool): """stochastic ("all or nothing") or deterministic (multiplicity altering) removal""" self.all_or_nothing = all_or_nothing @@ -20,7 +20,7 @@ def register(self, builder): def __call__(self): """see, e.g., the naive scheme in Algorithm 1 in [Curtis et al. 2016](https://doi.org/10.1016/j.jcp.2016.06.029)""" - self.particulator.deposition_removal( + self.particulator.sedimentation_removal( all_or_nothing=self.all_or_nothing, length_scale=self.length_scale, ) diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 710f889d66..fbf89de93a 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -566,17 +566,17 @@ def thaw_instantaneous(self): temperature=self.environment["T"], ) - def deposition_removal(self, *, all_or_nothing, length_scale): + def sedimentation_removal(self, *, all_or_nothing, length_scale): # TODO: to be moved into backend (and attributes?) prob_zero_events = self.formulae.trivia.poissonian_avoidance_function - def backend_deposition_removal_all_or_nothing( + def backend_sedimentation_removal_all_or_nothing( *, relative_fall_velocity, length_scale, timestep, u01 ): pass - def backend_deposition_removal_deterministic( + def backend_sedimentation_removal_deterministic( *, relative_fall_velocity, multiplicity, length_scale, timestep ): for i, velocity in enumerate(relative_fall_velocity): @@ -586,14 +586,14 @@ def backend_deposition_removal_deterministic( multiplicity[i] *= survive_prob if all_or_nothing: - backend_deposition_removal_all_or_nothing( + backend_sedimentation_removal_all_or_nothing( relative_fall_velocity=self.attributes["relative fall velocity"].data, length_scale=length_scale, timestep=self.dt, u01=None, ) else: - backend_deposition_removal_deterministic( + backend_sedimentation_removal_deterministic( relative_fall_velocity=self.attributes["relative fall velocity"].data, multiplicity=self.attributes["multiplicity"].data, length_scale=length_scale, diff --git a/tests/unit_tests/dynamics/test_deposition_removal.py b/tests/unit_tests/dynamics/test_sedimentation_removal.py similarity index 83% rename from tests/unit_tests/dynamics/test_deposition_removal.py rename to tests/unit_tests/dynamics/test_sedimentation_removal.py index 58d127c027..6bd32083f8 100644 --- a/tests/unit_tests/dynamics/test_deposition_removal.py +++ b/tests/unit_tests/dynamics/test_sedimentation_removal.py @@ -1,4 +1,4 @@ -from PySDM.dynamics import DepositionRemoval +from PySDM.dynamics import SedimentationRemoval from PySDM.physics import si from PySDM.environments import Box from PySDM.builder import Builder @@ -9,10 +9,10 @@ import numpy as np -class TestDepositionRemoval: +class TestSedimentationRemoval: @staticmethod @pytest.mark.parametrize("all_or_nothing", (True, False)) - def test_convergence_wrt_dt(all_or_nothing, plot=True): + def test_convergence_wrt_dt(all_or_nothing, plot=False): # arrange dt = 1 * si.s dv = 666 * si.m**3 @@ -26,7 +26,7 @@ def test_convergence_wrt_dt(all_or_nothing, plot=True): environment=Box(dv=dv, dt=dt), backend=backend_instance, ) - builder.add_dynamic(DepositionRemoval(all_or_nothing=all_or_nothing)) + builder.add_dynamic(SedimentationRemoval(all_or_nothing=all_or_nothing)) particulator = builder.build( attributes={ "multiplicity": np.asarray(multiplicity), @@ -47,10 +47,12 @@ def test_convergence_wrt_dt(all_or_nothing, plot=True): pyplot.title(f"{all_or_nothing=}") pyplot.xlabel("time [s]") pyplot.ylabel("particle concentration [m$^{-3}$]") - pyplot.plot(output["time"], output["particle concentration"]) + pyplot.semilogy(output["time"], output["particle concentration"]) + if plot: pyplot.show() else: pyplot.clf() # assert + # TODO! From f05ae89753586b8e1c9ad73cde262d6375b7f9da Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sat, 8 Nov 2025 15:34:46 +0100 Subject: [PATCH 3/3] extended test cases --- .../dynamics/test_sedimentation_removal.py | 75 ++++++++++++------- 1 file changed, 47 insertions(+), 28 deletions(-) diff --git a/tests/unit_tests/dynamics/test_sedimentation_removal.py b/tests/unit_tests/dynamics/test_sedimentation_removal.py index 6bd32083f8..a915f6bf39 100644 --- a/tests/unit_tests/dynamics/test_sedimentation_removal.py +++ b/tests/unit_tests/dynamics/test_sedimentation_removal.py @@ -11,43 +11,62 @@ class TestSedimentationRemoval: @staticmethod - @pytest.mark.parametrize("all_or_nothing", (True, False)) - def test_convergence_wrt_dt(all_or_nothing, plot=False): + @pytest.mark.parametrize("all_or_nothing", (False,)) + def test_convergence_wrt_dt(all_or_nothing, plot=True): # arrange - dt = 1 * si.s - dv = 666 * si.m**3 - n_steps = 100 - multiplicity = [1, 10, 100, 1000] - water_mass = [1 * si.ug, 2 * si.ug, 3 * si.ug, 4 * si.ug] + dts = 5 * si.s, 0.5 * si.s + dvs = 1e2 * si.m**3, 1e3 * si.m**3, 1e4 * si.m**3 + t_max = 500 * si.s + multiplicities = 1e5, 1e6, 1e7, 1e8 + water_masses = 1 * si.ug, 2 * si.ug, 3 * si.ug, 4 * si.ug backend_instance = CPU() - builder = Builder( - n_sd=len(multiplicity), - environment=Box(dv=dv, dt=dt), - backend=backend_instance, - ) - builder.add_dynamic(SedimentationRemoval(all_or_nothing=all_or_nothing)) - particulator = builder.build( - attributes={ - "multiplicity": np.asarray(multiplicity), - "signed water mass": np.asarray(water_mass), - }, - products=(ParticleConcentration(), SuperDropletCountPerGridbox(), Time()), - ) - # act - output = {name: [] for name in particulator.products} - for step in range(n_steps): - if step != 0: - particulator.run(steps=1) - for name, product in particulator.products.items(): - output[name].append(product.get() + 0) + output = {} + for dt in dts: + for dv in dvs: + builder = Builder( + n_sd=len(multiplicities), + environment=Box(dv=dv, dt=dt), + backend=backend_instance, + ) + builder.add_dynamic(SedimentationRemoval(all_or_nothing=all_or_nothing)) + particulator = builder.build( + attributes={ + "multiplicity": np.asarray(multiplicities), + "signed water mass": np.asarray(water_masses), + }, + products=( + ParticleConcentration(), + SuperDropletCountPerGridbox(), + Time(), + ), + ) + key = f"{dt=} {dv=}" + output[key] = {name: [] for name in particulator.products} + while particulator.n_steps * dt <= t_max: + if len(output[key]["time"]) != 0: + particulator.run(steps=1) + for name, product in particulator.products.items(): + output[key][name].append(product.get() + 0) # plot pyplot.title(f"{all_or_nothing=}") pyplot.xlabel("time [s]") pyplot.ylabel("particle concentration [m$^{-3}$]") - pyplot.semilogy(output["time"], output["particle concentration"]) + for dt in dts: + for dv in dvs: + key = f"{dt=} {dv=}" + pyplot.plot( + output[key]["time"], + output[key]["particle concentration"], + label=key, + linewidth=4 + 3 * np.log10(dt), + ) + pyplot.gca().set_yscale("log") + pyplot.gca().set_xlim(left=0, right=t_max) + pyplot.legend() + pyplot.grid() if plot: pyplot.show()