diff --git a/PySDM/dynamics/__init__.py b/PySDM/dynamics/__init__.py index d296bf472..45725a2da 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.sedimentation_removal import SedimentationRemoval diff --git a/PySDM/dynamics/sedimentation_removal.py b/PySDM/dynamics/sedimentation_removal.py new file mode 100644 index 000000000..bdb3414f7 --- /dev/null +++ b/PySDM/dynamics/sedimentation_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 SedimentationRemoval: + 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.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 5919feb80..fbf89de93 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 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_sedimentation_removal_all_or_nothing( + *, relative_fall_velocity, length_scale, timestep, u01 + ): + pass + + def backend_sedimentation_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_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_sedimentation_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 deb122ac9..0aad561a4 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_sedimentation_removal.py b/tests/unit_tests/dynamics/test_sedimentation_removal.py new file mode 100644 index 000000000..6bd32083f --- /dev/null +++ b/tests/unit_tests/dynamics/test_sedimentation_removal.py @@ -0,0 +1,58 @@ +from PySDM.dynamics import SedimentationRemoval +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 TestSedimentationRemoval: + @staticmethod + @pytest.mark.parametrize("all_or_nothing", (True, False)) + def test_convergence_wrt_dt(all_or_nothing, plot=False): + # 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(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) + + # plot + pyplot.title(f"{all_or_nothing=}") + pyplot.xlabel("time [s]") + pyplot.ylabel("particle concentration [m$^{-3}$]") + pyplot.semilogy(output["time"], output["particle concentration"]) + + if plot: + pyplot.show() + else: + pyplot.clf() + + # assert + # TODO!