From 10a05b1b21b747bb28be63bde9a2b2d1a13c378c Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Tue, 17 Oct 2023 09:36:55 +0200 Subject: [PATCH] mass instead of volume as base attribute (incl. introduction of physics.particle_shape_and_density) (#1147) Co-authored-by: Oleksii Bulenok Co-authored-by: Agnieszka Makulska Co-authored-by: derlk --- .gitignore | 5 +- PySDM/attributes/impl/mapper.py | 2 + PySDM/attributes/physics/__init__.py | 1 + PySDM/attributes/physics/mass.py | 11 + .../physics/relative_fall_velocity.py | 8 +- PySDM/attributes/physics/volume.py | 12 +- .../impl_common/freezing_attributes.py | 4 +- .../impl_numba/methods/collisions_methods.py | 197 +++++++++--------- .../methods/condensation_methods.py | 60 +++--- .../impl_numba/methods/freezing_methods.py | 31 ++- .../impl_numba/methods/physics_methods.py | 22 ++ .../scipy_ode_condensation_solver.py | 30 +-- .../methods/collisions_methods.py | 126 +++++------ .../methods/condensation_methods.py | 84 ++++---- .../methods/freezing_methods.py | 26 ++- .../methods/physics_methods.py | 30 +++ .../test_helpers/cpp2python.py | 10 +- PySDM/builder.py | 18 +- .../breakup_fragmentations/__init__.py | 7 +- .../breakup_fragmentations/always_n.py | 11 +- .../breakup_fragmentations/constant_mass.py | 17 ++ .../breakup_fragmentations/constant_size.py | 18 -- .../breakup_fragmentations/expon_frag.py | 12 ++ .../breakup_fragmentations/exponential.py | 15 +- .../breakup_fragmentations/feingold1988.py | 18 +- .../breakup_fragmentations/gaussian.py | 15 +- .../breakup_fragmentations/impl/__init__.py | 4 + .../impl/volume_based.py | 26 +++ .../breakup_fragmentations/lowlist82.py | 14 +- .../breakup_fragmentations/slams.py | 13 +- .../breakup_fragmentations/straub2010.py | 15 +- .../coalescence_efficiencies/lowlist1982.py | 26 +-- PySDM/dynamics/collisions/collision.py | 8 +- PySDM/dynamics/freezing.py | 8 +- PySDM/dynamics/relaxed_velocity.py | 12 +- PySDM/formulae.py | 4 + PySDM/initialisation/init_fall_momenta.py | 6 +- PySDM/particulator.py | 18 +- PySDM/physics/__init__.py | 1 + PySDM/physics/constants.py | 2 + .../fragmentation_function/__init__.py | 5 +- .../{constant_size.py => constant_mass.py} | 4 +- .../fragmentation_function/expon_frag.py | 11 +- .../fragmentation_function/exponential.py | 8 + .../{feingold1988frag.py => feingold1988.py} | 4 +- .../particle_shape_and_density/__init__.py | 6 + .../liquid_spheres.py | 25 +++ .../mixed_phase_spheres.py | 31 +++ .../porous_spheroids.py | 10 + PySDM/physics/trivia.py | 8 +- PySDM/products/aqueous_chemistry/acidity.py | 1 + PySDM/products/freezing/ice_water_content.py | 6 +- PySDM/products/impl/moment_product.py | 4 +- .../size_spectral/effective_radius.py | 2 + PySDM/products/size_spectral/mean_radius.py | 1 + .../size_spectral/particle_concentration.py | 10 +- .../total_particle_concentration.py | 2 +- .../Alpert_and_Knopf_2016/simulation.py | 2 + .../Arabas_et_al_2023/aida.ipynb | 1 + .../Arabas_et_al_2023/fig_11.ipynb | 1 + .../Arabas_et_al_2023/fig_convergence.ipynb | 1 + .../Arabas_et_al_2023/make_particulator.py | 1 + .../PySDM_examples/Berry_1967/settings.py | 2 +- .../Bieli_et_al_2022/settings.py | 6 +- .../Bulenok_2023_MasterThesis/setups.py | 4 +- .../Niedermeier_et_al_2014/fig_2.ipynb | 1 + .../Niedermeier_et_al_2014/settings.py | 4 +- .../PySDM_examples/Srivastava_1982/example.py | 4 +- .../Srivastava_1982/simulation.py | 2 +- .../deJong_Mackay_et_al_2023/figs_6_7_8.ipynb | 14 +- .../deJong_Mackay_et_al_2023/settings1D.py | 5 +- .../deJong_Mackay_et_al_2023/settings_0D.py | 12 +- .../deJong_Mackay_et_al_2023/simulation_ss.py | 3 +- .../box/bieli_et_al_2022/test_moments.py | 2 +- .../test_fig_6.py | 7 +- .../test_fig_7.py | 16 +- .../test_fig_8.py | 2 +- .../arabas_et_al_2015/test_freezing.py | 1 + .../test_conservation.py | 10 +- .../test_temperature_profile.py | 1 + .../parcel/yang_et_al_2018/test_just_do_it.py | 4 +- .../attributes/test_fall_velocity.py | 6 +- .../backends/test_freezing_methods.py | 5 +- .../backends/test_physics_methods.py | 70 ++++++- tests/unit_tests/dummy_particulator.py | 2 +- .../collisions/test_fragmentations.py | 182 +++++++++------- .../dynamics/collisions/test_sdm_breakup.py | 146 ++++++++----- .../collisions/test_sdm_single_cell.py | 34 ++- .../impl/test_particle_attributes.py | 10 +- .../test_particle_shape_and_density.py | 40 ++++ 90 files changed, 1082 insertions(+), 604 deletions(-) create mode 100644 PySDM/attributes/physics/mass.py create mode 100644 PySDM/dynamics/collisions/breakup_fragmentations/constant_mass.py delete mode 100644 PySDM/dynamics/collisions/breakup_fragmentations/constant_size.py create mode 100644 PySDM/dynamics/collisions/breakup_fragmentations/expon_frag.py create mode 100644 PySDM/dynamics/collisions/breakup_fragmentations/impl/__init__.py create mode 100644 PySDM/dynamics/collisions/breakup_fragmentations/impl/volume_based.py rename PySDM/physics/fragmentation_function/{constant_size.py => constant_mass.py} (57%) create mode 100644 PySDM/physics/fragmentation_function/exponential.py rename PySDM/physics/fragmentation_function/{feingold1988frag.py => feingold1988.py} (60%) create mode 100644 PySDM/physics/particle_shape_and_density/__init__.py create mode 100644 PySDM/physics/particle_shape_and_density/liquid_spheres.py create mode 100644 PySDM/physics/particle_shape_and_density/mixed_phase_spheres.py create mode 100644 PySDM/physics/particle_shape_and_density/porous_spheroids.py create mode 100644 tests/unit_tests/physics/test_particle_shape_and_density.py diff --git a/.gitignore b/.gitignore index e731acff4..3997eb71f 100644 --- a/.gitignore +++ b/.gitignore @@ -134,4 +134,7 @@ dmypy.json .code-workspace # Mac stuff -**.DS_Store \ No newline at end of file +**.DS_Store + +#Jetbrains +.idea/ \ No newline at end of file diff --git a/PySDM/attributes/impl/mapper.py b/PySDM/attributes/impl/mapper.py index 1590601ad..4d18c06a9 100644 --- a/PySDM/attributes/impl/mapper.py +++ b/PySDM/attributes/impl/mapper.py @@ -24,6 +24,7 @@ Temperature, TerminalVelocity, Volume, + WaterMass, WetToCriticalVolumeRatio, ) from PySDM.attributes.physics.critical_supersaturation import CriticalSupersaturation @@ -100,6 +101,7 @@ "critical supersaturation": lambda _, __: CriticalSupersaturation, "equilibrium supersaturation": lambda _, __: EquilibriumSupersaturation, "wet to critical volume ratio": lambda _, __: WetToCriticalVolumeRatio, + "water mass": lambda _, __: WaterMass, } diff --git a/PySDM/attributes/physics/__init__.py b/PySDM/attributes/physics/__init__.py index fec980a51..09d981571 100644 --- a/PySDM/attributes/physics/__init__.py +++ b/PySDM/attributes/physics/__init__.py @@ -8,6 +8,7 @@ from .dry_volume import DryVolume from .equilibrium_supersaturation import EquilibriumSupersaturation from .heat import Heat +from .mass import WaterMass from .multiplicities import Multiplicities from .radius import Radius, SquareRootOfRadius from .relative_fall_velocity import RelativeFallMomentum, RelativeFallVelocity diff --git a/PySDM/attributes/physics/mass.py b/PySDM/attributes/physics/mass.py new file mode 100644 index 000000000..44acc5e87 --- /dev/null +++ b/PySDM/attributes/physics/mass.py @@ -0,0 +1,11 @@ +""" +particle mass, derived attribute for coalescence +in simulation involving mixed-phase clouds, positive values correspond to +liquid water and negative values to ice +""" +from PySDM.attributes.impl import ExtensiveAttribute + + +class WaterMass(ExtensiveAttribute): + def __init__(self, builder): + super().__init__(builder, name="water mass") diff --git a/PySDM/attributes/physics/relative_fall_velocity.py b/PySDM/attributes/physics/relative_fall_velocity.py index 7a49f3296..b55fb80b8 100644 --- a/PySDM/attributes/physics/relative_fall_velocity.py +++ b/PySDM/attributes/physics/relative_fall_velocity.py @@ -14,15 +14,13 @@ def __init__(self, builder): class RelativeFallVelocity(DerivedAttribute): def __init__(self, builder): self.momentum = builder.get_attribute("relative fall momentum") - self.volume = builder.get_attribute("volume") - self.rho_w = builder.formulae.constants.rho_w # TODO #798 + self.water_mass = builder.get_attribute("water mass") super().__init__( builder, name="relative fall velocity", - dependencies=(self.momentum, self.volume), + dependencies=(self.momentum, self.water_mass), ) def recalculate(self): - self.data.ratio(self.momentum.get(), self.volume.get()) - self.data[:] *= 1 / self.rho_w # TODO #798 + self.data.ratio(self.momentum.get(), self.water_mass.get()) diff --git a/PySDM/attributes/physics/volume.py b/PySDM/attributes/physics/volume.py index 5110aa3f4..3c9c4f786 100644 --- a/PySDM/attributes/physics/volume.py +++ b/PySDM/attributes/physics/volume.py @@ -1,11 +1,15 @@ """ -particle (wet) volume, key attribute for coalescence +particle volume (derived from water mass); in simulation involving mixed-phase clouds, positive values correspond to liquid water and negative values to ice """ -from PySDM.attributes.impl.extensive_attribute import ExtensiveAttribute +from PySDM.attributes.impl import DerivedAttribute -class Volume(ExtensiveAttribute): +class Volume(DerivedAttribute): def __init__(self, builder): - super().__init__(builder, name="volume") + self.water_mass = builder.get_attribute("water mass") + super().__init__(builder, name="volume", dependencies=(self.water_mass,)) + + def recalculate(self): + self.particulator.backend.volume_of_water_mass(self.data, self.water_mass.get()) diff --git a/PySDM/backends/impl_common/freezing_attributes.py b/PySDM/backends/impl_common/freezing_attributes.py index f695b976a..71942eaf5 100644 --- a/PySDM/backends/impl_common/freezing_attributes.py +++ b/PySDM/backends/impl_common/freezing_attributes.py @@ -5,7 +5,7 @@ class SingularAttributes( - namedtuple("SingularAttributes", ("freezing_temperature", "wet_volume")) + namedtuple("SingularAttributes", ("freezing_temperature", "water_mass")) ): """groups attributes required in singular regime""" @@ -13,7 +13,7 @@ class SingularAttributes( class TimeDependentAttributes( - namedtuple("TimeDependentAttributes", ("immersed_surface_area", "wet_volume")) + namedtuple("TimeDependentAttributes", ("immersed_surface_area", "water_mass")) ): """groups attributes required in time-dependent regime""" diff --git a/PySDM/backends/impl_numba/methods/collisions_methods.py b/PySDM/backends/impl_numba/methods/collisions_methods.py index e1cbd0f07..91896c40c 100644 --- a/PySDM/backends/impl_numba/methods/collisions_methods.py +++ b/PySDM/backends/impl_numba/methods/collisions_methods.py @@ -60,13 +60,15 @@ def coalesce( # pylint: disable=too-many-arguments @numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) def compute_transfer_multiplicities( - gamma, j, k, multiplicity, volume, fragment_size_i, max_multiplicity + gamma, j, k, multiplicity, particle_mass, fragment_mass_i, max_multiplicity ): # pylint: disable=too-many-arguments overflow_flag = False gamma_j_k = 0 take_from_j_test = multiplicity[k] take_from_j = 0 - new_mult_k_test = ((volume[j] + volume[k]) / fragment_size_i) * multiplicity[k] + new_mult_k_test = ( + (particle_mass[j] + particle_mass[k]) / fragment_mass_i + ) * multiplicity[k] new_mult_k = multiplicity[k] for m in range(int(gamma)): # check for overflow of multiplicity @@ -84,7 +86,7 @@ def compute_transfer_multiplicities( take_from_j_test += new_mult_k_test new_mult_k_test = ( - new_mult_k_test * (volume[j] / fragment_size_i) + new_mult_k_test + new_mult_k_test * (particle_mass[j] / fragment_mass_i) + new_mult_k_test ) return take_from_j, new_mult_k, gamma_j_k, overflow_flag @@ -138,20 +140,20 @@ def break_up( # pylint: disable=too-many-arguments,c,too-many-locals multiplicity, gamma, attributes, - fragment_size, + fragment_mass, max_multiplicity, breakup_rate, breakup_rate_deficit, warn_overflows, - volume, + particle_mass, ): # breakup0 guarantees take_from_j <= multiplicity[j] take_from_j, new_mult_k, gamma_j_k, overflow_flag = compute_transfer_multiplicities( gamma[i], j, k, multiplicity, - volume, - fragment_size[i], + particle_mass, + fragment_mass[i], max_multiplicity, ) gamma_deficit = gamma[i] - gamma_j_k @@ -181,19 +183,23 @@ def break_up_while( multiplicity, gamma, attributes, - fragment_size, + fragment_mass, max_multiplicity, breakup_rate, breakup_rate_deficit, warn_overflows, - volume, + particle_mass, ): # pylint: disable=too-many-arguments,unused-argument,too-many-locals gamma_deficit = gamma[i] overflow_flag = False while gamma_deficit > 0: if multiplicity[k] == multiplicity[j]: take_from_j = multiplicity[j] - new_mult_k = (volume[j] + volume[k]) / fragment_size[i] * multiplicity[k] + new_mult_k = ( + (particle_mass[j] + particle_mass[k]) + / fragment_mass[i] + * multiplicity[k] + ) # check for overflow if new_mult_k > max_multiplicity: @@ -215,8 +221,8 @@ def break_up_while( j, k, multiplicity, - volume, - fragment_size[i], + particle_mass, + fragment_mass[i], max_multiplicity, ) @@ -317,7 +323,7 @@ def __collision_coalescence_breakup_body( rand, Ec, Eb, - fragment_size, + fragment_mass, healthy, cell_id, coalescence_rate, @@ -326,7 +332,7 @@ def __collision_coalescence_breakup_body( is_first_in_pair, max_multiplicity, warn_overflows, - volume, + particle_mass, ): # pylint: disable=not-an-iterable,too-many-nested-blocks,too-many-locals for i in numba.prange(length // 2): @@ -357,12 +363,12 @@ def __collision_coalescence_breakup_body( multiplicity, gamma, attributes, - fragment_size, + fragment_mass, max_multiplicity, breakup_rate, breakup_rate_deficit, warn_overflows, - volume, + particle_mass, ) flag_zero_multiplicity(j, k, multiplicity, healthy) @@ -387,10 +393,10 @@ def __ll82_coalescence_check_body(*, Ec, dl): @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) def __straub_fragmentation_body( - *, CW, gam, ds, v_max, frag_size, rand, Nr1, Nr2, Nr3, Nr4, Nrt, d34 + *, CW, gam, ds, v_max, frag_volume, rand, Nr1, Nr2, Nr3, Nr4, Nrt, d34 ): # pylint: disable=too-many-arguments,too-many-locals for i in numba.prange( # pylint: disable=not-an-iterable - len(frag_size) + len(frag_volume) ): straub_Nr(i, Nr1, Nr2, Nr3, Nr4, Nrt, CW, gam) sigma1 = straub_sigma1(CW[i]) @@ -416,23 +422,24 @@ def __straub_fragmentation_body( Nr4, ) Nrt[i] = Nr1[i] + Nr2[i] + Nr3[i] + Nr4[i] + if Nrt[i] == 0.0: - frag_size[i] = 0.0 + diameter = 0.0 else: if rand[i] < Nr1[i] / Nrt[i]: X = rand[i] * Nrt[i] / Nr1[i] lnarg = mu1 + np.sqrt(2) * sigma1 * straub_erfinv(X) - frag_size[i] = np.exp(lnarg) + diameter = np.exp(lnarg) elif rand[i] < (Nr2[i] + Nr1[i]) / Nrt[i]: X = (rand[i] * Nrt[i] - Nr1[i]) / Nr2[i] - frag_size[i] = mu2 + np.sqrt(2) * sigma2 * straub_erfinv(X) + diameter = mu2 + np.sqrt(2) * sigma2 * straub_erfinv(X) elif rand[i] < (Nr3[i] + Nr2[i] + Nr1[i]) / Nrt[i]: X = (rand[i] * Nrt[i] - Nr1[i] - Nr2[i]) / Nr3[i] - frag_size[i] = mu3 + np.sqrt(2) * sigma3 * straub_erfinv(X) + diameter = mu3 + np.sqrt(2) * sigma3 * straub_erfinv(X) else: - frag_size[i] = d34[i] + diameter = d34[i] - frag_size[i] = frag_size[i] ** 3 * const.PI / 6 + frag_volume[i] = diameter**3 * const.PI / 6 self.__straub_fragmentation_body = __straub_fragmentation_body elif self.formulae.fragmentation_function.__name__ == "LowList1982Nf": @@ -447,15 +454,15 @@ def __straub_fragmentation_body( @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) def __ll82_fragmentation_body( - *, CKE, W, W2, St, ds, dl, dcoal, frag_size, rand, Rf, Rs, Rd, tol + *, CKE, W, W2, St, ds, dl, dcoal, frag_volume, rand, Rf, Rs, Rd, tol ): # pylint: disable=too-many-branches,too-many-locals,too-many-statements for i in numba.prange( # pylint: disable=not-an-iterable - len(frag_size) + len(frag_volume) ): if dl[i] <= 0.4e-3: - frag_size[i] = dcoal[i] ** 3 * const.PI / 6 + frag_volume[i] = dcoal[i] ** 3 * const.PI / 6 elif ds[i] == 0.0 or dl[i] == 0.0: - frag_size[i] = 1e-18 + frag_volume[i] = 1e-18 else: ll82_Nr(i, Rf, Rs, Rd, CKE, W, W2) if rand[i] <= Rf[i]: # filament breakup @@ -469,20 +476,20 @@ def __ll82_fragmentation_body( rand[i] = rand[i] / Rf[i] if rand[i] <= H1 / Hsum: X = max(rand[i] * Hsum / H1, tol) - frag_size[i] = mu1 + np.sqrt(2) * sigma1 * ll82_erfinv( - 2 * X - 1 - ) + frag_volume[i] = mu1 + np.sqrt( + 2 + ) * sigma1 * ll82_erfinv(2 * X - 1) elif rand[i] <= (H1 + H2) / Hsum: X = (rand[i] * Hsum - H1) / H2 - frag_size[i] = mu2 + np.sqrt(2) * sigma2 * ll82_erfinv( - 2 * X - 1 - ) + frag_volume[i] = mu2 + np.sqrt( + 2 + ) * sigma2 * ll82_erfinv(2 * X - 1) else: X = min((rand[i] * Hsum - H1 - H2) / H3, 1.0 - tol) lnarg = mu3 + np.sqrt(2) * sigma3 * ll82_erfinv( 2 * X - 1 ) - frag_size[i] = np.exp(lnarg) + frag_volume[i] = np.exp(lnarg) elif rand[i] <= Rf[i] + Rs[i]: # sheet breakup (H1, mu1, sigma1) = ll82_params_s1(dl[i], ds[i], dcoal[i]) @@ -493,15 +500,15 @@ def __ll82_fragmentation_body( rand[i] = (rand[i] - Rf[i]) / (Rs[i]) if rand[i] <= H1 / Hsum: X = max(rand[i] * Hsum / H1, tol) - frag_size[i] = mu1 + np.sqrt(2) * sigma1 * ll82_erfinv( - 2 * X - 1 - ) + frag_volume[i] = mu1 + np.sqrt( + 2 + ) * sigma1 * ll82_erfinv(2 * X - 1) else: X = min((rand[i] * Hsum - H1) / H2, 1.0 - tol) lnarg = mu2 + np.sqrt(2) * sigma2 * ll82_erfinv( 2 * X - 1 ) - frag_size[i] = np.exp(lnarg) + frag_volume[i] = np.exp(lnarg) else: # disk breakup (H1, mu1, sigma1) = ll82_params_d1( @@ -513,20 +520,20 @@ def __ll82_fragmentation_body( rand[i] = (rand[i] - Rf[i] - Rs[i]) / Rd[i] if rand[i] <= H1 / Hsum: X = max(rand[i] * Hsum / H1, tol) - frag_size[i] = mu1 + np.sqrt(2) * sigma1 * ll82_erfinv( - 2 * X - 1 - ) + frag_volume[i] = mu1 + np.sqrt( + 2 + ) * sigma1 * ll82_erfinv(2 * X - 1) else: X = min((rand[i] * Hsum - H1) / H2, 1 - tol) lnarg = mu2 + np.sqrt(2) * sigma2 * ll82_erfinv( 2 * X - 1 ) - frag_size[i] = np.exp(lnarg) + frag_volume[i] = np.exp(lnarg) - frag_size[i] = ( - frag_size[i] * 0.01 + frag_volume[i] = ( + frag_volume[i] * 0.01 ) # diameter in cm; convert to m - frag_size[i] = frag_size[i] ** 3 * const.PI / 6 + frag_volume[i] = frag_volume[i] ** 3 * const.PI / 6 self.__ll82_fragmentation_body = __ll82_fragmentation_body elif self.formulae.fragmentation_function.__name__ == "Gaussian": @@ -534,26 +541,26 @@ def __ll82_fragmentation_body( @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) def __gauss_fragmentation_body( - *, mu, sigma, frag_size, rand + *, mu, sigma, frag_volume, rand ): # pylint: disable=too-many-arguments for i in numba.prange( # pylint: disable=not-an-iterable - len(frag_size) + len(frag_volume) ): - frag_size[i] = mu + sigma * erfinv_approx(rand[i]) + frag_volume[i] = mu + sigma * erfinv_approx(rand[i]) self.__gauss_fragmentation_body = __gauss_fragmentation_body - elif self.formulae.fragmentation_function.__name__ == "Feingold1988Frag": - feingold1988_frag_size = self.formulae.fragmentation_function.frag_size + elif self.formulae.fragmentation_function.__name__ == "Feingold1988": + feingold1988_frag_volume = self.formulae.fragmentation_function.frag_volume @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) # pylint: disable=too-many-arguments def __feingold1988_fragmentation_body( - *, scale, frag_size, x_plus_y, rand, fragtol + *, scale, frag_volume, x_plus_y, rand, fragtol ): for i in numba.prange( # pylint: disable=not-an-iterable - len(frag_size) + len(frag_volume) ): - frag_size[i] = feingold1988_frag_size( + frag_volume[i] = feingold1988_frag_volume( scale, rand[i], x_plus_y[i], fragtol ) @@ -706,7 +713,7 @@ def collision_coalescence_breakup( rand, Ec, Eb, - fragment_size, + fragment_mass, healthy, cell_id, coalescence_rate, @@ -714,7 +721,7 @@ def collision_coalescence_breakup( breakup_rate_deficit, is_first_in_pair, warn_overflows, - volume, + particle_mass, max_multiplicity, ): # pylint: disable=too-many-locals @@ -727,7 +734,7 @@ def collision_coalescence_breakup( rand=rand.data, Ec=Ec.data, Eb=Eb.data, - fragment_size=fragment_size.data, + fragment_mass=fragment_mass.data, healthy=healthy.data, cell_id=cell_id.data, coalescence_rate=coalescence_rate.data, @@ -736,31 +743,31 @@ def collision_coalescence_breakup( is_first_in_pair=is_first_in_pair.indicator.data, max_multiplicity=max_multiplicity, warn_overflows=warn_overflows, - volume=volume.data, + particle_mass=particle_mass.data, ) @staticmethod @numba.njit(**{**conf.JIT_FLAGS}) # pylint: disable=too-many-arguments - def __fragmentation_limiters(n_fragment, frag_size, vmin, nfmax, x_plus_y): - for i in numba.prange(len(frag_size)): # pylint: disable=not-an-iterable + def __fragmentation_limiters(n_fragment, frag_volume, vmin, nfmax, x_plus_y): + for i in numba.prange(len(frag_volume)): # pylint: disable=not-an-iterable if x_plus_y[i] == 0.0: - frag_size[i] = 0.0 + frag_volume[i] = 0.0 n_fragment[i] = 1.0 else: - if np.isnan(frag_size[i]) or frag_size[i] == 0.0: - frag_size[i] = x_plus_y[i] - frag_size[i] = min(frag_size[i], x_plus_y[i]) - if nfmax is not None and x_plus_y[i] / frag_size[i] > nfmax: - frag_size[i] = x_plus_y[i] / nfmax - elif frag_size[i] < vmin: - frag_size[i] = x_plus_y[i] - n_fragment[i] = x_plus_y[i] / frag_size[i] - - def fragmentation_limiters(self, *, n_fragment, frag_size, vmin, nfmax, x_plus_y): + if np.isnan(frag_volume[i]) or frag_volume[i] == 0.0: + frag_volume[i] = x_plus_y[i] + frag_volume[i] = min(frag_volume[i], x_plus_y[i]) + if nfmax is not None and x_plus_y[i] / frag_volume[i] > nfmax: + frag_volume[i] = x_plus_y[i] / nfmax + elif frag_volume[i] < vmin: + frag_volume[i] = x_plus_y[i] + n_fragment[i] = x_plus_y[i] / frag_volume[i] + + def fragmentation_limiters(self, *, n_fragment, frag_volume, vmin, nfmax, x_plus_y): self.__fragmentation_limiters( n_fragment=n_fragment.data, - frag_size=frag_size.data, + frag_volume=frag_volume.data, vmin=vmin, nfmax=nfmax, x_plus_y=x_plus_y.data, @@ -768,7 +775,7 @@ def fragmentation_limiters(self, *, n_fragment, frag_size, vmin, nfmax, x_plus_y @staticmethod @numba.njit(**{**conf.JIT_FLAGS}) - def __slams_fragmentation_body(n_fragment, frag_size, x_plus_y, probs, rand): + def __slams_fragmentation_body(n_fragment, frag_volume, x_plus_y, probs, rand): for i in numba.prange(len(n_fragment)): # pylint: disable=not-an-iterable probs[i] = 0.0 n_fragment[i] = 1 @@ -777,17 +784,17 @@ def __slams_fragmentation_body(n_fragment, frag_size, x_plus_y, probs, rand): if rand[i] < probs[i]: n_fragment[i] = n + 2 break - frag_size[i] = x_plus_y[i] / n_fragment[i] + frag_volume[i] = x_plus_y[i] / n_fragment[i] def slams_fragmentation( - self, n_fragment, frag_size, x_plus_y, probs, rand, vmin, nfmax + self, n_fragment, frag_volume, x_plus_y, probs, rand, vmin, nfmax ): # pylint: disable=too-many-arguments self.__slams_fragmentation_body( - n_fragment.data, frag_size.data, x_plus_y.data, probs.data, rand.data + n_fragment.data, frag_volume.data, x_plus_y.data, probs.data, rand.data ) self.__fragmentation_limiters( n_fragment=n_fragment.data, - frag_size=frag_size.data, + frag_volume=frag_volume.data, vmin=vmin, nfmax=nfmax, x_plus_y=x_plus_y.data, @@ -796,19 +803,19 @@ def slams_fragmentation( @staticmethod @numba.njit(**{**conf.JIT_FLAGS}) # pylint: disable=too-many-arguments - def __exp_fragmentation_body(*, scale, frag_size, rand, tol=1e-5): + def __exp_fragmentation_body(*, scale, frag_volume, rand, tol=1e-5): """ Exponential PDF """ - for i in numba.prange(len(frag_size)): # pylint: disable=not-an-iterable - frag_size[i] = -scale * np.log(max(1 - rand[i], tol)) + for i in numba.prange(len(frag_volume)): # pylint: disable=not-an-iterable + frag_volume[i] = -scale * np.log(max(1 - rand[i], tol)) def exp_fragmentation( self, *, n_fragment, scale, - frag_size, + frag_volume, x_plus_y, rand, vmin, @@ -817,13 +824,13 @@ def exp_fragmentation( ): self.__exp_fragmentation_body( scale=scale, - frag_size=frag_size.data, + frag_volume=frag_volume.data, rand=rand.data, tol=tol, ) self.__fragmentation_limiters( n_fragment=n_fragment.data, - frag_size=frag_size.data, + frag_volume=frag_volume.data, x_plus_y=x_plus_y.data, vmin=vmin, nfmax=nfmax, @@ -834,7 +841,7 @@ def feingold1988_fragmentation( *, n_fragment, scale, - frag_size, + frag_volume, x_plus_y, rand, fragtol, @@ -843,7 +850,7 @@ def feingold1988_fragmentation( ): self.__feingold1988_fragmentation_body( scale=scale, - frag_size=frag_size.data, + frag_volume=frag_volume.data, x_plus_y=x_plus_y.data, rand=rand.data, fragtol=fragtol, @@ -851,24 +858,24 @@ def feingold1988_fragmentation( self.__fragmentation_limiters( n_fragment=n_fragment.data, - frag_size=frag_size.data, + frag_volume=frag_volume.data, x_plus_y=x_plus_y.data, vmin=vmin, nfmax=nfmax, ) def gauss_fragmentation( - self, *, n_fragment, mu, sigma, frag_size, x_plus_y, rand, vmin, nfmax + self, *, n_fragment, mu, sigma, frag_volume, x_plus_y, rand, vmin, nfmax ): self.__gauss_fragmentation_body( mu=mu, sigma=sigma, - frag_size=frag_size.data, + frag_volume=frag_volume.data, rand=rand.data, ) self.__fragmentation_limiters( n_fragment=n_fragment.data, - frag_size=frag_size.data, + frag_volume=frag_volume.data, x_plus_y=x_plus_y.data, vmin=vmin, nfmax=nfmax, @@ -882,7 +889,7 @@ def straub_fragmentation( CW, gam, ds, - frag_size, + frag_volume, v_max, x_plus_y, rand, @@ -899,7 +906,7 @@ def straub_fragmentation( CW=CW.data, gam=gam.data, ds=ds.data, - frag_size=frag_size.data, + frag_volume=frag_volume.data, v_max=v_max.data, rand=rand.data, Nr1=Nr1.data, @@ -911,7 +918,7 @@ def straub_fragmentation( ) self.__fragmentation_limiters( n_fragment=n_fragment.data, - frag_size=frag_size.data, + frag_volume=frag_volume.data, x_plus_y=x_plus_y.data, vmin=vmin, nfmax=nfmax, @@ -929,7 +936,7 @@ def ll82_fragmentation( ds, dl, dcoal, - frag_size, + frag_volume, x_plus_y, rand, vmin, @@ -947,7 +954,7 @@ def ll82_fragmentation( ds=ds.data, dl=dl.data, dcoal=dcoal.data, - frag_size=frag_size.data, + frag_volume=frag_volume.data, rand=rand.data, Rf=Rf.data, Rs=Rs.data, @@ -956,7 +963,7 @@ def ll82_fragmentation( ) self.__fragmentation_limiters( n_fragment=n_fragment.data, - frag_size=frag_size.data, + frag_volume=frag_volume.data, x_plus_y=x_plus_y.data, vmin=vmin, nfmax=nfmax, diff --git a/PySDM/backends/impl_numba/methods/condensation_methods.py b/PySDM/backends/impl_numba/methods/condensation_methods.py index 490e66fd2..cf3927d0e 100644 --- a/PySDM/backends/impl_numba/methods/condensation_methods.py +++ b/PySDM/backends/impl_numba/methods/condensation_methods.py @@ -21,7 +21,7 @@ def condensation( solver, n_cell, cell_start_arg, - v, + water_mass, v_cr, multiplicity, vdry, @@ -50,7 +50,7 @@ def condensation( n_threads=n_threads, n_cell=n_cell, cell_start_arg=cell_start_arg.data, - v=v.data, + water_mass=water_mass.data, v_cr=v_cr.data, multiplicity=multiplicity.data, vdry=vdry.data, @@ -84,7 +84,7 @@ def _condensation( n_threads, n_cell, cell_start_arg, - v, + water_mass, v_cr, multiplicity, vdry, @@ -138,7 +138,7 @@ def _condensation( n_ripening, RH_max_in_cell, ) = solver( - v, + water_mass, v_cr, multiplicity, vdry, @@ -257,7 +257,7 @@ def make_step_impl( ): @numba.njit(**jit_flags) def step_impl( # pylint: disable=too-many-arguments,too-many-locals - v, + water_mass, v_cr, multiplicity, vdry, @@ -277,7 +277,7 @@ def step_impl( # pylint: disable=too-many-arguments,too-many-locals fake, ): timestep /= n_substeps - ml_old = calculate_ml_old(v, multiplicity, cell_idx) + ml_old = calculate_ml_old(water_mass, multiplicity, cell_idx) count_activating, count_deactivating, count_ripening = 0, 0, 0 RH_max = 0 success = True @@ -309,7 +309,7 @@ def step_impl( # pylint: disable=too-many-arguments,too-many-locals T, p, RH, - v, + water_mass, v_cr, multiplicity, vdry, @@ -357,13 +357,13 @@ def step_impl( # pylint: disable=too-many-arguments,too-many-locals return step_impl @staticmethod - def make_calculate_ml_old(jit_flags, const): + def make_calculate_ml_old(jit_flags): @numba.njit(**jit_flags) - def calculate_ml_old(volume, multiplicity, cell_idx): + def calculate_ml_old(water_mass, multiplicity, cell_idx): result = 0 for drop in cell_idx: - if volume[drop] > 0: - result += multiplicity[drop] * volume[drop] * const.rho_w + if water_mass[drop] > 0: + result += multiplicity[drop] * water_mass[drop] return result return calculate_ml_old @@ -383,6 +383,8 @@ def make_calculate_ml_new( # pylint: disable=too-many-statements,too-many-local phys_lambdaD, phys_dk_D, phys_dk_K, + volume_to_mass, + mass_to_volume, within_tolerance, max_iters, RH_rtol, @@ -410,7 +412,7 @@ def calculate_ml_new( # pylint: disable=too-many-arguments,too-many-statements, T, p, RH, - v, + water_mass, v_cr, multiplicity, vdry, @@ -431,13 +433,14 @@ def calculate_ml_new( # pylint: disable=too-many-arguments,too-many-statements, lambdaK = phys_lambdaK(T, p) lambdaD = phys_lambdaD(DTp, T) for drop in cell_idx: - if v[drop] < 0: + v_drop = mass_to_volume(water_mass[drop]) + if v_drop < 0: continue - x_old = x(v[drop]) - r_old = radius(v[drop]) + x_old = x(v_drop) + r_old = radius(v_drop) x_insane = x(vdry[drop] / 100) rd3 = vdry[drop] / const.PI_4_3 - sgm = phys_sigma(T, v[drop], vdry[drop], f_org[drop]) + sgm = phys_sigma(T, v_drop, vdry[drop], f_org[drop]) RH_eq = phys_RH_eq(r_old, T, kappa[drop], rd3, sgm) if not within_tolerance(np.abs(RH - RH_eq), RH, RH_rtol): Dr = phys_dk_D(DTp, r_old, lambdaD) @@ -523,16 +526,17 @@ def calculate_ml_new( # pylint: disable=too-many-arguments,too-many-statements, else: x_new = x_old - v_new = volume_of_x(x_new) - result += multiplicity[drop] * v_new * const.rho_w + mass_new = volume_to_mass(volume_of_x(x_new)) + mass_cr = volume_to_mass(v_cr[drop]) + result += multiplicity[drop] * mass_new if not fake: - if v_new > v_cr[drop] and v_new > v[drop]: + if mass_new > mass_cr and mass_new > water_mass[drop]: n_activated_and_growing += multiplicity[drop] - if v_new > v_cr[drop] > v[drop]: + if mass_new > mass_cr > water_mass[drop]: n_activating += multiplicity[drop] - if v_new < v_cr[drop] < v[drop]: + if mass_new < mass_cr < water_mass[drop]: n_deactivating += multiplicity[drop] - v[drop] = v_new + water_mass[drop] = mass_new n_ripening = n_activated_and_growing if n_deactivating > 0 else 0 return result, success, n_activating, n_deactivating, n_ripening @@ -573,6 +577,8 @@ def make_condensation_solver( dx_dt=self.formulae.condensation_coordinate.dx_dt, volume=self.formulae.condensation_coordinate.volume, x=self.formulae.condensation_coordinate.x, + mass_to_volume=self.formulae.particle_shape_and_density.mass_to_volume, + volume_to_mass=self.formulae.particle_shape_and_density.volume_to_mass, timestep=timestep, dt_range=dt_range, adaptive=adaptive, @@ -608,6 +614,8 @@ def make_condensation_solver_impl( dx_dt, volume, x, + volume_to_mass, + mass_to_volume, timestep, dt_range, adaptive, @@ -623,7 +631,7 @@ def make_condensation_solver_impl( **{"parallel": False, "cache": False, "fastmath": fastmath}, } - calculate_ml_old = CondensationMethods.make_calculate_ml_old(jit_flags, const) + calculate_ml_old = CondensationMethods.make_calculate_ml_old(jit_flags) calculate_ml_new = CondensationMethods.make_calculate_ml_new( jit_flags=jit_flags, dx_dt=dx_dt, @@ -637,6 +645,8 @@ def make_condensation_solver_impl( phys_lambdaD=phys_lambdaD, phys_dk_D=phys_dk_D, phys_dk_K=phys_dk_K, + volume_to_mass=volume_to_mass, + mass_to_volume=mass_to_volume, within_tolerance=within_tolerance, max_iters=max_iters, RH_rtol=RH_rtol, @@ -670,7 +680,7 @@ def make_condensation_solver_impl( @numba.njit(**jit_flags) def solve( # pylint: disable=too-many-arguments - v, + water_mass, v_cr, multiplicity, vdry, @@ -690,7 +700,7 @@ def solve( # pylint: disable=too-many-arguments n_substeps, ): args = ( - v, + water_mass, v_cr, multiplicity, vdry, diff --git a/PySDM/backends/impl_numba/methods/freezing_methods.py b/PySDM/backends/impl_numba/methods/freezing_methods.py index 9b14e5127..e65a8c2f9 100644 --- a/PySDM/backends/impl_numba/methods/freezing_methods.py +++ b/PySDM/backends/impl_numba/methods/freezing_methods.py @@ -16,7 +16,6 @@ class FreezingMethods(BackendMethods): def __init__(self): BackendMethods.__init__(self) - const = self.formulae.constants unfrozen_and_saturated = self.formulae.trivia.unfrozen_and_saturated frozen_and_above_freezing_point = ( self.formulae.trivia.frozen_and_above_freezing_point @@ -25,15 +24,15 @@ def __init__(self): @numba.njit( **{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath, "parallel": False} ) - def _freeze(volume, i): - volume[i] = -1 * volume[i] * const.rho_w / const.rho_i + def _freeze(water_mass, i): + water_mass[i] = -1 * water_mass[i] # TODO #599: change thd (latent heat)! @numba.njit( **{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath, "parallel": False} ) - def _thaw(volume, i): - volume[i] = -1 * volume[i] * const.rho_i / const.rho_w + def _thaw(water_mass, i): + water_mass[i] = -1 * water_mass[i] # TODO #599: change thd (latent heat)! @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) @@ -45,16 +44,16 @@ def freeze_singular_body( if attributes.freezing_temperature[i] == 0: continue if thaw and frozen_and_above_freezing_point( - attributes.wet_volume[i], temperature[cell[i]] + attributes.water_mass[i], temperature[cell[i]] ): - _thaw(attributes.wet_volume, i) + _thaw(attributes.water_mass, i) elif ( unfrozen_and_saturated( - attributes.wet_volume[i], relative_humidity[cell[i]] + attributes.water_mass[i], relative_humidity[cell[i]] ) and temperature[cell[i]] <= attributes.freezing_temperature[i] ): - _freeze(attributes.wet_volume, i) + _freeze(attributes.water_mass, i) self.freeze_singular_body = freeze_singular_body @@ -73,17 +72,17 @@ def freeze_time_dependent_body( # pylint: disable=unused-argument,too-many-argu freezing_temperature, thaw, ): - n_sd = len(attributes.wet_volume) + n_sd = len(attributes.water_mass) for i in numba.prange(n_sd): # pylint: disable=not-an-iterable if attributes.immersed_surface_area[i] == 0: continue cell_id = cell[i] if thaw and frozen_and_above_freezing_point( - attributes.wet_volume[i], temperature[cell_id] + attributes.water_mass[i], temperature[cell_id] ): - _thaw(attributes.wet_volume, i) + _thaw(attributes.water_mass, i) elif unfrozen_and_saturated( - attributes.wet_volume[i], relative_humidity[cell_id] + attributes.water_mass[i], relative_humidity[cell_id] ): rate = j_het(a_w_ice[cell_id]) # TODO #594: this assumes constant T throughout timestep, can we do better? @@ -91,7 +90,7 @@ def freeze_time_dependent_body( # pylint: disable=unused-argument,too-many-argu -rate * attributes.immersed_surface_area[i] * timestep ) if rand[i] < prob: - _freeze(attributes.wet_volume, i) + _freeze(attributes.water_mass, i) # if record_freezing_temperature: # freezing_temperature[i] = temperature[cell_id] @@ -103,7 +102,7 @@ def freeze_singular( self.freeze_singular_body( SingularAttributes( freezing_temperature=attributes.freezing_temperature.data, - wet_volume=attributes.wet_volume.data, + water_mass=attributes.water_mass.data, ), temperature.data, relative_humidity.data, @@ -129,7 +128,7 @@ def freeze_time_dependent( rand.data, TimeDependentAttributes( immersed_surface_area=attributes.immersed_surface_area.data, - wet_volume=attributes.wet_volume.data, + water_mass=attributes.water_mass.data, ), timestep, cell.data, diff --git a/PySDM/backends/impl_numba/methods/physics_methods.py b/PySDM/backends/impl_numba/methods/physics_methods.py index 845af1c9c..e036bfeb1 100644 --- a/PySDM/backends/impl_numba/methods/physics_methods.py +++ b/PySDM/backends/impl_numba/methods/physics_methods.py @@ -20,6 +20,8 @@ def __init__(self): # pylint: disable=too-many-locals phys_sigma = self.formulae.surface_tension.sigma phys_volume = self.formulae.trivia.volume phys_r_cr = self.formulae.hygroscopicity.r_cr + phys_mass_to_volume = self.formulae.particle_shape_and_density.mass_to_volume + phys_volume_to_mass = self.formulae.particle_shape_and_density.volume_to_mass const = self.formulae.constants @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) @@ -80,6 +82,20 @@ def a_w_ice_body( self.a_w_ice_body = a_w_ice_body + @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) + def volume_of_mass(volume, mass): + for i in prange(volume.shape[0]): # pylint: disable=not-an-iterable + volume[i] = phys_mass_to_volume(mass[i]) + + self.volume_of_mass = volume_of_mass + + @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) + def mass_of_volume(mass, volume): + for i in prange(volume.shape[0]): # pylint: disable=not-an-iterable + mass[i] = phys_volume_to_mass(volume[i]) + + self.mass_of_volume = mass_of_volume + def temperature_pressure_RH( self, *, rhod, thd, water_vapour_mixing_ratio, T, p, RH ): @@ -119,3 +135,9 @@ def a_w_ice(self, *, T, p, RH, water_vapour_mixing_ratio, a_w_ice): water_vapour_mixing_ratio_in=water_vapour_mixing_ratio.data, a_w_ice_out=a_w_ice.data, ) + + def volume_of_water_mass(self, volume, mass): + self.volume_of_mass(volume.data, mass.data) + + def mass_of_water_volume(self, mass, volume): + self.mass_of_volume(mass.data, volume.data) diff --git a/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py b/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py index bb72583b7..842c57f7c 100644 --- a/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py +++ b/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py @@ -11,7 +11,7 @@ from PySDM.backends import Numba from PySDM.backends.impl_numba.conf import JIT_FLAGS -from PySDM.physics.constants_defaults import PI_4_3, T0, rho_w +from PySDM.physics.constants_defaults import PI_4_3, T0 idx_thd = 0 idx_x = 1 @@ -34,7 +34,7 @@ def _condensation( n_threads=1, n_cell=particulator.mesh.n_cell, cell_start_arg=particulator.attributes.cell_start.data, - v=particulator.attributes["volume"].data, + water_mass=particulator.attributes["water mass"].data, v_cr=None, multiplicity=particulator.attributes["multiplicity"].data, vdry=particulator.attributes["dry volume"].data, @@ -63,12 +63,15 @@ def _condensation( RH_max=RH_max.data, success=success.data, ) + particulator.attributes.mark_updated("water mass") @lru_cache() def _make_solve(formulae): # pylint: disable=too-many-statements,too-many-locals x = formulae.condensation_coordinate.x - volume = formulae.condensation_coordinate.volume + volume_of_x = formulae.condensation_coordinate.volume + volume_to_mass = formulae.particle_shape_and_density.volume_to_mass + mass_to_volume = formulae.particle_shape_and_density.mass_to_volume dx_dt = formulae.condensation_coordinate.dx_dt pvs_C = formulae.saturation_vapour_pressure.pvs_Celsius lv = formulae.latent_heat.lv @@ -89,7 +92,7 @@ def _make_solve(formulae): # pylint: disable=too-many-statements,too-many-local @numba.njit(**{**JIT_FLAGS, **{"parallel": False}}) def _liquid_water_mixing_ratio(n, x, m_d_mean): - return np.sum(n * volume(x)) * rho_w / m_d_mean + return np.sum(n * volume_to_mass(volume_of_x(x))) / m_d_mean @numba.njit(**{**JIT_FLAGS, **{"parallel": False}}) def _impl( # pylint: disable=too-many-arguments,too-many-locals @@ -115,7 +118,7 @@ def _impl( # pylint: disable=too-many-arguments,too-many-locals lambdaD = phys_lambdaD(DTp, T) lambdaK = phys_lambdaK(T, p) for i, x_i in enumerate(x): - v = volume(x_i) + v = volume_of_x(x_i) r = phys_radius(v) Dr = phys_diff_kin_D(DTp, r, lambdaD) Kr = phys_diff_kin_K(KTp, r, lambdaK) @@ -134,7 +137,7 @@ def _impl( # pylint: disable=too-many-arguments,too-many-locals ) d_water_vapour_mixing_ratio__dt = ( dot_water_vapour_mixing_ratio - - np.sum(n * volume(x) * dy_dt[idx_x:]) * rho_w / m_d_mean + - np.sum(n * volume_to_mass(volume_of_x(x)) * dy_dt[idx_x:]) / m_d_mean ) dy_dt[idx_thd] = dot_thd + phys_dthd_dt( rhod, thd, T, d_water_vapour_mixing_ratio__dt, lv @@ -192,8 +195,8 @@ def _odesys( # pylint: disable=too-many-arguments,too-many-locals return dy_dt def solve( # pylint: disable=too-many-arguments,too-many-locals - v, - _, + water_mass, + __, multiplicity, vdry, cell_idx, @@ -206,15 +209,15 @@ def solve( # pylint: disable=too-many-arguments,too-many-locals d_water_vapour_mixing_ratio__dt, drhod_dt, m_d_mean, - __, ___, - dt, ____, + dt, + _____, ): n_sd_in_cell = len(cell_idx) y0 = np.empty(n_sd_in_cell + idx_x) y0[idx_thd] = thd - y0[idx_x:] = x(v[cell_idx]) + y0[idx_x:] = x(mass_to_volume(water_mass[cell_idx])) total_water_mixing_ratio = ( water_vapour_mixing_ratio + _liquid_water_mixing_ratio(multiplicity[cell_idx], y0[idx_x:], m_d_mean) @@ -254,9 +257,8 @@ def solve( # pylint: disable=too-many-arguments,too-many-locals m_new = 0 for i in range(n_sd_in_cell): - v_new = volume(y1[idx_x + i]) - m_new += multiplicity[cell_idx[i]] * v_new * rho_w - v[cell_idx[i]] = v_new + water_mass[cell_idx[i]] = volume_to_mass(volume_of_x(y1[idx_x + i])) + m_new += multiplicity[cell_idx[i]] * water_mass[cell_idx[i]] return ( integ.success, diff --git a/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py b/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py index f4599ad2b..40a2b99da 100644 --- a/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py +++ b/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py @@ -77,8 +77,8 @@ int64_t j, int64_t k, VectorView multiplicity, - VectorView volume, - real_type fragment_size_i, + VectorView particle_mass, + real_type fragment_mass_i, int64_t max_multiplicity, real_type *take_from_j, @@ -87,7 +87,7 @@ real_type gamma_j_k = 0; real_type take_from_j_test = multiplicity[k]; take_from_j[0] = 0; - real_type new_mult_k_test = ((volume[j] + volume[k]) / fragment_size_i * multiplicity[k]); + real_type new_mult_k_test = ((particle_mass[j] + particle_mass[k]) / fragment_mass_i * multiplicity[k]); new_mult_k[0] = multiplicity[k]; for (auto m = 0; m < (int64_t) (gamma); m += 1) { @@ -104,7 +104,7 @@ gamma_j_k = m + 1; take_from_j_test += new_mult_k_test; - new_mult_k_test = new_mult_k_test * (volume[j] / fragment_size_i) + new_mult_k_test; + new_mult_k_test = new_mult_k_test * (particle_mass[j] / fragment_mass_i) + new_mult_k_test; } return gamma_j_k; } @@ -168,11 +168,11 @@ VectorView multiplicity, VectorView gamma, VectorView attributes, - VectorView fragment_size, + VectorView fragment_mass, int64_t max_multiplicity, VectorView breakup_rate, VectorView breakup_rate_deficit, - VectorView volume, + VectorView particle_mass, int64_t n_sd, int64_t n_attr ) { @@ -183,8 +183,8 @@ j, k, multiplicity, - volume, - fragment_size[i], + particle_mass, + fragment_mass[i], max_multiplicity, take_from_j, new_mult_k @@ -359,7 +359,7 @@ def __collision_coalescence_breakup_body(self): "rand", "Ec", "Eb", - "fragment_size", + "fragment_mass", "healthy", "cell_id", "coalescence_rate", @@ -367,7 +367,7 @@ def __collision_coalescence_breakup_body(self): "breakup_rate_deficit", "is_first_in_pair", "max_multiplicity", - "volume", + "particle_mass", "n_sd", "n_attr", ), @@ -407,11 +407,11 @@ def __collision_coalescence_breakup_body(self): multiplicity, gamma, attributes, - fragment_size, + fragment_mass, max_multiplicity, breakup_rate, breakup_rate_deficit, - volume, + particle_mass, n_sd, n_attr ); @@ -524,10 +524,10 @@ def __cell_id_body(self): @cached_property def __exp_fragmentation_body(self): return trtc.For( - param_names=("scale", "frag_size", "rand", "tol"), + param_names=("scale", "frag_volume", "rand", "tol"), name_iter="i", body=""" - frag_size[i] = -scale * log(max(1 - rand[i], tol)); + frag_volume[i] = -scale * log(max(1 - rand[i], tol)); """, ) @@ -536,7 +536,7 @@ def __fragmentation_limiters_body(self): return trtc.For( param_names=( "n_fragment", - "frag_size", + "frag_volume", "x_plus_y", "vmin", "nfmax", @@ -544,30 +544,30 @@ def __fragmentation_limiters_body(self): ), name_iter="i", body=""" - frag_size[i] = min(frag_size[i], x_plus_y[i]); + frag_volume[i] = min(frag_volume[i], x_plus_y[i]); if (nfmax_is_not_none) { - if (x_plus_y[i] / frag_size[i] > nfmax) { - frag_size[i] = x_plus_y[i] / nfmax; + if (x_plus_y[i] / frag_volume[i] > nfmax) { + frag_volume[i] = x_plus_y[i] / nfmax; } } - else if (frag_size[i] < vmin) { - frag_size[i] = x_plus_y[i]; + else if (frag_volume[i] < vmin) { + frag_volume[i] = x_plus_y[i]; } - else if (frag_size[i] == 0.0) { - frag_size[i] = x_plus_y[i]; + else if (frag_volume[i] == 0.0) { + frag_volume[i] = x_plus_y[i]; } - n_fragment[i] = x_plus_y[i] / frag_size[i]; + n_fragment[i] = x_plus_y[i] / frag_volume[i]; """, ) @cached_property def __gauss_fragmentation_body(self): return trtc.For( - param_names=("mu", "sigma", "frag_size", "rand"), + param_names=("mu", "sigma", "frag_volume", "rand"), name_iter="i", body=f""" - frag_size[i] = mu + sigma * {self.formulae.trivia.erfinv_approx.c_inline( + frag_volume[i] = mu + sigma * {self.formulae.trivia.erfinv_approx.c_inline( c="rand[i]" )}; """.replace( @@ -578,7 +578,7 @@ def __gauss_fragmentation_body(self): @cached_property def __slams_fragmentation_body(self): return trtc.For( - param_names=("n_fragment", "frag_size", "x_plus_y", "probs", "rand"), + param_names=("n_fragment", "frag_volume", "x_plus_y", "probs", "rand"), name_iter="i", body=""" probs[i] = 0.0; @@ -591,17 +591,17 @@ def __slams_fragmentation_body(self): break; } } - frag_size[i] = x_plus_y[i] / n_fragment[i]; + frag_volume[i] = x_plus_y[i] / n_fragment[i]; """, ) @cached_property def __feingold1988_fragmentation_body(self): return trtc.For( - param_names=("scale", "frag_size", "x_plus_y", "rand", "fragtol"), + param_names=("scale", "frag_volume", "x_plus_y", "rand", "fragtol"), name_iter="i", body=f""" - frag_size[i] = {self.formulae.fragmentation_function.frag_size.c_inline( + frag_volume[i] = {self.formulae.fragmentation_function.frag_volume.c_inline( scale="scale", rand="rand[i]", x_plus_y="x_plus_y[i]", @@ -656,7 +656,7 @@ def __straub_fragmentation_body(self): "CW", "gam", "ds", - "frag_size", + "frag_volume", "v_max", "rand", "Nr1", @@ -683,25 +683,25 @@ def __straub_fragmentation_body(self): auto lnarg = mu1 + sqrt(2.0) * sigma1 * {self.formulae.trivia.erfinv_approx.c_inline( c="X" )}; - frag_size[i] = exp(lnarg); + frag_volume[i] = exp(lnarg); }} else if (rand[i] < (Nr2[i] + Nr1[i]) / Nrt[i]) {{ auto X = (rand[i] * Nrt[i] - Nr1[i]) / Nr2[i]; - frag_size[i] = mu2 + sqrt(2.0) * sigma2 * {self.formulae.trivia.erfinv_approx.c_inline( + frag_volume[i] = mu2 + sqrt(2.0) * sigma2 * {self.formulae.trivia.erfinv_approx.c_inline( c="X" )}; }} else if (rand[i] < (Nr3[i] + Nr2[i] + Nr1[i]) / Nrt[i]) {{ auto X = (rand[i] * Nrt[i] - Nr1[i] - Nr2[i]) / Nr3[i]; - frag_size[i] = mu3 + sqrt(2.0) * sigma3 * {self.formulae.trivia.erfinv_approx.c_inline( + frag_volume[i] = mu3 + sqrt(2.0) * sigma3 * {self.formulae.trivia.erfinv_approx.c_inline( c="X" )}; }} else {{ - frag_size[i] = d34[i]; + frag_volume[i] = d34[i]; }} - frag_size[i] = pow(frag_size[i], 3) * {const.PI} / 6; + frag_volume[i] = pow(frag_volume[i], 3) * {const.PI} / 6; """.replace( "real_type", self._get_c_type() ), @@ -824,7 +824,7 @@ def collision_coalescence_breakup( # pylint: disable=unused-argument,too-many-l rand, Ec, Eb, - fragment_size, + fragment_mass, healthy, cell_id, coalescence_rate, @@ -832,7 +832,7 @@ def collision_coalescence_breakup( # pylint: disable=unused-argument,too-many-l breakup_rate_deficit, is_first_in_pair, warn_overflows, - volume, + particle_mass, max_multiplicity, ): if warn_overflows: @@ -851,7 +851,7 @@ def collision_coalescence_breakup( # pylint: disable=unused-argument,too-many-l rand.data, Ec.data, Eb.data, - fragment_size.data, + fragment_mass.data, healthy.data, cell_id.data, coalescence_rate.data, @@ -859,7 +859,7 @@ def collision_coalescence_breakup( # pylint: disable=unused-argument,too-many-l breakup_rate_deficit.data, is_first_in_pair.indicator.data, trtc.DVInt64(max_multiplicity), - volume.data, + particle_mass.data, n_sd, n_attr, ), @@ -950,7 +950,7 @@ def exp_fragmentation( *, n_fragment, scale, - frag_size, + frag_volume, x_plus_y, rand, vmin, @@ -958,20 +958,20 @@ def exp_fragmentation( tol=1e-5, ): self.__exp_fragmentation_body.launch_n( - n=len(frag_size), + n=len(frag_volume), args=( self._get_floating_point(scale), - frag_size.data, + frag_volume.data, rand.data, self._get_floating_point(tol), ), ) self.__fragmentation_limiters_body.launch_n( - n=len(frag_size), + n=len(frag_volume), args=( n_fragment.data, - frag_size.data, + frag_volume.data, x_plus_y.data, self._get_floating_point(vmin), self._get_floating_point(nfmax if nfmax else -1), @@ -980,23 +980,23 @@ def exp_fragmentation( ) def gauss_fragmentation( - self, *, n_fragment, mu, sigma, frag_size, x_plus_y, rand, vmin, nfmax + self, *, n_fragment, mu, sigma, frag_volume, x_plus_y, rand, vmin, nfmax ): self.__gauss_fragmentation_body.launch_n( - n=len(frag_size), + n=len(frag_volume), args=( self._get_floating_point(mu), self._get_floating_point(sigma), - frag_size.data, + frag_volume.data, rand.data, ), ) self.__fragmentation_limiters_body.launch_n( - n=len(frag_size), + n=len(frag_volume), args=( n_fragment.data, - frag_size.data, + frag_volume.data, x_plus_y.data, self._get_floating_point(vmin), self._get_floating_point(nfmax if nfmax else -1), @@ -1005,13 +1005,13 @@ def gauss_fragmentation( ) def slams_fragmentation( - self, n_fragment, frag_size, x_plus_y, probs, rand, vmin, nfmax + self, n_fragment, frag_volume, x_plus_y, probs, rand, vmin, nfmax ): # pylint: disable=too-many-arguments self.__slams_fragmentation_body.launch_n( n=(len(n_fragment)), args=( n_fragment.data, - frag_size.data, + frag_volume.data, x_plus_y.data, probs.data, rand.data, @@ -1019,10 +1019,10 @@ def slams_fragmentation( ) self.__fragmentation_limiters_body.launch_n( - n=len(frag_size), + n=len(frag_volume), args=( n_fragment.data, - frag_size.data, + frag_volume.data, x_plus_y.data, self._get_floating_point(vmin), self._get_floating_point(nfmax if nfmax else -1), @@ -1035,7 +1035,7 @@ def feingold1988_fragmentation( *, n_fragment, scale, - frag_size, + frag_volume, x_plus_y, rand, fragtol, @@ -1043,10 +1043,10 @@ def feingold1988_fragmentation( nfmax, ): self.__feingold1988_fragmentation_body.launch_n( - n=len(frag_size), + n=len(frag_volume), args=( self._get_floating_point(scale), - frag_size.data, + frag_volume.data, x_plus_y.data, rand.data, self._get_floating_point(fragtol), @@ -1054,10 +1054,10 @@ def feingold1988_fragmentation( ) self.__fragmentation_limiters_body.launch_n( - n=len(frag_size), + n=len(frag_volume), args=( n_fragment.data, - frag_size.data, + frag_volume.data, x_plus_y.data, self._get_floating_point(vmin), self._get_floating_point(nfmax if nfmax else -1), @@ -1073,7 +1073,7 @@ def straub_fragmentation( CW, gam, ds, - frag_size, + frag_volume, v_max, x_plus_y, rand, @@ -1087,12 +1087,12 @@ def straub_fragmentation( d34, ): self.__straub_fragmentation_body.launch_n( - n=len(frag_size), + n=len(frag_volume), args=( CW.data, gam.data, ds.data, - frag_size.data, + frag_volume.data, v_max.data, rand.data, Nr1.data, @@ -1105,10 +1105,10 @@ def straub_fragmentation( ) self.__fragmentation_limiters_body.launch_n( - n=len(frag_size), + n=len(frag_volume), args=( n_fragment.data, - frag_size.data, + frag_volume.data, x_plus_y.data, self._get_floating_point(vmin), self._get_floating_point(nfmax if nfmax else -1), diff --git a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py index bb586398f..96be42f74 100644 --- a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py +++ b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py @@ -38,13 +38,12 @@ class CondensationMethods( @cached_property def __calculate_m_l(self): - const = self.formulae.constants return trtc.For( - ("ml", "v", "multiplicity", "cell_id"), - "i", - f""" - atomicAdd((real_type*) &ml[cell_id[i]], multiplicity[i] * v[i] * {const.rho_w}); - """.replace( + param_names=("ml", "water_mass", "multiplicity", "cell_id"), + name_iter="i", + body=""" + atomicAdd((real_type*) &ml[cell_id[i]], multiplicity[i] * water_mass[i]); + """.replace( "real_type", self._get_c_type() ), ) @@ -66,12 +65,12 @@ def __init__(self): self.vars_data: Optional[Dict] = None @cached_property - def __update_volume(self): + def __update_drop_masses(self): phys = self.formulae const = phys.constants return trtc.For( - ( - "v", + param_names=( + "water_mass", "vdry", *CondensationMethods.keys, "_kappa", @@ -82,8 +81,8 @@ def __update_volume(self): "max_iters", "cell_id", ), - "i", - f""" + name_iter="i", + body=f""" struct Minfun {{ static __device__ real_type value(real_type x_new, void* args_p) {{ auto args = static_cast(args_p); @@ -127,12 +126,15 @@ def __update_volume(self): auto _lambdaK = lambdaK[cell_id[i]]; auto _lambdaD = lambdaD[cell_id[i]]; - auto x_old = {phys.condensation_coordinate.x.c_inline(volume="v[i]")}; - auto r_old = {phys.trivia.radius.c_inline(volume="v[i]")}; + auto v_old = {phys.particle_shape_and_density.mass_to_volume.c_inline( + mass="water_mass[i]" + )}; + auto x_old = {phys.condensation_coordinate.x.c_inline(volume="v_old")}; + auto r_old = {phys.trivia.radius.c_inline(volume="v_old")}; auto x_insane = {phys.condensation_coordinate.x.c_inline(volume="vdry[i]/100")}; auto rd3 = vdry[i] / {const.PI_4_3}; auto sgm = {phys.surface_tension.sigma.c_inline( - T="_T", v_wet="v[i]", v_dry="vdry[i]", f_org="_f_org[i]" + T="_T", v_wet="v", v_dry="vdry[i]", f_org="_f_org[i]" )}; auto RH_eq = {phys.hygroscopicity.RH_eq.c_inline( r="r_old", T="_T", kp="_kappa[i]", rd3="rd3", sgm="sgm" @@ -211,7 +213,10 @@ def __update_volume(self): x_new = x_old; }} }} - v[i] = {phys.condensation_coordinate.volume.c_inline(x="x_new")}; + auto v_new = {phys.condensation_coordinate.volume.c_inline(x="x_new")}; + water_mass[i] = {phys.particle_shape_and_density.volume_to_mass.c_inline( + volume="v_new" + )}; """.replace( "real_type", self._get_c_type() ), @@ -220,7 +225,7 @@ def __update_volume(self): @cached_property def __pre_for(self): return trtc.For( - ( + param_names=( "dthd_dt_pred", "d_water_vapour_mixing_ratio__dt_predicted", "drhod_dt_pred", @@ -235,8 +240,8 @@ def __pre_for(self): "RH_max", "dv", ), - "i", - """ + name_iter="i", + body=""" dthd_dt_pred[i] = (pthd[i] - thd[i]) / dt; d_water_vapour_mixing_ratio__dt_predicted[i] = ( predicted_water_vapour_mixing_ratio[i] - water_vapour_mixing_ratio[i] @@ -256,7 +261,7 @@ def __pre_for(self): def __pre(self): phys = self.formulae return trtc.For( - ( + param_names=( *CondensationMethods.keys, "dthd_dt_pred", "d_water_vapour_mixing_ratio__dt_predicted", @@ -267,8 +272,8 @@ def __pre(self): "dt", "RH_max", ), - "i", - f""" + name_iter="i", + body=f""" pthd[i] += dt * dthd_dt_pred[i] / 2; predicted_water_vapour_mixing_ratio[i] += ( dt * d_water_vapour_mixing_ratio__dt_predicted[i] / 2 @@ -304,7 +309,7 @@ def __pre(self): def __post(self): phys = self.formulae return trtc.For( - ( + param_names=( "dthd_dt_pred", "d_water_vapour_mixing_ratio__dt_predicted", "drhod_dt_pred", @@ -318,8 +323,8 @@ def __post(self): "T", "lv", ), - "i", - f""" + name_iter="i", + body=f""" auto dml_dt = (ml_new[i] - ml_old[i]) / dt; auto d_water_vapour_mixing_ratio__dt_corrected = - dml_dt / m_d[i]; auto dthd_dt_corr = {phys.state_variable_triplet.dthd_dt.c_inline( @@ -339,10 +344,11 @@ def __post(self): ), ) - def calculate_m_l(self, ml, v, multiplicity, cell_id): + def calculate_m_l(self, ml, water_mass, multiplicity, cell_id): ml[:] = 0 self.__calculate_m_l.launch_n( - len(multiplicity), (ml.data, v.data, multiplicity.data, cell_id.data) + n=len(multiplicity), + args=(ml.data, water_mass.data, multiplicity.data, cell_id.data), ) # pylint: disable=unused-argument,too-many-locals @@ -353,7 +359,7 @@ def condensation( solver, n_cell, cell_start_arg, - v, + water_mass, v_cr, multiplicity, vdry, @@ -388,8 +394,8 @@ def condensation( self.rhod_copy.fill(rhod) self.__pre_for.launch_n( - n_cell, - ( + n=n_cell, + args=( self.dthd_dt_pred.data, self.d_water_vapour_mixing_ratio__dt_predicted.data, self.drhod_dt_pred.data, @@ -406,12 +412,12 @@ def condensation( ), ) timestep /= n_substeps - self.calculate_m_l(self.ml_old, v, multiplicity, cell_id) + self.calculate_m_l(self.ml_old, water_mass, multiplicity, cell_id) for _ in range(n_substeps): self.__pre.launch_n( - n_cell, - ( + n=n_cell, + args=( *self.vars_data.values(), self.dthd_dt_pred.data, self.d_water_vapour_mixing_ratio__dt_predicted.data, @@ -423,10 +429,10 @@ def condensation( RH_max.data, ), ) - self.__update_volume.launch_n( - len(multiplicity), - ( - v.data, + self.__update_drop_masses.launch_n( + n=len(multiplicity), + args=( + water_mass.data, vdry.data, *self.vars_data.values(), kappa.data, @@ -438,10 +444,10 @@ def condensation( cell_id.data, ), ) - self.calculate_m_l(self.ml_new, v, multiplicity, cell_id) + self.calculate_m_l(self.ml_new, water_mass, multiplicity, cell_id) self.__post.launch_n( - n_cell, - ( + n=n_cell, + args=( self.dthd_dt_pred.data, self.d_water_vapour_mixing_ratio__dt_predicted.data, self.drhod_dt_pred.data, diff --git a/PySDM/backends/impl_thrust_rtc/methods/freezing_methods.py b/PySDM/backends/impl_thrust_rtc/methods/freezing_methods.py index 148851b74..3e4c20bb7 100644 --- a/PySDM/backends/impl_thrust_rtc/methods/freezing_methods.py +++ b/PySDM/backends/impl_thrust_rtc/methods/freezing_methods.py @@ -13,12 +13,11 @@ class FreezingMethods(ThrustRTCBackendMethods): @cached_property def freeze_time_dependent_body(self): - const = self.formulae.constants return trtc.For( param_names=( "rand", "immersed_surface_area", - "wet_volume", + "water_mass", "timestep", "cell", "a_w_ice", @@ -32,12 +31,12 @@ def freeze_time_dependent_body(self): return; }} if (thaw && {self.formulae.trivia.frozen_and_above_freezing_point.c_inline( - volume="wet_volume[i]", + water_mass="water_mass[i]", temperature="temperature[cell[i]]" )}) {{ - wet_volume[i] = -1 * wet_volume[i] * {const.rho_i} / {const.rho_w}; + water_mass[i] = -1 * water_mass[i]; }} else if ({self.formulae.trivia.unfrozen_and_saturated.c_inline( - volume="wet_volume[i]", + water_mass="water_mass[i]", relative_humidity="relative_humidity[cell[i]]" )}) {{ auto rate = {self.formulae.heterogeneous_ice_nucleation_rate.j_het.c_inline( @@ -47,7 +46,7 @@ def freeze_time_dependent_body(self): -rate * immersed_surface_area[i] * timestep ); if (rand[i] < prob) {{ - wet_volume[i] = -1 * wet_volume[i] * {const.rho_w} / {const.rho_i}; + water_mass[i] = -1 * water_mass[i]; }} }} """.replace( @@ -57,11 +56,10 @@ def freeze_time_dependent_body(self): @cached_property def freeze_singular_body(self): - const = self.formulae.constants return trtc.For( param_names=( "freezing_temperature", - "wet_volume", + "water_mass", "temperature", "relative_humidity", "cell", @@ -73,17 +71,17 @@ def freeze_singular_body(self): return; }} if (thaw && {self.formulae.trivia.frozen_and_above_freezing_point.c_inline( - volume="wet_volume[i]", + water_mass="water_mass[i]", temperature="temperature[cell[i]]" )}) {{ - wet_volume[i] = -1 * wet_volume[i] * {const.rho_i} / {const.rho_w}; + water_mass[i] = -1 * water_mass[i]; }} else if ( {self.formulae.trivia.unfrozen_and_saturated.c_inline( - volume="wet_volume[i]", + water_mass="water_mass[i]", relative_humidity="relative_humidity[cell[i]]" )} && temperature[cell[i]] <= freezing_temperature[i] ) {{ - wet_volume[i] = -1 * wet_volume[i] * {const.rho_w} / {const.rho_i}; + water_mass[i] = -1 * water_mass[i]; }} """.replace( "real_type", self._get_c_type() @@ -99,7 +97,7 @@ def freeze_singular( n=n_sd, args=( attributes.freezing_temperature.data, - attributes.wet_volume.data, + attributes.water_mass.data, temperature.data, relative_humidity.data, cell.data, @@ -128,7 +126,7 @@ def freeze_time_dependent( # pylint: disable=unused-argument args=( rand.data, attributes.immersed_surface_area.data, - attributes.wet_volume.data, + attributes.water_mass.data, self._get_floating_point(timestep), cell.data, a_w_ice.data, diff --git a/PySDM/backends/impl_thrust_rtc/methods/physics_methods.py b/PySDM/backends/impl_thrust_rtc/methods/physics_methods.py index 6e13d949f..3308c96dc 100644 --- a/PySDM/backends/impl_thrust_rtc/methods/physics_methods.py +++ b/PySDM/backends/impl_thrust_rtc/methods/physics_methods.py @@ -89,6 +89,30 @@ def __terminal_velocity_body(self): ), ) + @cached_property + def __volume_of_mass_body(self): + return trtc.For( + param_names=("volume", "mass"), + name_iter="i", + body=f""" + volume[i] = {self.formulae.particle_shape_and_density.mass_to_volume.c_inline(mass="mass[i]")}; + """.replace( + "real_type", self._get_c_type() + ), + ) + + @cached_property + def __mass_of_volume_body(self): + return trtc.For( + param_names=("mass", "volume"), + name_iter="i", + body=f""" + mass[i] = {self.formulae.particle_shape_and_density.volume_to_mass.c_inline(volume="volume[i]")}; + """.replace( + "real_type", self._get_c_type() + ), + ) + @nice_thrust(**NICE_THRUST_FLAGS) def critical_volume(self, *, v_cr, kappa, f_org, v_dry, v_wet, T, cell): self.__critical_volume_body.launch_n( @@ -136,3 +160,9 @@ def explicit_euler(self, y, dt, dy_dt): dt = self._get_floating_point(dt) dy_dt = self._get_floating_point(dy_dt) self.__explicit_euler_body.launch_n(y.shape[0], (y.data, dt, dy_dt)) + + def volume_of_water_mass(self, volume, mass): + self.__volume_of_mass_body.launch_n(volume.shape[0], (volume.data, mass.data)) + + def mass_of_water_volume(self, mass, volume): + self.__mass_of_volume_body.launch_n(mass.shape[0], (mass.data, volume.data)) diff --git a/PySDM/backends/impl_thrust_rtc/test_helpers/cpp2python.py b/PySDM/backends/impl_thrust_rtc/test_helpers/cpp2python.py index 8d2dd1fb9..12ce47b09 100644 --- a/PySDM/backends/impl_thrust_rtc/test_helpers/cpp2python.py +++ b/PySDM/backends/impl_thrust_rtc/test_helpers/cpp2python.py @@ -190,7 +190,15 @@ def to_numba(name, args, iter_var, body): f""" def make(self): import numpy as np - from numpy import floor, ceil, exp, log, power, sqrt, arctanh as atanh, sinh, arcsinh as asinh + from numpy import ( + floor, ceil, + exp, log, + power, sqrt, + arctanh as atanh, + arcsinh as asinh, + sinh, + maximum, minimum + ) import numba @numba.njit(parallel=False, {JIT_OPTS}) diff --git a/PySDM/builder.py b/PySDM/builder.py index 4f9f54e03..310a2894e 100644 --- a/PySDM/builder.py +++ b/PySDM/builder.py @@ -8,14 +8,15 @@ from PySDM.attributes.impl.mapper import get_class as attr_class from PySDM.attributes.numerics.cell_id import CellID +from PySDM.attributes.physics import WaterMass from PySDM.attributes.physics.multiplicities import Multiplicities -from PySDM.attributes.physics.volume import Volume from PySDM.impl.particle_attributes_factory import ParticleAttributesFactory from PySDM.impl.wall_timer import WallTimer from PySDM.initialisation.discretise_multiplicities import ( # TODO #324 discretise_multiplicities, ) from PySDM.particulator import Particulator +from PySDM.physics.particle_shape_and_density import LiquidSpheres, MixedPhaseSpheres class Builder: @@ -25,7 +26,7 @@ def __init__(self, n_sd, backend): self.particulator = Particulator(n_sd, backend) self.req_attr = { "multiplicity": Multiplicities(self), - "volume": Volume(self), + "water mass": WaterMass(self), "cell id": CellID(self), } self.aerosol_radius_threshold = 0 @@ -87,6 +88,19 @@ def build( DeprecationWarning, ) + if "volume" in attributes and "water mass" not in attributes: + assert self.particulator.formulae.particle_shape_and_density.__name__ in ( + LiquidSpheres.__name__, + MixedPhaseSpheres.__name__, + ), "implied volume-to-mass conversion is only supported for spherical particles" + attributes[ + "water mass" + ] = self.particulator.formulae.particle_shape_and_density.volume_to_mass( + attributes["volume"] + ) + del attributes["volume"] + self.request_attribute("volume") + for dynamic in self.particulator.dynamics.values(): dynamic.register(self) diff --git a/PySDM/dynamics/collisions/breakup_fragmentations/__init__.py b/PySDM/dynamics/collisions/breakup_fragmentations/__init__.py index 413025948..85cf4239b 100644 --- a/PySDM/dynamics/collisions/breakup_fragmentations/__init__.py +++ b/PySDM/dynamics/collisions/breakup_fragmentations/__init__.py @@ -2,9 +2,10 @@ TODO #744 """ from .always_n import AlwaysN -from .constant_size import ConstantSize -from .exponential import ExponFrag -from .feingold1988 import Feingold1988Frag +from .constant_mass import ConstantMass +from .expon_frag import ExponFrag +from .exponential import Exponential +from .feingold1988 import Feingold1988 from .gaussian import Gaussian from .lowlist82 import LowList1982Nf from .slams import SLAMS diff --git a/PySDM/dynamics/collisions/breakup_fragmentations/always_n.py b/PySDM/dynamics/collisions/breakup_fragmentations/always_n.py index ab02a47b6..7c76a4148 100644 --- a/PySDM/dynamics/collisions/breakup_fragmentations/always_n.py +++ b/PySDM/dynamics/collisions/breakup_fragmentations/always_n.py @@ -4,17 +4,14 @@ class AlwaysN: # pylint: disable=too-many-instance-attributes - def __init__(self, n, vmin=0.0, nfmax=None): + def __init__(self, n): self.particulator = None self.N = n - self.vmin = vmin - self.nfmax = nfmax - def __call__(self, nf, frag_size, u01, is_first_in_pair): + def __call__(self, nf, frag_mass, u01, is_first_in_pair): nf.fill(self.N) - frag_size.sum(self.particulator.attributes["volume"], is_first_in_pair) - frag_size /= self.N + frag_mass.sum(self.particulator.attributes["water mass"], is_first_in_pair) + frag_mass /= self.N def register(self, builder): self.particulator = builder.particulator - builder.request_attribute("volume") diff --git a/PySDM/dynamics/collisions/breakup_fragmentations/constant_mass.py b/PySDM/dynamics/collisions/breakup_fragmentations/constant_mass.py new file mode 100644 index 000000000..d7951d6c6 --- /dev/null +++ b/PySDM/dynamics/collisions/breakup_fragmentations/constant_mass.py @@ -0,0 +1,17 @@ +""" +Always produces fragments of mass c in a given collisional breakup +""" + + +class ConstantMass: # pylint: disable=too-many-instance-attributes + def __init__(self, c): + self.particulator = None + self.C = c + + def __call__(self, nf, frag_mass, u01, is_first_in_pair): + frag_mass[:] = self.C + nf.sum(self.particulator.attributes["water mass"], is_first_in_pair) + nf /= self.C + + def register(self, builder): + self.particulator = builder.particulator diff --git a/PySDM/dynamics/collisions/breakup_fragmentations/constant_size.py b/PySDM/dynamics/collisions/breakup_fragmentations/constant_size.py deleted file mode 100644 index 4e63bb4d9..000000000 --- a/PySDM/dynamics/collisions/breakup_fragmentations/constant_size.py +++ /dev/null @@ -1,18 +0,0 @@ -""" -Always produces fragments of size c in a given collisional breakup -""" - - -class ConstantSize: # pylint: disable=too-many-instance-attributes - def __init__(self, c): - self.particulator = None - self.C = c - - def __call__(self, nf, frag_size, u01, is_first_in_pair): - frag_size[:] = self.C - nf.sum(self.particulator.attributes["volume"], is_first_in_pair) - nf /= self.C - - def register(self, builder): - self.particulator = builder.particulator - builder.request_attribute("volume") diff --git a/PySDM/dynamics/collisions/breakup_fragmentations/expon_frag.py b/PySDM/dynamics/collisions/breakup_fragmentations/expon_frag.py new file mode 100644 index 000000000..210fa376c --- /dev/null +++ b/PySDM/dynamics/collisions/breakup_fragmentations/expon_frag.py @@ -0,0 +1,12 @@ +""" +DEPRECATED +P(x) = exp(-x / lambda); lambda specified in volume units +""" +import warnings + +from .exponential import Exponential + + +class ExponFrag(Exponential): # pylint: disable=too-few-public-methods + def __init_subclass__(cls): + warnings.warn("Class has been renamed", DeprecationWarning) diff --git a/PySDM/dynamics/collisions/breakup_fragmentations/exponential.py b/PySDM/dynamics/collisions/breakup_fragmentations/exponential.py index 3b64bfccc..d41e04a48 100644 --- a/PySDM/dynamics/collisions/breakup_fragmentations/exponential.py +++ b/PySDM/dynamics/collisions/breakup_fragmentations/exponential.py @@ -1,31 +1,34 @@ """ -P(x) = exp(-x / lambda); specified in volume units +P(x) = exp(-x / lambda); lambda specified in volume units """ # TODO #796: introduce common code with Feingold fragmentation, including possible limiter +from .impl import VolumeBasedFragmentationFunction -class ExponFrag: +class Exponential(VolumeBasedFragmentationFunction): def __init__(self, scale, vmin=0.0, nfmax=None): - self.particulator = None + super().__init__() self.scale = scale self.vmin = vmin self.nfmax = nfmax self.sum_of_volumes = None def register(self, builder): - self.particulator = builder.particulator + super().register(builder) self.sum_of_volumes = self.particulator.PairwiseStorage.empty( self.particulator.n_sd // 2, dtype=float ) - def __call__(self, nf, frag_size, u01, is_first_in_pair): + def compute_fragment_number_and_volumes( + self, nf, frag_volume, u01, is_first_in_pair + ): self.sum_of_volumes.sum( self.particulator.attributes["volume"], is_first_in_pair ) self.particulator.backend.exp_fragmentation( n_fragment=nf, scale=self.scale, - frag_size=frag_size, + frag_volume=frag_volume, x_plus_y=self.sum_of_volumes, rand=u01, vmin=self.vmin, diff --git a/PySDM/dynamics/collisions/breakup_fragmentations/feingold1988.py b/PySDM/dynamics/collisions/breakup_fragmentations/feingold1988.py index ce11351df..9c56669e8 100644 --- a/PySDM/dynamics/collisions/breakup_fragmentations/feingold1988.py +++ b/PySDM/dynamics/collisions/breakup_fragmentations/feingold1988.py @@ -1,13 +1,16 @@ """ P(m; x, y) = nu^2 * (x+y) exp(-m * nu) -nu = 1/m* where m* is a scaling factor for fragment size dist. +nu = 1/m* where m* is a scaling factor for fragment volume dist. see [Feingold et al. 1999](https://doi.org/10.1175/1520-0469(1999)056<4100:TIOGCC>2.0.CO;2) """ +from .impl import VolumeBasedFragmentationFunction -class Feingold1988Frag: # pylint: disable=too-many-instance-attributes +class Feingold1988( + VolumeBasedFragmentationFunction +): # pylint: disable=too-many-instance-attributes def __init__(self, scale, fragtol=1e-3, vmin=0.0, nfmax=None): - self.particulator = None + super().__init__() self.scale = scale self.fragtol = fragtol self.vmin = vmin @@ -15,20 +18,21 @@ def __init__(self, scale, fragtol=1e-3, vmin=0.0, nfmax=None): self.sum_of_volumes = None def register(self, builder): - self.particulator = builder.particulator - builder.request_attribute("volume") + super().register(builder) self.sum_of_volumes = self.particulator.PairwiseStorage.empty( self.particulator.n_sd // 2, dtype=float ) - def __call__(self, nf, frag_size, u01, is_first_in_pair): + def compute_fragment_number_and_volumes( + self, nf, frag_volume, u01, is_first_in_pair + ): self.sum_of_volumes.sum( self.particulator.attributes["volume"], is_first_in_pair ) self.particulator.backend.feingold1988_fragmentation( n_fragment=nf, scale=self.scale, - frag_size=frag_size, + frag_volume=frag_volume, x_plus_y=self.sum_of_volumes, rand=u01, fragtol=self.fragtol, diff --git a/PySDM/dynamics/collisions/breakup_fragmentations/gaussian.py b/PySDM/dynamics/collisions/breakup_fragmentations/gaussian.py index e65b223e7..bd6c02db8 100644 --- a/PySDM/dynamics/collisions/breakup_fragmentations/gaussian.py +++ b/PySDM/dynamics/collisions/breakup_fragmentations/gaussian.py @@ -1,11 +1,14 @@ """ P(x) = exp(-(x-mu)^2 / 2 sigma^2); mu and sigma are volumes """ +from .impl import VolumeBasedFragmentationFunction -class Gaussian: # pylint: disable=too-many-instance-attributes +class Gaussian( + VolumeBasedFragmentationFunction +): # pylint: disable=too-many-instance-attributes def __init__(self, mu, sigma, vmin=0.0, nfmax=None): - self.particulator = None + super().__init__() self.mu = mu self.sigma = sigma self.vmin = vmin @@ -13,12 +16,14 @@ def __init__(self, mu, sigma, vmin=0.0, nfmax=None): self.sum_of_volumes = None def register(self, builder): - self.particulator = builder.particulator + super().register(builder) self.sum_of_volumes = self.particulator.PairwiseStorage.empty( self.particulator.n_sd // 2, dtype=float ) - def __call__(self, nf, frag_size, u01, is_first_in_pair): + def compute_fragment_number_and_volumes( + self, nf, frag_volume, u01, is_first_in_pair + ): self.sum_of_volumes.sum( self.particulator.attributes["volume"], is_first_in_pair ) @@ -26,7 +31,7 @@ def __call__(self, nf, frag_size, u01, is_first_in_pair): n_fragment=nf, mu=self.mu, sigma=self.sigma, - frag_size=frag_size, + frag_volume=frag_volume, x_plus_y=self.sum_of_volumes, rand=u01, vmin=self.vmin, diff --git a/PySDM/dynamics/collisions/breakup_fragmentations/impl/__init__.py b/PySDM/dynamics/collisions/breakup_fragmentations/impl/__init__.py new file mode 100644 index 000000000..3b6014a62 --- /dev/null +++ b/PySDM/dynamics/collisions/breakup_fragmentations/impl/__init__.py @@ -0,0 +1,4 @@ +""" +Abstractions and common code for fragmentation functions +""" +from .volume_based import VolumeBasedFragmentationFunction diff --git a/PySDM/dynamics/collisions/breakup_fragmentations/impl/volume_based.py b/PySDM/dynamics/collisions/breakup_fragmentations/impl/volume_based.py new file mode 100644 index 000000000..62bcb8c08 --- /dev/null +++ b/PySDM/dynamics/collisions/breakup_fragmentations/impl/volume_based.py @@ -0,0 +1,26 @@ +""" +Base class for volume-based fragmentation functions +""" + + +class VolumeBasedFragmentationFunction: + def __init__(self): + self.particulator = None + + def __call__(self, nf, frag_mass, u01, is_first_in_pair): + frag_volume_aliased_to_mass = frag_mass + self.compute_fragment_number_and_volumes( + nf, frag_volume_aliased_to_mass, u01, is_first_in_pair + ) + self.particulator.backend.mass_of_water_volume( + frag_mass, frag_volume_aliased_to_mass + ) + + def compute_fragment_number_and_volumes( + self, nf, frag_volume, u01, is_first_in_pair + ): + raise NotImplementedError() + + def register(self, builder): + self.particulator = builder.particulator + builder.request_attribute("volume") diff --git a/PySDM/dynamics/collisions/breakup_fragmentations/lowlist82.py b/PySDM/dynamics/collisions/breakup_fragmentations/lowlist82.py index de287dc3a..7a7a681a8 100644 --- a/PySDM/dynamics/collisions/breakup_fragmentations/lowlist82.py +++ b/PySDM/dynamics/collisions/breakup_fragmentations/lowlist82.py @@ -1,12 +1,13 @@ """ See [Low & List 1982](https://doi.org/10.1175/1520-0469(1982)039<1607:CCABOR>2.0.CO;2) """ +from .impl import VolumeBasedFragmentationFunction -class LowList1982Nf: +class LowList1982Nf(VolumeBasedFragmentationFunction): # pylint: disable=too-many-instance-attributes def __init__(self, vmin=0.0, nfmax=None): - self.particulator = None + super().__init__() self.vmin = vmin self.nfmax = nfmax self.arrays = {} @@ -15,13 +16,12 @@ def __init__(self, vmin=0.0, nfmax=None): self.const = None def register(self, builder): - self.particulator = builder.particulator + super().register(builder) self.sum_of_volumes = self.particulator.PairwiseStorage.empty( self.particulator.n_sd // 2, dtype=float ) self.const = self.particulator.formulae.constants builder.request_attribute("radius") - builder.request_attribute("volume") builder.request_attribute("relative fall velocity") for key in ("Sc", "St", "tmp", "tmp2", "CKE", "We", "W2", "ds", "dl", "dcoal"): self.arrays[key] = self.particulator.PairwiseStorage.empty( @@ -32,7 +32,9 @@ def register(self, builder): self.particulator.n_sd // 2, dtype=float ) - def __call__(self, nf, frag_size, u01, is_first_in_pair): + def compute_fragment_number_and_volumes( + self, nf, frag_volume, u01, is_first_in_pair + ): self.sum_of_volumes.sum( self.particulator.attributes["volume"], is_first_in_pair ) @@ -91,7 +93,7 @@ def __call__(self, nf, frag_size, u01, is_first_in_pair): ds=self.arrays["ds"], dl=self.arrays["dl"], dcoal=self.arrays["dcoal"], - frag_size=frag_size, + frag_volume=frag_volume, x_plus_y=self.sum_of_volumes, rand=u01, vmin=self.vmin, diff --git a/PySDM/dynamics/collisions/breakup_fragmentations/slams.py b/PySDM/dynamics/collisions/breakup_fragmentations/slams.py index ae95c9a83..fb5e78b29 100644 --- a/PySDM/dynamics/collisions/breakup_fragmentations/slams.py +++ b/PySDM/dynamics/collisions/breakup_fragmentations/slams.py @@ -2,23 +2,26 @@ Based on [Jokulsdottir & Archer 2016 (GMD)](https://doi.org/10.5194/gmd-9-1455-2016) for ocean particles """ +from .impl import VolumeBasedFragmentationFunction -class SLAMS: +class SLAMS(VolumeBasedFragmentationFunction): def __init__(self, vmin=0.0, nfmax=None): - self.particulator = None + super().__init__() self.p_vec = None self.sum_of_volumes = None self.vmin = vmin self.nfmax = nfmax - def __call__(self, nf, frag_size, u01, is_first_in_pair): + def compute_fragment_number_and_volumes( + self, nf, frag_volume, u01, is_first_in_pair + ): self.sum_of_volumes.sum( self.particulator.attributes["volume"], is_first_in_pair ) self.particulator.backend.slams_fragmentation( n_fragment=nf, - frag_size=frag_size, + frag_volume=frag_volume, x_plus_y=self.sum_of_volumes, probs=self.p_vec, rand=u01, @@ -27,7 +30,7 @@ def __call__(self, nf, frag_size, u01, is_first_in_pair): ) def register(self, builder): - self.particulator = builder.particulator + super().register(builder) self.p_vec = self.particulator.PairwiseStorage.empty( self.particulator.n_sd // 2, dtype=float ) diff --git a/PySDM/dynamics/collisions/breakup_fragmentations/straub2010.py b/PySDM/dynamics/collisions/breakup_fragmentations/straub2010.py index 16f1a82ec..83eae6f8b 100644 --- a/PySDM/dynamics/collisions/breakup_fragmentations/straub2010.py +++ b/PySDM/dynamics/collisions/breakup_fragmentations/straub2010.py @@ -3,11 +3,13 @@ """ from PySDM.physics.constants import si +from .impl import VolumeBasedFragmentationFunction -class Straub2010Nf: + +class Straub2010Nf(VolumeBasedFragmentationFunction): # pylint: disable=too-many-instance-attributes def __init__(self, vmin=0.0, nfmax=None): - self.particulator = None + super().__init__() self.vmin = vmin self.nfmax = nfmax self.arrays = {} @@ -17,7 +19,7 @@ def __init__(self, vmin=0.0, nfmax=None): self.const = None def register(self, builder): - self.particulator = builder.particulator + super().register(builder) self.max_size = self.particulator.PairwiseStorage.empty( self.particulator.n_sd // 2, dtype=float ) @@ -26,7 +28,6 @@ def register(self, builder): ) self.const = self.particulator.formulae.constants builder.request_attribute("radius") - builder.request_attribute("volume") builder.request_attribute("relative fall velocity") for key in ("Sc", "tmp", "tmp2", "CKE", "We", "gam", "CW", "ds"): self.arrays[key] = self.particulator.PairwiseStorage.empty( @@ -37,7 +38,9 @@ def register(self, builder): self.particulator.n_sd // 2, dtype=float ) - def __call__(self, nf, frag_size, u01, is_first_in_pair): + def compute_fragment_number_and_volumes( + self, nf, frag_volume, u01, is_first_in_pair + ): self.max_size.max(self.particulator.attributes["volume"], is_first_in_pair) self.sum_of_volumes.sum( self.particulator.attributes["volume"], is_first_in_pair @@ -82,7 +85,7 @@ def __call__(self, nf, frag_size, u01, is_first_in_pair): CW=self.arrays["CW"], gam=self.arrays["gam"], ds=self.arrays["ds"], - frag_size=frag_size, + frag_volume=frag_volume, v_max=self.max_size, x_plus_y=self.sum_of_volumes, rand=u01, diff --git a/PySDM/dynamics/collisions/coalescence_efficiencies/lowlist1982.py b/PySDM/dynamics/collisions/coalescence_efficiencies/lowlist1982.py index 30e403232..545949196 100644 --- a/PySDM/dynamics/collisions/coalescence_efficiencies/lowlist1982.py +++ b/PySDM/dynamics/collisions/coalescence_efficiencies/lowlist1982.py @@ -8,14 +8,12 @@ class LowList1982Ec: # pylint: disable=too-many-instance-attributes - def __init__(self, vmin=0.0, nfmax=None): + def __init__(self): self.particulator = None - self.vmin = vmin - self.nfmax = nfmax self.arrays = {} self.ll82_tmp = {} self.max_size = None - self.sum_of_volumes = None + self.sum_of_masses = None self.const = None def register(self, builder): @@ -23,12 +21,12 @@ def register(self, builder): self.max_size = self.particulator.PairwiseStorage.empty( self.particulator.n_sd // 2, dtype=float ) - self.sum_of_volumes = self.particulator.PairwiseStorage.empty( + self.sum_of_masses = self.particulator.PairwiseStorage.empty( self.particulator.n_sd // 2, dtype=float ) self.const = self.particulator.formulae.constants builder.request_attribute("radius") - builder.request_attribute("volume") + builder.request_attribute("water mass") builder.request_attribute("relative fall velocity") for key in ("Sc", "St", "dS", "tmp", "tmp2", "CKE", "Et", "ds", "dl"): self.arrays[key] = self.particulator.PairwiseStorage.empty( @@ -36,9 +34,9 @@ def register(self, builder): ) def __call__(self, output, is_first_in_pair): - self.max_size.max(self.particulator.attributes["volume"], is_first_in_pair) - self.sum_of_volumes.sum( - self.particulator.attributes["volume"], is_first_in_pair + self.max_size.max(self.particulator.attributes["water mass"], is_first_in_pair) + self.sum_of_masses.sum( + self.particulator.attributes["water mass"], is_first_in_pair ) self.arrays["ds"].min(self.particulator.attributes["radius"], is_first_in_pair) self.arrays["ds"] *= 2 @@ -46,7 +44,9 @@ def __call__(self, output, is_first_in_pair): self.arrays["dl"] *= 2 # compute the surface energy, CKE - self.arrays["Sc"].sum(self.particulator.attributes["volume"], is_first_in_pair) + self.arrays["Sc"].sum( + self.particulator.attributes["water mass"], is_first_in_pair + ) self.arrays["Sc"] **= 2 / 3 self.arrays["Sc"] *= ( self.const.PI * self.const.sgm_w * (6 / self.const.PI) ** (2 / 3) @@ -62,13 +62,15 @@ def __call__(self, output, is_first_in_pair): self.arrays["dS"].fill(self.arrays["St"]) self.arrays["dS"] -= self.arrays["Sc"] - self.arrays["tmp"].sum(self.particulator.attributes["volume"], is_first_in_pair) + self.arrays["tmp"].sum( + self.particulator.attributes["water mass"], is_first_in_pair + ) self.arrays["tmp2"].distance( self.particulator.attributes["relative fall velocity"], is_first_in_pair ) self.arrays["tmp2"] **= 2 self.arrays["CKE"].multiply( - self.particulator.attributes["volume"], is_first_in_pair + self.particulator.attributes["water mass"], is_first_in_pair ) self.arrays["CKE"].divide_if_not_zero(self.arrays["tmp"]) self.arrays["CKE"] *= self.arrays["tmp2"] diff --git a/PySDM/dynamics/collisions/collision.py b/PySDM/dynamics/collisions/collision.py index da3d2d863..a0773f28c 100644 --- a/PySDM/dynamics/collisions/collision.py +++ b/PySDM/dynamics/collisions/collision.py @@ -81,7 +81,7 @@ def __init__( self.kernel_temp = None self.n_fragment = None - self.fragment_size = None + self.fragment_mass = None self.Ec_temp = None self.Eb_temp = None self.norm_factor_temp = None @@ -149,7 +149,7 @@ def register(self, builder): self.n_fragment = self.particulator.PairwiseStorage.empty( **empty_args_pairwise ) - self.fragment_size = self.particulator.PairwiseStorage.empty( + self.fragment_mass = self.particulator.PairwiseStorage.empty( **empty_args_pairwise ) self.Ec_temp = self.particulator.PairwiseStorage.empty( @@ -206,7 +206,7 @@ def step(self): self.compute_coalescence_efficiency(self.Ec_temp, self.is_first_in_pair) self.compute_breakup_efficiency(self.Eb_temp, self.is_first_in_pair) self.compute_number_of_fragments( - self.n_fragment, self.fragment_size, rand_frag, self.is_first_in_pair + self.n_fragment, self.fragment_mass, rand_frag, self.is_first_in_pair ) else: proc_rand = None @@ -221,7 +221,7 @@ def step(self): rand=proc_rand, Ec=self.Ec_temp, Eb=self.Eb_temp, - fragment_size=self.fragment_size, + fragment_mass=self.fragment_mass, coalescence_rate=self.coalescence_rate, breakup_rate=self.breakup_rate, breakup_rate_deficit=self.breakup_rate_deficit, diff --git a/PySDM/dynamics/freezing.py b/PySDM/dynamics/freezing.py index 4f8ec0761..736ac2780 100644 --- a/PySDM/dynamics/freezing.py +++ b/PySDM/dynamics/freezing.py @@ -22,6 +22,8 @@ def __init__(self, *, singular=True, record_freezing_temperature=False, thaw=Fal def register(self, builder): self.particulator = builder.particulator + assert builder.formulae.particle_shape_and_density.supports_mixed_phase() + builder.request_attribute("volume") if self.singular or self.record_freezing_temperature: builder.request_attribute("freezing temperature") @@ -54,7 +56,7 @@ def __call__(self): freezing_temperature=self.particulator.attributes[ "freezing temperature" ], - wet_volume=self.particulator.attributes["volume"], + water_mass=self.particulator.attributes["water mass"], ), temperature=self.particulator.environment["T"], relative_humidity=self.particulator.environment["RH"], @@ -69,7 +71,7 @@ def __call__(self): immersed_surface_area=self.particulator.attributes[ "immersed surface area" ], - wet_volume=self.particulator.attributes["volume"], + water_mass=self.particulator.attributes["water mass"], ), timestep=self.particulator.dt, cell=self.particulator.attributes["cell id"], @@ -85,6 +87,6 @@ def __call__(self): thaw=self.thaw, ) - self.particulator.attributes.mark_updated("volume") + self.particulator.attributes.mark_updated("water mass") if self.record_freezing_temperature: self.particulator.attributes.mark_updated("freezing temperature") diff --git a/PySDM/dynamics/relaxed_velocity.py b/PySDM/dynamics/relaxed_velocity.py index bffd92765..b68c2cac9 100644 --- a/PySDM/dynamics/relaxed_velocity.py +++ b/PySDM/dynamics/relaxed_velocity.py @@ -28,9 +28,8 @@ def __init__(self, c: float = 8, constant: bool = False): self.particulator = None self.fall_momentum_attr = None self.terminal_vel_attr = None - self.volume_attr = None + self.water_mass_attr = None self.sqrt_radius_attr = None - self.rho_w = None # TODO #798 - we plan to use masses instead of volumes soon self.tmp_momentum_diff = None self.tmp_tau = None @@ -63,13 +62,11 @@ def register(self, builder): "relative fall momentum" ) self.terminal_vel_attr: Attribute = builder.get_attribute("terminal velocity") - self.volume_attr: Attribute = builder.get_attribute("volume") + self.water_mass_attr: Attribute = builder.get_attribute("water mass") self.sqrt_radius_attr: Attribute = builder.get_attribute( "square root of radius" ) - self.rho_w: float = builder.formulae.constants.rho_w # TODO #798 - self.tmp_momentum_diff = self.create_storage(self.particulator.n_sd) self.tmp_tau = self.create_storage(self.particulator.n_sd) self.tmp_scale = self.create_storage(self.particulator.n_sd) @@ -77,11 +74,8 @@ def register(self, builder): def __call__(self): # calculate momentum difference self.tmp_momentum_diff.product( - self.terminal_vel_attr.get(), self.volume_attr.get() + self.terminal_vel_attr.get(), self.water_mass_attr.get() ) - self.tmp_momentum_diff *= ( - self.rho_w - ) # TODO #798 - we plan to use masses instead of volumes soon self.tmp_momentum_diff -= self.fall_momentum_attr.get() if not self.tmp_tau_init or not self.constant: diff --git a/PySDM/formulae.py b/PySDM/formulae.py index eb95e3f69..6e8625712 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -42,6 +42,7 @@ def __init__( # pylint: disable=too-many-locals freezing_temperature_spectrum: str = "Null", heterogeneous_ice_nucleation_rate: str = "Null", fragmentation_function: str = "AlwaysN", + particle_shape_and_density: str = "LiquidSpheres", handle_all_breakups: bool = False, ): # initialisation of the fields below is just to silence pylint and to enable code hints @@ -61,6 +62,7 @@ def __init__( # pylint: disable=too-many-locals self.freezing_temperature_spectrum = freezing_temperature_spectrum self.heterogeneous_ice_nucleation_rate = heterogeneous_ice_nucleation_rate self.fragmentation_function = fragmentation_function + self.particle_shape_and_density = particle_shape_and_density components = tuple(i for i in dir(self) if not i.startswith("__")) constants_defaults = { @@ -208,6 +210,8 @@ def _c_inline(fun, return_type=None, constants=None, **args): source = source.replace("np.power(", "np.pow(") source = source.replace("np.arctanh(", "atanh(") source = source.replace("np.arcsinh(", "asinh(") + source = source.replace("np.minimum(", "min(") + source = source.replace("np.maximum(", "max(") for pkg in ("np", "math"): source = source.replace(f"{pkg}.", "") source = source.replace(", )", ")") diff --git a/PySDM/initialisation/init_fall_momenta.py b/PySDM/initialisation/init_fall_momenta.py index 861679bae..5d44dc9f9 100644 --- a/PySDM/initialisation/init_fall_momenta.py +++ b/PySDM/initialisation/init_fall_momenta.py @@ -15,7 +15,7 @@ def init_fall_momenta( volume: np.ndarray, rho_w: float, # TODO #798 - we plan to use masses instead of volumes soon zero: bool = False, - terminal_velocity_approx=Interpolation, + terminal_velocity_approx=Interpolation, # TODO #1155 ): """ Calculate default values of the @@ -33,7 +33,7 @@ def init_fall_momenta( if zero: return np.zeros_like(volume) - particulator = Particulator(0, CPU(Formulae())) + particulator = Particulator(0, CPU(Formulae())) # TODO #1155 approximation = terminal_velocity_approx(particulator=particulator) @@ -44,4 +44,4 @@ def init_fall_momenta( approximation(output=output, radius=radii) - return output.to_ndarray() * volume * rho_w # TODO #798 + return output.to_ndarray() * volume * rho_w # TODO #798 this assumes no ice diff --git a/PySDM/particulator.py b/PySDM/particulator.py index beea22090..17d1bb1d5 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -118,7 +118,7 @@ def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_orde solver=self.condensation_solver, n_cell=self.mesh.n_cell, cell_start_arg=self.attributes.cell_start, - v=self.attributes["volume"], + water_mass=self.attributes["water mass"], multiplicity=self.attributes["multiplicity"], vdry=self.attributes["dry volume"], idx=self.attributes._ParticleAttributes__idx, @@ -143,7 +143,7 @@ def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_orde success=success, cell_id=self.attributes["cell id"], ) - self.attributes.mark_updated("volume") + self.attributes.mark_updated("water mass") def collision_coalescence_breakup( self, @@ -153,7 +153,7 @@ def collision_coalescence_breakup( rand, Ec, Eb, - fragment_size, + fragment_mass, coalescence_rate, breakup_rate, breakup_rate_deficit, @@ -176,7 +176,7 @@ def collision_coalescence_breakup( rand=rand, Ec=Ec, Eb=Eb, - fragment_size=fragment_size, + fragment_mass=fragment_mass, healthy=healthy, cell_id=cell_id, coalescence_rate=coalescence_rate, @@ -184,7 +184,7 @@ def collision_coalescence_breakup( breakup_rate_deficit=breakup_rate_deficit, is_first_in_pair=is_first_in_pair, warn_overflows=warn_overflows, - volume=self.attributes["volume"], + particle_mass=self.attributes["water mass"], max_multiplicity=max_multiplicity, ) else: @@ -312,9 +312,9 @@ def moments( moment_0, moments, specs: dict, - attr_name="volume", + attr_name="water mass", attr_range=(-np.inf, np.inf), - weighting_attribute="volume", + weighting_attribute="water mass", weighting_rank=0, skip_division_by_m0=False, ): @@ -369,8 +369,8 @@ def spectrum_moments( attr, rank, attr_bins, - attr_name="volume", - weighting_attribute="volume", + attr_name="water mass", + weighting_attribute="water mass", weighting_rank=0, ): attr_data = self.attributes[attr] diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py index 30d0f2d4e..a5a3e3931 100644 --- a/PySDM/physics/__init__.py +++ b/PySDM/physics/__init__.py @@ -17,6 +17,7 @@ impl, latent_heat, particle_advection, + particle_shape_and_density, saturation_vapour_pressure, state_variable_triplet, surface_tension, diff --git a/PySDM/physics/constants.py b/PySDM/physics/constants.py index a21c2b1d6..ab9c42a6f 100644 --- a/PySDM/physics/constants.py +++ b/PySDM/physics/constants.py @@ -27,6 +27,8 @@ def convert_to(value, unit): PI = sci.pi PI_4_3 = PI * 4 / 3 LN_2 = np.log(2) +ZERO_MASS = 0 * si.kg +ZERO_VOLUME = 0 * si.m**3 TWO = 2 THREE = 3 FOUR = 4 diff --git a/PySDM/physics/fragmentation_function/__init__.py b/PySDM/physics/fragmentation_function/__init__.py index 2d558d448..88a0bb6e6 100644 --- a/PySDM/physics/fragmentation_function/__init__.py +++ b/PySDM/physics/fragmentation_function/__init__.py @@ -2,9 +2,10 @@ fragmentation functions for use with breakup """ from .always_n import AlwaysN -from .constant_size import ConstantSize +from .constant_mass import ConstantMass from .expon_frag import ExponFrag -from .feingold1988frag import Feingold1988Frag +from .exponential import Exponential +from .feingold1988 import Feingold1988 from .gaussian import Gaussian from .lowlist82 import LowList1982Nf from .slams import SLAMS diff --git a/PySDM/physics/fragmentation_function/constant_size.py b/PySDM/physics/fragmentation_function/constant_mass.py similarity index 57% rename from PySDM/physics/fragmentation_function/constant_size.py rename to PySDM/physics/fragmentation_function/constant_mass.py index 2c8a88d7c..b02507d5d 100644 --- a/PySDM/physics/fragmentation_function/constant_size.py +++ b/PySDM/physics/fragmentation_function/constant_mass.py @@ -1,8 +1,8 @@ """ -Formulae supporting `PySDM.dynamics.collisions.breakup_fragmentations.constant_size` +Formulae supporting `PySDM.dynamics.collisions.breakup_fragmentations.constant_mass` """ -class ConstantSize: # pylint: disable=too-few-public-methods +class ConstantMass: # pylint: disable=too-few-public-methods def __init__(self, _): pass diff --git a/PySDM/physics/fragmentation_function/expon_frag.py b/PySDM/physics/fragmentation_function/expon_frag.py index ac4970de2..7094ccc4c 100644 --- a/PySDM/physics/fragmentation_function/expon_frag.py +++ b/PySDM/physics/fragmentation_function/expon_frag.py @@ -1,8 +1,11 @@ """ -Formulae supporting `PySDM.dynamics.collisions.breakup_fragmentations.exponential` +Formulae supporting `PySDM.dynamics.collisions.breakup_fragmentations.expon_frag` """ +import warnings +from .exponential import Exponential -class ExponFrag: # pylint: disable=too-few-public-methods - def __init__(self, _): - pass + +class ExponFrag(Exponential): # pylint: disable=too-few-public-methods + def __init_subclass__(cls): + warnings.warn("Class has been renamed", DeprecationWarning) diff --git a/PySDM/physics/fragmentation_function/exponential.py b/PySDM/physics/fragmentation_function/exponential.py new file mode 100644 index 000000000..f9323771a --- /dev/null +++ b/PySDM/physics/fragmentation_function/exponential.py @@ -0,0 +1,8 @@ +""" +Formulae supporting `PySDM.dynamics.collisions.breakup_fragmentations.exponential` +""" + + +class Exponential: # pylint: disable=too-few-public-methods + def __init__(self, _): + pass diff --git a/PySDM/physics/fragmentation_function/feingold1988frag.py b/PySDM/physics/fragmentation_function/feingold1988.py similarity index 60% rename from PySDM/physics/fragmentation_function/feingold1988frag.py rename to PySDM/physics/fragmentation_function/feingold1988.py index 220b5285a..f6e454751 100644 --- a/PySDM/physics/fragmentation_function/feingold1988frag.py +++ b/PySDM/physics/fragmentation_function/feingold1988.py @@ -4,10 +4,10 @@ import numpy as np -class Feingold1988Frag: # pylint: disable=too-few-public-methods +class Feingold1988: # pylint: disable=too-few-public-methods def __init__(self, _): pass @staticmethod - def frag_size(_, scale, rand, x_plus_y, fragtol): + def frag_volume(_, scale, rand, x_plus_y, fragtol): return -scale * np.log(max(1 - rand * scale / x_plus_y, fragtol)) diff --git a/PySDM/physics/particle_shape_and_density/__init__.py b/PySDM/physics/particle_shape_and_density/__init__.py new file mode 100644 index 000000000..8dd4350fd --- /dev/null +++ b/PySDM/physics/particle_shape_and_density/__init__.py @@ -0,0 +1,6 @@ +""" +particle shape and density relationships +""" +from .liquid_spheres import LiquidSpheres +from .mixed_phase_spheres import MixedPhaseSpheres +from .porous_spheroids import PorousSpheroid diff --git a/PySDM/physics/particle_shape_and_density/liquid_spheres.py b/PySDM/physics/particle_shape_and_density/liquid_spheres.py new file mode 100644 index 000000000..67e6747f0 --- /dev/null +++ b/PySDM/physics/particle_shape_and_density/liquid_spheres.py @@ -0,0 +1,25 @@ +""" +spherical particles with constant density of water +""" +import numpy as np + + +class LiquidSpheres: + def __init__(self, _): + pass + + @staticmethod + def supports_mixed_phase(_=None): + return False + + @staticmethod + def mass_to_volume(const, mass): + return mass / const.rho_w + + @staticmethod + def volume_to_mass(const, volume): + return const.rho_w * volume + + @staticmethod + def radius_to_mass(const, radius): + return const.rho_w * const.PI_4_3 * np.power(radius, const.THREE) diff --git a/PySDM/physics/particle_shape_and_density/mixed_phase_spheres.py b/PySDM/physics/particle_shape_and_density/mixed_phase_spheres.py new file mode 100644 index 000000000..b6de0f06b --- /dev/null +++ b/PySDM/physics/particle_shape_and_density/mixed_phase_spheres.py @@ -0,0 +1,31 @@ +""" +spherical particles with constant density of water or ice +""" +import numpy as np + + +class MixedPhaseSpheres: + def __init__(self, _): + pass + + @staticmethod + def supports_mixed_phase(_=None): + return True + + @staticmethod + def mass_to_volume(const, mass): + return ( + max(const.ZERO_MASS, mass) / const.rho_w + + min(const.ZERO_MASS, mass) / const.rho_i + ) + + @staticmethod + def volume_to_mass(const, volume): + return ( + np.maximum(const.ZERO_VOLUME, volume) * const.rho_w + + np.minimum(const.ZERO_VOLUME, volume) * const.rho_i + ) + + @staticmethod + def radius_to_mass(const, radius): + raise NotImplementedError() diff --git a/PySDM/physics/particle_shape_and_density/porous_spheroids.py b/PySDM/physics/particle_shape_and_density/porous_spheroids.py new file mode 100644 index 000000000..94d6cf3d6 --- /dev/null +++ b/PySDM/physics/particle_shape_and_density/porous_spheroids.py @@ -0,0 +1,10 @@ +""" +for mixed-phase microphysics as in +[Shima et al. 2020](https://doi.org/10.5194/gmd-13-4107-2020) +""" + + +class PorousSpheroid: # pylint: disable=too-few-public-methods + @staticmethod + def supports_mixed_phase(_=None): + return True diff --git a/PySDM/physics/trivia.py b/PySDM/physics/trivia.py index a1c96e4a9..4615487eb 100644 --- a/PySDM/physics/trivia.py +++ b/PySDM/physics/trivia.py @@ -75,12 +75,12 @@ def th_std(const, p, T): return T * np.power(const.p1000 / p, const.Rd_over_c_pd) @staticmethod - def unfrozen_and_saturated(_, volume, relative_humidity): - return volume > 0 and relative_humidity > 1 + def unfrozen_and_saturated(_, water_mass, relative_humidity): + return water_mass > 0 and relative_humidity > 1 @staticmethod - def frozen_and_above_freezing_point(const, volume, temperature): - return volume < 0 and temperature > const.T0 + def frozen_and_above_freezing_point(const, water_mass, temperature): + return water_mass < 0 and temperature > const.T0 @staticmethod def erfinv_approx(const, c): diff --git a/PySDM/products/aqueous_chemistry/acidity.py b/PySDM/products/aqueous_chemistry/acidity.py index 88141480d..b46e717fd 100644 --- a/PySDM/products/aqueous_chemistry/acidity.py +++ b/PySDM/products/aqueous_chemistry/acidity.py @@ -42,6 +42,7 @@ def _impl(self, **kwargs): self.formulae.trivia.volume(self.radius_range[0]), self.formulae.trivia.volume(self.radius_range[1]), ), + filter_attr="volume", weighting_attribute="volume", weighting_rank=self.weighting_rank, ) diff --git a/PySDM/products/freezing/ice_water_content.py b/PySDM/products/freezing/ice_water_content.py index 6d76dc2f9..10181bd42 100644 --- a/PySDM/products/freezing/ice_water_content.py +++ b/PySDM/products/freezing/ice_water_content.py @@ -13,16 +13,16 @@ def __init__(self, unit="kg/m^3", name=None, specific=False): def _impl(self, **kwargs): self._download_moment_to_buffer( - attr="volume", rank=1, filter_range=(-np.inf, 0) + attr="water mass", rank=1, filter_range=(-np.inf, 0) ) result = self.buffer.copy() self._download_moment_to_buffer( - attr="volume", rank=0, filter_range=(-np.inf, 0) + attr="water mass", rank=0, filter_range=(-np.inf, 0) ) conc = self.buffer - result[:] *= -self.formulae.constants.rho_i * conc / self.particulator.mesh.dv + result[:] *= -1 * conc / self.particulator.mesh.dv if self.specific: self._download_to_buffer(self.particulator.environment["rhod"]) diff --git a/PySDM/products/impl/moment_product.py b/PySDM/products/impl/moment_product.py index bd4dd21da..002f016bf 100644 --- a/PySDM/products/impl/moment_product.py +++ b/PySDM/products/impl/moment_product.py @@ -26,9 +26,9 @@ def _download_moment_to_buffer( *, attr, rank, - filter_attr="volume", + filter_attr="water mass", filter_range=(-np.inf, np.inf), - weighting_attribute="volume", + weighting_attribute="water mass", weighting_rank=0, skip_division_by_m0=False, ): diff --git a/PySDM/products/size_spectral/effective_radius.py b/PySDM/products/size_spectral/effective_radius.py index 760d17b94..dc3fc9cf1 100644 --- a/PySDM/products/size_spectral/effective_radius.py +++ b/PySDM/products/size_spectral/effective_radius.py @@ -39,12 +39,14 @@ def _impl(self, **kwargs): attr="volume", rank=2 / 3, filter_range=self.volume_range, + filter_attr="volume", ) tmp[:] = self.buffer[:] self._download_moment_to_buffer( attr="volume", rank=1, filter_range=self.volume_range, + filter_attr="volume", ) EffectiveRadius.__get_impl(self.buffer, tmp) return self.buffer diff --git a/PySDM/products/size_spectral/mean_radius.py b/PySDM/products/size_spectral/mean_radius.py index 5576a3f69..901bdc556 100644 --- a/PySDM/products/size_spectral/mean_radius.py +++ b/PySDM/products/size_spectral/mean_radius.py @@ -24,6 +24,7 @@ def _impl(self, **kwargs): self.formulae.trivia.volume(self.radius_range[0]), self.formulae.trivia.volume(self.radius_range[1]), ), + filter_attr="volume", ) self.buffer[:] /= self.formulae.constants.PI_4_3 ** (1 / 3) return self.buffer diff --git a/PySDM/products/size_spectral/particle_concentration.py b/PySDM/products/size_spectral/particle_concentration.py index b24369485..c9b859a62 100644 --- a/PySDM/products/size_spectral/particle_concentration.py +++ b/PySDM/products/size_spectral/particle_concentration.py @@ -22,11 +22,15 @@ def __init__( def _impl(self, **kwargs): self._download_moment_to_buffer( - attr="volume", + attr="water mass", rank=0, filter_range=( - self.formulae.trivia.volume(self.radius_range[0]), - self.formulae.trivia.volume(self.radius_range[1]), + self.formulae.particle_shape_and_density.volume_to_mass( + self.formulae.trivia.volume(radius=self.radius_range[0]) + ), + self.formulae.particle_shape_and_density.volume_to_mass( + self.formulae.trivia.volume(self.radius_range[1]) + ), ), ) return super()._impl(**kwargs) diff --git a/PySDM/products/size_spectral/total_particle_concentration.py b/PySDM/products/size_spectral/total_particle_concentration.py index ba66cd629..b5ac589f5 100644 --- a/PySDM/products/size_spectral/total_particle_concentration.py +++ b/PySDM/products/size_spectral/total_particle_concentration.py @@ -9,5 +9,5 @@ def __init__(self, name=None, unit="m^-3", stp=False): super().__init__(name=name, unit=unit, specific=False, stp=stp) def _impl(self, **kwargs): - self._download_moment_to_buffer(attr="volume", rank=0) + self._download_moment_to_buffer(attr="water mass", rank=0) return super()._impl(**kwargs) diff --git a/examples/PySDM_examples/Alpert_and_Knopf_2016/simulation.py b/examples/PySDM_examples/Alpert_and_Knopf_2016/simulation.py index 0867a2cc5..26574fea5 100644 --- a/examples/PySDM_examples/Alpert_and_Knopf_2016/simulation.py +++ b/examples/PySDM_examples/Alpert_and_Knopf_2016/simulation.py @@ -131,6 +131,7 @@ def plot_j_het(self, variant: str, abifm_params_case: str, ylim=None): assert variant in ("apparent", "actual") formulae = Formulae( + particle_shape_and_density="MixedPhaseSpheres", heterogeneous_ice_nucleation_rate="ABIFM", constants={ "ABIFM_M": self.cases[abifm_params_case]["ABIFM_m"], @@ -212,6 +213,7 @@ def simulation( seed=seed, heterogeneous_ice_nucleation_rate=heterogeneous_ice_nucleation_rate, constants=constants, + particle_shape_and_density="MixedPhaseSpheres", ) builder = Builder(n_sd=n_sd, backend=CPU(formulae=formulae)) env = Box(dt=time_step, dv=volume) diff --git a/examples/PySDM_examples/Arabas_et_al_2023/aida.ipynb b/examples/PySDM_examples/Arabas_et_al_2023/aida.ipynb index 4f99be360..14b2cd07a 100644 --- a/examples/PySDM_examples/Arabas_et_al_2023/aida.ipynb +++ b/examples/PySDM_examples/Arabas_et_al_2023/aida.ipynb @@ -65,6 +65,7 @@ " 'volume': 84.5 * si.m**3,\n", " 'kappa': .1, # hygroscopicity\n", " 'formulae': Formulae(\n", + " particle_shape_and_density=\"MixedPhaseSpheres\",\n", " saturation_vapour_pressure='MurphyKoop2005',\n", " heterogeneous_ice_nucleation_rate='ABIFM',\n", " freezing_temperature_spectrum='Niemand_et_al_2012',\n", diff --git a/examples/PySDM_examples/Arabas_et_al_2023/fig_11.ipynb b/examples/PySDM_examples/Arabas_et_al_2023/fig_11.ipynb index 150b8be9f..a83113a39 100644 --- a/examples/PySDM_examples/Arabas_et_al_2023/fig_11.ipynb +++ b/examples/PySDM_examples/Arabas_et_al_2023/fig_11.ipynb @@ -64,6 +64,7 @@ "outputs": [], "source": [ "formulae = Formulae(\n", + " particle_shape_and_density=\"MixedPhaseSpheres\",\n", " freezing_temperature_spectrum='Niemand_et_al_2012',\n", " heterogeneous_ice_nucleation_rate='ABIFM',\n", " constants=FREEZING_CONSTANTS[\"dust\"]\n", diff --git a/examples/PySDM_examples/Arabas_et_al_2023/fig_convergence.ipynb b/examples/PySDM_examples/Arabas_et_al_2023/fig_convergence.ipynb index 95b457cf2..9195d72c5 100644 --- a/examples/PySDM_examples/Arabas_et_al_2023/fig_convergence.ipynb +++ b/examples/PySDM_examples/Arabas_et_al_2023/fig_convergence.ipynb @@ -90,6 +90,7 @@ "products = (IceWaterContent(name=\"qi\"),)\n", "\n", "formulae = Formulae(\n", + " particle_shape_and_density=\"MixedPhaseSpheres\",\n", " heterogeneous_ice_nucleation_rate=\"Constant\",\n", " constants={\"J_HET\": rate / immersed_surface_area},\n", " seed=seed,\n", diff --git a/examples/PySDM_examples/Arabas_et_al_2023/make_particulator.py b/examples/PySDM_examples/Arabas_et_al_2023/make_particulator.py index 86f7d9ded..60cf453a8 100644 --- a/examples/PySDM_examples/Arabas_et_al_2023/make_particulator.py +++ b/examples/PySDM_examples/Arabas_et_al_2023/make_particulator.py @@ -34,6 +34,7 @@ def make_particulator( "constants": constants, "freezing_temperature_spectrum": shima_T_fz, "heterogeneous_ice_nucleation_rate": "ABIFM", + "particle_shape_and_density": "MixedPhaseSpheres", } formulae = Formulae(**formulae_ctor_args) backend = CPU(formulae) diff --git a/examples/PySDM_examples/Berry_1967/settings.py b/examples/PySDM_examples/Berry_1967/settings.py index a587c1444..6acc7f16c 100644 --- a/examples/PySDM_examples/Berry_1967/settings.py +++ b/examples/PySDM_examples/Berry_1967/settings.py @@ -25,7 +25,7 @@ def __init__(self, steps: Optional[list] = None): 1e1 * si.metres**3 ) # 1e6 -> overflows on ThrustRTC (32-bit int multiplicities) self.norm_factor = self.n_part * self.dv - self.rho = 1000 * si.kilogram / si.metre**3 + self.rho = self.formulae.constants.rho_w self.dt = 1 * si.seconds self.adaptive = False self.seed = 44 diff --git a/examples/PySDM_examples/Bieli_et_al_2022/settings.py b/examples/PySDM_examples/Bieli_et_al_2022/settings.py index de907a904..ec2c8e301 100644 --- a/examples/PySDM_examples/Bieli_et_al_2022/settings.py +++ b/examples/PySDM_examples/Bieli_et_al_2022/settings.py @@ -1,12 +1,12 @@ from pystrict import strict from PySDM.dynamics.collisions.breakup_efficiencies import ConstEb -from PySDM.dynamics.collisions.breakup_fragmentations import Feingold1988Frag +from PySDM.dynamics.collisions.breakup_fragmentations import Feingold1988 from PySDM.dynamics.collisions.coalescence_efficiencies import ConstEc from PySDM.dynamics.collisions.collision_kernels import Golovin from PySDM.formulae import Formulae from PySDM.initialisation.spectra import Gamma -from PySDM.physics.constants import si +from PySDM.physics import si from PySDM.physics.constants_defaults import rho_w @@ -28,7 +28,7 @@ def __init__(self, formulae: Formulae = None): self.vmin = 1.0 * si.um**3 self.nfmax = 10 self.fragtol = 1e-3 - self.fragmentation = Feingold1988Frag( + self.fragmentation = Feingold1988( scale=self.k * self.theta, fragtol=self.fragtol, vmin=self.vmin, diff --git a/examples/PySDM_examples/Bulenok_2023_MasterThesis/setups.py b/examples/PySDM_examples/Bulenok_2023_MasterThesis/setups.py index 53ee378e6..4d7ad3afc 100644 --- a/examples/PySDM_examples/Bulenok_2023_MasterThesis/setups.py +++ b/examples/PySDM_examples/Bulenok_2023_MasterThesis/setups.py @@ -4,7 +4,7 @@ from PySDM.dynamics import Collision from PySDM.dynamics.collisions.breakup_efficiencies import ConstEb -from PySDM.dynamics.collisions.breakup_fragmentations import ConstantSize +from PySDM.dynamics.collisions.breakup_fragmentations import ConstantMass from PySDM.dynamics.collisions.coalescence_efficiencies import ConstEc from PySDM.dynamics.collisions.collision_kernels import ConstantK from PySDM.physics import si @@ -65,7 +65,7 @@ def setup_simulation(settings, n_sd, seed, double_precision=True): collision_kernel=ConstantK(a=collision_rate), coalescence_efficiency=ConstEc(settings.srivastava_c / collision_rate), breakup_efficiency=NO_BOUNCE, - fragmentation_function=ConstantSize(c=settings.frag_mass / settings.rho), + fragmentation_function=ConstantMass(c=settings.frag_mass / settings.rho), warn_overflows=False, adaptive=False, ), diff --git a/examples/PySDM_examples/Niedermeier_et_al_2014/fig_2.ipynb b/examples/PySDM_examples/Niedermeier_et_al_2014/fig_2.ipynb index 89e9c8814..773ee0445 100644 --- a/examples/PySDM_examples/Niedermeier_et_al_2014/fig_2.ipynb +++ b/examples/PySDM_examples/Niedermeier_et_al_2014/fig_2.ipynb @@ -42,6 +42,7 @@ "source": [ "formulae = Formulae(\n", " heterogeneous_ice_nucleation_rate='ABIFM',\n", + " particle_shape_and_density=\"MixedPhaseSpheres\",\n", " constants = {\n", " 'ABIFM_M': 70,\n", " 'ABIFM_C': -10\n", diff --git a/examples/PySDM_examples/Niedermeier_et_al_2014/settings.py b/examples/PySDM_examples/Niedermeier_et_al_2014/settings.py index 824f46ed7..0dda9eda4 100644 --- a/examples/PySDM_examples/Niedermeier_et_al_2014/settings.py +++ b/examples/PySDM_examples/Niedermeier_et_al_2014/settings.py @@ -12,7 +12,7 @@ class Settings: def __init__( self, *, - formulae: Optional[Formulae] = None, + formulae: Optional[Formulae], ccn_sampling_n: int = 11, in_sampling_n: int = 20, initial_temperature: float, @@ -23,7 +23,7 @@ def __init__( self.timestep = timestep self.initial_temperature = initial_temperature - self.formulae = formulae or Formulae() + self.formulae = formulae self.initial_relative_humidity = 0.985 self.vertical_velocity = 20 * si.cm / si.s diff --git a/examples/PySDM_examples/Srivastava_1982/example.py b/examples/PySDM_examples/Srivastava_1982/example.py index a6192022a..d5b01cd3e 100644 --- a/examples/PySDM_examples/Srivastava_1982/example.py +++ b/examples/PySDM_examples/Srivastava_1982/example.py @@ -8,7 +8,7 @@ from PySDM.dynamics import Collision from PySDM.dynamics.collisions.breakup_efficiencies import ConstEb -from PySDM.dynamics.collisions.breakup_fragmentations import ConstantSize +from PySDM.dynamics.collisions.breakup_fragmentations import ConstantMass from PySDM.dynamics.collisions.coalescence_efficiencies import ConstEc from PySDM.dynamics.collisions.collision_kernels import ConstantK @@ -29,7 +29,7 @@ def coalescence_and_breakup_eq13( collision_kernel=ConstantK(a=collision_rate), coalescence_efficiency=ConstEc(settings.srivastava_c / collision_rate), breakup_efficiency=NO_BOUNCE, - fragmentation_function=ConstantSize(c=settings.frag_mass / settings.rho), + fragmentation_function=ConstantMass(c=settings.frag_mass), warn_overflows=warn_overflows, ), ) diff --git a/examples/PySDM_examples/Srivastava_1982/simulation.py b/examples/PySDM_examples/Srivastava_1982/simulation.py index 8a9c999f5..3bc076adc 100644 --- a/examples/PySDM_examples/Srivastava_1982/simulation.py +++ b/examples/PySDM_examples/Srivastava_1982/simulation.py @@ -26,7 +26,7 @@ def build(self, n_sd, seed, products): backend=self.settings.backend_class( formulae=Formulae( constants={"rho_w": self.settings.rho}, - fragmentation_function="ConstantSize", + fragmentation_function="ConstantMass", seed=seed, ), double_precision=self.double_precision, diff --git a/examples/PySDM_examples/deJong_Mackay_et_al_2023/figs_6_7_8.ipynb b/examples/PySDM_examples/deJong_Mackay_et_al_2023/figs_6_7_8.ipynb index afa77d086..daee04b17 100644 --- a/examples/PySDM_examples/deJong_Mackay_et_al_2023/figs_6_7_8.ipynb +++ b/examples/PySDM_examples/deJong_Mackay_et_al_2023/figs_6_7_8.ipynb @@ -23,8 +23,8 @@ "import numpy as np\n", "\n", "from PySDM.dynamics.collisions.coalescence_efficiencies import ConstEc, Straub2010Ec\n", - "from PySDM.dynamics.collisions.breakup_fragmentations import ExponFrag, AlwaysN, Straub2010Nf\n", - "from PySDM.physics.constants import si" + "from PySDM.dynamics.collisions.breakup_fragmentations import Exponential, AlwaysN, Straub2010Nf\n", + "from PySDM.physics import si" ] }, { @@ -42,7 +42,7 @@ "source": [ "rmin = 0.1 * si.um\n", "vmin = 4/3 * np.pi * rmin**3\n", - "settings = Settings0D(fragmentation=AlwaysN(n=8, vmin=vmin))\n", + "settings = Settings0D(fragmentation=AlwaysN(n=8))\n", "\n", "settings.n_sd = 2**13\n", "settings.radius_bins_edges = np.logspace(\n", @@ -170,7 +170,7 @@ "res = run_box_breakup(settings0, [0])\n", "ax[1].step(res.x, res.y[0]*settings0.rho, color='k', linestyle='--', label='initial (mean = X$_0$)')\n", "for (i, mu) in enumerate(mu_vals):\n", - " settings = Settings0D(fragmentation=ExponFrag(scale=mu, vmin=vmin, nfmax=nfmax))\n", + " settings = Settings0D(fragmentation=Exponential(scale=mu, vmin=vmin, nfmax=nfmax))\n", " settings._steps = [0, 120] # pylint: disable=protected-access\n", " settings.dt = 1 * si.s\n", " settings.warn_overflows = False\n", @@ -217,8 +217,8 @@ "settings.n_sd = 2**11\n", "settings.warn_overflows = False\n", "settings.radius_bins_edges = np.logspace(\n", - " np.log10(4.0 * si.um), np.log10(1e4 * si.um), num=64, endpoint=True\n", - " )\n", + " np.log10(4.0 * si.um), np.log10(1e4 * si.um), num=64, endpoint=True\n", + ")\n", "vmin = X0 * 1e-3\n", "nfmax = 10\n", "settings.coal_eff=Straub2010Ec()\n", @@ -255,7 +255,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3.9.13 ('edjPySDM')", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, diff --git a/examples/PySDM_examples/deJong_Mackay_et_al_2023/settings1D.py b/examples/PySDM_examples/deJong_Mackay_et_al_2023/settings1D.py index 2ebd3a4bb..1d64b1919 100644 --- a/examples/PySDM_examples/deJong_Mackay_et_al_2023/settings1D.py +++ b/examples/PySDM_examples/deJong_Mackay_et_al_2023/settings1D.py @@ -59,7 +59,10 @@ def __init__( frag_scale_r = 30 * si.um frag_scale_v = frag_scale_r**3 * 4 / 3 * np.pi self.fragmentation_function = Gaussian( - mu=frag_scale_v, sigma=frag_scale_v / 2, vmin=(1 * si.um) ** 3, nfmax=20 + mu=frag_scale_v, + sigma=frag_scale_v / 2, + vmin=(1 * si.um) ** 3, + nfmax=20, ) super().__init__( diff --git a/examples/PySDM_examples/deJong_Mackay_et_al_2023/settings_0D.py b/examples/PySDM_examples/deJong_Mackay_et_al_2023/settings_0D.py index b8c7c26b9..ada1bed10 100644 --- a/examples/PySDM_examples/deJong_Mackay_et_al_2023/settings_0D.py +++ b/examples/PySDM_examples/deJong_Mackay_et_al_2023/settings_0D.py @@ -3,13 +3,13 @@ import numpy as np from pystrict import strict +from PySDM.dynamics.collisions import breakup_fragmentations from PySDM.dynamics.collisions.breakup_efficiencies import ConstEb -from PySDM.dynamics.collisions.breakup_fragmentations import ExponFrag from PySDM.dynamics.collisions.coalescence_efficiencies import Berry1967 from PySDM.dynamics.collisions.collision_kernels import Geometric from PySDM.formulae import Formulae -from PySDM.initialisation.spectra import Exponential -from PySDM.physics.constants import si +from PySDM.initialisation import spectra +from PySDM.physics import si TRIVIA = Formulae().trivia @@ -37,10 +37,12 @@ def __init__( self._steps = [0] self.kernel = Geometric() self.coal_eff = Berry1967() - self.fragmentation = fragmentation or ExponFrag(scale=self.frag_scale) + self.fragmentation = fragmentation or breakup_fragmentations.Exponential( + scale=self.frag_scale + ) self.vmin = 0.0 self.break_eff = ConstEb(1.0) # no "bouncing" - self.spectrum = Exponential(norm_factor=self.norm_factor, scale=self.X0) + self.spectrum = spectra.Exponential(norm_factor=self.norm_factor, scale=self.X0) self.radius_bins_edges = np.logspace( np.log10(0.01 * si.um), np.log10(5000 * si.um), num=64, endpoint=True ) diff --git a/examples/PySDM_examples/deJong_Mackay_et_al_2023/simulation_ss.py b/examples/PySDM_examples/deJong_Mackay_et_al_2023/simulation_ss.py index 70099b231..78d1054dc 100644 --- a/examples/PySDM_examples/deJong_Mackay_et_al_2023/simulation_ss.py +++ b/examples/PySDM_examples/deJong_Mackay_et_al_2023/simulation_ss.py @@ -27,7 +27,8 @@ def run_to_steady_state(parameterization, n_sd, steps, nruns=1, dt=1 * si.s): settings = Settings0D( seed=7 ** (irun + 1), fragmentation=Straub2010Nf( - vmin=(0.01 * si.mm) ** 3 * np.pi / 6, nfmax=10000 + vmin=(0.01 * si.mm) ** 3 * np.pi / 6, + nfmax=10000, ), ) settings.coal_eff = Straub2010Ec() diff --git a/tests/smoke_tests/box/bieli_et_al_2022/test_moments.py b/tests/smoke_tests/box/bieli_et_al_2022/test_moments.py index 32d66e7cc..79fa5b824 100644 --- a/tests/smoke_tests/box/bieli_et_al_2022/test_moments.py +++ b/tests/smoke_tests/box/bieli_et_al_2022/test_moments.py @@ -12,7 +12,7 @@ def test_moments(plot=False): # arrange - settings = Settings(Formulae(fragmentation_function="Feingold1988Frag")) + settings = Settings(Formulae(fragmentation_function="Feingold1988")) if "CI" in os.environ: settings.n_sd = 10 else: diff --git a/tests/smoke_tests/box/dejong_and_mackay_et_al_2023/test_fig_6.py b/tests/smoke_tests/box/dejong_and_mackay_et_al_2023/test_fig_6.py index 70bb7df5f..381420146 100644 --- a/tests/smoke_tests/box/dejong_and_mackay_et_al_2023/test_fig_6.py +++ b/tests/smoke_tests/box/dejong_and_mackay_et_al_2023/test_fig_6.py @@ -8,10 +8,9 @@ from PySDM.backends import CPU, GPU from PySDM.dynamics.collisions.breakup_fragmentations import AlwaysN from PySDM.dynamics.collisions.coalescence_efficiencies import ConstEc, Straub2010Ec -from PySDM.physics.constants import si +from PySDM.physics import si R_MIN = 0.1 * si.um -V_MIN = 4 / 3 * np.pi * R_MIN**3 EC_VALS = [1.0, 0.9, 0.8] BINS = 32 N_SD = 2**10 @@ -21,9 +20,9 @@ "backend_class", (CPU, pytest.param(GPU, marks=pytest.mark.xfail(strict=True))), # TODO #987 ) -def test_fig_3_reduced_resolution(backend_class, plot=False): +def test_fig_6_reduced_resolution(backend_class, plot=False): # arrange - settings = Settings0D(fragmentation=AlwaysN(n=8, vmin=V_MIN), seed=44) + settings = Settings0D(fragmentation=AlwaysN(n=8), seed=44) settings.n_sd = N_SD settings.radius_bins_edges = np.logspace( np.log10(10 * si.um), np.log10(10000 * si.um), num=BINS, endpoint=True diff --git a/tests/smoke_tests/box/dejong_and_mackay_et_al_2023/test_fig_7.py b/tests/smoke_tests/box/dejong_and_mackay_et_al_2023/test_fig_7.py index a963b3b55..799c8e1be 100644 --- a/tests/smoke_tests/box/dejong_and_mackay_et_al_2023/test_fig_7.py +++ b/tests/smoke_tests/box/dejong_and_mackay_et_al_2023/test_fig_7.py @@ -6,9 +6,9 @@ from PySDM_examples.deJong_Mackay_et_al_2023 import Settings0D, run_box_breakup from PySDM.backends import CPU, GPU -from PySDM.dynamics.collisions.breakup_fragmentations import AlwaysN, ExponFrag +from PySDM.dynamics.collisions.breakup_fragmentations import AlwaysN, Exponential from PySDM.dynamics.collisions.coalescence_efficiencies import ConstEc -from PySDM.physics.constants import si +from PySDM.physics import si CMAP = matplotlib.cm.get_cmap("viridis") N_SD = 2**12 @@ -21,13 +21,13 @@ def bins_edges(num): ) -class TestFig4: +class TestFig7: @staticmethod @pytest.mark.parametrize( "backend_class", (CPU, pytest.param(GPU, marks=pytest.mark.xfail(strict=True))), # TODO #987 ) - def test_fig_4a(backend_class, plot=False): + def test_fig_7a(backend_class, plot=False): # arrange settings0 = Settings0D(seed=44) settings0.n_sd = N_SD @@ -102,7 +102,7 @@ def test_fig_4a(backend_class, plot=False): ) @staticmethod - def test_fig_4b(backend_class, plot=False): # pylint: disable=too-many-locals + def test_fig_7b(backend_class, plot=False): # pylint: disable=too-many-locals # arrange settings0 = Settings0D() settings0.n_sd = N_SD @@ -118,8 +118,10 @@ def test_fig_4b(backend_class, plot=False): # pylint: disable=too-many-locals data_x[lbl], data_y[lbl] = res.x, res.y for i, mu_val in enumerate(mu_vals): settings = Settings0D( - fragmentation=ExponFrag( - scale=mu_val, vmin=(1 * si.um) ** 3, nfmax=None + fragmentation=Exponential( + scale=mu_val, + vmin=(1 * si.um) ** 3, + nfmax=None, ), warn_overflows=False, seed=44, diff --git a/tests/smoke_tests/box/dejong_and_mackay_et_al_2023/test_fig_8.py b/tests/smoke_tests/box/dejong_and_mackay_et_al_2023/test_fig_8.py index 483e5015d..310740b05 100644 --- a/tests/smoke_tests/box/dejong_and_mackay_et_al_2023/test_fig_8.py +++ b/tests/smoke_tests/box/dejong_and_mackay_et_al_2023/test_fig_8.py @@ -15,7 +15,7 @@ "backend_class", (CPU, pytest.param(GPU, marks=pytest.mark.xfail(strict=True))), # TODO #987 ) -def test_fig_5(backend_class, plot=False): +def test_fig_8(backend_class, plot=False): # arrange settings = Settings0D( fragmentation=Straub2010Nf(vmin=Settings0D.X0 * 1e-3, nfmax=10), diff --git a/tests/smoke_tests/kinematic_2d/arabas_et_al_2015/test_freezing.py b/tests/smoke_tests/kinematic_2d/arabas_et_al_2015/test_freezing.py index 744fd645c..7dcc189d2 100644 --- a/tests/smoke_tests/kinematic_2d/arabas_et_al_2015/test_freezing.py +++ b/tests/smoke_tests/kinematic_2d/arabas_et_al_2015/test_freezing.py @@ -21,6 +21,7 @@ def test_freezing(singular): # Arrange settings = Settings( Formulae( + particle_shape_and_density="MixedPhaseSpheres", seed=44, condensation_coordinate="VolumeLogarithm", fastmath=True, diff --git a/tests/smoke_tests/parcel/arabas_and_shima_2017/test_conservation.py b/tests/smoke_tests/parcel/arabas_and_shima_2017/test_conservation.py index 5ba8ecbc3..b8cd0a0b0 100644 --- a/tests/smoke_tests/parcel/arabas_and_shima_2017/test_conservation.py +++ b/tests/smoke_tests/parcel/arabas_and_shima_2017/test_conservation.py @@ -6,15 +6,14 @@ from PySDM.backends import CPU, GPU from PySDM.backends.impl_numba.test_helpers import scipy_ode_condensation_solver -from PySDM.physics.constants_defaults import rho_w def liquid_water_mixing_ratio(simulation: Simulation): - droplet_volume = simulation.particulator.attributes["volume"].to_ndarray()[0] + droplet_mass = simulation.particulator.attributes["water mass"].to_ndarray()[0] droplet_number = simulation.particulator.attributes["multiplicity"].to_ndarray()[0] - droplet_mass = droplet_number * droplet_volume * rho_w + mass_of_all_droplets = droplet_number * droplet_mass env = simulation.particulator.environment - return droplet_mass / env.mass_of_dry_air + return mass_of_all_droplets / env.mass_of_dry_air @pytest.mark.parametrize("settings_idx", range(len(w_avgs))) @@ -51,9 +50,8 @@ def test_water_mass_conservation(settings_idx, mass_of_dry_air, scheme, coord): total_water_mixing_ratio = simulation.particulator.environment[ "water_vapour_mixing_ratio" ].to_ndarray() + liquid_water_mixing_ratio(simulation) - significant = 6 if scheme == "GPU" else 14 # TODO #540 np.testing.assert_approx_equal( - total_water_mixing_ratio, initial_total_water_mixing_ratio, significant + total_water_mixing_ratio, initial_total_water_mixing_ratio, significant=6 ) if scheme != "SciPy": assert simulation.particulator.products["S_max"].get() >= output["S"][-1] diff --git a/tests/smoke_tests/parcel/niedermeier_et_al_2013/test_temperature_profile.py b/tests/smoke_tests/parcel/niedermeier_et_al_2013/test_temperature_profile.py index 22d63cab9..bcc58535c 100644 --- a/tests/smoke_tests/parcel/niedermeier_et_al_2013/test_temperature_profile.py +++ b/tests/smoke_tests/parcel/niedermeier_et_al_2013/test_temperature_profile.py @@ -11,6 +11,7 @@ def test_temperature_profile(initial_temperature, plot=False): # arrange formulae = Formulae( + particle_shape_and_density="MixedPhaseSpheres", heterogeneous_ice_nucleation_rate="ABIFM", constants={"ABIFM_M": 54.48, "ABIFM_C": -10.67}, ) diff --git a/tests/smoke_tests/parcel/yang_et_al_2018/test_just_do_it.py b/tests/smoke_tests/parcel/yang_et_al_2018/test_just_do_it.py index f98662348..72918b580 100644 --- a/tests/smoke_tests/parcel/yang_et_al_2018/test_just_do_it.py +++ b/tests/smoke_tests/parcel/yang_et_al_2018/test_just_do_it.py @@ -5,7 +5,7 @@ from PySDM.backends import CPU, GPU from PySDM.backends.impl_numba.test_helpers import scipy_ode_condensation_solver -from PySDM.physics.constants import si +from PySDM.physics import si # TODO #527 @@ -15,7 +15,7 @@ def test_just_do_it(scheme, adaptive, backend_class=CPU): # Arrange if scheme == "SciPy" and (not adaptive or backend_class is GPU): - return + pytest.skip() settings = Settings(dt_output=10 * si.second) settings.adaptive = adaptive diff --git a/tests/unit_tests/attributes/test_fall_velocity.py b/tests/unit_tests/attributes/test_fall_velocity.py index a1811d281..d39184059 100644 --- a/tests/unit_tests/attributes/test_fall_velocity.py +++ b/tests/unit_tests/attributes/test_fall_velocity.py @@ -73,11 +73,7 @@ def test_fall_velocity_calculation(default_attributes, backend_class): assert np.allclose( particulator.attributes["relative fall velocity"].to_ndarray(), particulator.attributes["relative fall momentum"].to_ndarray() - / ( - # TODO #798 - we plan to use masses instead of volumes soon - particulator.formulae.constants.rho_w - * particulator.attributes["volume"].to_ndarray() - ), + / (particulator.attributes["water mass"].to_ndarray()), ) diff --git a/tests/unit_tests/backends/test_freezing_methods.py b/tests/unit_tests/backends/test_freezing_methods.py index 93749f908..2fdcd5a0e 100644 --- a/tests/unit_tests/backends/test_freezing_methods.py +++ b/tests/unit_tests/backends/test_freezing_methods.py @@ -26,7 +26,7 @@ def test_no_subsaturated_freezing(self): @pytest.mark.parametrize("epsilon", (0, 1e-5)) def test_thaw(backend_class, singular, thaw, epsilon): # arrange - formulae = Formulae() + formulae = Formulae(particle_shape_and_density="MixedPhaseSpheres") builder = Builder(n_sd=1, backend=backend_class(formulae=formulae)) env = Box(dt=1 * si.s, dv=1 * si.m**3) builder.set_environment(env) @@ -71,7 +71,7 @@ def test_freeze_singular(backend_class): multiplicity = 1e10 steps = 1 - formulae = Formulae() + formulae = Formulae(particle_shape_and_density="MixedPhaseSpheres") builder = Builder(n_sd=n_sd, backend=backend_class(formulae=formulae)) env = Box(dt=dt, dv=dv) builder.set_environment(env) @@ -136,6 +136,7 @@ def low(t): output = {} formulae = Formulae( + particle_shape_and_density="MixedPhaseSpheres", heterogeneous_ice_nucleation_rate="Constant", constants={"J_HET": rate / immersed_surface_area}, seed=seed, diff --git a/tests/unit_tests/backends/test_physics_methods.py b/tests/unit_tests/backends/test_physics_methods.py index 7a9135cc3..09fcd5d08 100644 --- a/tests/unit_tests/backends/test_physics_methods.py +++ b/tests/unit_tests/backends/test_physics_methods.py @@ -1,10 +1,12 @@ # pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring import numpy as np +import pytest +from PySDM import Formulae from PySDM.physics import si -class TestPhysicsMethods: # pylint: disable=too-few-public-methods +class TestPhysicsMethods: @staticmethod def test_temperature_pressure_RH(backend_class): # Arrange @@ -34,3 +36,69 @@ def test_temperature_pressure_RH(backend_class): assert 282 * si.K < T.amin() < 283 * si.K assert 820 * si.hPa < p.amin() < 830 * si.hPa assert 1.10 < RH.amin() < 1.11 + + @staticmethod + @pytest.mark.parametrize("variant", ("LiquidSpheres", "MixedPhaseSpheres")) + def test_mass_to_volume(backend_class, variant): + # Arrange + formulae = Formulae(particle_shape_and_density=variant) + backend = backend_class(formulae, double_precision=True) + sut = backend.volume_of_water_mass + mass = np.asarray([1.0, -1.0]) + mass_in = backend.Storage.from_ndarray(mass) + volume_out = backend.Storage.from_ndarray(np.zeros_like(mass_in)) + + # Act + sut(volume=volume_out, mass=mass_in) + + # Assert + assert (mass_in.to_ndarray() == mass).all() + if variant == "LiquidSpheres": + assert ( + volume_out.to_ndarray() + == mass_in.to_ndarray() / formulae.constants.rho_w + ).all() + elif variant == "MixedPhaseSpheres": + assert ( + volume_out.to_ndarray() + == np.where( + mass < 0, + mass / formulae.constants.rho_i, + mass / formulae.constants.rho_w, + ) + ).all() + else: + raise NotImplementedError() + + @staticmethod + @pytest.mark.parametrize("variant", ("LiquidSpheres", "MixedPhaseSpheres")) + def test_volume_to_mass(backend_class, variant): + # Arrange + formulae = Formulae(particle_shape_and_density=variant) + backend = backend_class(formulae, double_precision=True) + sut = backend.mass_of_water_volume + volume = np.asarray([1.0, -1.0]) + volume_in = backend.Storage.from_ndarray(volume) + mass_out = backend.Storage.from_ndarray(np.zeros_like(volume_in)) + + # Act + sut(volume=volume_in, mass=mass_out) + + # Assert + assert (volume_in.to_ndarray() == volume).all() + if variant == "LiquidSpheres": + assert ( + mass_out.to_ndarray() + == volume_in.to_ndarray() * formulae.constants.rho_w + ).all() + elif variant == "MixedPhaseSpheres": + assert ( + mass_out.to_ndarray() + == np.where( + volume < 0, + volume * formulae.constants.rho_i, + volume * formulae.constants.rho_w, + ) + ).all() + else: + raise NotImplementedError() diff --git a/tests/unit_tests/dummy_particulator.py b/tests/unit_tests/dummy_particulator.py index 2072dd4f5..53db9258d 100644 --- a/tests/unit_tests/dummy_particulator.py +++ b/tests/unit_tests/dummy_particulator.py @@ -8,7 +8,7 @@ class DummyParticulator(Builder, Particulator): def __init__(self, backend_class, n_sd=0, formulae=None, grid=None): - backend = backend_class(formulae) + backend = backend_class(formulae, double_precision=True) Builder.__init__(self, n_sd, backend) Particulator.__init__(self, n_sd, backend) self.particulator = self diff --git a/tests/unit_tests/dynamics/collisions/test_fragmentations.py b/tests/unit_tests/dynamics/collisions/test_fragmentations.py index fba6facfa..7c6331f12 100644 --- a/tests/unit_tests/dynamics/collisions/test_fragmentations.py +++ b/tests/unit_tests/dynamics/collisions/test_fragmentations.py @@ -1,6 +1,4 @@ # pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring -from timeit import default_timer as timer - import matplotlib.pyplot as plt import numpy as np import pytest @@ -10,14 +8,14 @@ from PySDM.dynamics.collisions.breakup_fragmentations import ( SLAMS, AlwaysN, - ConstantSize, - ExponFrag, - Feingold1988Frag, + ConstantMass, + Exponential, + Feingold1988, Gaussian, Straub2010Nf, ) from PySDM.environments import Box -from PySDM.physics import si +from PySDM.physics import constants_defaults, si ARBITRARY_VALUE_BETWEEN_0_AND_1 = 0.5 @@ -34,9 +32,12 @@ class TestFragmentations: # pylint: disable=too-few-public-methods "fragmentation_fn", ( AlwaysN(n=2), - ExponFrag(scale=1e6 * si.um**3), - Feingold1988Frag(scale=1e6 * si.um**3), - Gaussian(mu=2e6 * si.um**3, sigma=1e6 * si.um**3), + Exponential(scale=1e6 * si.um**3), + Feingold1988(scale=1e6 * si.um**3), + Gaussian( + mu=2e6 * si.um**3, + sigma=1e6 * si.um**3, + ), SLAMS(), Straub2010Nf(), ), @@ -62,7 +63,7 @@ def test_fragmentation_fn_call(fragmentation_fn, backend_class): _PairwiseStorage = builder.particulator.PairwiseStorage _Indicator = builder.particulator.PairIndicator nf = _PairwiseStorage.from_ndarray(np.zeros_like(fragments)) - frag_size = _PairwiseStorage.from_ndarray(np.zeros_like(fragments)) + frag_mass = _PairwiseStorage.from_ndarray(np.zeros_like(fragments)) is_first_in_pair = _Indicator(length=volume.size) is_first_in_pair.indicator = builder.particulator.Storage.from_ndarray( np.asarray([True, False]) @@ -70,19 +71,33 @@ def test_fragmentation_fn_call(fragmentation_fn, backend_class): u01 = dummy_u01(builder, fragments.size) # act - sut(nf, frag_size, u01, is_first_in_pair) + sut(nf, frag_mass, u01, is_first_in_pair) # Assert np.testing.assert_array_less([0.99], nf.to_ndarray()) - np.testing.assert_array_less([0.0], frag_size.to_ndarray()) + np.testing.assert_array_less([0.0], frag_mass.to_ndarray()) + + np.testing.assert_approx_equal( + nf[0] * frag_mass[0], np.sum(volume) * constants_defaults.rho_w + ) @staticmethod @pytest.mark.parametrize( "fragmentation_fn", [ - ExponFrag(scale=1 * si.um**3, vmin=6660.0 * si.um**3), - Feingold1988Frag(scale=1 * si.um**3, vmin=6660.0 * si.um**3), - Gaussian(mu=2 * si.um**3, sigma=1 * si.um**3, vmin=6660.0 * si.um**3), + Exponential( + scale=1 * si.um**3, + vmin=6660.0 * si.um**3, + ), + Feingold1988( + scale=1 * si.um**3, + vmin=6660.0 * si.um**3, + ), + Gaussian( + mu=2 * si.um**3, + sigma=1 * si.um**3, + vmin=6660.0 * si.um**3, + ), SLAMS(vmin=6660.0 * si.um**3), Straub2010Nf(vmin=6660.0 * si.um**3), pytest.param(AlwaysN(n=10), marks=pytest.mark.xfail(strict=True)), @@ -109,7 +124,7 @@ def test_fragmentation_limiters_vmin(fragmentation_fn, backend_class): _PairwiseStorage = builder.particulator.PairwiseStorage _Indicator = builder.particulator.PairIndicator nf = _PairwiseStorage.from_ndarray(np.zeros_like(fragments)) - frag_size = _PairwiseStorage.from_ndarray(np.zeros_like(fragments)) + frag_mass = _PairwiseStorage.from_ndarray(np.zeros_like(fragments)) is_first_in_pair = _Indicator(length=volume.size) is_first_in_pair.indicator = builder.particulator.Storage.from_ndarray( np.asarray([True, False]) @@ -117,21 +132,25 @@ def test_fragmentation_limiters_vmin(fragmentation_fn, backend_class): u01 = dummy_u01(builder, fragments.size) # act - sut(nf, frag_size, u01, is_first_in_pair) + sut(nf, frag_mass, u01, is_first_in_pair) # Assert np.testing.assert_array_equal([1.0], nf.to_ndarray()) np.testing.assert_array_equal( - [(6660.0 + 440.0) * si.um**3], frag_size.to_ndarray() + [(6660.0 + 440.0) * si.um**3 * constants_defaults.rho_w], + frag_mass.to_ndarray(), ) @staticmethod @pytest.mark.parametrize( "fragmentation_fn", [ - ExponFrag(scale=1.0 * si.cm**3), - Feingold1988Frag(scale=1.0 * si.cm**3), - Gaussian(mu=1.0 * si.cm**3, sigma=1e6 * si.um**3), + Exponential(scale=1.0 * si.cm**3), + Feingold1988(scale=1.0 * si.cm**3), + Gaussian( + mu=1.0 * si.cm**3, + sigma=1e6 * si.um**3, + ), SLAMS(), Straub2010Nf(), pytest.param(AlwaysN(n=0.01), marks=pytest.mark.xfail(strict=True)), @@ -158,7 +177,7 @@ def test_fragmentation_limiters_vmax(fragmentation_fn, backend_class): _PairwiseStorage = builder.particulator.PairwiseStorage _Indicator = builder.particulator.PairIndicator nf = _PairwiseStorage.from_ndarray(np.zeros_like(fragments)) - frag_size = _PairwiseStorage.from_ndarray(np.zeros_like(fragments)) + frag_mass = _PairwiseStorage.from_ndarray(np.zeros_like(fragments)) is_first_in_pair = _Indicator(length=volume.size) is_first_in_pair.indicator = builder.particulator.Storage.from_ndarray( np.asarray([True, False]) @@ -166,21 +185,30 @@ def test_fragmentation_limiters_vmax(fragmentation_fn, backend_class): u01 = dummy_u01(builder, fragments.size) # act - sut(nf, frag_size, u01, is_first_in_pair) + sut(nf, frag_mass, u01, is_first_in_pair) # Assert np.testing.assert_array_less([0.999], nf.to_ndarray()) np.testing.assert_array_less( - frag_size.to_ndarray(), [(6661.0 + 440.0) * si.um**3] + frag_mass.to_ndarray(), + [(6661.0 + 440.0) * si.um**3 * constants_defaults.rho_w], + ) + + np.testing.assert_approx_equal( + nf[0] * frag_mass[0], np.sum(volume) * constants_defaults.rho_w ) @staticmethod @pytest.mark.parametrize( "fragmentation_fn", [ - ExponFrag(scale=1.0 * si.um**3, nfmax=2), - Feingold1988Frag(scale=1.0 * si.um**3, nfmax=2), - Gaussian(mu=1.0 * si.um**3, sigma=1e6 * si.um**3, nfmax=2), + Exponential(scale=1.0 * si.um**3, nfmax=2), + Feingold1988(scale=1.0 * si.um**3, nfmax=2), + Gaussian( + mu=1.0 * si.um**3, + sigma=1e6 * si.um**3, + nfmax=2, + ), SLAMS(nfmax=2), Straub2010Nf(nfmax=2), pytest.param(AlwaysN(n=10), marks=pytest.mark.xfail(strict=True)), @@ -207,7 +235,7 @@ def test_fragmentation_limiters_nfmax(fragmentation_fn, backend_class): _PairwiseStorage = builder.particulator.PairwiseStorage _Indicator = builder.particulator.PairIndicator nf = _PairwiseStorage.from_ndarray(np.zeros_like(fragments)) - frag_size = _PairwiseStorage.from_ndarray(np.zeros_like(fragments)) + frag_mass = _PairwiseStorage.from_ndarray(np.zeros_like(fragments)) is_first_in_pair = _Indicator(length=volume.size) is_first_in_pair.indicator = builder.particulator.Storage.from_ndarray( np.asarray([True, False]) @@ -215,12 +243,16 @@ def test_fragmentation_limiters_nfmax(fragmentation_fn, backend_class): u01 = dummy_u01(builder, fragments.size) # act - sut(nf, frag_size, u01, is_first_in_pair) + sut(nf, frag_mass, u01, is_first_in_pair) # Assert np.testing.assert_array_less(nf.to_ndarray(), [2.0 + 1e-6]) np.testing.assert_array_less( - [((6660.0 + 440.0) / 2 - 1) * si.um**3], frag_size.to_ndarray() + [((6660.0 + 440.0) / 2 - 1) * si.um**3], frag_mass.to_ndarray() + ) + + np.testing.assert_approx_equal( + nf[0] * frag_mass[0], np.sum(volume) * constants_defaults.rho_w ) @staticmethod @@ -228,9 +260,12 @@ def test_fragmentation_limiters_nfmax(fragmentation_fn, backend_class): "fragmentation_fn", ( AlwaysN(n=2), - ExponFrag(scale=1e6 * si.um**3), - Feingold1988Frag(scale=1e6 * si.um**3), - Gaussian(mu=2e6 * si.um**3, sigma=1e6 * si.um**3), + Exponential(scale=1e6 * si.um**3), + Feingold1988(scale=1e6 * si.um**3), + Gaussian( + mu=2e6 * si.um**3, + sigma=1e6 * si.um**3, + ), SLAMS(), Straub2010Nf(), ), @@ -269,51 +304,35 @@ def test_fragmentation_fn_distribution( rns = np.linspace(1e-6, 1 - 1e-6, n) for i, rn in enumerate(rns): - print("i", i) - start = timer() - _PairwiseStorage = builder.particulator.PairwiseStorage _Indicator = builder.particulator.PairIndicator nf = _PairwiseStorage.from_ndarray( np.zeros_like(fragments, dtype=np.double) ) - frag_size = _PairwiseStorage.from_ndarray( + frag_mass = _PairwiseStorage.from_ndarray( np.zeros_like(fragments, dtype=np.double) ) is_first_in_pair = _Indicator(length=volume.size) is_first_in_pair.indicator = builder.particulator.Storage.from_ndarray( np.asarray([True, False]) ) - u01 = _PairwiseStorage.from_ndarray( - np.asarray([rn]) - ) # (np.random.rand(*fragments.shape)) - print("u01", u01.data) - - end = timer() - print("elapsed time setup", end - start) - - start = timer() + u01 = _PairwiseStorage.from_ndarray(np.asarray([rn])) # act - sut(nf, frag_size, u01, is_first_in_pair) + sut(nf, frag_mass, u01, is_first_in_pair) - end = timer() - print("elapsed time sut", end - start) - print(nf.data) - print(frag_size.data) res[i][0] = nf[0] - res[i][1] = frag_size[0] + res[i][1] = frag_mass[0] # Assert np.testing.assert_array_less([0.99], nf.to_ndarray()) - np.testing.assert_array_less([0.0], frag_size.to_ndarray()) + np.testing.assert_array_less([0.0], frag_mass.to_ndarray()) + + np.testing.assert_approx_equal( + nf[0] * frag_mass[0], np.sum(volume) * constants_defaults.rho_w + ) res = np.asarray(sorted(res, key=lambda x: x[1], reverse=True)) - print(res[:, 0]) - unique_nfs, nfs_counts = np.unique(res[:, 0], return_counts=True) - unique_frag_size, frag_sizes_counts = np.unique(res[:, 1], return_counts=True) - print("nfs", unique_nfs, nfs_counts) - print("frag_sizes", unique_frag_size, frag_sizes_counts) plt.hist(res[:, 0]) if plot: @@ -321,28 +340,42 @@ def test_fragmentation_fn_distribution( @staticmethod @pytest.mark.parametrize( - "fragmentation_fn, volume, expected_nf", + "fragmentation_fn, water_mass, expected_nf", ( ( - ConstantSize(c=4 * si.um**3), - np.asarray([400.0 * si.um**3, 600.0 * si.um**3]), + ConstantMass(c=4 * si.um**3), + np.asarray( + [ + 400.0 * si.um**3, + 600.0 * si.um**3, + ] + ), + 250, + ), + ( + AlwaysN(n=250), + np.asarray( + [ + 400.0 * si.um**3, + 600.0 * si.um**3, + ] + ), 250, ), - (AlwaysN(n=250), np.asarray([400.0 * si.um**3, 600.0 * si.um**3]), 250), ), ) - def test_fragmentation_nf_and_frag_size_equals( # TODO #987 + def test_fragmentation_nf_and_frag_mass_equals( # TODO #987 fragmentation_fn, - volume, + water_mass, expected_nf, backend_class=CPU, ): # arrange - expected_frag_size = np.sum(volume) / expected_nf + expected_frag_mass = np.sum(water_mass) / expected_nf fragments = np.asarray([-1.0]) builder = Builder( - volume.size, + water_mass.size, backend_class( Formulae(fragmentation_function=fragmentation_fn.__class__.__name__) ), @@ -352,22 +385,27 @@ def test_fragmentation_nf_and_frag_size_equals( # TODO #987 sut.register(builder) builder.set_environment(Box(dv=None, dt=None)) _ = builder.build( - attributes={"volume": volume, "multiplicity": np.ones_like(volume)} + attributes={ + "water mass": water_mass, + "multiplicity": np.ones_like(water_mass), + } ) _PairwiseStorage = builder.particulator.PairwiseStorage _Indicator = builder.particulator.PairIndicator nf = _PairwiseStorage.from_ndarray(np.zeros_like(fragments)) - frag_size = _PairwiseStorage.from_ndarray(np.zeros_like(fragments)) - is_first_in_pair = _Indicator(length=volume.size) + frag_mass = _PairwiseStorage.from_ndarray(np.zeros_like(fragments)) + is_first_in_pair = _Indicator(length=water_mass.size) is_first_in_pair.indicator = builder.particulator.Storage.from_ndarray( np.asarray([True, False]) ) u01 = _PairwiseStorage.from_ndarray(np.ones_like(fragments)) # act - sut(nf, frag_size, u01, is_first_in_pair) + sut(nf, frag_mass, u01, is_first_in_pair) # Assert - np.testing.assert_array_equal(nf.to_ndarray(), expected_nf) - np.testing.assert_array_equal([expected_frag_size], frag_size.to_ndarray()) + np.testing.assert_approx_equal(nf.to_ndarray(), expected_nf) + np.testing.assert_array_almost_equal( + [expected_frag_mass], frag_mass.to_ndarray() + ) diff --git a/tests/unit_tests/dynamics/collisions/test_sdm_breakup.py b/tests/unit_tests/dynamics/collisions/test_sdm_breakup.py index dacc7fa4c..b01d713ae 100644 --- a/tests/unit_tests/dynamics/collisions/test_sdm_breakup.py +++ b/tests/unit_tests/dynamics/collisions/test_sdm_breakup.py @@ -7,19 +7,29 @@ from PySDM.backends import CPU from PySDM.backends.impl_common.pair_indicator import make_PairIndicator from PySDM.dynamics import Breakup +from PySDM.dynamics.collisions import breakup_fragmentations from PySDM.dynamics.collisions.breakup_efficiencies import ConstEb -from PySDM.dynamics.collisions.breakup_fragmentations import AlwaysN, ExponFrag +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 -from PySDM.initialisation.spectra import Exponential from PySDM.physics import si from PySDM.physics.trivia import Trivia from PySDM.products.size_spectral import ParticleVolumeVersusRadiusLogarithmSpectrum +def volume_to_mass(particulator, volume): + if isinstance(volume, (list, tuple)): + return list( + map(particulator.formulae.particle_shape_and_density.volume_to_mass, volume) + ) + + raise NotImplementedError + + class TestSDMBreakup: @staticmethod @pytest.mark.parametrize( @@ -102,7 +112,7 @@ def test_single_collision_bounce(params, backend_class): gamma = particulator.PairwiseStorage.from_ndarray(np.asarray([params["gamma"]])) rand = particulator.PairwiseStorage.from_ndarray(np.asarray([params["rand"]])) - fragment_size = particulator.PairwiseStorage.from_ndarray( + fragment_mass = particulator.PairwiseStorage.from_ndarray( np.array([50 * si.um**3], dtype=float) ) is_first_in_pair = make_PairIndicator(backend)(n_sd) @@ -114,7 +124,7 @@ def test_single_collision_bounce(params, backend_class): rand=rand, Ec=pairwise_zeros, Eb=pairwise_zeros, - fragment_size=fragment_size, + fragment_mass=fragment_mass, coalescence_rate=general_zeros, breakup_rate=general_zeros, breakup_rate_deficit=general_zeros, @@ -190,7 +200,9 @@ def test_breakup_counters(params, backend_class): # pylint: disable=too-many-lo np.array([params["Eb"]] * n_pairs) ) breakup_rate = particulator.Storage.from_ndarray(np.array([0])) - frag_size = particulator.PairwiseStorage.from_ndarray(np.array([2.0] * n_pairs)) + frag_mass = particulator.PairwiseStorage.from_ndarray( + np.array(volume_to_mass(particulator, [2 * si.m**3]) * n_pairs) + ) is_first_in_pair = particulator.PairIndicator(n_sd) is_first_in_pair.indicator[:] = particulator.Storage.from_ndarray( np.asarray(params["is_first_in_pair"]) @@ -203,7 +215,7 @@ def test_breakup_counters(params, backend_class): # pylint: disable=too-many-lo rand=rand, Ec=pairwise_zeros, Eb=Eb, - fragment_size=frag_size, + fragment_mass=frag_mass, coalescence_rate=general_zeros, breakup_rate=breakup_rate, breakup_rate_deficit=general_zeros, @@ -230,7 +242,7 @@ def test_breakup_counters(params, backend_class): # pylint: disable=too-many-lo "v_expected": [0.5, 0.5], "expected_deficit": [0.0], "is_first_in_pair": [True, False], - "frag_size": [0.5], + "frag_volume": [0.5], }, { "gamma": [2.0], @@ -240,7 +252,7 @@ def test_breakup_counters(params, backend_class): # pylint: disable=too-many-lo "v_expected": [1, 1], "expected_deficit": [0.0], "is_first_in_pair": [True, False], - "frag_size": [1.0], + "frag_volume": [1.0], }, { "gamma": [2.0], @@ -250,7 +262,7 @@ def test_breakup_counters(params, backend_class): # pylint: disable=too-many-lo "v_expected": [0.5, 0.5], "expected_deficit": [1.0], "is_first_in_pair": [True, False], - "frag_size": [0.5], + "frag_volume": [0.5], }, { "gamma": [2.0], @@ -260,7 +272,27 @@ def test_breakup_counters(params, backend_class): # pylint: disable=too-many-lo "v_expected": [1.0, 0.5], "expected_deficit": [1.0], "is_first_in_pair": [True, False], - "frag_size": [0.5], + "frag_volume": [0.5], + }, + { + "gamma": [2.0], + "n_init": [9, 2], + "v_init": [1, 2], + "n_expected": [1, 12], + "v_expected": [1, 1], + "expected_deficit": [0.0], + "is_first_in_pair": [True, False], + "frag_volume": [1.0], + }, + { + "gamma": [1.0], + "n_init": [12, 1], + "v_init": [1, 1], + "n_expected": [11, 2], + "v_expected": [1, 1], + "expected_deficit": [0.0], + "is_first_in_pair": [True, False], + "frag_volume": [1.0], }, { "gamma": [1.0], @@ -270,7 +302,7 @@ def test_breakup_counters(params, backend_class): # pylint: disable=too-many-lo "v_expected": [2, 4], "expected_deficit": [0.0], "is_first_in_pair": [True, False], - "frag_size": [4], + "frag_volume": [4], }, { "gamma": [1.0], @@ -280,7 +312,7 @@ def test_breakup_counters(params, backend_class): # pylint: disable=too-many-lo "v_expected": [2, 4], "expected_deficit": [0.0], "is_first_in_pair": [True, False], - "frag_size": [4], + "frag_volume": [4], }, { "gamma": [3.0], @@ -290,7 +322,7 @@ def test_breakup_counters(params, backend_class): # pylint: disable=too-many-lo "v_expected": [2, 4], "expected_deficit": [0.0], "is_first_in_pair": [True, False], - "frag_size": [4], + "frag_volume": [4], }, { "gamma": [0.0], @@ -300,7 +332,7 @@ def test_breakup_counters(params, backend_class): # pylint: disable=too-many-lo "v_expected": [2, 6], "expected_deficit": [0.0], "is_first_in_pair": [True, False], - "frag_size": [4], + "frag_volume": [4], }, ], ) @@ -311,7 +343,7 @@ def test_attribute_update_single_breakup( # Arrange n_init = params["n_init"] n_sd = len(n_init) - builder = Builder(n_sd, backend_class()) + builder = Builder(n_sd, backend_class(double_precision=True)) builder.set_environment(Box(dv=np.NaN, dt=np.NaN)) particulator = builder.build( attributes={ @@ -334,8 +366,10 @@ def test_attribute_update_single_breakup( Eb = particulator.PairwiseStorage.from_ndarray(np.array(Eb)) breakup_rate = particulator.Storage.from_ndarray(np.array([0])) breakup_rate_deficit = particulator.Storage.from_ndarray(np.array([0])) - frag_size = particulator.PairwiseStorage.from_ndarray( - np.asarray(params["frag_size"], dtype=float) + + frag_mass = volume_to_mass(particulator, params["frag_volume"]) + frag_mass = particulator.PairwiseStorage.from_ndarray( + np.asarray(frag_mass, dtype=float) ) is_first_in_pair = particulator.PairIndicator(n_sd) is_first_in_pair.indicator[:] = particulator.Storage.from_ndarray( @@ -349,7 +383,7 @@ def test_attribute_update_single_breakup( rand=rand, Ec=pairwise_zeros, Eb=Eb, - fragment_size=frag_size, + fragment_mass=frag_mass, coalescence_rate=general_zeros, breakup_rate=breakup_rate, breakup_rate_deficit=breakup_rate_deficit, @@ -364,18 +398,18 @@ def test_attribute_update_single_breakup( particulator.attributes["multiplicity"].to_ndarray(), np.array(params["n_expected"]), ), - "v": lambda: np.testing.assert_array_equal( + "v": lambda: np.testing.assert_array_almost_equal( particulator.attributes["volume"].to_ndarray(), np.array(params["v_expected"]), ), - "conserve": lambda: np.testing.assert_equal( + "conserve": lambda: np.testing.assert_almost_equal( np.sum( particulator.attributes["multiplicity"].to_ndarray() * particulator.attributes["volume"].to_ndarray() ), np.sum(np.array(params["n_init"]) * np.array(params["v_init"])), ), - "deficit": lambda: np.testing.assert_equal( + "deficit": lambda: np.testing.assert_almost_equal( breakup_rate_deficit.to_ndarray(), np.array(params["expected_deficit"]) ), }[flag]() @@ -389,33 +423,33 @@ def test_attribute_update_single_breakup( "n_init": [64, 2], "v_init": [128, 128], "is_first_in_pair": [True, False], - "frag_size": [128], + "frag_volume": [128], }, { "n_init": [20, 4], "v_init": [1, 2], "is_first_in_pair": [True, False], - "frag_size": [1.0], + "frag_volume": [1.0], }, { "n_init": [3, 1], "v_init": [1, 1], "is_first_in_pair": [True, False], - "frag_size": [0.5], + "frag_volume": [0.5], }, { "n_init": [64, 2], "v_init": [8, 16], "is_first_in_pair": [True, False], "n_fragment": [6], - "frag_size": [4.0], + "frag_volume": [4.0], }, { "n_init": [64, 2], "v_init": [6, 2], "is_first_in_pair": [True, False], "n_fragment": [2], - "frag_size": [4.0], + "frag_volume": [4.0], }, ], ) @@ -424,7 +458,7 @@ def test_attribute_update_n_breakups( ): # pylint: disable=too-many-locals # Arrange - assert len(params["frag_size"]) == 1 + assert len(params["frag_volume"]) == 1 def run_simulation(_n_times, _gamma): n_init = params["n_init"] @@ -452,8 +486,10 @@ def run_simulation(_n_times, _gamma): Eb = particulator.PairwiseStorage.from_ndarray(np.array(Eb)) breakup_rate = particulator.Storage.from_ndarray(np.array([0])) breakup_rate_deficit = particulator.Storage.from_ndarray(np.array([0])) - frag_size = particulator.PairwiseStorage.from_ndarray( - np.asarray(params["frag_size"], dtype=float) + frag_mass = particulator.PairwiseStorage.from_ndarray( + np.asarray( + volume_to_mass(particulator, params["frag_volume"]), dtype=float + ) ) is_first_in_pair = particulator.PairIndicator(n_sd) is_first_in_pair.indicator[:] = particulator.Storage.from_ndarray( @@ -468,7 +504,7 @@ def run_simulation(_n_times, _gamma): rand=rand, Ec=pairwise_zeros, Eb=Eb, - fragment_size=frag_size, + fragment_mass=frag_mass, coalescence_rate=general_zeros, breakup_rate=breakup_rate, breakup_rate_deficit=breakup_rate_deficit, @@ -485,11 +521,11 @@ def run_simulation(_n_times, _gamma): run2 = run_simulation(_n_times=_n, _gamma=[1]) # Assert - np.testing.assert_array_equal( + np.testing.assert_array_almost_equal( run1[0], run2[0], ) - np.testing.assert_array_equal( + np.testing.assert_array_almost_equal( run1[1], run2[1], ) @@ -502,7 +538,7 @@ def test_multiplicity_overflow(backend=CPU()): # pylint: disable=too-many-local "n_init": [1, 3], "v_init": [1, 1], "is_first_in_pair": [True, False], - "frag_size": [2e-10], + "frag_volume": [2e-10], } n_init = params["n_init"] n_sd = len(n_init) @@ -529,8 +565,8 @@ def test_multiplicity_overflow(backend=CPU()): # pylint: disable=too-many-local Eb = particulator.PairwiseStorage.from_ndarray(np.array(Eb)) breakup_rate = particulator.Storage.from_ndarray(np.array([0])) breakup_rate_deficit = particulator.Storage.from_ndarray(np.array([0])) - frag_size = particulator.PairwiseStorage.from_ndarray( - np.array(params["frag_size"], dtype=float) + frag_mass = particulator.PairwiseStorage.from_ndarray( + np.array(volume_to_mass(particulator, params["frag_volume"]), dtype=float) ) is_first_in_pair = particulator.PairIndicator(n_sd) is_first_in_pair.indicator[:] = particulator.Storage.from_ndarray( @@ -544,7 +580,7 @@ def test_multiplicity_overflow(backend=CPU()): # pylint: disable=too-many-local rand=rand, Ec=pairwise_zeros, Eb=Eb, - fragment_size=frag_size, + fragment_mass=frag_mass, coalescence_rate=general_zeros, breakup_rate=breakup_rate, breakup_rate_deficit=breakup_rate_deficit, @@ -553,7 +589,7 @@ def test_multiplicity_overflow(backend=CPU()): # pylint: disable=too-many-local max_multiplicity=DEFAULTS.max_multiplicity, ) assert breakup_rate_deficit[0] > 0 - np.testing.assert_equal( + np.testing.assert_almost_equal( np.sum( particulator.attributes["multiplicity"].to_ndarray() * particulator.attributes["volume"].to_ndarray() @@ -571,7 +607,7 @@ def test_same_multiplicity_overflow_no_substeps( "n_init": [1, 1], "v_init": [1, 1], "is_first_in_pair": [True, False], - "frag_size": [0.5], + "frag_volume": [0.5], } n_init = params["n_init"] n_sd = len(n_init) @@ -598,8 +634,8 @@ def test_same_multiplicity_overflow_no_substeps( Eb = particulator.PairwiseStorage.from_ndarray(np.array(Eb)) breakup_rate = particulator.Storage.from_ndarray(np.array([0])) breakup_rate_deficit = particulator.Storage.from_ndarray(np.array([0])) - frag_size = particulator.PairwiseStorage.from_ndarray( - np.array(params["frag_size"], dtype=float) + frag_mass = particulator.PairwiseStorage.from_ndarray( + np.array(volume_to_mass(particulator, params["frag_volume"]), dtype=float) ) is_first_in_pair = particulator.PairIndicator(n_sd) is_first_in_pair.indicator[:] = particulator.Storage.from_ndarray( @@ -613,7 +649,7 @@ def test_same_multiplicity_overflow_no_substeps( rand=rand, Ec=pairwise_zeros, Eb=Eb, - fragment_size=frag_size, + fragment_mass=frag_mass, coalescence_rate=general_zeros, breakup_rate=breakup_rate, breakup_rate_deficit=breakup_rate_deficit, @@ -642,7 +678,7 @@ def test_same_multiplicity_overflow_no_substeps( "v_expected": [1, 1], "expected_deficit": [0.0], "is_first_in_pair": [True, False], - "frag_size": [1.25], + "frag_volume": [1.25], }, { "gamma": [1.0], @@ -652,7 +688,7 @@ def test_same_multiplicity_overflow_no_substeps( "v_expected": [1, 1], "expected_deficit": [0.0], "is_first_in_pair": [True, False], - "frag_size": [1 / 1.3], + "frag_volume": [1 / 1.3], }, { "gamma": [2.0], @@ -662,7 +698,7 @@ def test_same_multiplicity_overflow_no_substeps( "v_expected": [1, 2 / 3], "expected_deficit": [1.0], "is_first_in_pair": [True, False], - "frag_size": [1 / 1.4], + "frag_volume": [1 / 1.4], }, ], ) @@ -696,8 +732,8 @@ def test_noninteger_fragments( Eb = particulator.PairwiseStorage.from_ndarray(np.array(Eb)) breakup_rate = particulator.Storage.from_ndarray(np.array([0])) breakup_rate_deficit = particulator.Storage.from_ndarray(np.array([0])) - frag_size = particulator.PairwiseStorage.from_ndarray( - np.asarray(params["frag_size"], dtype=float) + frag_mass = particulator.PairwiseStorage.from_ndarray( + np.asarray(volume_to_mass(particulator, params["frag_volume"]), dtype=float) ) is_first_in_pair = particulator.PairIndicator(n_sd) is_first_in_pair.indicator[:] = particulator.Storage.from_ndarray( @@ -711,7 +747,7 @@ def test_noninteger_fragments( rand=rand, Ec=pairwise_zeros, Eb=Eb, - fragment_size=frag_size, + fragment_mass=frag_mass, coalescence_rate=general_zeros, breakup_rate=breakup_rate, breakup_rate_deficit=breakup_rate_deficit, @@ -759,14 +795,14 @@ def test_nonnegative_even_if_overflow( norm_factor = 100 / si.cm**3 * si.m**3 X0 = Trivia.volume(const, radius=30.531 * si.micrometres) - spectrum = Exponential(norm_factor=norm_factor, scale=X0) + spectrum = spectra.Exponential(norm_factor=norm_factor, scale=X0) attributes = {} attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity( spectrum ).sample(n_sd) mu = Trivia.volume(const, radius=100 * si.um) - fragmentation = ExponFrag(scale=mu) + fragmentation = breakup_fragmentations.Exponential(scale=mu) kernel = Geometric() coal_eff = ConstEc(Ec=0.01) break_eff = ConstEb(Eb=1.0) @@ -806,7 +842,7 @@ def test_nonnegative_even_if_overflow( "v_expected": [0.5, 0.5], "expected_deficit": [0.0], "is_first_in_pair": [True, False], - "frag_size": [0.5], + "frag_volume": [0.5], }, { "gamma": [3.0], @@ -816,7 +852,7 @@ def test_nonnegative_even_if_overflow( "v_expected": [1, 1], "expected_deficit": [0.0], "is_first_in_pair": [True, False], - "frag_size": [1.0], + "frag_volume": [1.0], }, ], ) @@ -850,8 +886,10 @@ def test_while_loop_multi_breakup( Eb = particulator.PairwiseStorage.from_ndarray(np.array(Eb)) breakup_rate = particulator.Storage.from_ndarray(np.array([0])) breakup_rate_deficit = particulator.Storage.from_ndarray(np.array([0])) - frag_size = particulator.PairwiseStorage.from_ndarray( - np.array(params["frag_size"], dtype=float) + + frag_mass = volume_to_mass(particulator, params["frag_volume"]) + frag_mass = particulator.PairwiseStorage.from_ndarray( + np.array(frag_mass, dtype=float) ) is_first_in_pair = particulator.PairIndicator(n_sd) is_first_in_pair.indicator[:] = particulator.Storage.from_ndarray( @@ -865,7 +903,7 @@ def test_while_loop_multi_breakup( rand=rand, Ec=pairwise_zeros, Eb=Eb, - fragment_size=frag_size, + fragment_mass=frag_mass, coalescence_rate=general_zeros, breakup_rate=breakup_rate, breakup_rate_deficit=breakup_rate_deficit, diff --git a/tests/unit_tests/dynamics/collisions/test_sdm_single_cell.py b/tests/unit_tests/dynamics/collisions/test_sdm_single_cell.py index 3bf824e43..5e85423c4 100644 --- a/tests/unit_tests/dynamics/collisions/test_sdm_single_cell.py +++ b/tests/unit_tests/dynamics/collisions/test_sdm_single_cell.py @@ -7,6 +7,7 @@ from PySDM.backends.impl_common.indexed_storage import make_IndexedStorage from PySDM.backends.impl_common.pair_indicator import make_PairIndicator from PySDM.dynamics import Coalescence +from PySDM.physics import constants_defaults from .conftest import backend_fill, get_dummy_particulator_and_coalescence @@ -51,15 +52,19 @@ def test_single_collision(backend_class, v_2, T_2, n_2): np.round(particles["temperature"].to_ndarray().astype(float), 7), ) - assert np.sum( - particles["multiplicity"].to_ndarray() * particles["volume"].to_ndarray() - ) == np.sum(n_2 * v_2) + np.testing.assert_approx_equal( + np.sum( + particles["multiplicity"].to_ndarray() + * particles["volume"].to_ndarray() + ), + np.sum(n_2 * v_2), + ) assert np.sum(particulator.attributes["multiplicity"].to_ndarray()) == np.sum( n_2 ) - np.amin(n_2) if np.amin(n_2) > 0: - assert np.amax(particulator.attributes["volume"].to_ndarray()) == np.sum( - v_2 + np.testing.assert_approx_equal( + np.amax(particulator.attributes["volume"].to_ndarray()), np.sum(v_2) ) assert np.amax(particulator.attributes["multiplicity"].to_ndarray()) == max( np.amax(n_2) - np.amin(n_2), np.amin(n_2) @@ -124,15 +129,22 @@ def _compute_gamma(prob, rand, is_first_in_pair, out): state = particulator.attributes gamma = min(p, max(n_2[0] // n_2[1], n_2[1] // n_2[1])) assert np.amin(state["multiplicity"]) >= 0 - assert np.sum( - state["multiplicity"].to_ndarray() * state["volume"].to_ndarray() - ) == np.sum(n_2 * v_2) + np.testing.assert_approx_equal( + np.sum( + state["multiplicity"].to_ndarray() * state["water mass"].to_ndarray() + ), + np.sum(n_2 * v_2 * constants_defaults.rho_w), + ) + np.testing.assert_approx_equal( + np.sum(state["multiplicity"].to_ndarray() * state["volume"].to_ndarray()), + np.sum(n_2 * v_2), + ) assert np.sum(state["multiplicity"].to_ndarray()) == np.sum( n_2 ) - gamma * np.amin(n_2) - assert ( - np.amax(state["volume"].to_ndarray()) - == gamma * v_2[np.argmax(n_2)] + v_2[np.argmax(n_2) - 1] + np.testing.assert_approx_equal( + np.amax(state["volume"].to_ndarray()), + gamma * v_2[np.argmax(n_2)] + v_2[np.argmax(n_2) - 1], ) assert np.amax(state["multiplicity"].to_ndarray()) == max( np.amax(n_2) - gamma * np.amin(n_2), np.amin(n_2) diff --git a/tests/unit_tests/impl/test_particle_attributes.py b/tests/unit_tests/impl/test_particle_attributes.py index a87649b82..7c020a620 100644 --- a/tests/unit_tests/impl/test_particle_attributes.py +++ b/tests/unit_tests/impl/test_particle_attributes.py @@ -24,17 +24,17 @@ def make_indexed_storage(backend, iterable, idx=None): class TestParticleAttributes: @staticmethod @pytest.mark.parametrize( - "volume, multiplicity", + "water_mass, multiplicity", [ pytest.param(np.array([1.0, 1, 1, 1]), np.array([1, 1, 1, 1])), pytest.param(np.array([1.0, 2, 1, 1]), np.array([2, 0, 2, 0])), pytest.param(np.array([1.0, 1, 4]), np.array([5, 0, 0])), ], ) - def test_housekeeping(backend_class, volume, multiplicity): + def test_housekeeping(backend_class, water_mass, multiplicity): # Arrange particulator = DummyParticulator(backend_class, n_sd=len(multiplicity)) - attributes = {"multiplicity": multiplicity, "volume": volume} + attributes = {"multiplicity": multiplicity, "water mass": water_mass} particulator.build(attributes, int_caster=np.int64) sut = particulator.attributes sut.healthy = False @@ -47,8 +47,8 @@ def test_housekeeping(backend_class, volume, multiplicity): assert sut.super_droplet_count == (multiplicity != 0).sum() assert sut["multiplicity"].to_ndarray().sum() == multiplicity.sum() assert ( - sut["volume"].to_ndarray() * sut["multiplicity"].to_ndarray() - ).sum() == (volume * multiplicity).sum() + sut["water mass"].to_ndarray() * sut["multiplicity"].to_ndarray() + ).sum() == (water_mass * multiplicity).sum() @staticmethod @pytest.mark.parametrize( diff --git a/tests/unit_tests/physics/test_particle_shape_and_density.py b/tests/unit_tests/physics/test_particle_shape_and_density.py new file mode 100644 index 000000000..eac9ab6a3 --- /dev/null +++ b/tests/unit_tests/physics/test_particle_shape_and_density.py @@ -0,0 +1,40 @@ +# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring +import pytest + +from PySDM.formulae import Formulae +from PySDM.physics import constants_defaults +from PySDM.physics.dimensional_analysis import DimensionalAnalysis + + +class TestParticleShapeAndDensity: + @staticmethod + @pytest.mark.parametrize("variant", ("LiquidSpheres", "MixedPhaseSpheres")) + def test_mass_to_volume_units(variant): + with DimensionalAnalysis(): + # Arrange + formulae = Formulae(particle_shape_and_density=variant) + si = constants_defaults.si + sut = formulae.particle_shape_and_density.mass_to_volume + mass = 1 * si.gram + + # Act + volume = sut(mass) + + # Assert + assert volume.check("[volume]") + + @staticmethod + @pytest.mark.parametrize("variant", ("LiquidSpheres", "MixedPhaseSpheres")) + def test_volume_to_mass_units(variant): + with DimensionalAnalysis(): + # Arrange + formulae = Formulae(particle_shape_and_density=variant) + si = constants_defaults.si + sut = formulae.particle_shape_and_density.volume_to_mass + volume = 1 * si.micrometre**3 + + # Act + mass = sut(volume) + + # Assert + assert mass.check("[mass]")