diff --git a/PySDM/backends/impl_numba/methods/__init__.py b/PySDM/backends/impl_numba/methods/__init__.py
index ca35c0e652..73d7ceed33 100644
--- a/PySDM/backends/impl_numba/methods/__init__.py
+++ b/PySDM/backends/impl_numba/methods/__init__.py
@@ -12,3 +12,4 @@
from .pair_methods import PairMethods
from .physics_methods import PhysicsMethods
from .terminal_velocity_methods import TerminalVelocityMethods
+from .seeding_methods import SeedingMethods
diff --git a/PySDM/backends/impl_numba/methods/condensation_methods.py b/PySDM/backends/impl_numba/methods/condensation_methods.py
index 59d3440bbd..5c7f71e94d 100644
--- a/PySDM/backends/impl_numba/methods/condensation_methods.py
+++ b/PySDM/backends/impl_numba/methods/condensation_methods.py
@@ -445,7 +445,7 @@ def calculate_ml_new( # pylint: disable=too-many-branches,too-many-arguments,to
v_drop = formulae.particle_shape_and_density__mass_to_volume(
attributes.water_mass[drop]
)
- if v_drop < 0:
+ if v_drop <= 0:
continue
x_old = formulae.condensation_coordinate__x(v_drop)
r_old = formulae.trivia__radius(v_drop)
diff --git a/PySDM/backends/impl_numba/methods/seeding_methods.py b/PySDM/backends/impl_numba/methods/seeding_methods.py
new file mode 100644
index 0000000000..a1e6372ba6
--- /dev/null
+++ b/PySDM/backends/impl_numba/methods/seeding_methods.py
@@ -0,0 +1,68 @@
+""" CPU implementation of backend methods for particle injections """
+
+from functools import cached_property
+
+import numba
+
+from PySDM.backends.impl_common.backend_methods import BackendMethods
+
+
+class SeedingMethods(BackendMethods): # pylint: disable=too-few-public-methods
+ @cached_property
+ def _seeding(self):
+ @numba.njit(**{**self.default_jit_flags, "parallel": False})
+ def body( # pylint: disable=too-many-arguments
+ idx,
+ multiplicity,
+ extensive_attributes,
+ seeded_particle_index,
+ seeded_particle_multiplicity,
+ seeded_particle_extensive_attributes,
+ number_of_super_particles_to_inject: int,
+ ):
+ number_of_super_particles_already_injected = 0
+ # TODO #1387 start enumerating from the end of valid particle set
+ for i, mult in enumerate(multiplicity):
+ if (
+ number_of_super_particles_to_inject
+ == number_of_super_particles_already_injected
+ ):
+ break
+ if mult == 0:
+ idx[i] = -1
+ s = seeded_particle_index[
+ number_of_super_particles_already_injected
+ ]
+ number_of_super_particles_already_injected += 1
+ multiplicity[i] = seeded_particle_multiplicity[s]
+ for a in range(len(extensive_attributes)):
+ extensive_attributes[a, i] = (
+ seeded_particle_extensive_attributes[a, s]
+ )
+ assert (
+ number_of_super_particles_to_inject
+ == number_of_super_particles_already_injected
+ )
+
+ return body
+
+ def seeding(
+ self,
+ *,
+ idx,
+ multiplicity,
+ extensive_attributes,
+ seeded_particle_index,
+ seeded_particle_multiplicity,
+ seeded_particle_extensive_attributes,
+ number_of_super_particles_to_inject: int,
+ ):
+ self._seeding(
+ idx=idx.data,
+ multiplicity=multiplicity.data,
+ extensive_attributes=extensive_attributes.data,
+ seeded_particle_index=seeded_particle_index.data,
+ seeded_particle_multiplicity=seeded_particle_multiplicity.data,
+ seeded_particle_extensive_attributes=seeded_particle_extensive_attributes.data,
+ number_of_super_particles_to_inject=number_of_super_particles_to_inject,
+ )
diff --git a/PySDM/backends/numba.py b/PySDM/backends/numba.py
index 69833b161c..a3b558d112 100644
--- a/PySDM/backends/numba.py
+++ b/PySDM/backends/numba.py
@@ -28,6 +28,7 @@ class Numba( # pylint: disable=too-many-ancestors,duplicate-code
methods.DisplacementMethods,
methods.TerminalVelocityMethods,
methods.IsotopeMethods,
+ methods.SeedingMethods,
):
Storage = ImportedStorage
Random = ImportedRandom
@@ -75,3 +76,4 @@ def __init__(self, formulae=None, double_precision=True, override_jit_flags=None
methods.DisplacementMethods.__init__(self)
methods.TerminalVelocityMethods.__init__(self)
methods.IsotopeMethods.__init__(self)
+ methods.SeedingMethods.__init__(self)
diff --git a/PySDM/builder.py b/PySDM/builder.py
index 4ab070ad6b..81dad5e1a6 100644
--- a/PySDM/builder.py
+++ b/PySDM/builder.py
@@ -150,4 +150,8 @@ def build(
for key in self.particulator.dynamics:
self.particulator.timers[key] = WallTimer()
+ if (attributes["multiplicity"] == 0).any():
+ self.particulator.attributes.healthy = False
+ self.particulator.attributes.sanitize()
+
return self.particulator
diff --git a/PySDM/dynamics/__init__.py b/PySDM/dynamics/__init__.py
index cf2def8da4..66056cf0d4 100644
--- a/PySDM/dynamics/__init__.py
+++ b/PySDM/dynamics/__init__.py
@@ -15,3 +15,4 @@
from PySDM.dynamics.eulerian_advection import EulerianAdvection
from PySDM.dynamics.freezing import Freezing
from PySDM.dynamics.relaxed_velocity import RelaxedVelocity
+from PySDM.dynamics.seeding import Seeding
diff --git a/PySDM/dynamics/seeding.py b/PySDM/dynamics/seeding.py
new file mode 100644
index 0000000000..3655e62a66
--- /dev/null
+++ b/PySDM/dynamics/seeding.py
@@ -0,0 +1,94 @@
+""" particle injection handling, requires initalising a simulation with
+enough particles flagged with NaN multiplicity (translated to zeros
+at multiplicity discretisation """
+
+from collections.abc import Sized
+
+import numpy as np
+
+from PySDM.dynamics.impl import register_dynamic
+from PySDM.initialisation import discretise_multiplicities
+
+
+@register_dynamic()
+class Seeding:
+ def __init__(
+ self,
+ *,
+ super_droplet_injection_rate: callable,
+ seeded_particle_extensive_attributes: dict,
+ seeded_particle_multiplicity: Sized,
+ ):
+ for attr in seeded_particle_extensive_attributes.values():
+ assert len(seeded_particle_multiplicity) == len(attr)
+ self.particulator = None
+ self.super_droplet_injection_rate = super_droplet_injection_rate
+ self.seeded_particle_extensive_attributes = seeded_particle_extensive_attributes
+ self.seeded_particle_multiplicity = seeded_particle_multiplicity
+ self.rnd = None
+ self.u01 = None
+ self.index = None
+
+ def register(self, builder):
+ self.particulator = builder.particulator
+
+ def post_register_setup_when_attributes_are_known(self):
+ if tuple(self.particulator.attributes.get_extensive_attribute_keys()) != tuple(
+ self.seeded_particle_extensive_attributes.keys()
+ ):
+ raise ValueError(
+ f"extensive attributes ({self.seeded_particle_extensive_attributes.keys()})"
+ " do not match those used in particulator"
+ f" ({self.particulator.attributes.get_extensive_attribute_keys()})"
+ )
+
+ self.index = self.particulator.Index.identity_index(
+ len(self.seeded_particle_multiplicity)
+ )
+ if len(self.seeded_particle_multiplicity) > 1:
+ self.rnd = self.particulator.Random(
+ len(self.seeded_particle_multiplicity), self.particulator.formulae.seed
+ )
+ self.u01 = self.particulator.Storage.empty(
+ len(self.seeded_particle_multiplicity), dtype=float
+ )
+ self.seeded_particle_multiplicity = (
+ self.particulator.IndexedStorage.from_ndarray(
+ self.index,
+ discretise_multiplicities(
+ np.asarray(self.seeded_particle_multiplicity)
+ ),
+ )
+ )
+ self.seeded_particle_extensive_attributes = (
+ self.particulator.IndexedStorage.from_ndarray(
+ self.index,
+ np.asarray(list(self.seeded_particle_extensive_attributes.values())),
+ )
+ )
+
+ def __call__(self):
+ if self.particulator.n_steps == 0:
+ self.post_register_setup_when_attributes_are_known()
+
+ time = self.particulator.n_steps * self.particulator.dt
+ number_of_super_particles_to_inject = self.super_droplet_injection_rate(time)
+
+ if number_of_super_particles_to_inject > 0:
+ assert number_of_super_particles_to_inject <= len(
+ self.seeded_particle_multiplicity
+ )
+
+ if self.rnd is not None:
+ self.u01.urand(self.rnd)
+ # TODO #1387 make shuffle smarter
+ # e.g. don't need to shuffle if only one type of seed particle
+ # or if the number of super particles to inject
+ # is equal to the number of possible seeds
+ self.index.shuffle(self.u01)
+ self.particulator.seeding(
+ seeded_particle_index=self.index,
+ number_of_super_particles_to_inject=number_of_super_particles_to_inject,
+ seeded_particle_multiplicity=self.seeded_particle_multiplicity,
+ seeded_particle_extensive_attributes=self.seeded_particle_extensive_attributes,
+ )
diff --git a/PySDM/impl/particle_attributes.py b/PySDM/impl/particle_attributes.py
index f12bb6fabc..5e486db29f 100644
--- a/PySDM/impl/particle_attributes.py
+++ b/PySDM/impl/particle_attributes.py
@@ -118,3 +118,8 @@ def get_extensive_attribute_keys(self):
def has_attribute(self, attr):
return attr in self.__attributes
+
+ def reset_idx(self):
+ self.__valid_n_sd = self.__idx.shape[0]
+ self.__idx.reset_index()
+ self.healthy = False
diff --git a/PySDM/initialisation/discretise_multiplicities.py b/PySDM/initialisation/discretise_multiplicities.py
index b085d6e438..1daf978cf0 100644
--- a/PySDM/initialisation/discretise_multiplicities.py
+++ b/PySDM/initialisation/discretise_multiplicities.py
@@ -6,11 +6,22 @@
def discretise_multiplicities(values_arg):
- values_int = values_arg.round().astype(np.int64)
+ """any NaN values in the input array are ignored and flagged
+ with zero multiplicities in the output array"""
+
+ values_int = np.where(np.isnan(values_arg), 0, values_arg).round().astype(np.int64)
if np.issubdtype(values_arg.dtype, np.floating):
+ if np.isnan(values_arg).all():
+ return values_int
+
+ if not np.logical_or(values_int > 0, np.isnan(values_arg)).all():
+ raise ValueError(
+ f"int-casting resulted in multiplicity of zero (min(y_float)={min(values_arg)})"
+ )
+
percent_diff = 100 * abs(
- 1 - np.sum(values_arg) / np.sum(values_int.astype(float))
+ 1 - np.nansum(values_arg) / np.sum(values_int.astype(float))
)
if percent_diff > 1:
raise ValueError(
@@ -18,9 +29,4 @@ def discretise_multiplicities(values_arg):
f" due to casting multiplicities to ints"
)
- if not (values_int > 0).all():
- raise ValueError(
- f"int-casting resulted in multiplicity of zero (min(y_float)={min(values_arg)})"
- )
-
return values_int
diff --git a/PySDM/initialisation/init_fall_momenta.py b/PySDM/initialisation/init_fall_momenta.py
index 57342e2cbd..36f523894a 100644
--- a/PySDM/initialisation/init_fall_momenta.py
+++ b/PySDM/initialisation/init_fall_momenta.py
@@ -5,7 +5,6 @@
import numpy as np
-from PySDM.backends import CPU
from PySDM.dynamics.terminal_velocity import GunnKinzer1949
from PySDM.formulae import Formulae
from PySDM.particulator import Particulator
@@ -31,6 +30,8 @@ def init_fall_momenta(
if zero:
return np.zeros_like(water_mass)
+ from PySDM.backends import CPU # pylint: disable=import-outside-toplevel
+
particulator = Particulator(0, CPU(Formulae())) # TODO #1155
approximation = terminal_velocity_approx(particulator=particulator)
diff --git a/PySDM/particulator.py b/PySDM/particulator.py
index 08e933d813..7ab1098ce6 100644
--- a/PySDM/particulator.py
+++ b/PySDM/particulator.py
@@ -13,7 +13,7 @@
class Particulator: # pylint: disable=too-many-public-methods,too-many-instance-attributes
- def __init__(self, n_sd, backend: BackendMethods):
+ def __init__(self, n_sd, backend):
assert isinstance(backend, BackendMethods)
self.__n_sd = n_sd
@@ -438,3 +438,44 @@ def isotopic_fractionation(self, heavy_isotopes: tuple):
self.backend.isotopic_fractionation()
for isotope in heavy_isotopes:
self.attributes.mark_updated(f"moles_{isotope}")
+
+ def seeding(
+ self,
+ *,
+ seeded_particle_index,
+ seeded_particle_multiplicity,
+ seeded_particle_extensive_attributes,
+ number_of_super_particles_to_inject,
+ ):
+ n_null = self.n_sd - self.attributes.super_droplet_count
+ if n_null == 0:
+ raise ValueError(
+ "No available seeds to inject. Please provide particles with nan filled attributes."
+ )
+
+ if number_of_super_particles_to_inject > n_null:
+ raise ValueError(
+ "Trying to inject more super particles than space available."
+ )
+
+ if number_of_super_particles_to_inject > len(seeded_particle_multiplicity):
+ raise ValueError(
+ "Trying to inject multiple super particles with the same attributes. \
+ Instead increase multiplicity of injected particles."
+ )
+
+ self.backend.seeding(
+ idx=self.attributes._ParticleAttributes__idx,
+ multiplicity=self.attributes["multiplicity"],
+ extensive_attributes=self.attributes.get_extensive_attribute_storage(),
+ seeded_particle_index=seeded_particle_index,
+ seeded_particle_multiplicity=seeded_particle_multiplicity,
+ seeded_particle_extensive_attributes=seeded_particle_extensive_attributes,
+ number_of_super_particles_to_inject=number_of_super_particles_to_inject,
+ )
+ self.attributes.reset_idx()
+ self.attributes.sanitize()
+
+ self.attributes.mark_updated("multiplicity")
+ for key in self.attributes.get_extensive_attribute_keys():
+ self.attributes.mark_updated(key)
diff --git a/examples/PySDM_examples/seeding/__init__.py b/examples/PySDM_examples/seeding/__init__.py
new file mode 100644
index 0000000000..e4cb5d00e7
--- /dev/null
+++ b/examples/PySDM_examples/seeding/__init__.py
@@ -0,0 +1,2 @@
+from .settings import Settings
+from .simulation import Simulation
diff --git a/examples/PySDM_examples/seeding/hello_world.ipynb b/examples/PySDM_examples/seeding/hello_world.ipynb
new file mode 100644
index 0000000000..9f868178a4
--- /dev/null
+++ b/examples/PySDM_examples/seeding/hello_world.ipynb
@@ -0,0 +1,221 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "b1c0a002-7f76-4a4c-9e77-6aae747d7fab",
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2024-08-20T22:17:43.753049Z",
+ "start_time": "2024-08-20T22:17:40.315223Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "from matplotlib import pyplot\n",
+ "from PySDM import Formulae\n",
+ "from PySDM.physics import in_unit, si\n",
+ "from open_atmos_jupyter_utils import show_plot\n",
+ "from PySDM_examples.seeding import Settings, Simulation"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "1695dd83-b78b-46f3-a50c-47cc601c5920",
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2024-08-20T22:18:08.299964Z",
+ "start_time": "2024-08-20T22:17:43.765032Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "n_sd_initial = 100\n",
+ "n_sd_seeding = 100\n",
+ "rain_water_radius_threshold = .1 * si.mm\n",
+ "formulae = Formulae(seed=100)\n",
+ "\n",
+ "simulations = {\n",
+ " case: Simulation(\n",
+ " Settings(\n",
+ " n_sd_initial=n_sd_initial,\n",
+ " n_sd_seeding=n_sd_seeding,\n",
+ " rain_water_radius_threshold=rain_water_radius_threshold,\n",
+ " super_droplet_injection_rate={\n",
+ " 'seeding': lambda time: 1 if 5 * si.min < time < 10 * si.min else 0,\n",
+ " 'no seeding': lambda _: 0,\n",
+ " }[case],\n",
+ " formulae=formulae,\n",
+ " )\n",
+ " )\n",
+ " for case in ('seeding', 'no seeding')\n",
+ "} "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "af543973-a177-4fe9-9ae5-b0fb8511323e",
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2024-08-20T22:19:04.834137Z",
+ "start_time": "2024-08-20T22:18:08.367114Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "outputs = {case: simulations[case].run() for case in simulations}"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "45bbaff8-d6ea-4685-a10b-2c4af167b89a",
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2024-08-20T22:19:14.817471Z",
+ "start_time": "2024-08-20T22:19:04.867861Z"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": "\n\n\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "eb2e804a4ee64c3483bbc526299d3471",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./hello_world_seeding.pdf\n\n\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "e1257a7fbd6d41dcb5887fb3e4fca118",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./hello_world_no_seeding.pdf…"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "for case, output in outputs.items():\n",
+ " time = output['products']['time']\n",
+ " water_mass = output['attributes']['water mass']\n",
+ " \n",
+ " fig, axs = pyplot.subplot_mosaic(\n",
+ " [['a', 'b', 'c']],\n",
+ " sharey=True,\n",
+ " figsize=(12, 4),\n",
+ " tight_layout=True\n",
+ " )\n",
+ " \n",
+ " for drop_id in range(water_mass.shape[1]):\n",
+ " axs['a'].plot(\n",
+ " in_unit(water_mass[:, drop_id], si.ng),\n",
+ " in_unit(time, si.min),\n",
+ " color=\"navy\" if np.isfinite(water_mass[0, drop_id]) else \"red\",\n",
+ " linewidth=0.333,\n",
+ " )\n",
+ " axs['a'].set_ylabel(\"time [min]\")\n",
+ " axs['a'].set_xlabel(\"drop mass [ng]\")\n",
+ " axs['a'].grid()\n",
+ " axs['a'].set_xscale(\"log\")\n",
+ " axs['a'].set_xlim(1e-6, 1e8)\n",
+ "\n",
+ " axs['b'].plot(\n",
+ " output['products']['sd_count'],\n",
+ " in_unit(time, si.min),\n",
+ " marker='.',\n",
+ " color='green',\n",
+ " )\n",
+ " axs['b'].set_xlabel(\"super droplet count\")\n",
+ " axs['b'].grid()\n",
+ " axs['b'].set_xlim(95, 125)\n",
+ "\n",
+ " axs['c'].plot(\n",
+ " in_unit(output['products']['rain water mixing ratio'], si.g/si.kg),\n",
+ " in_unit(time, si.min),\n",
+ " marker='.',\n",
+ " color='green', \n",
+ " )\n",
+ " axs['c'].set_xlabel(f\"rain water mixing ratio [g/kg] (radius > {in_unit(rain_water_radius_threshold, si.mm)} mm)\")\n",
+ " axs['c'].grid()\n",
+ " axs['c'].set_xlim(0, 3)\n",
+ "\n",
+ " fig.suptitle(case)\n",
+ " show_plot(f\"hello_world_{case.replace(' ', '_')}.pdf\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a4742eaf-dccb-48b7-9138-74ae5c579cd8",
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2024-08-20T22:19:14.881814Z",
+ "start_time": "2024-08-20T22:19:14.872933Z"
+ }
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3.10.9 ('pysdm')",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.9"
+ },
+ "vscode": {
+ "interpreter": {
+ "hash": "b14f34a08619f4a218d80d7380beed3f0c712c89ff93e7183219752d640ed427"
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/examples/PySDM_examples/seeding/seeding_no_collisions.ipynb b/examples/PySDM_examples/seeding/seeding_no_collisions.ipynb
new file mode 100644
index 0000000000..7670746837
--- /dev/null
+++ b/examples/PySDM_examples/seeding/seeding_no_collisions.ipynb
@@ -0,0 +1,164 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "from matplotlib import pyplot\n",
+ "from PySDM import Formulae\n",
+ "from PySDM.physics import in_unit, si\n",
+ "from open_atmos_jupyter_utils import show_plot\n",
+ "from PySDM_examples.seeding import Settings, Simulation"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n_sd_initial = 100\n",
+ "n_sd_seeding = 100\n",
+ "formulae = Formulae(seed=100)\n",
+ "\n",
+ "simulations = {\n",
+ " case: Simulation(\n",
+ " Settings(\n",
+ " n_sd_initial=n_sd_initial,\n",
+ " n_sd_seeding=n_sd_seeding,\n",
+ " rain_water_radius_threshold=0 * si.mm,\n",
+ " super_droplet_injection_rate={\n",
+ " 'seeding': lambda time: 1 if 15 * si.min < time < 20 * si.min else 0,\n",
+ " 'no seeding': lambda _: 0,\n",
+ " }[case],\n",
+ " formulae=formulae,\n",
+ " enable_collisions=False,\n",
+ " )\n",
+ " )\n",
+ " for case in ('seeding', 'no seeding')\n",
+ "} "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "outputs = {case: simulations[case].run() for case in simulations}"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": "\n\n\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "65e46cc48f3f4c28bcf9fba385cb721f",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HTML(value=\"./seeding_no_collisions.pdf
\")"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "fig, axs = pyplot.subplot_mosaic(\n",
+ " [['a', 'b', 'c']],\n",
+ " sharey=True,\n",
+ " figsize=(12, 4),\n",
+ " tight_layout=True\n",
+ " )\n",
+ "\n",
+ "for case, output in outputs.items():\n",
+ " time = output['products']['time']\n",
+ " idx = np.where(time > 10 * si.min)[0]\n",
+ "\n",
+ " axs['a'].plot(\n",
+ " output['products']['sd_count'][idx],\n",
+ " in_unit(time[idx], si.min),\n",
+ " marker='.',\n",
+ " label=case,\n",
+ " )\n",
+ " axs['a'].set_xlabel(\"super droplet count\")\n",
+ " axs['a'].set_ylabel(\"height [m]\")\n",
+ " axs['a'].axhspan(15, 20, color=\"grey\", alpha=0.2)\n",
+ "\n",
+ " axs['b'].plot(\n",
+ " output['products']['r_eff'][idx],\n",
+ " in_unit(time[idx], si.min),\n",
+ " marker='.',\n",
+ " )\n",
+ " axs['b'].set_xlabel(\"r_eff [um]\")\n",
+ " axs['b'].axhspan(15, 20, color=\"grey\", alpha=0.2)\n",
+ "\n",
+ " axs['c'].plot(\n",
+ " output['products']['n_drop'][idx],\n",
+ " in_unit(time[idx], si.min),\n",
+ " marker='.',\n",
+ " )\n",
+ " axs['c'].set_xlabel(\"n_drop [cm-3]\")\n",
+ " axs['c'].axhspan(15, 20, color=\"grey\", alpha=0.2)\n",
+ "\n",
+ "axs['a'].legend()\n",
+ "fig.suptitle(\"parcel with no collisions\")\n",
+ "show_plot(\"seeding_no_collisions.pdf\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3.10.9 ('pysdm')",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.9"
+ },
+ "orig_nbformat": 4,
+ "vscode": {
+ "interpreter": {
+ "hash": "b14f34a08619f4a218d80d7380beed3f0c712c89ff93e7183219752d640ed427"
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/examples/PySDM_examples/seeding/settings.py b/examples/PySDM_examples/seeding/settings.py
new file mode 100644
index 0000000000..051595e454
--- /dev/null
+++ b/examples/PySDM_examples/seeding/settings.py
@@ -0,0 +1,65 @@
+import numpy as np
+from pystrict import strict
+from PySDM.initialisation.spectra import Lognormal
+from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity
+from PySDM.physics import si
+from PySDM import Formulae
+
+
+@strict
+class Settings:
+ def __init__(
+ self,
+ *,
+ super_droplet_injection_rate: callable,
+ n_sd_initial: int,
+ n_sd_seeding: int,
+ rain_water_radius_threshold: float,
+ formulae: Formulae,
+ enable_collisions: bool = True,
+ ):
+ self.enable_collisions = enable_collisions
+ self.formulae = formulae
+ self.n_sd_initial = n_sd_initial
+ self.n_sd_seeding = n_sd_seeding
+ self.rain_water_radius_threshold = rain_water_radius_threshold
+
+ self.t_max = 25 * si.min
+ self.w_max = 3 * si.m / si.s
+ self.w_min = 0.025 * si.m / si.s
+
+ self.timestep = 15 * si.s
+ self.mass_of_dry_air = 666 * si.kg
+
+ self.updraft = (
+ lambda t: self.w_min
+ + (self.w_max - self.w_min)
+ * np.maximum(0, np.sin(t / self.t_max * 2 * np.pi)) ** 2
+ )
+ self.initial_water_vapour_mixing_ratio = 666 / 30 * si.g / si.kg
+ self.initial_total_pressure = 1000 * si.hPa
+ self.initial_temperature = 300 * si.K
+ self.initial_aerosol_kappa = 0.5
+ self.initial_aerosol_dry_radii = Lognormal(
+ norm_factor=200 / si.mg * self.mass_of_dry_air,
+ m_mode=75 * si.nm,
+ s_geom=1.6,
+ )
+ self.super_droplet_injection_rate = super_droplet_injection_rate
+
+ r_dry, n_in_dv = ConstantMultiplicity(
+ Lognormal(
+ norm_factor=10 / si.mg * self.mass_of_dry_air,
+ m_mode=1 * si.um,
+ s_geom=1.1,
+ )
+ ).sample(
+ n_sd=n_sd_seeding
+ ) # TODO #1367: does not to be the same?
+ v_dry = self.formulae.trivia.volume(radius=r_dry)
+ self.seeded_particle_multiplicity = n_in_dv
+ self.seeded_particle_extensive_attributes = {
+ "water mass": [0.0001 * si.ng] * n_sd_seeding,
+ "dry volume": v_dry,
+ "kappa times dry volume": 0.8 * v_dry,
+ }
diff --git a/examples/PySDM_examples/seeding/simulation.py b/examples/PySDM_examples/seeding/simulation.py
new file mode 100644
index 0000000000..668fb58e18
--- /dev/null
+++ b/examples/PySDM_examples/seeding/simulation.py
@@ -0,0 +1,97 @@
+import numpy as np
+
+from PySDM_examples.seeding.settings import Settings
+
+from PySDM import Builder
+from PySDM.backends import CPU
+from PySDM.environments import Parcel
+from PySDM.dynamics import Condensation, AmbientThermodynamics, Coalescence, Seeding
+from PySDM.dynamics.collisions.collision_kernels import Geometric
+from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity
+from PySDM import products
+from PySDM.physics import si
+
+
+class Simulation:
+ def __init__(self, settings: Settings):
+ builder = Builder(
+ n_sd=settings.n_sd_seeding + settings.n_sd_initial,
+ backend=CPU(
+ formulae=settings.formulae, override_jit_flags={"parallel": False}
+ ),
+ environment=Parcel(
+ dt=settings.timestep,
+ mass_of_dry_air=settings.mass_of_dry_air,
+ w=settings.updraft,
+ initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio,
+ p0=settings.initial_total_pressure,
+ T0=settings.initial_temperature,
+ ),
+ )
+ builder.add_dynamic(AmbientThermodynamics())
+ builder.add_dynamic(Condensation())
+ if settings.enable_collisions:
+ builder.add_dynamic(Coalescence(collision_kernel=Geometric()))
+ builder.add_dynamic(
+ Seeding(
+ super_droplet_injection_rate=settings.super_droplet_injection_rate,
+ seeded_particle_multiplicity=settings.seeded_particle_multiplicity,
+ seeded_particle_extensive_attributes=settings.seeded_particle_extensive_attributes,
+ )
+ )
+
+ r_dry, n_in_dv = ConstantMultiplicity(
+ settings.initial_aerosol_dry_radii
+ ).sample(n_sd=settings.n_sd_initial, backend=builder.particulator.backend)
+ attributes = builder.particulator.environment.init_attributes(
+ n_in_dv=n_in_dv, kappa=settings.initial_aerosol_kappa, r_dry=r_dry
+ )
+ self.particulator = builder.build(
+ attributes={
+ k: np.pad(
+ array=v,
+ pad_width=(0, settings.n_sd_seeding),
+ mode="constant",
+ constant_values=np.nan if k == "multiplicity" else 0,
+ )
+ for k, v in attributes.items()
+ },
+ products=(
+ products.SuperDropletCountPerGridbox(name="sd_count"),
+ products.Time(),
+ products.WaterMixingRatio(
+ radius_range=(settings.rain_water_radius_threshold, np.inf),
+ name="rain water mixing ratio",
+ ),
+ products.EffectiveRadius(
+ name="r_eff",
+ unit="um",
+ radius_range=(0.5 * si.um, 25 * si.um),
+ ),
+ products.ParticleConcentration(
+ name="n_drop",
+ unit="cm^-3",
+ radius_range=(0.5 * si.um, 25 * si.um),
+ ),
+ ),
+ )
+ self.n_steps = int(settings.t_max // settings.timestep)
+
+ def run(self):
+ output = {
+ "attributes": {"water mass": []},
+ "products": {key: [] for key in self.particulator.products},
+ }
+ for step in range(self.n_steps + 1):
+ if step != 0:
+ self.particulator.run(steps=1)
+ for key, attr in output["attributes"].items():
+ data = self.particulator.attributes[key].to_ndarray(raw=True)
+ data[data == 0] = np.nan
+ attr.append(data)
+ for key, prod in output["products"].items():
+ prod.append(float(self.particulator.products[key].get()))
+ for out in ("attributes", "products"):
+ for key, val in output[out].items():
+ output[out][key] = np.array(val)
+ return output
diff --git a/tests/examples_tests/conftest.py b/tests/examples_tests/conftest.py
index a9de5ac523..969bae90ea 100644
--- a/tests/examples_tests/conftest.py
+++ b/tests/examples_tests/conftest.py
@@ -53,6 +53,7 @@ def findfiles(path, regex):
"deJong_Azimi",
"Bulenok_2023_MasterThesis",
"Shipway_and_Hill_2012",
+ "seeding",
],
"multi-process_b": [
"Bartman_2020_MasterThesis",
diff --git a/tests/smoke_tests/parcel_a/seeding/__init__.py b/tests/smoke_tests/parcel_a/seeding/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/tests/smoke_tests/parcel_a/seeding/test_hello_world.py b/tests/smoke_tests/parcel_a/seeding/test_hello_world.py
new file mode 100644
index 0000000000..607490eed8
--- /dev/null
+++ b/tests/smoke_tests/parcel_a/seeding/test_hello_world.py
@@ -0,0 +1,52 @@
+""" tests ensuring values on plots match those in the paper """
+
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from PySDM_examples.utils import notebook_vars
+from PySDM_examples import seeding
+
+PLOT = False
+
+
+@pytest.fixture(scope="session", name="variables")
+def variables_fixture():
+ return notebook_vars(
+ file=Path(seeding.__file__).parent / "hello_world.ipynb",
+ plot=PLOT,
+ )
+
+
+class TestHelloWorld:
+ @staticmethod
+ def test_sd_count(variables):
+ minimum = variables["n_sd_initial"]
+ maximum = minimum + variables["n_sd_seeding"]
+ assert variables["outputs"]["seeding"]["products"]["sd_count"][0] == minimum
+ assert (
+ minimum
+ < variables["outputs"]["seeding"]["products"]["sd_count"][-1]
+ < maximum
+ )
+ np.testing.assert_equal(
+ variables["outputs"]["no seeding"]["products"]["sd_count"], minimum
+ )
+
+ @staticmethod
+ def test_final_rain_water_mixing_ratio_larger_with_seeding(variables):
+ assert (
+ variables["outputs"]["seeding"]["products"]["rain water mixing ratio"][-1]
+ > variables["outputs"]["no seeding"]["products"]["rain water mixing ratio"][
+ -1
+ ]
+ )
+
+ @staticmethod
+ def test_rain_water_earlier_with_seeding(variables):
+ assert np.count_nonzero(
+ variables["outputs"]["seeding"]["products"]["rain water mixing ratio"]
+ ) > np.count_nonzero(
+ variables["outputs"]["no seeding"]["products"]["rain water mixing ratio"]
+ )
diff --git a/tests/smoke_tests/parcel_a/seeding/test_seeding_no_collisions.py b/tests/smoke_tests/parcel_a/seeding/test_seeding_no_collisions.py
new file mode 100644
index 0000000000..436a933dbc
--- /dev/null
+++ b/tests/smoke_tests/parcel_a/seeding/test_seeding_no_collisions.py
@@ -0,0 +1,37 @@
+""" tests ensuring values on plots match those in the paper """
+
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from PySDM_examples.utils import notebook_vars
+from PySDM_examples import seeding
+
+PLOT = False
+
+
+@pytest.fixture(scope="session", name="variables")
+def variables_fixture():
+ return notebook_vars(
+ file=Path(seeding.__file__).parent / "seeding_no_collisions.ipynb",
+ plot=PLOT,
+ )
+
+
+class TestSeedingNoCollisions:
+ @staticmethod
+ # seeding has smaller cloud drops than no seeding
+ def test_reff(variables):
+ np.testing.assert_array_less(
+ variables["outputs"]["seeding"]["products"]["r_eff"],
+ variables["outputs"]["no seeding"]["products"]["r_eff"] + 1e-6,
+ )
+
+ @staticmethod
+ # seeding has more cloud drops than no seeding
+ def test_n_drop(variables):
+ np.testing.assert_array_less(
+ variables["outputs"]["no seeding"]["products"]["n_drop"],
+ variables["outputs"]["seeding"]["products"]["n_drop"] + 1e-6,
+ )
diff --git a/tests/unit_tests/backends/test_seeding_methods.py b/tests/unit_tests/backends/test_seeding_methods.py
new file mode 100644
index 0000000000..d7002785b5
--- /dev/null
+++ b/tests/unit_tests/backends/test_seeding_methods.py
@@ -0,0 +1,158 @@
+""" Seeding backend tests of injection logic """
+
+from contextlib import nullcontext
+
+import numpy as np
+import pytest
+
+from PySDM import Builder
+from PySDM.environments import Box
+from PySDM.backends import CPU
+from PySDM.physics import si
+
+
+class TestSeeding:
+ max_number_to_inject = 4
+
+ @staticmethod
+ @pytest.mark.parametrize(
+ "n_sd, number_of_super_particles_to_inject, context",
+ (
+ (1, 1, nullcontext()),
+ (
+ 1,
+ 2,
+ pytest.raises(
+ ValueError, match="inject more super particles than space available"
+ ),
+ ),
+ (max_number_to_inject, max_number_to_inject - 1, nullcontext()),
+ (max_number_to_inject, max_number_to_inject, nullcontext()),
+ (
+ max_number_to_inject + 2,
+ max_number_to_inject + 1,
+ pytest.raises(
+ ValueError,
+ match="inject multiple super particles with the same attributes",
+ ),
+ ),
+ ),
+ )
+ def test_number_of_super_particles_to_inject(
+ n_sd,
+ number_of_super_particles_to_inject,
+ context,
+ dt=1,
+ dv=1,
+ ):
+ # arrange
+ builder = Builder(n_sd, CPU(), Box(dt, dv))
+ particulator = builder.build(
+ attributes={
+ "multiplicity": np.full(n_sd, np.nan),
+ "water mass": np.zeros(n_sd),
+ },
+ )
+
+ seeded_particle_extensive_attributes = {
+ "water mass": [0.0001 * si.ng] * TestSeeding.max_number_to_inject,
+ }
+ seeded_particle_multiplicity = [1] * TestSeeding.max_number_to_inject
+
+ seeded_particle_index = particulator.Index.identity_index(
+ len(seeded_particle_multiplicity)
+ )
+ seeded_particle_multiplicity = particulator.IndexedStorage.from_ndarray(
+ seeded_particle_index,
+ np.asarray(seeded_particle_multiplicity),
+ )
+ seeded_particle_extensive_attributes = particulator.IndexedStorage.from_ndarray(
+ seeded_particle_index,
+ np.asarray(list(seeded_particle_extensive_attributes.values())),
+ )
+
+ # act
+ with context:
+ particulator.seeding(
+ seeded_particle_index=seeded_particle_index,
+ seeded_particle_multiplicity=seeded_particle_multiplicity,
+ seeded_particle_extensive_attributes=seeded_particle_extensive_attributes,
+ number_of_super_particles_to_inject=number_of_super_particles_to_inject,
+ )
+
+ # assert
+ assert (
+ number_of_super_particles_to_inject
+ == particulator.attributes.super_droplet_count
+ )
+
+ @staticmethod
+ @pytest.mark.parametrize(
+ "seeded_particle_index, context",
+ (
+ ([0, 0, 0], nullcontext()),
+ ([0, 1, 2], nullcontext()),
+ ([2, 1, 0], nullcontext()),
+ (
+ [0],
+ pytest.raises(
+ ValueError,
+ match=" multiple super particles with the same attributes",
+ ),
+ ),
+ ),
+ )
+ def test_seeded_particle_index_multiplicity_extensive_attributes(
+ seeded_particle_index,
+ context,
+ n_sd=3,
+ number_of_super_particles_to_inject=3,
+ ):
+ # arrange
+ builder = Builder(n_sd, CPU(), Box(dt=np.nan, dv=np.nan))
+ particulator = builder.build(
+ attributes={
+ "multiplicity": np.full(n_sd, np.nan),
+ "water mass": np.zeros(n_sd),
+ },
+ )
+
+ seeded_particle_extensive_attributes = {
+ "water mass": [0.0001, 0.0003, 0.0002],
+ }
+ seeded_particle_multiplicity = [1, 2, 3]
+
+ seeded_particle_index_impl = particulator.Index.from_ndarray(
+ np.asarray(seeded_particle_index)
+ )
+ seeded_particle_multiplicity_impl = particulator.IndexedStorage.from_ndarray(
+ seeded_particle_index_impl,
+ np.asarray(seeded_particle_multiplicity),
+ )
+ seeded_particle_extensive_attributes_impl = (
+ particulator.IndexedStorage.from_ndarray(
+ seeded_particle_index_impl,
+ np.asarray(list(seeded_particle_extensive_attributes.values())),
+ )
+ )
+
+ # act
+ with context:
+ particulator.seeding(
+ seeded_particle_index=seeded_particle_index_impl,
+ seeded_particle_multiplicity=seeded_particle_multiplicity_impl,
+ seeded_particle_extensive_attributes=seeded_particle_extensive_attributes_impl,
+ number_of_super_particles_to_inject=number_of_super_particles_to_inject,
+ )
+
+ # assert
+ np.testing.assert_array_equal(
+ particulator.attributes["multiplicity"].to_ndarray(),
+ np.asarray(seeded_particle_multiplicity)[seeded_particle_index],
+ )
+ np.testing.assert_array_equal(
+ particulator.attributes["water mass"].to_ndarray(),
+ np.asarray(seeded_particle_extensive_attributes["water mass"])[
+ seeded_particle_index
+ ],
+ )
diff --git a/tests/unit_tests/dynamics/displacement/test_sedimentation.py b/tests/unit_tests/dynamics/displacement/test_sedimentation.py
index 8028625220..776230fdc4 100644
--- a/tests/unit_tests/dynamics/displacement/test_sedimentation.py
+++ b/tests/unit_tests/dynamics/displacement/test_sedimentation.py
@@ -1,4 +1,4 @@
-# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring
+# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring,no-member
import numpy as np
from .displacement_settings import DisplacementSettings
diff --git a/tests/unit_tests/dynamics/test_seeding.py b/tests/unit_tests/dynamics/test_seeding.py
new file mode 100644
index 0000000000..b7dad9fa6a
--- /dev/null
+++ b/tests/unit_tests/dynamics/test_seeding.py
@@ -0,0 +1,251 @@
+""" Seeding dynamic tests """
+
+from collections import namedtuple
+
+import numpy as np
+import pytest
+
+from matplotlib import pyplot
+
+from PySDM import Builder
+from PySDM.products import SuperDropletCountPerGridbox, Time
+from PySDM.backends import CPU
+from PySDM.backends.impl_common.index import make_Index
+from PySDM.backends.impl_common.indexed_storage import make_IndexedStorage
+from PySDM.dynamics import Seeding
+from PySDM.environments import Box
+from PySDM.physics import si
+
+
+class TestSeeding:
+ @staticmethod
+ def test_zero_injection_rate_same_as_no_seeding_monodisperse(
+ plot=False, backend_instance=CPU()
+ ):
+ """a not-so-unit test checking that results of a box simulation
+ are the same without seeding as with a zero injection rate"""
+ # arrange
+ n_sd_seeding = 100
+ n_sd_initial = 100
+
+ def simulation(*, dynamics):
+ t_max = 20 * si.min
+ timestep = 15 * si.s
+ dv = 1 * si.cm**3
+
+ builder = Builder(
+ backend=backend_instance,
+ n_sd=n_sd_seeding + n_sd_initial,
+ environment=Box(dt=timestep, dv=dv),
+ )
+ for dynamic in dynamics:
+ builder.add_dynamic(dynamic)
+
+ particulator = builder.build(
+ attributes={
+ k: np.pad(
+ array=v,
+ pad_width=(0, n_sd_seeding),
+ mode="constant",
+ constant_values=np.nan if k == "multiplicity" else 0,
+ )
+ for k, v in {
+ "volume": np.ones(n_sd_initial) * si.um**3,
+ "multiplicity": np.ones(n_sd_initial),
+ }.items()
+ },
+ products=(
+ SuperDropletCountPerGridbox(name="sd_count"),
+ Time(),
+ ),
+ )
+ products = {"sd_count": [], "time": []}
+ for step in range(int(t_max // timestep) + 1):
+ if step != 0:
+ particulator.run(steps=1)
+ for key, val in products.items():
+ val.append(float(particulator.products[key].get()))
+ for key in products:
+ products[key] = np.array(products[key])
+ return products
+
+ # act
+ common_seeding_ctor_args = {
+ "seeded_particle_multiplicity": [1],
+ "seeded_particle_extensive_attributes": {
+ "water mass": [0.001 * si.ng],
+ },
+ }
+ output = {
+ "zero_injection_rate": simulation(
+ dynamics=(
+ Seeding(
+ super_droplet_injection_rate=lambda time: 0,
+ **common_seeding_ctor_args,
+ ),
+ )
+ ),
+ "positive_injection_rate": simulation(
+ dynamics=(
+ Seeding(
+ super_droplet_injection_rate=lambda time: (
+ 1 if 10 * si.min <= time < 15 * si.min else 0
+ ),
+ **common_seeding_ctor_args,
+ ),
+ )
+ ),
+ "no_seeding": simulation(dynamics=()),
+ }
+
+ # plot
+ pyplot.plot(
+ output["zero_injection_rate"]["sd_count"],
+ output["zero_injection_rate"]["time"],
+ lw=2,
+ label="with seeding dynamic, zero injection rate",
+ )
+ pyplot.plot(
+ output["positive_injection_rate"]["sd_count"],
+ output["positive_injection_rate"]["time"],
+ lw=2,
+ ls="--",
+ label="with seeding dynamic, positive injection rate",
+ )
+ pyplot.plot(
+ output["no_seeding"]["sd_count"],
+ output["no_seeding"]["time"],
+ lw=2,
+ ls=":",
+ label="without seeding dynamic",
+ )
+ pyplot.xlabel("sd_count")
+ pyplot.ylabel("time")
+ pyplot.grid()
+ pyplot.legend()
+ if plot:
+ pyplot.show()
+ else:
+ pyplot.clf()
+
+ # assert
+ np.testing.assert_array_equal(
+ output["zero_injection_rate"]["sd_count"], output["no_seeding"]["sd_count"]
+ )
+ np.testing.assert_array_equal(np.diff(output["no_seeding"]["sd_count"]), 0)
+ assert np.amax(np.diff(output["positive_injection_rate"]["sd_count"])) >= 0
+ assert output["positive_injection_rate"]["sd_count"][0] == n_sd_initial
+ assert (
+ n_sd_initial
+ < output["positive_injection_rate"]["sd_count"][-1]
+ < n_sd_initial + n_sd_seeding
+ )
+
+ @staticmethod
+ def test_attribute_set_match():
+ # arrange
+ extensive_attributes = ["a", "b", "c"]
+ seeding_attributes = ["d", "e", "f"]
+
+ builder = namedtuple(typename="MockBuilder", field_names=("particulator",))(
+ particulator=namedtuple(
+ typename="MockParticulator", field_names=("n_steps", "attributes")
+ )(
+ n_steps=0,
+ attributes=namedtuple(
+ typename="MockAttributes",
+ field_names=("get_extensive_attribute_keys",),
+ )(get_extensive_attribute_keys=lambda: extensive_attributes),
+ )
+ )
+
+ dynamic = Seeding(
+ super_droplet_injection_rate=lambda time: 0,
+ seeded_particle_extensive_attributes={
+ k: [np.nan] for k in seeding_attributes
+ },
+ seeded_particle_multiplicity=[0],
+ )
+ dynamic.register(builder)
+
+ # act
+ with pytest.raises(ValueError) as excinfo:
+ dynamic()
+
+ # assert
+ assert "do not match those used in particulator" in str(excinfo)
+
+ @staticmethod
+ @pytest.mark.parametrize(
+ "super_droplet_injection_rate, reservoir_length",
+ (
+ (0, 0), # shuffle not called
+ (0, 2),
+ (1, 10),
+ (2, 10),
+ ),
+ )
+ def test_seeded_particle_shuffle(
+ super_droplet_injection_rate, reservoir_length, backend=CPU()
+ ):
+ # arrange
+ extensive_attributes = ["a"]
+ seeding_attributes = ["a"]
+
+ class MockParticulator: # pylint: disable=too-few-public-methods
+ def __init__(self):
+ self.seeding_call_count = 0
+ self.indices = []
+
+ def seeding(
+ self,
+ *,
+ seeded_particle_index,
+ number_of_super_particles_to_inject, # pylint: disable=unused-argument
+ seeded_particle_multiplicity, # pylint: disable=unused-argument
+ seeded_particle_extensive_attributes, # pylint: disable=unused-argument
+ ):
+ self.seeding_call_count += 1
+ self.indices.append(seeded_particle_index.to_ndarray())
+
+ Index = make_Index(backend)
+ IndexedStorage = make_IndexedStorage(backend)
+ Random = None if reservoir_length == 0 else backend.Random
+ formulae = None if reservoir_length == 0 else backend.formulae
+ Storage = None if reservoir_length == 0 else backend.Storage
+ n_steps = 0
+ dt = np.nan
+ attributes = namedtuple(
+ typename="MockAttributes",
+ field_names=("get_extensive_attribute_keys",),
+ )(get_extensive_attribute_keys=lambda: extensive_attributes)
+
+ builder = namedtuple(typename="MockBuilder", field_names=("particulator",))(
+ particulator=MockParticulator()
+ )
+
+ dynamic = Seeding(
+ super_droplet_injection_rate=lambda t: super_droplet_injection_rate,
+ seeded_particle_extensive_attributes={
+ k: [np.nan] * reservoir_length for k in seeding_attributes
+ },
+ seeded_particle_multiplicity=[1] * reservoir_length,
+ )
+ dynamic.register(builder)
+
+ # act
+ dynamic()
+ dynamic.particulator.n_steps += 1
+ dynamic()
+
+ # assert
+ assert dynamic.particulator.seeding_call_count == (
+ 2 if super_droplet_injection_rate > 0 else 0
+ )
+ if super_droplet_injection_rate > 0:
+ assert (
+ dynamic.particulator.indices[0] != dynamic.particulator.indices[1]
+ ).any()
+ assert sorted(dynamic.particulator.indices[0]) == sorted(
+ dynamic.particulator.indices[1]
+ )
diff --git a/tests/unit_tests/impl/test_particle_attributes.py b/tests/unit_tests/impl/test_particle_attributes.py
index 4f83be1776..4a8f687399 100644
--- a/tests/unit_tests/impl/test_particle_attributes.py
+++ b/tests/unit_tests/impl/test_particle_attributes.py
@@ -1,4 +1,4 @@
-# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring
+# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring,no-member
import numpy as np
import pytest
diff --git a/tests/unit_tests/initialisation/test_discretise_multiplicities.py b/tests/unit_tests/initialisation/test_discretise_multiplicities.py
new file mode 100644
index 0000000000..daf765e31b
--- /dev/null
+++ b/tests/unit_tests/initialisation/test_discretise_multiplicities.py
@@ -0,0 +1,79 @@
+""" checks verifying imeplementation of `discretise_multiplicities` function """
+
+import pytest
+import numpy as np
+
+from PySDM.initialisation import discretise_multiplicities
+
+
+class TestDiscretiseMultiplicities:
+ @staticmethod
+ def test_reporting_zero_multiplicity():
+ """should raise an ValueError if int-casting results in zeros"""
+ # arrange
+ values = np.asarray([0.1], dtype=np.float64)
+
+ # act
+ with pytest.raises(ValueError) as excinfo:
+ discretise_multiplicities(values)
+
+ # assert
+ assert "int-casting resulted in multiplicity of zero" in str(excinfo.value)
+
+ @staticmethod
+ def test_reporting_sum_error():
+ """should err if sum of int-casted values differs a lot from the sum of floats"""
+ # arrange
+ values = np.asarray(
+ [
+ 1.49,
+ ],
+ dtype=np.float64,
+ )
+
+ # act
+ with pytest.raises(ValueError) as excinfo:
+ discretise_multiplicities(values)
+
+ # assert
+ assert (
+ "error in total real-droplet number due to casting multiplicities to ints"
+ in str(excinfo.value)
+ )
+
+ @staticmethod
+ def test_nans_converted_to_zeros():
+ """checks for how NaN values are treated"""
+ # arrange
+ values = np.asarray(
+ [
+ np.nan,
+ 1,
+ ],
+ dtype=np.float64,
+ )
+
+ # act
+ ints = discretise_multiplicities(values)
+
+ # assert
+ assert ints[0] == 0
+
+ @staticmethod
+ def test_all_nans_to_all_zeros():
+ """checks for how NaN values are treated"""
+ # arrange
+ values = np.asarray(
+ [
+ np.nan,
+ np.nan,
+ np.nan,
+ ],
+ dtype=np.float64,
+ )
+
+ # act
+ ints = discretise_multiplicities(values)
+
+ # assert
+ assert (ints == 0).all()
diff --git a/tests/unit_tests/test_particulator.py b/tests/unit_tests/test_particulator.py
index 8c94688ad8..db2d19f00d 100644
--- a/tests/unit_tests/test_particulator.py
+++ b/tests/unit_tests/test_particulator.py
@@ -1,4 +1,6 @@
# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring
+from collections import namedtuple
+
import pytest
from .dummy_particulator import DummyParticulator
@@ -53,3 +55,88 @@ class DP(DummyParticulator):
assert particulator.attributes.updated == [
f"moles_{isotope}" for isotope in isotopes
]
+
+ @staticmethod
+ def test_seeding_marks_modified_attributes_as_updated(backend_class):
+ # arrange
+ storage = backend_class().Storage.empty(1, dtype=int)
+ abc = ["a", "b", "c"]
+
+ class ParticleAttributes:
+ def __init__(self):
+ self.updated = []
+ self.super_droplet_count = -1
+ self.__idx = storage # pylint: disable=unused-private-member
+ self.idx_reset = False
+ self.sane = False
+
+ def get_extensive_attribute_storage(self):
+ return storage
+
+ def get_extensive_attribute_keys(self):
+ return abc
+
+ def __getitem__(self, item):
+ return storage
+
+ def mark_updated(self, attr):
+ self.updated += [attr]
+
+ def reset_idx(self):
+ self.idx_reset = True
+
+ def sanitize(self):
+ self.sane = True
+
+ class DP(DummyParticulator):
+ pass
+
+ particulator = DP(backend_class, 44)
+ particulator.attributes = ParticleAttributes()
+ # fmt: off
+ particulator.backend.seeding = (
+ lambda
+ idx,
+ multiplicity,
+ extensive_attributes,
+ seeded_particle_index,
+ seeded_particle_multiplicity,
+ seeded_particle_extensive_attributes,
+ number_of_super_particles_to_inject:
+ None
+ )
+ # fmt: on
+
+ # act
+ particulator.seeding(
+ seeded_particle_index=storage,
+ seeded_particle_multiplicity=storage,
+ seeded_particle_extensive_attributes=storage,
+ number_of_super_particles_to_inject=0,
+ )
+
+ # assert
+ assert particulator.attributes.updated == ["multiplicity"] + abc
+ assert particulator.attributes.idx_reset
+ assert particulator.attributes.sane
+
+ @staticmethod
+ def test_seeding_fails_if_no_null_super_droplets_availale(backend_class):
+ # arrange
+ a_number = 44
+
+ particulator = DummyParticulator(backend_class, n_sd=a_number)
+ storage = backend_class().Storage.empty(1, dtype=int)
+
+ particulator.attributes = namedtuple(
+ typename="_", field_names=("super_droplet_count",)
+ )(super_droplet_count=a_number)
+
+ # act
+ with pytest.raises(ValueError, match="No available seeds to inject"):
+ particulator.seeding(
+ seeded_particle_index=storage,
+ seeded_particle_multiplicity=storage,
+ seeded_particle_extensive_attributes=storage,
+ number_of_super_particles_to_inject=0,
+ )