diff --git a/PySDM/dynamics/collisions/collision.py b/PySDM/dynamics/collisions/collision.py index 2a917911ce..9f4e940572 100644 --- a/PySDM/dynamics/collisions/collision.py +++ b/PySDM/dynamics/collisions/collision.py @@ -18,6 +18,7 @@ from PySDM.dynamics.collisions.breakup_efficiencies import ConstEb from PySDM.dynamics.collisions.breakup_fragmentations import AlwaysN from PySDM.dynamics.collisions.coalescence_efficiencies import ConstEc +from PySDM.dynamics.collisions import collision_kernels from PySDM.dynamics.impl.random_generator_optimizer import RandomGeneratorOptimizer from PySDM.dynamics.impl.random_generator_optimizer_nopair import ( RandomGeneratorOptimizerNoPair, @@ -42,7 +43,6 @@ class Collision: # pylint: disable=too-many-instance-attributes def __init__( self, *, - collision_kernel, coalescence_efficiency, breakup_efficiency, fragmentation_function, @@ -63,7 +63,6 @@ def __init__( self.warn_overflows = warn_overflows self.max_multiplicity = DEFAULTS.max_multiplicity - self.collision_kernel = collision_kernel self.compute_coalescence_efficiency = coalescence_efficiency self.compute_breakup_efficiency = breakup_efficiency self.compute_number_of_fragments = fragmentation_function @@ -98,6 +97,8 @@ def __init__( self.breakup_rate = None self.breakup_rate_deficit = None + self.collision_kernel = None + def register(self, builder): self.particulator = builder.particulator rnd_args = { @@ -136,6 +137,10 @@ def register(self, builder): self.stats_dt_min.fill(np.nan) self.rnd_opt_coll.register(builder) + self.collision_kernel = getattr( + collision_kernels, + self.particulator.formulae.collision_kernel_liquid_liquid.__name__, + )() self.collision_kernel.register(builder) if self.croupier is None: @@ -295,7 +300,6 @@ class Coalescence(Collision): def __init__( self, *, - collision_kernel, coalescence_efficiency=ConstEc(Ec=1), croupier=None, optimized_random=False, @@ -306,7 +310,6 @@ def __init__( breakup_efficiency = ConstEb(Eb=0) fragmentation_function = AlwaysN(n=1) super().__init__( - collision_kernel=collision_kernel, coalescence_efficiency=coalescence_efficiency, breakup_efficiency=breakup_efficiency, fragmentation_function=fragmentation_function, @@ -324,7 +327,6 @@ class Breakup(Collision): def __init__( self, *, - collision_kernel, fragmentation_function, croupier=None, optimized_random=False, @@ -336,7 +338,6 @@ def __init__( coalescence_efficiency = ConstEc(Ec=0.0) breakup_efficiency = ConstEb(Eb=1.0) super().__init__( - collision_kernel=collision_kernel, coalescence_efficiency=coalescence_efficiency, breakup_efficiency=breakup_efficiency, fragmentation_function=fragmentation_function, diff --git a/PySDM/dynamics/collisions/collision_kernels/__init__.py b/PySDM/dynamics/collisions/collision_kernels/__init__.py index 26ddc443e9..ef7c1bc3d7 100644 --- a/PySDM/dynamics/collisions/collision_kernels/__init__.py +++ b/PySDM/dynamics/collisions/collision_kernels/__init__.py @@ -13,3 +13,4 @@ from .linear import Linear from .simple_geometric import SimpleGeometric from .long1974 import Long1974 +from .neglect import Neglect diff --git a/PySDM/dynamics/collisions/collision_kernels/constantK.py b/PySDM/dynamics/collisions/collision_kernels/constantK.py index d906be151a..a87c9a3197 100644 --- a/PySDM/dynamics/collisions/collision_kernels/constantK.py +++ b/PySDM/dynamics/collisions/collision_kernels/constantK.py @@ -4,12 +4,11 @@ class ConstantK: - def __init__(self, a): - self.a = a + def __init__(self): self.particulator = None def __call__(self, output, is_first_in_pair): - output.fill(self.a) + output.fill(self.particulator.formulae.constants.CONSTANTK_a) def register(self, builder): self.particulator = builder.particulator diff --git a/PySDM/dynamics/collisions/collision_kernels/geometric.py b/PySDM/dynamics/collisions/collision_kernels/geometric.py index 365cfa96b5..df8b7ede25 100644 --- a/PySDM/dynamics/collisions/collision_kernels/geometric.py +++ b/PySDM/dynamics/collisions/collision_kernels/geometric.py @@ -7,15 +7,13 @@ class Geometric(Gravitational): - def __init__(self, collection_efficiency=1.0, x="volume"): - super().__init__() - self.collection_efficiency = collection_efficiency - self.x = x - def __call__(self, output, is_first_in_pair): output.sum(self.particulator.attributes["radius"], is_first_in_pair) output **= 2 - output *= const.PI * self.collection_efficiency + output *= ( + const.PI + * self.particulator.formulae.constants.GEOMETRIC_collection_efficiency + ) self.pair_tmp.distance( self.particulator.attributes["relative fall velocity"], is_first_in_pair ) diff --git a/PySDM/dynamics/collisions/collision_kernels/golovin.py b/PySDM/dynamics/collisions/collision_kernels/golovin.py index 581fd6a75e..cafa7fc2aa 100644 --- a/PySDM/dynamics/collisions/collision_kernels/golovin.py +++ b/PySDM/dynamics/collisions/collision_kernels/golovin.py @@ -4,43 +4,15 @@ [Golovin 1963](http://mi.mathnet.ru/dan27630)) """ -import numpy as np -from scipy import special - class Golovin: - def __init__(self, b): - self.b = b + def __init__(self): self.particulator = None def __call__(self, output, is_first_in_pair): output.sum(self.particulator.attributes["volume"], is_first_in_pair) - output *= self.b + output *= self.particulator.formulae.constants.GOLOVIN_b def register(self, builder): self.particulator = builder.particulator builder.request_attribute("volume") - - def analytic_solution(self, x, t, x_0, N_0): - tau = 1 - np.exp(-N_0 * self.b * x_0 * t) - - if isinstance(x, np.ndarray): - func = np.vectorize( - lambda i: Golovin.analytic_solution_helper(x[int(i)], tau, x_0) - ) - result = np.fromfunction(func, x.shape, dtype=float) - return result - - return Golovin.analytic_solution_helper(x, tau, x_0) - - @staticmethod - def analytic_solution_helper(x, tau, x_0): - sqrt_tau = np.sqrt(tau) - result = float( - (1 - tau) - * 1 - / (x * np.sqrt(tau)) - * special.ive(1, 2 * x / x_0 * sqrt_tau) # pylint: disable=no-member - * np.exp(-(1 + tau - 2 * sqrt_tau) * x / x_0) - ) - return result diff --git a/PySDM/dynamics/collisions/collision_kernels/linear.py b/PySDM/dynamics/collisions/collision_kernels/linear.py index a5e9d219a0..1011d1189c 100644 --- a/PySDM/dynamics/collisions/collision_kernels/linear.py +++ b/PySDM/dynamics/collisions/collision_kernels/linear.py @@ -4,15 +4,13 @@ class Linear: - def __init__(self, a, b): - self.a = a - self.b = b + def __init__(self): self.particulator = None def __call__(self, output, is_first_in_pair): output.sum_pair(self.particulator.attributes["volume"], is_first_in_pair) - output *= self.b - output += self.a + output *= self.particulator.formulae.constants.LINEAR_b + output += self.particulator.formulae.constants.LINEAR_a def register(self, builder): self.particulator = builder.particulator diff --git a/PySDM/dynamics/collisions/collision_kernels/neglect.py b/PySDM/dynamics/collisions/collision_kernels/neglect.py new file mode 100644 index 0000000000..64f5e60996 --- /dev/null +++ b/PySDM/dynamics/collisions/collision_kernels/neglect.py @@ -0,0 +1,14 @@ +""" +kernel for disabling collision +""" + + +class Neglect: + def __init__(self): + self.particulator = None + + def __call__(self, output, is_first_in_pair): + output.fill(0) + + def register(self, builder): + self.particulator = builder.particulator diff --git a/PySDM/dynamics/collisions/collision_kernels/simple_geometric.py b/PySDM/dynamics/collisions/collision_kernels/simple_geometric.py index 005c61927b..774cdb4476 100644 --- a/PySDM/dynamics/collisions/collision_kernels/simple_geometric.py +++ b/PySDM/dynamics/collisions/collision_kernels/simple_geometric.py @@ -4,10 +4,9 @@ class SimpleGeometric: - def __init__(self, C): + def __init__(self): self.particulator = None self.pair_tmp = None - self.C = C def register(self, builder): self.particulator = builder.particulator @@ -18,7 +17,7 @@ def register(self, builder): ) def __call__(self, output, is_first_in_pair): - output[:] = self.C + output[:] = self.particulator.formulae.constants.SIMPLE_GEOMETRIC_c self.pair_tmp.sum(self.particulator.attributes["radius"], is_first_in_pair) self.pair_tmp **= 2 output *= self.pair_tmp diff --git a/PySDM/formulae.py b/PySDM/formulae.py index 6e53819b8a..7d814a7963 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -55,6 +55,7 @@ def __init__( # pylint: disable=too-many-locals,too-many-arguments freezing_temperature_spectrum: str = "Null", heterogeneous_ice_nucleation_rate: str = "Null", homogeneous_ice_nucleation_rate: str = "Null", + collision_kernel: str = "Neglect+Neglect+Neglect", fragmentation_function: str = "AlwaysN", isotope_equilibrium_fractionation_factors: str = "Null", isotope_kinetic_fractionation_factors: str = "Null", @@ -96,6 +97,9 @@ def __init__( # pylint: disable=too-many-locals,too-many-arguments self.freezing_temperature_spectrum = freezing_temperature_spectrum self.heterogeneous_ice_nucleation_rate = heterogeneous_ice_nucleation_rate self.homogeneous_ice_nucleation_rate = homogeneous_ice_nucleation_rate + self.collision_kernel_liquid_liquid = collision_kernel.split("+")[0] + self.collision_kernel_liquid_ice = collision_kernel.split("+")[1] + self.collision_kernel_ice_ice = collision_kernel.split("+")[2] self.fragmentation_function = fragmentation_function self.isotope_equilibrium_fractionation_factors = ( isotope_equilibrium_fractionation_factors diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py index 82cabe3bd5..a431369a58 100644 --- a/PySDM/physics/__init__.py +++ b/PySDM/physics/__init__.py @@ -54,5 +54,6 @@ terminal_velocity, terminal_velocity_ice, bulk_phase_partitioning, + collision_kernel, ) from .constants import convert_to, in_unit, si diff --git a/PySDM/physics/collision_kernel/__init__.py b/PySDM/physics/collision_kernel/__init__.py new file mode 100644 index 0000000000..c43dfa1941 --- /dev/null +++ b/PySDM/physics/collision_kernel/__init__.py @@ -0,0 +1,13 @@ +""" +Collision kernels physics +""" + +from .neglect import Neglect +from .constantK import ConstantK +from .electric import Electric +from .geometric import Geometric +from .golovin import Golovin +from .hydrodynamic import Hydrodynamic +from .linear import Linear +from .simple_geometric import SimpleGeometric +from .long1974 import Long1974 diff --git a/PySDM/physics/collision_kernel/constantK.py b/PySDM/physics/collision_kernel/constantK.py new file mode 100644 index 0000000000..efc03f737f --- /dev/null +++ b/PySDM/physics/collision_kernel/constantK.py @@ -0,0 +1,8 @@ +""" +Formulae supporting `PySDM.dynamics.collisions.collision_kernels.ConstantK` +""" + + +class ConstantK: + def __init__(self, _): + pass diff --git a/PySDM/physics/collision_kernel/electric.py b/PySDM/physics/collision_kernel/electric.py new file mode 100644 index 0000000000..b74c071e9e --- /dev/null +++ b/PySDM/physics/collision_kernel/electric.py @@ -0,0 +1,8 @@ +""" +Formulae supporting `PySDM.dynamics.collisions.collision_kernels.Electric` +""" + + +class Electric: + def __init__(self, _): + pass diff --git a/PySDM/physics/collision_kernel/geometric.py b/PySDM/physics/collision_kernel/geometric.py new file mode 100644 index 0000000000..ebb4544dda --- /dev/null +++ b/PySDM/physics/collision_kernel/geometric.py @@ -0,0 +1,8 @@ +""" +Formulae supporting `PySDM.dynamics.collisions.collision_kernels.Geometric` +""" + + +class Geometric: + def __init__(self, _): + pass diff --git a/PySDM/physics/collision_kernel/golovin.py b/PySDM/physics/collision_kernel/golovin.py new file mode 100644 index 0000000000..0ec3e26acf --- /dev/null +++ b/PySDM/physics/collision_kernel/golovin.py @@ -0,0 +1,36 @@ +""" +Formulae supporting `PySDM.dynamics.collisions.collision_kernels.golovin` +""" + +import numpy as np +from scipy import special + + +class Golovin: + def __init__(self, _): + pass + + +@staticmethod +def analytic_solution(b, x, t, x_0, N_0): + + def analytic_solution_helper(x, tau, x_0): + sqrt_tau = np.sqrt(tau) + result = float( + (1 - tau) + * 1 + / (x * np.sqrt(tau)) + * special.ive(1, 2 * x / x_0 * sqrt_tau) # pylint: disable=no-member + * np.exp(-(1 + tau - 2 * sqrt_tau) * x / x_0) + ) + return result + + tau = 1 - np.exp(-N_0 * b * x_0 * t) + + if isinstance(x, np.ndarray): + func = np.vectorize(lambda i: analytic_solution_helper(x[int(i)], tau, x_0)) + result = np.fromfunction(func, x.shape, dtype=float) + else: + result = analytic_solution_helper(x, tau, x_0) + + return result diff --git a/PySDM/physics/collision_kernel/hydrodynamic.py b/PySDM/physics/collision_kernel/hydrodynamic.py new file mode 100644 index 0000000000..e86164ffea --- /dev/null +++ b/PySDM/physics/collision_kernel/hydrodynamic.py @@ -0,0 +1,8 @@ +""" +Formulae supporting `PySDM.dynamics.collisions.collision_kernels.Hydrodynamic` +""" + + +class Hydrodynamic: + def __init__(self, _): + pass diff --git a/PySDM/physics/collision_kernel/linear.py b/PySDM/physics/collision_kernel/linear.py new file mode 100644 index 0000000000..04b62665f0 --- /dev/null +++ b/PySDM/physics/collision_kernel/linear.py @@ -0,0 +1,8 @@ +""" +Formulae supporting `PySDM.dynamics.collisions.collision_kernels.Linear` +""" + + +class Linear: + def __init__(self, _): + pass diff --git a/PySDM/physics/collision_kernel/long1974.py b/PySDM/physics/collision_kernel/long1974.py new file mode 100644 index 0000000000..86770a2180 --- /dev/null +++ b/PySDM/physics/collision_kernel/long1974.py @@ -0,0 +1,8 @@ +""" +Formulae supporting `PySDM.dynamics.collisions.collision_kernels.Long1974` +""" + + +class Long1974: + def __init__(self, _): + pass diff --git a/PySDM/physics/collision_kernel/neglect.py b/PySDM/physics/collision_kernel/neglect.py new file mode 100644 index 0000000000..1e5621aa96 --- /dev/null +++ b/PySDM/physics/collision_kernel/neglect.py @@ -0,0 +1,12 @@ +""" +disable liquid liquid collisions +""" + + +class Neglect: + def __init__(self, _): + pass + + @staticmethod + def collision_kernel(const, output, is_first_in_pair): + pass diff --git a/PySDM/physics/collision_kernel/simple_geometric.py b/PySDM/physics/collision_kernel/simple_geometric.py new file mode 100644 index 0000000000..f2473f3402 --- /dev/null +++ b/PySDM/physics/collision_kernel/simple_geometric.py @@ -0,0 +1,8 @@ +""" +Formulae supporting `PySDM.dynamics.collisions.collision_kernels.SimpleGeometric` +""" + + +class SimpleGeometric: + def __init__(self, _): + pass diff --git a/PySDM/physics/constants_defaults.py b/PySDM/physics/constants_defaults.py index 70e29479c3..06900546fd 100644 --- a/PySDM/physics/constants_defaults.py +++ b/PySDM/physics/constants_defaults.py @@ -370,6 +370,22 @@ J_HOM = np.nan """ constant ice nucleation rates """ +GOLOVIN_b = 5e3 * si.s +""" Golovin collision kernel constant """ + +CONSTANTK_a = np.nan +""" CONSTANTK collision kernel constant """ + +LINEAR_a = np.nan +LINEAR_b = np.nan +""" LINEAR collision kernel constants """ + +GEOMETRIC_collection_efficiency = 1 +""" GEOMETRIC collision kernel constants """ + +SIMPLE_GEOMETRIC_c = 1 +""" SIMPLE_GEOMETRIC collision kernel constants """ + STRAUB_E_D1 = 0.04 * si.cm """ [Straub et al. 2010](https://doi.org/10.1175/2009JAS3175.1) """ STRAUB_MU2 = 0.095 * si.cm diff --git a/docs/markdown/pysdm_landing.md b/docs/markdown/pysdm_landing.md index 0d88caacc1..b684e66d33 100644 --- a/docs/markdown/pysdm_landing.md +++ b/docs/markdown/pysdm_landing.md @@ -146,16 +146,17 @@ Instantiation of the [``Particulator``](https://open-atmos.github.io/PySDM/PySDM ```Julia Builder = pyimport("PySDM").Builder +Formulae = pyimport("PySDM").Formulae Box = pyimport("PySDM.environments").Box Coalescence = pyimport("PySDM.dynamics").Coalescence -Golovin = pyimport("PySDM.dynamics.collisions.collision_kernels").Golovin CPU = pyimport("PySDM.backends").CPU ParticleVolumeVersusRadiusLogarithmSpectrum = pyimport("PySDM.products").ParticleVolumeVersusRadiusLogarithmSpectrum radius_bins_edges = 10 .^ range(log10(10*si.um), log10(5e3*si.um), length=32) env = Box(dt=1 * si.s, dv=1e6 * si.m^3) -builder = Builder(n_sd=n_sd, backend=CPU(), environment=env, dynamics=(Coalescence(collision_kernel=Golovin(b=1.5e3 / si.s)),)) +formulae = Formulae(collision_kernel_liquid_liquid="Golovin", constants=Dict("GOLOVIN_b" => 1.5e3 / si.s)) +builder = Builder(n_sd=n_sd, backend=CPU(formulae), environment=env, dynamics=(Coalescence(),)) products = [ParticleVolumeVersusRadiusLogarithmSpectrum(radius_bins_edges=radius_bins_edges, name="dv/dlnr")] particulator = builder.build(attributes, products) ``` @@ -165,16 +166,17 @@ particulator = builder.build(attributes, products) ```Matlab Builder = py.importlib.import_module('PySDM').Builder; +Formulae = py.importlib.import_module('PySDM').Formulae; Box = py.importlib.import_module('PySDM.environments').Box; Coalescence = py.importlib.import_module('PySDM.dynamics').Coalescence; -Golovin = py.importlib.import_module('PySDM.dynamics.collisions.collision_kernels').Golovin; CPU = py.importlib.import_module('PySDM.backends').CPU; ParticleVolumeVersusRadiusLogarithmSpectrum = py.importlib.import_module('PySDM.products').ParticleVolumeVersusRadiusLogarithmSpectrum; radius_bins_edges = logspace(log10(10 * si.um), log10(5e3 * si.um), 32); env = Box(pyargs('dt', 1 * si.s, 'dv', 1e6 * si.m ^ 3)); -builder = Builder(pyargs('n_sd', int32(n_sd), 'backend', CPU(), 'environment', env, 'dynamics', py.tuple({Coalescence(pyargs('collision_kernel', Golovin(1.5e3 / si.s)))}))); +formulae = Formulae(pyargs('collision_kernel_liquid_liquid', 'Golovin', 'constants', py.dict(pyargs('GOLOVIN_b', 1.5e3 / si.s)))); +builder = Builder(pyargs('n_sd', int32(n_sd), 'backend', CPU(formulae), 'environment', env, 'dynamics', py.tuple({Coalescence()}))); products = py.list({ ParticleVolumeVersusRadiusLogarithmSpectrum(pyargs( ... 'radius_bins_edges', py.numpy.array(radius_bins_edges), ... 'name', 'dv/dlnr' ... @@ -188,16 +190,17 @@ particulator = builder.build(attributes, products); ```Python import numpy as np from PySDM import Builder +from PySDM import Formulae from PySDM.environments import Box from PySDM.dynamics import Coalescence -from PySDM.dynamics.collisions.collision_kernels import Golovin from PySDM.backends import CPU from PySDM.products import ParticleVolumeVersusRadiusLogarithmSpectrum radius_bins_edges = np.logspace(np.log10(10 * si.um), np.log10(5e3 * si.um), num=32) env = Box(dt=1 * si.s, dv=1e6 * si.m ** 3) -builder = Builder(n_sd=n_sd, backend=CPU(), environment=env, dynamics=(Coalescence(collision_kernel=Golovin(b=1.5e3 / si.s)),)) +formulae = Formulae(collision_kernel_liquid_liquid='Golovin', constants={'GOLOVIN_b': 1.5e3 / si.s}) +builder = Builder(n_sd=n_sd, backend=CPU(), environment=env, dynamics=(Coalescence(),)) products = [ParticleVolumeVersusRadiusLogarithmSpectrum(radius_bins_edges=radius_bins_edges, name='dv/dlnr')] particulator = builder.build(attributes, products) ``` @@ -308,9 +311,6 @@ The component submodules used to create this simulation are visualized below: DT[dt :float] -->|passed as arg to| ENV_INIT DV[dv :float] -->|passed as arg to| ENV_INIT ENV[":Box"] -->|passed as arg to| BUILDER_INIT - B["b: float"] --->|passed as arg to| KERNEL_INIT(["Golovin.__init__()"]) - KERNEL_INIT -->|instantiates| KERNEL - KERNEL[collision_kernel: Golovin] -->|passed as arg to| COAL_INIT(["Coalesncence.__init__()"]) COAL_INIT -->|instantiates| COAL PRODUCTS[products: list] ----->|passed as arg to| BUILDER_BUILD NORM_FACTOR[norm_factor: float]-->|passed as arg to| EXP_INIT @@ -468,7 +468,6 @@ AmbientThermodynamics = py.importlib.import_module('PySDM.dynamics').AmbientTher Condensation = py.importlib.import_module('PySDM.dynamics').Condensation; Parcel = py.importlib.import_module('PySDM.environments').Parcel; Builder = py.importlib.import_module('PySDM').Builder; -Formulae = py.importlib.import_module('PySDM').Formulae; products = py.importlib.import_module('PySDM.products'); env = Parcel(pyargs( ... @@ -486,9 +485,8 @@ output_interval = 4; output_points = 40; n_sd = 256; -formulae = Formulae(); builder = Builder(pyargs( ... - 'backend', CPU(formulae), ... + 'backend', CPU(), ... 'n_sd', int32(n_sd), ... 'environment', env, ... 'dynamics', py.tuple({AmbientThermodynamics(), Condensation()}) ... diff --git a/examples/PySDM_examples/Berry_1967/settings.py b/examples/PySDM_examples/Berry_1967/settings.py index bd52b857a4..23a98bfae1 100644 --- a/examples/PySDM_examples/Berry_1967/settings.py +++ b/examples/PySDM_examples/Berry_1967/settings.py @@ -4,7 +4,6 @@ from pystrict import strict from PySDM import Formulae -from PySDM.dynamics.collisions import collision_kernels from PySDM.initialisation import spectra from PySDM.physics import si @@ -15,10 +14,16 @@ def __init__( self, steps: Optional[list] = None, terminal_velocity_variant: str = "GunnKinzer1949", + kernel: str = None, + constants: dict = None, ): steps = steps or [200 * i for i in range(10)] - self.formulae = Formulae(terminal_velocity=terminal_velocity_variant) + self.formulae = Formulae( + terminal_velocity=terminal_velocity_variant, + collision_kernel_liquid_liquid=kernel, + constants=constants or {}, + ) self.init_x_min = self.formulae.trivia.volume(radius=3.94 * si.micrometre) self.init_x_max = self.formulae.trivia.volume(radius=25 * si.micrometres) @@ -34,7 +39,6 @@ def __init__( self.adaptive = False self.seed = 44 self._steps = steps - self.kernel = collision_kernels.Geometric(collection_efficiency=1) self.spectrum = spectra.Exponential(norm_factor=self.norm_factor, scale=self.X0) # Note 220 instead of 200 for smoothing diff --git a/examples/PySDM_examples/Shima_et_al_2009/example.py b/examples/PySDM_examples/Shima_et_al_2009/example.py index 795b112774..92f0c2d4c8 100644 --- a/examples/PySDM_examples/Shima_et_al_2009/example.py +++ b/examples/PySDM_examples/Shima_et_al_2009/example.py @@ -19,9 +19,7 @@ def run(settings, backend=CPU, observers=()): n_sd=settings.n_sd, backend=backend(formulae=settings.formulae), environment=env, - dynamics=( - Coalescence(collision_kernel=settings.kernel, adaptive=settings.adaptive), - ), + dynamics=(Coalescence(adaptive=settings.adaptive),), ) attributes = {} sampling = ConstantMultiplicity(settings.spectrum) diff --git a/examples/PySDM_examples/Shima_et_al_2009/settings.py b/examples/PySDM_examples/Shima_et_al_2009/settings.py index 5e0d524bae..0a85419bfd 100644 --- a/examples/PySDM_examples/Shima_et_al_2009/settings.py +++ b/examples/PySDM_examples/Shima_et_al_2009/settings.py @@ -4,7 +4,6 @@ from pystrict import strict from PySDM import Formulae -from PySDM.dynamics.collisions.collision_kernels import Golovin from PySDM.initialisation import spectra from PySDM.physics import si @@ -13,7 +12,10 @@ class Settings: def __init__(self, steps: Optional[list] = None): steps = steps or [0, 1200, 2400, 3600] - self.formulae = Formulae() + self.formulae = Formulae( + collision_kernel_liquid_liquid="Golovin", + constants={"GOLOVIN_b": 1.5e3 / si.s}, + ) self.n_sd = 2**13 self.n_part = 2**23 / si.metre**3 self.X0 = self.formulae.trivia.volume(radius=30.531 * si.micrometres) @@ -24,7 +26,6 @@ def __init__(self, steps: Optional[list] = None): self.adaptive = False self.seed = 44 self.steps = steps - self.kernel = Golovin(b=1.5e3 / si.second) self.spectrum = spectra.Exponential(norm_factor=self.norm_factor, scale=self.X0) self.radius_bins_edges = np.logspace( np.log10(10 * si.um), np.log10(5e3 * si.um), num=128, endpoint=True diff --git a/examples/PySDM_examples/Shima_et_al_2009/spectrum_plotter.py b/examples/PySDM_examples/Shima_et_al_2009/spectrum_plotter.py index f989cec287..5f51c4cc9f 100644 --- a/examples/PySDM_examples/Shima_et_al_2009/spectrum_plotter.py +++ b/examples/PySDM_examples/Shima_et_al_2009/spectrum_plotter.py @@ -6,6 +6,7 @@ from PySDM_examples.Shima_et_al_2009.error_measure import error_measure from PySDM.physics.constants import si +from PySDM.physics.collision_kernel import golovin _matplotlib_version_3_3_3 = version.parse("3.3.0") _matplotlib_version_actual = version.parse(matplotlib.__version__) @@ -96,8 +97,12 @@ def plot_analytic_solution(self, settings, t, spectrum, title): else: def analytic_solution(x): - return settings.norm_factor * settings.kernel.analytic_solution( - x=x, t=t, x_0=settings.X0, N_0=settings.n_part + return settings.norm_factor * golovin.analytic_solution( + x=x, + t=t, + x_0=settings.X0, + N_0=settings.n_part, + b=self.settings.formulae.constants.GOLOVIN_b, ) volume_bins_edges = self.settings.formulae.trivia.volume( diff --git a/examples/PySDM_examples/Srivastava_1982/simulation.py b/examples/PySDM_examples/Srivastava_1982/simulation.py index e2e643d45c..73b5529cd6 100644 --- a/examples/PySDM_examples/Srivastava_1982/simulation.py +++ b/examples/PySDM_examples/Srivastava_1982/simulation.py @@ -8,9 +8,19 @@ class Simulation: def __init__( - self, n_steps, settings, collision_dynamic=None, double_precision=True + self, + *, + n_steps, + settings, + collision_kernel, + constants_overrides, + collision_dynamic=None, + double_precision=True, ): self.collision_dynamic = collision_dynamic + self.collision_kernel = collision_kernel + self.constants_overrides = constants_overrides + self.settings = settings self.n_steps = n_steps @@ -26,8 +36,9 @@ def build(self, n_sd, seed, products): builder = Builder( backend=self.settings.backend_class( formulae=Formulae( - constants={"rho_w": self.settings.rho}, + constants={"rho_w": self.settings.rho, **self.constants_overrides}, fragmentation_function="ConstantMass", + collision_kernel_liquid_liquid=self.collision_kernel, seed=seed, ), double_precision=self.double_precision, diff --git a/tests/smoke_tests/box/berry_1967/test_coalescence.py b/tests/smoke_tests/box/berry_1967/test_coalescence.py index 285022f3cb..01d741f5f5 100644 --- a/tests/smoke_tests/box/berry_1967/test_coalescence.py +++ b/tests/smoke_tests/box/berry_1967/test_coalescence.py @@ -8,26 +8,20 @@ from PySDM.backends import ThrustRTC from PySDM.builder import Builder from PySDM.dynamics import Coalescence -from PySDM.dynamics.collisions.collision_kernels import ( - Electric, - Geometric, - Golovin, - Hydrodynamic, -) from PySDM.environments import Box from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity @pytest.mark.parametrize("croupier", ("local", "global")) @pytest.mark.parametrize("adaptive", (True, False)) -@pytest.mark.parametrize("kernel", (Geometric(), Electric(), Hydrodynamic())) +@pytest.mark.parametrize("kernel", ("Geometric", "Electric", "Hydrodynamic")) def test_coalescence(backend_class, kernel, croupier, adaptive): if backend_class == ThrustRTC and croupier == "local": pytest.skip("TODO #358") if backend_class == ThrustRTC and adaptive and croupier == "global": pytest.skip("TODO #329") # Arrange - s = Settings() + s = Settings(kernel=kernel) s.formulae.seed = 0 steps = [0, 800] @@ -36,9 +30,7 @@ def test_coalescence(backend_class, kernel, croupier, adaptive): n_sd=s.n_sd, backend=backend_class(formulae=s.formulae), environment=env, - dynamics=( - Coalescence(collision_kernel=kernel, croupier=croupier, adaptive=adaptive), - ), + dynamics=(Coalescence(croupier=croupier, adaptive=adaptive),), ) attributes = {} attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity( @@ -63,8 +55,7 @@ def test_coalescence(backend_class, kernel, croupier, adaptive): @pytest.mark.xfail(struct.calcsize("P") * 8 == 32, reason="32 bit", strict=False) def test_coalescence_2_sd(backend_class): # Arrange - s = Settings() - s.kernel = Golovin(b=1.5e12) + s = Settings(kernel="Golovin", constants={"GOLOVIN_b": 1.5e12}) s.formulae.seed = 0 steps = [0, 200] s.n_sd = 2 @@ -74,7 +65,7 @@ def test_coalescence_2_sd(backend_class): n_sd=s.n_sd, backend=backend_class(formulae=s.formulae), environment=env, - dynamics=(Coalescence(collision_kernel=s.kernel, adaptive=False),), + dynamics=(Coalescence(adaptive=False),), ) attributes = {} attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity( diff --git a/tests/smoke_tests/box/shima_et_al_2009/test_lwc_constant.py b/tests/smoke_tests/box/shima_et_al_2009/test_lwc_constant.py index 1763e269e2..ccf13896fb 100644 --- a/tests/smoke_tests/box/shima_et_al_2009/test_lwc_constant.py +++ b/tests/smoke_tests/box/shima_et_al_2009/test_lwc_constant.py @@ -8,7 +8,6 @@ from PySDM.backends import ThrustRTC from PySDM.builder import Builder from PySDM.dynamics import Coalescence -from PySDM.dynamics.collisions.collision_kernels import Golovin from PySDM.environments import Box from PySDM.formulae import Formulae from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity @@ -41,7 +40,11 @@ def test_lwc_constant(backend_class, croupier, adaptive): if backend_class == ThrustRTC and adaptive and croupier == "global": # TODO #329 pytest.skip() # Arrange - formulae = Formulae(seed=256) + formulae = Formulae( + seed=256, + collision_kernel_liquid_liquid="Golovin", + constants={"GOLOVIN_b": 1.5e3 / si.s}, + ) n_sd = 2**14 steps = [0, 100, 200] X0 = formulae.trivia.volume(radius=30.531e-6) @@ -51,7 +54,6 @@ def test_lwc_constant(backend_class, croupier, adaptive): norm_factor = n_part * dv rho = 1000 * si.kilogram / si.metre**3 - kernel = Golovin(b=1.5e3) # [s-1] spectrum = Exponential(norm_factor=norm_factor, scale=X0) env = Box(dt=dt, dv=dv) @@ -59,9 +61,7 @@ def test_lwc_constant(backend_class, croupier, adaptive): n_sd=n_sd, backend=backend_class(formulae=formulae), environment=env, - dynamics=( - Coalescence(collision_kernel=kernel, croupier=croupier, adaptive=adaptive), - ), + dynamics=(Coalescence(croupier=croupier, adaptive=adaptive),), ) attributes = {} diff --git a/tests/smoke_tests/box/srivastava_1982/test_eq_10.py b/tests/smoke_tests/box/srivastava_1982/test_eq_10.py index eff8629407..3d34b4b3ab 100644 --- a/tests/smoke_tests/box/srivastava_1982/test_eq_10.py +++ b/tests/smoke_tests/box/srivastava_1982/test_eq_10.py @@ -13,7 +13,6 @@ from PySDM_examples.Srivastava_1982.simulation import Simulation from PySDM.dynamics import Coalescence -from PySDM.dynamics.collisions.collision_kernels import ConstantK from PySDM.physics import si ASSERT_PROD = SimProducts.Computed.mean_drop_volume_total_volume_ratio.name @@ -39,9 +38,9 @@ def test_pysdm_coalescence_is_close_to_analytic_coalescence( simulation = Simulation( n_steps=N_STEPS, settings=settings, - collision_dynamic=Coalescence( - collision_kernel=ConstantK(a=settings.srivastava_c) - ), + collision_dynamic=Coalescence(), + collision_kernel="ConstantK", + constants_overrides={"CONSTANTK_a": settings.srivastava_c}, ) x = np.arange(N_STEPS + 1, dtype=float) diff --git a/tests/unit_tests/attributes/test_diffusional_growth_mass_change.py b/tests/unit_tests/attributes/test_diffusional_growth_mass_change.py index 754fdcbc31..0dff42b3a1 100644 --- a/tests/unit_tests/attributes/test_diffusional_growth_mass_change.py +++ b/tests/unit_tests/attributes/test_diffusional_growth_mass_change.py @@ -40,7 +40,6 @@ def test_if_collision(backend_class): backend_class, dynamics=( Collision( - collision_kernel=np.nan, breakup_efficiency=np.nan, coalescence_efficiency=np.nan, fragmentation_function=np.nan, diff --git a/tests/unit_tests/attributes/test_fall_velocity.py b/tests/unit_tests/attributes/test_fall_velocity.py index 7846b17329..e1a37ff60b 100644 --- a/tests/unit_tests/attributes/test_fall_velocity.py +++ b/tests/unit_tests/attributes/test_fall_velocity.py @@ -6,10 +6,10 @@ import numpy as np import pytest +from PySDM import Formulae from PySDM.attributes.physics import RelativeFallVelocity, TerminalVelocity from PySDM.builder import Builder from PySDM.dynamics import Coalescence, RelaxedVelocity -from PySDM.dynamics.collisions.collision_kernels.constantK import ConstantK from PySDM.environments.box import Box from PySDM.physics import si @@ -80,18 +80,23 @@ def test_fall_velocity_calculation(default_attributes, backend_instance): ) @staticmethod - def test_conservation_of_momentum(default_attributes, backend_instance): + def test_conservation_of_momentum(default_attributes, backend_class): """ Test that conservation of momentum holds when many super-droplets coalesce """ env = Box(dt=1, dv=1) builder = Builder( n_sd=len(default_attributes["multiplicity"]), - backend=backend_instance, + backend=backend_class( + Formulae( + collision_kernel_liquid_liquid="ConstantK", + constants={"CONSTANTK_a": 1}, + ) + ), environment=env, dynamics=( RelaxedVelocity(), - Coalescence(collision_kernel=ConstantK(a=1), adaptive=False), + Coalescence(adaptive=False), ), ) diff --git a/tests/unit_tests/dynamics/collisions/conftest.py b/tests/unit_tests/dynamics/collisions/conftest.py index 710e950d75..0f38677875 100644 --- a/tests/unit_tests/dynamics/collisions/conftest.py +++ b/tests/unit_tests/dynamics/collisions/conftest.py @@ -48,7 +48,6 @@ def get_dummy_particulator_and_coalescence( particulator = DummyParticulator(backend, n_sd=n_length) particulator.environment = environment or Box(dv=1, dt=DEFAULTS.dt_coal_range[1]) coalescence = Coalescence( - collision_kernel=StubKernel(particulator.backend), optimized_random=optimized_random, substeps=substeps, adaptive=False, diff --git a/tests/unit_tests/dynamics/collisions/test_collision_kernel_call.py b/tests/unit_tests/dynamics/collisions/test_collision_kernel_call.py new file mode 100644 index 0000000000..b54498bedd --- /dev/null +++ b/tests/unit_tests/dynamics/collisions/test_collision_kernel_call.py @@ -0,0 +1,61 @@ +""" +test for collision kernel basics +""" + +import numpy as np +import pytest +from PySDM import Builder +from PySDM.formulae import Formulae, _choices +from PySDM.physics import collision_kernel +from PySDM.environments import Box +from PySDM.dynamics import Coalescence +from PySDM.products import ( + LiquidWaterContent, + ParticleConcentration, +) +from PySDM.physics import si + + +class TestCollisionKernel: # pylint: disable=too-few-public-methods + @staticmethod + @pytest.mark.parametrize("variant", _choices(collision_kernel)) + def test_collision_kernel_call(backend_class, variant): + if ( + variant == "Linear" + or variant == "ConstantK" + and backend_class.__name__ == "ThrustRTC" + ): + pytest.skip() + + formulae = Formulae( + collision_kernel=variant + "+Neglect+Neglect", + seed=666, + constants={ + "CONSTANTK_a": 1, + "LINEAR_a": 1, + "LINEAR_b": 1, + }, + ) + env = Box(dt=1 * si.s, dv=1 * si.m**3) + builder = Builder( + n_sd=2, + backend=backend_class(formulae=formulae), + environment=env, + dynamics=([Coalescence()]), + ) + particulator = builder.build( + products=(LiquidWaterContent(name="qc"), ParticleConcentration(name="nc")), + attributes={ + "multiplicity": np.ones(builder.particulator.n_sd), + "volume": np.ones(builder.particulator.n_sd) * si.mm**3, + }, + ) + + # act + nc_initial = particulator.products["nc"].get() + particulator.run(steps=1) + + # assert + qc, nc = (particulator.products[k].get() for k in ("qc", "nc")) + assert np.isfinite([qc, nc]).all() and qc > 0 and nc > 0 + assert nc <= nc_initial diff --git a/tests/unit_tests/dynamics/collisions/test_kernels.py b/tests/unit_tests/dynamics/collisions/test_kernels.py index 6010f1540a..24f4ad48dc 100644 --- a/tests/unit_tests/dynamics/collisions/test_kernels.py +++ b/tests/unit_tests/dynamics/collisions/test_kernels.py @@ -5,12 +5,12 @@ from PySDM import Builder from PySDM.backends import CPU from PySDM.dynamics.collisions.collision_kernels import ( - Golovin, SimpleGeometric, Long1974, ) from PySDM.environments import Box from PySDM.formulae import Formulae +from PySDM.physics.collision_kernel.golovin import analytic_solution class TestKernels: @@ -20,14 +20,13 @@ class TestKernels: ) def test_golovin_analytic_solution_underflow(x): # Arrange - formulae = Formulae() + formulae = Formulae(collision_kernel_liquid_liquid="Golovin") b = 1.5e3 x_0 = formulae.trivia.volume(radius=30.531e-6) N_0 = 2**23 - sut = Golovin(b) # Act - value = sut.analytic_solution(x=x, t=1200, x_0=x_0, N_0=N_0) + value = analytic_solution(b=b, x=x, t=1200, x_0=x_0, N_0=N_0) # Assert assert np.all(np.isfinite(value)) diff --git a/tests/unit_tests/dynamics/collisions/test_sdm_breakup.py b/tests/unit_tests/dynamics/collisions/test_sdm_breakup.py index ef041b31cc..bd00d93aa8 100644 --- a/tests/unit_tests/dynamics/collisions/test_sdm_breakup.py +++ b/tests/unit_tests/dynamics/collisions/test_sdm_breakup.py @@ -12,7 +12,6 @@ from PySDM.dynamics.collisions.breakup_fragmentations import AlwaysN from PySDM.dynamics.collisions.coalescence_efficiencies import ConstEc from PySDM.dynamics.collisions.collision import DEFAULTS, Collision -from PySDM.dynamics.collisions.collision_kernels import ConstantK, Geometric from PySDM.environments import Box from PySDM.initialisation import spectra from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity @@ -46,7 +45,6 @@ def test_nonadaptive_same_results_regardless_of_dt(dt, backend_class): "volume": np.asarray([100 * si.um**3, 100 * si.um**3]), } breakup = Breakup( - collision_kernel=ConstantK(1 * si.cm**3 / si.s), fragmentation_function=AlwaysN(4), adaptive=False, warn_overflows=False, @@ -57,7 +55,13 @@ def test_nonadaptive_same_results_regardless_of_dt(dt, backend_class): env = Box(dv=1 * si.cm**3, dt=dt) builder = Builder( n_sd, - backend_class(Formulae(fragmentation_function="AlwaysN")), + backend_class( + Formulae( + fragmentation_function="AlwaysN", + collision_kernel_liquid_liquid="ConstantK", + constants={"CONSTANTK_a": 1 * si.cm**3 / si.s}, + ) + ), environment=env, dynamics=(breakup,), ) @@ -785,9 +789,7 @@ def test_noninteger_fragments( }[flag]() @staticmethod - def test_nonnegative_even_if_overflow( - backend=CPU(), - ): # pylint: disable=too-many-locals + def test_nonnegative_even_if_overflow(): # pylint: disable=too-many-locals n_sd = 2**5 dv = 1 * si.m**3 @@ -804,16 +806,15 @@ def test_nonnegative_even_if_overflow( mu = Trivia.volume(const, radius=100 * si.um) fragmentation = breakup_fragmentations.Exponential(scale=mu) - kernel = Geometric() coal_eff = ConstEc(Ec=0.01) break_eff = ConstEb(Eb=1.0) breakup = Collision( - collision_kernel=kernel, breakup_efficiency=break_eff, coalescence_efficiency=coal_eff, fragmentation_function=fragmentation, warn_overflows=False, ) + backend = CPU(Formulae(collision_kernel_liquid_liquid="Geometric")) builder = Builder( n_sd=n_sd, backend=backend, environment=env, dynamics=(breakup,) ) diff --git a/tests/unit_tests/products/test_collision_rates.py b/tests/unit_tests/products/test_collision_rates.py index 4ff1248b9b..d262129e04 100644 --- a/tests/unit_tests/products/test_collision_rates.py +++ b/tests/unit_tests/products/test_collision_rates.py @@ -8,7 +8,6 @@ from PySDM.dynamics.collisions.breakup_efficiencies import ConstEb from PySDM.dynamics.collisions.breakup_fragmentations import AlwaysN from PySDM.dynamics.collisions.coalescence_efficiencies import ConstEc -from PySDM.dynamics.collisions.collision_kernels import ConstantK from PySDM.environments import Box from PySDM.physics import si from PySDM.products import ( @@ -69,11 +68,8 @@ class TestCollisionProducts: }, ], ) - def test_individual_dynamics_rates_nonadaptive(params, backend_instance): - if ( - backend_instance.__class__.__name__ == "ThrustRTC" - and params["enable_breakup"] - ): + def test_individual_dynamics_rates_nonadaptive(params, backend_class): + if backend_class.__name__ == "ThrustRTC" and params["enable_breakup"]: pytest.skip("# TODO #744") # Arrange @@ -82,8 +78,12 @@ def test_individual_dynamics_rates_nonadaptive(params, backend_instance): env = Box(**ENV_ARGS) - dynamic, products = _get_dynamics_and_products(params, adaptive=False) - builder = Builder(n_sd, backend_instance, environment=env, dynamics=(dynamic,)) + formulae, dynamic, products = _get_formulae_dynamics_and_products( + params, adaptive=False + ) + builder = Builder( + n_sd, backend_class(formulae), environment=env, dynamics=(dynamic,) + ) particulator = builder.build( attributes={ @@ -138,11 +138,13 @@ def test_no_collision_deficits_when_adaptive(params, n_init, backend_class=CPU): n_sd = len(n_init) env = Box(**ENV_ARGS) - dynamic, _ = _get_dynamics_and_products( + formulae, dynamic, _ = _get_formulae_dynamics_and_products( params, adaptive=True, kernel_a=1e4 * si.cm**3 / si.s ) - builder = Builder(n_sd, backend_class(), environment=env, dynamics=(dynamic,)) + builder = Builder( + n_sd, backend_class(formulae), environment=env, dynamics=(dynamic,) + ) particulator = builder.build( attributes={ @@ -185,8 +187,12 @@ def test_breakup_deficits_when_adaptive(params, backend_class=CPU): n_sd = len(n_init) env = Box(**ENV_ARGS) - dynamic, _ = _get_dynamics_and_products(params, adaptive=True) - builder = Builder(n_sd, backend_class(), environment=env, dynamics=(dynamic,)) + formulae, dynamic, _ = _get_formulae_dynamics_and_products( + params, adaptive=True + ) + builder = Builder( + n_sd, backend_class(formulae), environment=env, dynamics=(dynamic,) + ) particulator = builder.build( attributes={ @@ -232,13 +238,14 @@ def test_no_breakup_deficits_when_while_loop(params, backend_class=CPU): n_sd = len(n_init) env = Box(**ENV_ARGS) - dynamic, _ = _get_dynamics_and_products( + formulae, dynamic, _ = _get_formulae_dynamics_and_products( params, adaptive=True, kernel_a=1e4 * si.cm**3 / si.s ) + formulae.handle_all_breakups = True builder = Builder( n_sd, - backend_class(Formulae(handle_all_breakups=True)), + backend_class(formulae), environment=env, dynamics=(dynamic,), ) @@ -293,8 +300,12 @@ def test_rate_sums_single_cell(params, backend_class=CPU): n_sd = len(n_init) env = Box(**ENV_ARGS) - dynamic, products = _get_dynamics_and_products(params, adaptive=False) - builder = Builder(n_sd, backend_class(), environment=env, dynamics=(dynamic,)) + formulae, dynamic, products = _get_formulae_dynamics_and_products( + params, adaptive=False + ) + builder = Builder( + n_sd, backend_class(formulae), environment=env, dynamics=(dynamic,) + ) particulator = builder.build( attributes={ @@ -313,12 +324,15 @@ def test_rate_sums_single_cell(params, backend_class=CPU): assert (particulator.products["cr"].get()[0] == rhs_sum).all() -def _get_dynamics_and_products(params, adaptive, kernel_a=1e6 * si.cm**3 / si.s): - kernel = ConstantK(a=kernel_a) +def _get_formulae_dynamics_and_products( + params, adaptive, kernel_a=1e6 * si.cm**3 / si.s +): + formulae = Formulae( + collision_kernel_liquid_liquid="ConstantK", constants={"CONSTANTK_a": kernel_a} + ) if params["enable_breakup"]: if params["enable_coalescence"]: dynamic = Collision( - collision_kernel=kernel, coalescence_efficiency=ConstEc(Ec=params["Ec"]), breakup_efficiency=ConstEb(Eb=params["Eb"]), fragmentation_function=AlwaysN(n=params["nf"]), @@ -333,7 +347,6 @@ def _get_dynamics_and_products(params, adaptive, kernel_a=1e6 * si.cm**3 / si.s) ) else: dynamic = Breakup( - collision_kernel=kernel, fragmentation_function=AlwaysN(n=params["nf"]), adaptive=adaptive, ) @@ -345,7 +358,6 @@ def _get_dynamics_and_products(params, adaptive, kernel_a=1e6 * si.cm**3 / si.s) ) else: dynamic = Coalescence( - collision_kernel=kernel, coalescence_efficiency=ConstEc(Ec=1.0), adaptive=adaptive, ) @@ -354,7 +366,7 @@ def _get_dynamics_and_products(params, adaptive, kernel_a=1e6 * si.cm**3 / si.s) CollisionRateDeficitPerGridbox(name="crd"), CoalescenceRatePerGridbox(name="cor"), ) - return (dynamic, products) + return formulae, dynamic, products def _get_product_component_sums(params, products):