diff --git a/.github/workflows/julia.yml b/.github/workflows/julia.yml index 88a0c8c0a..c0315dc2a 100644 --- a/.github/workflows/julia.yml +++ b/.github/workflows/julia.yml @@ -21,8 +21,8 @@ jobs: with: python-version: 3.9 - run: pip install -e . - - run: pip install exdown pytest - - run: python -c "import exdown; code=exdown.extract('README.md'); f=open('readme.jl', 'w'); f.writelines(block.code for block in code if block.syntax=='Julia'); f.close()" + - run: pip install pytest-codeblocks pytest + - run: python -c "import pytest_codeblocks; code=pytest_codeblocks.extract_from_file('README.md'); f=open('readme.jl', 'w'); f.writelines(block.code for block in code if block.syntax=='Julia'); f.close()" - uses: julia-actions/setup-julia@v1 - run: cat -n readme.jl - run: julia readme.jl diff --git a/.github/workflows/matlab.yml b/.github/workflows/matlab.yml index fa6aba883..ac68da5a8 100644 --- a/.github/workflows/matlab.yml +++ b/.github/workflows/matlab.yml @@ -21,8 +21,8 @@ jobs: with: python-version: 3.8 - run: pip install -e . - - run: pip install exdown pytest - - run: python -c "import exdown; code=exdown.extract('README.md'); f=open('readme.m', 'w'); f.writelines(block.code for block in code if block.syntax=='Matlab'); f.close()" + - run: pip install pytest-codeblocks pytest + - run: python -c "import pytest_codeblocks; code=pytest_codeblocks.extract_from_file('README.md'); f=open('readme.m', 'w'); f.writelines(block.code for block in code if block.syntax=='Matlab'); f.close()" - run: cat -n readme.m - uses: matlab-actions/setup-matlab@v0 with: diff --git a/.github/workflows/python.yml b/.github/workflows/python.yml index 38b1d969b..55ed28b41 100644 --- a/.github/workflows/python.yml +++ b/.github/workflows/python.yml @@ -21,8 +21,8 @@ jobs: with: python-version: 3.9 - run: pip install -e . - - run: pip install exdown pytest - - run: python -c "import exdown; code=exdown.extract('README.md'); f=open('readme.py', 'w'); f.writelines(block.code for block in code if block.syntax=='Python'); f.close()" + - run: pip install pytest-codeblocks pytest + - run: python -c "import pytest_codeblocks; code=pytest_codeblocks.extract_from_file('README.md'); f=open('readme.py', 'w'); f.writelines(block.code for block in code if block.syntax=='Python'); f.close()" - run: python -We readme.py - run: sed -i 's/CPU/GPU/g' readme.py - run: python -We readme.py diff --git a/PySDM/attributes/physics/volume.py b/PySDM/attributes/physics/volume.py index 00b95ae9b..ae8192278 100644 --- a/PySDM/attributes/physics/volume.py +++ b/PySDM/attributes/physics/volume.py @@ -1,5 +1,7 @@ """ particle (wet) volume, key attribute for coalescence +in simmulation involving mixed-phase clouds, positive values correspond to +liquid water and nagative values to ice """ from PySDM.attributes.impl.extensive_attribute import ExtensiveAttribute diff --git a/PySDM/backends/impl_numba/methods/collisions_methods.py b/PySDM/backends/impl_numba/methods/collisions_methods.py index fb522ca32..b42e6b43d 100644 --- a/PySDM/backends/impl_numba/methods/collisions_methods.py +++ b/PySDM/backends/impl_numba/methods/collisions_methods.py @@ -16,7 +16,7 @@ def pair_indices(i, idx, is_first_in_pair): return j, k -class AlgorithmicMethods(BackendMethods): +class CollisionsMethods(BackendMethods): @staticmethod @numba.njit(**{**conf.JIT_FLAGS, **{'parallel': False}}) def adaptive_sdm_end_body(dt_left, n_cell, cell_start): @@ -141,51 +141,6 @@ def compute_gamma(self, gamma, rand, multiplicity, cell_id, cell_id.data, collision_rate_deficit.data, collision_rate.data, is_first_in_pair.indicator.data) - @staticmethod - @numba.njit(**conf.JIT_FLAGS) - # pylint: disable=too-many-arguments,too-many-locals - def linear_collection_efficiency_body( - params, output, radii, is_first_in_pair, idx, length, unit - ): - A, B, D1, D2, E1, E2, F1, F2, G1, G2, G3, Mf, Mg = params - output[:] = 0 - for i in numba.prange(length - 1): # pylint: disable=not-an-iterable - if is_first_in_pair[i]: - if radii[idx[i]] > radii[idx[i + 1]]: - r = radii[idx[i]] / unit - r_s = radii[idx[i + 1]] / unit - else: - r = radii[idx[i + 1]] / unit - r_s = radii[idx[i]] / unit - p = r_s / r - if p not in (0, 1): - G = (G1 / r) ** Mg + G2 + G3 * r - Gp = (1 - p) ** G - if Gp != 0: - D = D1 / r ** D2 - E = E1 / r ** E2 - F = (F1 / r) ** Mf + F2 - output[i // 2] = A + B * p + D / p ** F + E / Gp - output[i // 2] = max(0, output[i // 2]) - - def linear_collection_efficiency(self, params, output, radii, is_first_in_pair, unit): - return self.linear_collection_efficiency_body( - params, output.data, radii.data, is_first_in_pair.indicator.data, - radii.idx.data, len(is_first_in_pair), unit) - - @staticmethod - @numba.njit(**conf.JIT_FLAGS) - def interpolation_body(output, radius, factor, b, c): - for i in numba.prange(len(radius)): # pylint: disable=not-an-iterable - r_id = int(factor * radius[i]) - r_rest = ((factor * radius[i]) % 1) / factor - output[i] = b[r_id] + r_rest * c[r_id] - - def interpolation(self, output, radius, factor, b, c): - return self.interpolation_body( - output.data, radius.data, factor, b.data, c.data - ) - @staticmethod def make_cell_caretaker(idx, cell_start, scheme="default"): class CellCaretaker: # pylint: disable=too-few-public-methods @@ -207,11 +162,11 @@ def __init__(self, idx, cell_start, scheme): def __call__(self, cell_id, cell_idx, cell_start, idx): length = len(idx) if self.scheme == "counting_sort": - AlgorithmicMethods._counting_sort_by_cell_id_and_update_cell_start( + CollisionsMethods._counting_sort_by_cell_id_and_update_cell_start( self.tmp_idx.data, idx.data, cell_id.data, cell_idx.data, length, cell_start.data) elif self.scheme == "counting_sort_parallel": - AlgorithmicMethods._parallel_counting_sort_by_cell_id_and_update_cell_start( + CollisionsMethods._parallel_counting_sort_by_cell_id_and_update_cell_start( self.tmp_idx.data, idx.data, cell_id.data, cell_idx.data, length, cell_start.data, self.cell_starts.data) idx.data, self.tmp_idx.data = self.tmp_idx.data, idx.data diff --git a/PySDM/backends/impl_numba/methods/freezing_methods.py b/PySDM/backends/impl_numba/methods/freezing_methods.py index 7def58270..19a718a18 100644 --- a/PySDM/backends/impl_numba/methods/freezing_methods.py +++ b/PySDM/backends/impl_numba/methods/freezing_methods.py @@ -25,10 +25,11 @@ def _freeze(volume, i): @numba.njit(**{**conf.JIT_FLAGS, 'fastmath': self.formulae.fastmath}) def freeze_singular_body(attributes, temperature, relative_humidity, cell): - for i in numba.prange(len(attributes.freezing_temperature)): # pylint: disable=not-an-iterable + n_sd = len(attributes.freezing_temperature) + for i in numba.prange(n_sd): # pylint: disable=not-an-iterable if ( _unfrozen(attributes.wet_volume, i) and - relative_humidity[cell[i]] > 1 and + relative_humidity[cell[i]] > 1 and # TODO #599 as in Shima, but is it needed? temperature[cell[i]] <= attributes.freezing_temperature[i] ): _freeze(attributes.wet_volume, i) diff --git a/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py b/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py new file mode 100644 index 000000000..aee4fb463 --- /dev/null +++ b/PySDM/backends/impl_numba/methods/terminal_velocity_methods.py @@ -0,0 +1,56 @@ +""" +CPU implementation of backend methods for terminal velocities +""" +import numba +from PySDM.backends.impl_numba import conf +from PySDM.backends.impl_common.backend_methods import BackendMethods + + +class TerminalVelocityMethods(BackendMethods): + @staticmethod + @numba.njit(**conf.JIT_FLAGS) + # pylint: disable=too-many-arguments,too-many-locals + def linear_collection_efficiency_body( + params, output, radii, is_first_in_pair, idx, length, unit + ): + A, B, D1, D2, E1, E2, F1, F2, G1, G2, G3, Mf, Mg = params + output[:] = 0 + for i in numba.prange(length - 1): # pylint: disable=not-an-iterable + if is_first_in_pair[i]: + if radii[idx[i]] > radii[idx[i + 1]]: + r = radii[idx[i]] / unit + r_s = radii[idx[i + 1]] / unit + else: + r = radii[idx[i + 1]] / unit + r_s = radii[idx[i]] / unit + p = r_s / r + if p not in (0, 1): + G = (G1 / r) ** Mg + G2 + G3 * r + Gp = (1 - p) ** G + if Gp != 0: + D = D1 / r ** D2 + E = E1 / r ** E2 + F = (F1 / r) ** Mf + F2 + output[i // 2] = A + B * p + D / p ** F + E / Gp + output[i // 2] = max(0, output[i // 2]) + + def linear_collection_efficiency(self, params, output, radii, is_first_in_pair, unit): + return self.linear_collection_efficiency_body( + params, output.data, radii.data, is_first_in_pair.indicator.data, + radii.idx.data, len(is_first_in_pair), unit) + + @staticmethod + @numba.njit(**conf.JIT_FLAGS) + def interpolation_body(output, radius, factor, b, c): + for i in numba.prange(len(radius)): # pylint: disable=not-an-iterable + if radius[i] < 0: + output[i] = 0 + else: + r_id = int(factor * radius[i]) + r_rest = ((factor * radius[i]) % 1) / factor + output[i] = b[r_id] + r_rest * c[r_id] + + def interpolation(self, output, radius, factor, b, c): + return self.interpolation_body( + output.data, radius.data, factor, b.data, c.data + ) diff --git a/PySDM/backends/impl_numba/storage_impl.py b/PySDM/backends/impl_numba/storage_impl.py index be6c439e4..77a7f494f 100644 --- a/PySDM/backends/impl_numba/storage_impl.py +++ b/PySDM/backends/impl_numba/storage_impl.py @@ -55,7 +55,8 @@ def sum_out_of_place(output, a, b): @numba.njit(**conf.JIT_FLAGS) def power(output, exponent): - output[:] = np.power(output, exponent) + # TODO #599 (was: output[:] = np.power(output, exponent)) + output[:] = np.sign(output) * np.power(np.abs(output), exponent) @numba.njit(**{**conf.JIT_FLAGS, **{'parallel': False}}) diff --git a/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py b/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py index 1b45149b0..4e4bcecfa 100644 --- a/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py +++ b/PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py @@ -10,40 +10,6 @@ class CollisionsMethods(ThrustRTCBackendMethods): def __init__(self): super().__init__() - self.__linear_collection_efficiency_body = trtc.For( - ( - 'A', 'B', 'D1', 'D2', 'E1', 'E2', 'F1', 'F2', 'G1', 'G2', 'G3', 'Mf', 'Mg', - 'output', 'radii', 'is_first_in_pair', 'idx', 'unit' - ), - "i", - ''' - if (is_first_in_pair[i]) { - real_type r = 0; - real_type r_s = 0; - if (radii[idx[i]] > radii[idx[i + 1]]) { - r = radii[idx[i]] / unit; - r_s = radii[idx[i + 1]] / unit; - } - else { - r = radii[idx[i + 1]] / unit; - r_s = radii[idx[i]] / unit; - } - real_type p = r_s / r; - if (p != 0 && p != 1) { - real_type G = pow((G1 / r), Mg) + G2 + G3 * r; - real_type Gp = pow((1 - p), G); - if (Gp != 0) { - real_type D = D1 / pow(r, D2); - real_type E = E1 / pow(r, E2); - real_type F = pow((F1 / r), Mf) + F2; - output[int(i / 2)] = A + B * p + D / pow(p, F) + E / Gp; - if (output[int(i / 2)] < 0) { - output[int(i / 2)] = 0; - } - } - } - } - '''.replace("real_type", self._get_c_type())) self.__adaptive_sdm_gamma_body_1 = trtc.For( ('dt_todo', 'dt_left', 'dt_max'), @@ -164,16 +130,6 @@ def __init__(self): ''' ) - self.__interpolation_body = trtc.For( - ('output', 'radius', 'factor', 'a', 'b'), - 'i', - ''' - auto r_id = (int64_t)(factor * radius[i]); - auto r_rest = (factor * radius[i] - r_id) / factor; - output[i] = a[r_id] + r_rest * b[r_id]; - ''' - ) - self.__normalize_body_0 = trtc.For( ('cell_start', 'norm_factor', 'dt_div_dv'), "i", @@ -287,36 +243,6 @@ def compute_gamma(self, gamma, rand, multiplicity, cell_id, (gamma.data, rand.data, multiplicity.idx.data, multiplicity.data, cell_id.data, collision_rate_deficit.data, collision_rate.data)) - def linear_collection_efficiency(self, params, output, radii, is_first_in_pair, unit): - A, B, D1, D2, E1, E2, F1, F2, G1, G2, G3, Mf, Mg = params - dA = self._get_floating_point(A) - dB = self._get_floating_point(B) - dD1 = self._get_floating_point(D1) - dD2 = self._get_floating_point(D2) - dE1 = self._get_floating_point(E1) - dE2 = self._get_floating_point(E2) - dF1 = self._get_floating_point(F1) - dF2 = self._get_floating_point(F2) - dG1 = self._get_floating_point(G1) - dG2 = self._get_floating_point(G2) - dG3 = self._get_floating_point(G3) - dMf = self._get_floating_point(Mf) - dMg = self._get_floating_point(Mg) - dunit = self._get_floating_point(unit) - trtc.Fill(output.data, trtc.DVDouble(0)) - self.__linear_collection_efficiency_body.launch_n( - len(is_first_in_pair) - 1, - [dA, dB, dD1, dD2, dE1, dE2, dF1, dF2, dG1, dG2, dG3, dMf, dMg, - output.data, radii.data, is_first_in_pair.indicator.data, radii.idx.data, dunit]) - - @nice_thrust(**NICE_THRUST_FLAGS) - def interpolation(self, output, radius, factor, b, c): - factor_device = trtc.DVInt64(factor) - self.__interpolation_body.launch_n( - len(radius), - (output.data, radius.data, factor_device, b.data, c.data) - ) - # pylint: disable=unused-argument def make_cell_caretaker(self, idx, cell_start, scheme=None): return self._sort_by_cell_id_and_update_cell_start diff --git a/PySDM/backends/impl_thrust_rtc/methods/terminal_velocity_methods.py b/PySDM/backends/impl_thrust_rtc/methods/terminal_velocity_methods.py new file mode 100644 index 000000000..b23520618 --- /dev/null +++ b/PySDM/backends/impl_thrust_rtc/methods/terminal_velocity_methods.py @@ -0,0 +1,87 @@ +""" +GPU implementation of backend methods for terminal velocities +""" +from PySDM.backends.impl_thrust_rtc.conf import NICE_THRUST_FLAGS +from PySDM.backends.impl_thrust_rtc.nice_thrust import nice_thrust +from ..conf import trtc +from ..methods.thrust_rtc_backend_methods import ThrustRTCBackendMethods + + +class TerminalVelocityMethods(ThrustRTCBackendMethods): + def __init__(self): + super().__init__() + self.__linear_collection_efficiency_body = trtc.For( + ( + 'A', 'B', 'D1', 'D2', 'E1', 'E2', 'F1', 'F2', 'G1', 'G2', 'G3', 'Mf', 'Mg', + 'output', 'radii', 'is_first_in_pair', 'idx', 'unit' + ), + "i", + ''' + if (is_first_in_pair[i]) { + real_type r = 0; + real_type r_s = 0; + if (radii[idx[i]] > radii[idx[i + 1]]) { + r = radii[idx[i]] / unit; + r_s = radii[idx[i + 1]] / unit; + } + else { + r = radii[idx[i + 1]] / unit; + r_s = radii[idx[i]] / unit; + } + real_type p = r_s / r; + if (p != 0 && p != 1) { + real_type G = pow((G1 / r), Mg) + G2 + G3 * r; + real_type Gp = pow((1 - p), G); + if (Gp != 0) { + real_type D = D1 / pow(r, D2); + real_type E = E1 / pow(r, E2); + real_type F = pow((F1 / r), Mf) + F2; + output[int(i / 2)] = A + B * p + D / pow(p, F) + E / Gp; + if (output[int(i / 2)] < 0) { + output[int(i / 2)] = 0; + } + } + } + } + '''.replace("real_type", self._get_c_type())) + + # TODO #599 r<0 + self.__interpolation_body = trtc.For( + ('output', 'radius', 'factor', 'a', 'b'), + 'i', + ''' + auto r_id = (int64_t)(factor * radius[i]); + auto r_rest = (factor * radius[i] - r_id) / factor; + output[i] = a[r_id] + r_rest * b[r_id]; + ''' + ) + + def linear_collection_efficiency(self, params, output, radii, is_first_in_pair, unit): + A, B, D1, D2, E1, E2, F1, F2, G1, G2, G3, Mf, Mg = params + dA = self._get_floating_point(A) + dB = self._get_floating_point(B) + dD1 = self._get_floating_point(D1) + dD2 = self._get_floating_point(D2) + dE1 = self._get_floating_point(E1) + dE2 = self._get_floating_point(E2) + dF1 = self._get_floating_point(F1) + dF2 = self._get_floating_point(F2) + dG1 = self._get_floating_point(G1) + dG2 = self._get_floating_point(G2) + dG3 = self._get_floating_point(G3) + dMf = self._get_floating_point(Mf) + dMg = self._get_floating_point(Mg) + dunit = self._get_floating_point(unit) + trtc.Fill(output.data, trtc.DVDouble(0)) + self.__linear_collection_efficiency_body.launch_n( + len(is_first_in_pair) - 1, + [dA, dB, dD1, dD2, dE1, dE2, dF1, dF2, dG1, dG2, dG3, dMf, dMg, + output.data, radii.data, is_first_in_pair.indicator.data, radii.idx.data, dunit]) + + @nice_thrust(**NICE_THRUST_FLAGS) + def interpolation(self, output, radius, factor, b, c): + factor_device = trtc.DVInt64(factor) + self.__interpolation_body.launch_n( + len(radius), + (output.data, radius.data, factor_device, b.data, c.data) + ) diff --git a/PySDM/backends/numba.py b/PySDM/backends/numba.py index 8015e4d37..b6ff67151 100644 --- a/PySDM/backends/numba.py +++ b/PySDM/backends/numba.py @@ -2,7 +2,7 @@ Multi-threaded CPU backend using LLVM-powered just-in-time compilation """ -from PySDM.backends.impl_numba.methods.collisions_methods import AlgorithmicMethods +from PySDM.backends.impl_numba.methods.collisions_methods import CollisionsMethods from PySDM.backends.impl_numba.methods.pair_methods import PairMethods from PySDM.backends.impl_numba.methods.physics_methods import PhysicsMethods from PySDM.backends.impl_numba.methods.index_methods import IndexMethods @@ -11,13 +11,14 @@ from PySDM.backends.impl_numba.methods.freezing_methods import FreezingMethods from PySDM.backends.impl_numba.methods.chemistry_methods import ChemistryMethods from PySDM.backends.impl_numba.methods.displacement_methods import DisplacementMethods +from PySDM.backends.impl_numba.methods.terminal_velocity_methods import TerminalVelocityMethods from PySDM.backends.impl_numba.random import Random as ImportedRandom from PySDM.backends.impl_numba.storage import Storage as ImportedStorage from PySDM.formulae import Formulae class Numba( # pylint: disable=too-many-ancestors,duplicate-code - AlgorithmicMethods, + CollisionsMethods, PairMethods, IndexMethods, PhysicsMethods, @@ -25,7 +26,8 @@ class Numba( # pylint: disable=too-many-ancestors,duplicate-code ChemistryMethods, MomentsMethods, FreezingMethods, - DisplacementMethods + DisplacementMethods, + TerminalVelocityMethods ): Storage = ImportedStorage Random = ImportedRandom @@ -34,7 +36,7 @@ class Numba( # pylint: disable=too-many-ancestors,duplicate-code def __init__(self, formulae=None): self.formulae = formulae or Formulae() - AlgorithmicMethods.__init__(self) + CollisionsMethods.__init__(self) PairMethods.__init__(self) IndexMethods.__init__(self) PhysicsMethods.__init__(self) @@ -43,3 +45,4 @@ def __init__(self, formulae=None): MomentsMethods.__init__(self) FreezingMethods.__init__(self) DisplacementMethods.__init__(self) + TerminalVelocityMethods.__init__(self) diff --git a/PySDM/backends/thrust_rtc.py b/PySDM/backends/thrust_rtc.py index d2c0b8fb1..9883d9dd7 100644 --- a/PySDM/backends/thrust_rtc.py +++ b/PySDM/backends/thrust_rtc.py @@ -12,6 +12,8 @@ from PySDM.backends.impl_thrust_rtc.methods.moments_methods import MomentsMethods from PySDM.backends.impl_thrust_rtc.methods.condensation_methods import CondensationMethods from PySDM.backends.impl_thrust_rtc.methods.displacement_methods import DisplacementMethods +from PySDM.backends.impl_thrust_rtc.methods.terminal_velocity_methods import \ + TerminalVelocityMethods from PySDM.backends.impl_thrust_rtc.storage import make_storage_class from PySDM.backends.impl_thrust_rtc.random import Random as ImportedRandom from PySDM.formulae import Formulae @@ -25,7 +27,8 @@ class ThrustRTC( # pylint: disable=duplicate-code,too-many-ancestors PhysicsMethods, CondensationMethods, MomentsMethods, - DisplacementMethods + DisplacementMethods, + TerminalVelocityMethods ): ENABLE = True Random = ImportedRandom @@ -48,6 +51,7 @@ def __init__(self, formulae=None, double_precision=False, debug=False, verbose=F CondensationMethods.__init__(self) MomentsMethods.__init__(self) DisplacementMethods.__init__(self) + TerminalVelocityMethods.__init__(self) trtc.Set_Kernel_Debug(debug) trtc.Set_Verbose(verbose) diff --git a/PySDM/environments/impl/moist.py b/PySDM/environments/impl/moist.py index 00f0df0a4..58ea3a0dc 100644 --- a/PySDM/environments/impl/moist.py +++ b/PySDM/environments/impl/moist.py @@ -7,8 +7,10 @@ class Moist: - def __init__(self, dt, mesh, variables): + def __init__(self, dt, mesh, variables, mixed_phase=False): variables += ['qv', 'thd', 'T', 'p', 'RH'] + if mixed_phase: + variables += ['a_w_ice'] self.particulator = None self.dt = dt self.mesh = mesh diff --git a/PySDM/environments/kinematic_2d.py b/PySDM/environments/kinematic_2d.py index 837a679b6..2af7a67ac 100644 --- a/PySDM/environments/kinematic_2d.py +++ b/PySDM/environments/kinematic_2d.py @@ -7,13 +7,14 @@ from PySDM.impl.mesh import Mesh from PySDM.initialisation.equilibrate_wet_radii import equilibrate_wet_radii, default_rtol from PySDM.initialisation.discretise_multiplicities import discretise_multiplicities +from PySDM.initialisation.sampling import spectral_sampling from PySDM.environments.impl.moist import Moist from ..impl import arakawa_c class Kinematic2D(Moist): - def __init__(self, dt, grid, size, rhod_of): - super().__init__(dt, Mesh(grid, size), []) + def __init__(self, dt, grid, size, rhod_of, mixed_phase=False): + super().__init__(dt, Mesh(grid, size), [], mixed_phase=mixed_phase) self.rhod_of = rhod_of self.formulae = None @@ -32,28 +33,21 @@ def dv(self): def init_attributes(self, *, spatial_discretisation, kappa, - spectral_discretisation = None, - spectro_glacial_discretisation = None, + dry_radius_spectrum, rtol=default_rtol ): super().sync() self.notify() - assert spectro_glacial_discretisation is None or spectral_discretisation is None - attributes = {} with np.errstate(all='raise'): positions = spatial_discretisation.sample(self.mesh.grid, self.particulator.n_sd) attributes['cell id'], attributes['cell origin'], attributes['position in cell'] = \ self.mesh.cellular_attributes(positions) - if spectral_discretisation: - r_dry, n_per_kg = spectral_discretisation.sample(self.particulator.n_sd) - elif spectro_glacial_discretisation: - r_dry, T_fz, n_per_kg = spectro_glacial_discretisation.sample( - self.particulator.n_sd) - attributes['freezing temperature'] = T_fz - else: - raise NotImplementedError() + + r_dry, n_per_kg = spectral_sampling.ConstantMultiplicity( + spectrum=dry_radius_spectrum + ).sample(n_sd=self.particulator.n_sd) attributes['dry volume'] = self.formulae.trivia.volume(radius=r_dry) attributes['kappa times dry volume'] = kappa * attributes['dry volume'] diff --git a/PySDM/environments/parcel.py b/PySDM/environments/parcel.py index d31e8ec87..a1cb70772 100644 --- a/PySDM/environments/parcel.py +++ b/PySDM/environments/parcel.py @@ -23,7 +23,8 @@ def __init__( super().__init__( dt, Mesh.mesh_0d(), - ['rhod', 'z', 't'] + (['a_w_ice'] if mixed_phase else []) + ['rhod', 'z', 't'], + mixed_phase=mixed_phase ) self.p0 = p0 diff --git a/PySDM/exporters/netcdf_exporter.py b/PySDM/exporters/netcdf_exporter.py index 532d54b49..b3ed5012b 100644 --- a/PySDM/exporters/netcdf_exporter.py +++ b/PySDM/exporters/netcdf_exporter.py @@ -3,8 +3,11 @@ """ import numpy as np from scipy.io.netcdf import netcdf_file +from PySDM.products.impl.spectrum_moment_product import SpectrumMomentProduct +DIM_SUFFIX = '_bin_left_edges' + class NetCDFExporter: def __init__(self, storage, settings, simulator, filename): self.storage = storage @@ -20,12 +23,16 @@ def _write_settings(self, ncdf): def _create_dimensions(self, ncdf): ncdf.createDimension("T", len(self.settings.output_steps)) + for index, label in enumerate(self.XZ): ncdf.createDimension(label, self.settings.grid[index]) - ncdf.createDimension( - "ParticleVolume", - len(self.settings.formulae.trivia.volume(self.settings.r_bins_edges)) - 1 - ) + + for name, instance in self.simulator.products.items(): + if isinstance(instance, SpectrumMomentProduct): + ncdf.createDimension( + f"{name}{DIM_SUFFIX}", + len(instance.attr_bins_edges) - 1 + ) def _create_variables(self, ncdf): self.vars = {} @@ -38,7 +45,12 @@ def _create_variables(self, ncdf): * (1 / 2 + np.arange(self.settings.grid[index]))) self.vars[label].units = "metres" - # TODO #340 ParticleVolume var + for name, instance in self.simulator.products.items(): + if isinstance(instance, SpectrumMomentProduct): + label = f"{name}{DIM_SUFFIX}" + self.vars[label] = ncdf.createVariable(label, 'f', (label,)) + self.vars[label][:] = instance.attr_bins_edges.to_ndarray()[:-1] + self.vars[label].units = instance.attr_unit for name, instance in self.simulator.products.items(): if name in self.vars: @@ -46,7 +58,7 @@ def _create_variables(self, ncdf): n_dimensions = len(instance.shape) if n_dimensions == 3: - dimensions = ("T", "X", "Z", "ParticleVolume") + dimensions = ("T", "X", "Z", f"{name}{DIM_SUFFIX}") elif n_dimensions == 2: dimensions = ("T", "X", "Z") elif n_dimensions == 0: diff --git a/PySDM/initialisation/__init__.py b/PySDM/initialisation/__init__.py index ec15d47e5..2b2ad9abc 100644 --- a/PySDM/initialisation/__init__.py +++ b/PySDM/initialisation/__init__.py @@ -1,5 +1,5 @@ """ -Initialisation logic particle size spectra, sampling methods and +initialisation logic, particle size spectra, sampling methods and wet radii equilibration """ from .discretise_multiplicities import discretise_multiplicities diff --git a/PySDM/physics/terminal_velocity/gunn_and_kinzer.py b/PySDM/physics/terminal_velocity/gunn_and_kinzer.py index 09a2d17e8..21617067b 100644 --- a/PySDM/physics/terminal_velocity/gunn_and_kinzer.py +++ b/PySDM/physics/terminal_velocity/gunn_and_kinzer.py @@ -1,6 +1,7 @@ """ [Gunn & Kinzer 1949](https://doi.org/10.1175/1520-0469(1949)006%3C0243:TTVOFF%3E2.0.CO;2) - terminal velocities + terminal velocities used for both coalescence kernel evaluation as well as for particle + displacement """ import numba import numpy as np @@ -43,29 +44,6 @@ def __call__(self, output, radius): self.particulator.backend.interpolation(output, radius, self.factor, self.a, self.b) -class RogersYau: - """ - Rogers & Yau, equations: (8.5), (8.6), (8.8) - """ - def __init__(self, particles, - small_k=None, medium_k=None, large_k=None, - small_r_limit=None, medium_r_limit=None): - si = const.si - self.particles = particles - self.small_k = small_k or 1.19e6 / si.cm / si.s - self.medium_k = medium_k or 8e3 / si.s - self.large_k = large_k or 2.01e3 * si.cm ** (1 / 2) / si.s - self.small_r_limit = small_r_limit or 35 * si.um - self.medium_r_limit = medium_r_limit or 600 * si.um - - def __call__(self, output, radius): - self.particles.backend.terminal_velocity( - output.data, radius.data, - self.small_k, self.medium_k, self.large_k, - self.small_r_limit, self.medium_r_limit - ) - - # TODO #348 implement in backend logic class TpDependent: def __init__(self, _, small_r_limit): diff --git a/PySDM/physics/terminal_velocity/rogers_and_yau.py b/PySDM/physics/terminal_velocity/rogers_and_yau.py new file mode 100644 index 000000000..90909fb74 --- /dev/null +++ b/PySDM/physics/terminal_velocity/rogers_and_yau.py @@ -0,0 +1,24 @@ +""" +Rogers & Yau, equations: (8.5), (8.6), (8.8) +""" +from PySDM.physics import constants as const + + +class RogersYau: + def __init__(self, particles, + small_k=None, medium_k=None, large_k=None, + small_r_limit=None, medium_r_limit=None): + si = const.si + self.particles = particles + self.small_k = small_k or 1.19e6 / si.cm / si.s + self.medium_k = medium_k or 8e3 / si.s + self.large_k = large_k or 2.01e3 * si.cm ** (1 / 2) / si.s + self.small_r_limit = small_r_limit or 35 * si.um + self.medium_r_limit = medium_r_limit or 600 * si.um + + def __call__(self, output, radius): + self.particles.backend.terminal_velocity( + output.data, radius.data, + self.small_k, self.medium_k, self.large_k, + self.small_r_limit, self.medium_r_limit + ) diff --git a/PySDM/products/aqueous_chemistry/aqueous_mass_spectrum.py b/PySDM/products/aqueous_chemistry/aqueous_mass_spectrum.py index 84935aeb7..b91086cf8 100644 --- a/PySDM/products/aqueous_chemistry/aqueous_mass_spectrum.py +++ b/PySDM/products/aqueous_chemistry/aqueous_mass_spectrum.py @@ -11,7 +11,7 @@ class AqueousMassSpectrum(SpectrumMomentProduct): def __init__(self, key, dry_radius_bins_edges, specific=False, name=None, unit='kg/m^3'): - super().__init__(name=name, unit=unit) + super().__init__(name=name, unit=unit, attr_unit='m') self.key = key self.dry_radius_bins_edges = dry_radius_bins_edges self.molar_mass = Substance.from_formula(AQUEOUS_COMPOUNDS[key][0]).mass * si.g / si.mole diff --git a/PySDM/products/freezing/freezable_specific_concentration.py b/PySDM/products/freezing/freezable_specific_concentration.py index 98a56013a..3e7f57ea4 100644 --- a/PySDM/products/freezing/freezable_specific_concentration.py +++ b/PySDM/products/freezing/freezable_specific_concentration.py @@ -7,7 +7,7 @@ class FreezableSpecificConcentration(SpectrumMomentProduct): def __init__(self, temperature_bins_edges, name=None, unit='kg^-1 K^-1'): - super().__init__(name=name, unit=unit) + super().__init__(name=name, unit=unit, attr_unit='K') self.attr_bins_edges = temperature_bins_edges def register(self, builder): diff --git a/PySDM/products/freezing/ice_water_content.py b/PySDM/products/freezing/ice_water_content.py index e5528454c..b9266e94a 100644 --- a/PySDM/products/freezing/ice_water_content.py +++ b/PySDM/products/freezing/ice_water_content.py @@ -12,11 +12,12 @@ def __init__(self, unit='kg/kg', name=None, __specific=False): self.specific = __specific def _impl(self, **kwargs): - self._download_moment_to_buffer('volume', rank=0, filter_range=(-np.inf, 0)) - conc = self.buffer.copy() - self._download_moment_to_buffer('volume', rank=1, filter_range=(-np.inf, 0)) result = self.buffer.copy() + + self._download_moment_to_buffer('volume', rank=0, filter_range=(-np.inf, 0)) + conc = self.buffer + result[:] *= -const.rho_i * conc / self.particulator.mesh.dv if self.specific: diff --git a/PySDM/products/freezing/total_unfrozen_immersed_surface_area.py b/PySDM/products/freezing/total_unfrozen_immersed_surface_area.py index 008372f7f..e9b55e1df 100644 --- a/PySDM/products/freezing/total_unfrozen_immersed_surface_area.py +++ b/PySDM/products/freezing/total_unfrozen_immersed_surface_area.py @@ -19,4 +19,5 @@ def _impl(self, **kwargs): result = np.copy(self.buffer) self._download_moment_to_buffer(**params, rank=0) result[:] *= self.buffer + # TODO #599 per volume / per gridbox ? return result diff --git a/PySDM/products/impl/product.py b/PySDM/products/impl/product.py index 385a09c5b..8975728b5 100644 --- a/PySDM/products/impl/product.py +++ b/PySDM/products/impl/product.py @@ -17,7 +17,7 @@ class Product: def __init__(self, *, unit: str, name: Optional[str] = None): self.name = name or self._camel_case_to_words(self.__class__.__name__) - self._unit = self.__parse_unit(unit) + self._unit = self._parse_unit(unit) self.unit_magnitude_in_base_units = self._unit.to_base_units().magnitude self.__check_unit() @@ -36,7 +36,7 @@ def _download_to_buffer(self, storage): storage.download(self.buffer.ravel()) @staticmethod - def __parse_unit(unit: str): + def _parse_unit(unit: str): if unit in ('%', 'percent'): return .01 * _UNIT_REGISTRY.dimensionless if unit in ('PPB', 'ppb'): @@ -65,7 +65,7 @@ def __check_unit(self): raise AssertionError(f"unit parameter of {type(self).__name__}.__init__" f" is expected to have a non-empty default value") - default_unit = self.__parse_unit(default_unit_arg) + default_unit = self._parse_unit(default_unit_arg) if default_unit.to_base_units().magnitude != 1: raise AssertionError(f'default value "{default_unit_arg}"' diff --git a/PySDM/products/impl/spectrum_moment_product.py b/PySDM/products/impl/spectrum_moment_product.py index 1296a32f2..201d49214 100644 --- a/PySDM/products/impl/spectrum_moment_product.py +++ b/PySDM/products/impl/spectrum_moment_product.py @@ -7,9 +7,11 @@ class SpectrumMomentProduct(ABC, Product): - def __init__(self, name, unit): + def __init__(self, name, unit, attr_unit): super().__init__(name=name, unit=unit) + _ = self._parse_unit(unit) self.attr_bins_edges = None + self.attr_unit = attr_unit self.moment_0 = None self.moments = None diff --git a/PySDM/products/size_spectral/__init__.py b/PySDM/products/size_spectral/__init__.py index 208b35f23..4b1a20110 100644 --- a/PySDM/products/size_spectral/__init__.py +++ b/PySDM/products/size_spectral/__init__.py @@ -3,7 +3,10 @@ from .mean_radius import MeanRadius from .particle_size_spectrum import ParticleSizeSpectrumPerMass, ParticleSizeSpectrumPerVolume from .particles_concentration import ParticleConcentration, ParticleSpecificConcentration -from .particles_volume_spectrum import ParticlesVolumeSpectrum +from .particle_volume_versus_radius_logarithm_spectrum import \ + ParticleVolumeVersusRadiusLogarithmSpectrum from .total_particle_concentration import TotalParticleConcentration from .total_particle_specific_concentration import TotalParticleSpecificConcentration from .water_mixing_ratio import WaterMixingRatio +from .radius_binned_number_averaged_terminal_velocity \ + import RadiusBinnedNumberAveragedTerminalVelocity diff --git a/PySDM/products/size_spectral/particle_size_spectrum.py b/PySDM/products/size_spectral/particle_size_spectrum.py index cf7685cd5..3d25f9e17 100644 --- a/PySDM/products/size_spectral/particle_size_spectrum.py +++ b/PySDM/products/size_spectral/particle_size_spectrum.py @@ -11,7 +11,7 @@ def __init__(self, radius_bins_edges, name, unit, dry=False, normalise_by_dv=Fal self.volume_attr = 'dry volume' if dry else 'volume' self.radius_bins_edges = radius_bins_edges self.normalise_by_dv = normalise_by_dv - super().__init__(name=name, unit=unit) + super().__init__(name=name, unit=unit, attr_unit='m') def register(self, builder): builder.request_attribute(self.volume_attr) diff --git a/PySDM/products/size_spectral/particles_volume_spectrum.py b/PySDM/products/size_spectral/particle_volume_versus_radius_logarithm_spectrum.py similarity index 80% rename from PySDM/products/size_spectral/particles_volume_spectrum.py rename to PySDM/products/size_spectral/particle_volume_versus_radius_logarithm_spectrum.py index 13c7a56c0..c9f681dc3 100644 --- a/PySDM/products/size_spectral/particles_volume_spectrum.py +++ b/PySDM/products/size_spectral/particle_volume_versus_radius_logarithm_spectrum.py @@ -1,13 +1,15 @@ """ -the dn-dlnr particle volume spectrum per volume of air (uses natural logarithm) +n_V^e(ln(r)) particle volume spectrum per volume of air (uses natural logarithm), +i.e. volume of particles per volume of air having in the size range ln(r) to +ln(r) + dln(r) """ import numpy as np from PySDM.products.impl.spectrum_moment_product import SpectrumMomentProduct -class ParticlesVolumeSpectrum(SpectrumMomentProduct): +class ParticleVolumeVersusRadiusLogarithmSpectrum(SpectrumMomentProduct): def __init__(self, radius_bins_edges, name=None, unit='dimensionless'): - super().__init__(name=name, unit=unit) + super().__init__(name=name, unit=unit, attr_unit='m') self.radius_bins_edges = radius_bins_edges self.moment_0 = None self.moments = None diff --git a/PySDM/products/size_spectral/radius_binned_number_averaged_terminal_velocity.py b/PySDM/products/size_spectral/radius_binned_number_averaged_terminal_velocity.py new file mode 100644 index 000000000..da47177f3 --- /dev/null +++ b/PySDM/products/size_spectral/radius_binned_number_averaged_terminal_velocity.py @@ -0,0 +1,38 @@ +""" +Provides radius bin-resolved average terminal velocity (average is particle-number weighted) +""" +import numpy as np +from PySDM.products.impl.spectrum_moment_product import SpectrumMomentProduct + +ATTR = 'terminal velocity' +RANK = 1 + +class RadiusBinnedNumberAveragedTerminalVelocity(SpectrumMomentProduct): + def __init__(self, radius_bin_edges, name=None, unit='m/s'): + super().__init__(name=name, unit=unit, attr_unit='m') + self.radius_bin_edges = radius_bin_edges + + def register(self, builder): + builder.request_attribute(ATTR) + + volume_bin_edges = builder.particulator.formulae.trivia.volume(self.radius_bin_edges) + self.attr_bins_edges = builder.particulator.backend.Storage.from_ndarray(volume_bin_edges) + + super().register(builder) + + self.shape = (*builder.particulator.mesh.grid, len(self.attr_bins_edges) - 1) + + + def _impl(self, **kwargs): + vals = np.empty([self.particulator.mesh.n_cell, len(self.attr_bins_edges) - 1]) + + self._recalculate_spectrum_moment( + attr=ATTR, + rank=RANK, + ) + + for i in range(vals.shape[1]): + self._download_spectrum_moment_to_buffer(rank=RANK, bin_number=i) + vals[:, i] = self.buffer.ravel() + + return np.squeeze(vals.reshape(self.shape)) diff --git a/README.md b/README.md index 3f2d930bb..1be872931 100644 --- a/README.md +++ b/README.md @@ -169,14 +169,14 @@ Box = pyimport("PySDM.environments").Box Coalescence = pyimport("PySDM.dynamics").Coalescence Golovin = pyimport("PySDM.physics.coalescence_kernels").Golovin CPU = pyimport("PySDM.backends").CPU -ParticlesVolumeSpectrum = pyimport("PySDM.products").ParticlesVolumeSpectrum +ParticleVolumeVersusRadiusLogarithmSpectrum = pyimport("PySDM.products").ParticleVolumeVersusRadiusLogarithmSpectrum radius_bins_edges = 10 .^ range(log10(10*si.um), log10(5e3*si.um), length=32) builder = Builder(n_sd=n_sd, backend=CPU()) builder.set_environment(Box(dt=1 * si.s, dv=1e6 * si.m^3)) builder.add_dynamic(Coalescence(kernel=Golovin(b=1.5e3 / si.s))) -products = [ParticlesVolumeSpectrum(radius_bins_edges=radius_bins_edges, name="dv/dlnr")] +products = [ParticleVolumeVersusRadiusLogarithmSpectrum(radius_bins_edges=radius_bins_edges, name="dv/dlnr")] particulator = builder.build(attributes, products) ``` @@ -189,14 +189,14 @@ Box = py.importlib.import_module('PySDM.environments').Box; Coalescence = py.importlib.import_module('PySDM.dynamics').Coalescence; Golovin = py.importlib.import_module('PySDM.physics.coalescence_kernels').Golovin; CPU = py.importlib.import_module('PySDM.backends').CPU; -ParticlesVolumeSpectrum = py.importlib.import_module('PySDM.products').ParticlesVolumeSpectrum; +ParticleVolumeVersusRadiusLogarithmSpectrum = py.importlib.import_module('PySDM.products').ParticleVolumeVersusRadiusLogarithmSpectrum; radius_bins_edges = logspace(log10(10 * si.um), log10(5e3 * si.um), 32); builder = Builder(pyargs('n_sd', int32(n_sd), 'backend', CPU())); builder.set_environment(Box(pyargs('dt', 1 * si.s, 'dv', 1e6 * si.m ^ 3))); builder.add_dynamic(Coalescence(pyargs('kernel', Golovin(1.5e3 / si.s)))); -products = py.list({ ParticlesVolumeSpectrum(pyargs( ... +products = py.list({ ParticleVolumeVersusRadiusLogarithmSpectrum(pyargs( ... 'radius_bins_edges', py.numpy.array(radius_bins_edges), ... 'name', 'dv/dlnr' ... )) }); @@ -213,14 +213,14 @@ from PySDM.environments import Box from PySDM.dynamics import Coalescence from PySDM.physics.coalescence_kernels import Golovin from PySDM.backends import CPU -from PySDM.products import ParticlesVolumeSpectrum +from PySDM.products import ParticleVolumeVersusRadiusLogarithmSpectrum radius_bins_edges = np.logspace(np.log10(10 * si.um), np.log10(5e3 * si.um), num=32) builder = Builder(n_sd=n_sd, backend=CPU()) -builder.set_environment(Box(dt=1 * si.s, dv=1e6 * si.m**3)) +builder.set_environment(Box(dt=1 * si.s, dv=1e6 * si.m ** 3)) builder.add_dynamic(Coalescence(kernel=Golovin(b=1.5e3 / si.s))) -products = [ParticlesVolumeSpectrum(radius_bins_edges=radius_bins_edges, name='dv/dlnr')] +products = [ParticleVolumeVersusRadiusLogarithmSpectrum(radius_bins_edges=radius_bins_edges, name='dv/dlnr')] particulator = builder.build(attributes, products) ``` diff --git a/test-time-requirements.txt b/test-time-requirements.txt index 8af10bf0b..cf49cf5c5 100644 --- a/test-time-requirements.txt +++ b/test-time-requirements.txt @@ -4,4 +4,4 @@ ghapi pytest # note: if cloning both PySDM and PySDM examples, consider "pip install -e" -PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@8267691#egg=PySDM-examples +PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@ed59c6e#egg=PySDM-examples diff --git a/tests/smoke_tests/arabas_et_al_2015/dummy_storage.py b/tests/smoke_tests/arabas_et_al_2015/dummy_storage.py new file mode 100644 index 000000000..df9665ffc --- /dev/null +++ b/tests/smoke_tests/arabas_et_al_2015/dummy_storage.py @@ -0,0 +1,14 @@ +# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring +import numpy as np + + +class DummyStorage: + def __init__(self): + self.profiles = [] + + def init(*_): # pylint: disable=no-method-argument + pass + + def save(self, data: np.ndarray, step: int, name: str): # pylint: disable=unused-argument + if name == "qv_env": + self.profiles.append({"qv_env": np.mean(data, axis=0)}) diff --git a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py new file mode 100644 index 000000000..60f40b9fc --- /dev/null +++ b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py @@ -0,0 +1,57 @@ +# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring +import pytest +from PySDM_examples.Szumowski_et_al_1998 import Simulation +from PySDM_examples.Arabas_et_al_2015 import Settings, SpinUp +from PySDM import Formulae +from PySDM.physics import si +from PySDM.backends import CPU +from PySDM.physics.freezing_temperature_spectrum import niemand_et_al_2012 +from PySDM.physics.heterogeneous_ice_nucleation_rate import abifm +from .dummy_storage import DummyStorage + +# TODO #599 +niemand_et_al_2012.a = -0.517 +niemand_et_al_2012.b = 8.934 +abifm.m = 28.13797 +abifm.c = -2.92414 + + +@pytest.mark.parametrize("singular", ( + pytest.param(False, id="singular: False"), + pytest.param(True, id="singular: True") +)) +# pylint: disable=redefined-outer-name +def test_freezing(singular): + # Arrange + settings = Settings(Formulae( + seed=44, + condensation_coordinate='VolumeLogarithm', + fastmath=True, + freezing_temperature_spectrum='Niemand_et_al_2012', + heterogeneous_ice_nucleation_rate='ABIFM' + )) + settings.dt = .5 * si.second + settings.grid = (5, 15) + settings.n_sd_per_gridbox = 64 + + settings.simulation_time = 100 * settings.dt + settings.spin_up_time = 10 * settings.dt + + settings.output_interval = settings.dt # settings.simulation_time + + settings.processes['freezing'] = True + settings.processes['coalescence'] = False + + settings.freezing_singular = singular + settings.th_std0 -= 35 * si.K + settings.qv0 -= 7.15 * si.g/si.kg + + storage = DummyStorage() + simulation = Simulation(settings, storage, SpinUp=SpinUp, backend_class=CPU) + simulation.reinit() + + # Act + simulation.run() + + # Assert + assert (simulation.products['ice water content'].get() > 0).any() diff --git a/tests/smoke_tests/arabas_et_al_2015/test_gui_settings.py b/tests/smoke_tests/arabas_et_al_2015/test_gui_settings.py new file mode 100644 index 000000000..a474a9497 --- /dev/null +++ b/tests/smoke_tests/arabas_et_al_2015/test_gui_settings.py @@ -0,0 +1,25 @@ +# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring +from PySDM_examples.Arabas_et_al_2015 import Settings +from PySDM_examples.Szumowski_et_al_1998.gui_settings import GUISettings + + +class TestGUISettings: + @staticmethod + def test_instantiate(): + _ = GUISettings(Settings()) + + @staticmethod + def test_stream_function(): + # arrange + gui_settings = GUISettings(Settings()) + gui_settings.ui_rhod_w_max = None + failed = False + + # act + try: + _ = gui_settings.stream_function(0, 0, 0) + except AttributeError: + failed = True + + # assert + assert failed diff --git a/tests/smoke_tests/arabas_et_al_2015/test_spin_up.py b/tests/smoke_tests/arabas_et_al_2015/test_spin_up.py index fec14d217..10c905991 100644 --- a/tests/smoke_tests/arabas_et_al_2015/test_spin_up.py +++ b/tests/smoke_tests/arabas_et_al_2015/test_spin_up.py @@ -5,23 +5,12 @@ from PySDM_examples.Szumowski_et_al_1998 import Simulation from PySDM_examples.Arabas_et_al_2015 import Settings, SpinUp from PySDM.physics import si - +from PySDM.formulae import Formulae +from .dummy_storage import DummyStorage from ...backends_fixture import backend_class assert hasattr(backend_class, '_pytestfixturefunction') -class DummyStorage: - def __init__(self): - self.profiles = [] - - def init(*_): # pylint: disable=no-method-argument - pass - - def save(self, data: np.ndarray, step: int, name: str): # pylint: disable=unused-argument - if name == "qv_env": - self.profiles.append({"qv_env": np.mean(data, axis=0)}) - - @pytest.mark.parametrize("fastmath", ( pytest.param(False, id="fastmath: False"), pytest.param(True, id="fastmath: True") @@ -29,7 +18,7 @@ def save(self, data: np.ndarray, step: int, name: str): # pylint: disable=unuse # pylint: disable=redefined-outer-name def test_spin_up(backend_class, fastmath, plot=False): # Arrange - settings = Settings(fastmath=fastmath) + settings = Settings(Formulae(fastmath=fastmath)) settings.dt = .5 * si.second settings.grid = (3, 25) settings.simulation_time = 20 * settings.dt diff --git a/tests/smoke_tests/morrison_and_grabowski_2008/__init__.py b/tests/smoke_tests/morrison_and_grabowski_2008/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/tests/smoke_tests/morrison_and_grabowski_2008/test_just_do_it.py b/tests/smoke_tests/morrison_and_grabowski_2008/test_just_do_it.py deleted file mode 100644 index 7b2827097..000000000 --- a/tests/smoke_tests/morrison_and_grabowski_2008/test_just_do_it.py +++ /dev/null @@ -1,37 +0,0 @@ -# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring -from contextlib import AbstractContextManager -import numpy as np -from PySDM_examples.Szumowski_et_al_1998 import Simulation, Storage -from PySDM_examples.Morrison_and_Grabowski_2008 import ColdCumulus - - -class FlowFieldAsserts(AbstractContextManager): - def __init__(self, simulation): - self.particulator = simulation.particulator - self.panic = None - - def set_percent(self, percent): - if percent == 0: - return - advector = None - solvers = self.particulator.dynamics['EulerianAdvection'].solvers.mpdatas.values() - for solver in solvers: - assert advector is None or advector is solver.advector - advector = solver.advector - np.testing.assert_allclose(advector.get_component(1)[:, 0], 0) - np.testing.assert_allclose(advector.get_component(1)[:, -1], 0, atol=1e-15) - - def __exit__(self, exc_type, exc_val, exc_tb): - pass - - -def test_just_do_it(): - settings = ColdCumulus() - settings.kappa = 0 - for process in settings.processes: - settings.processes[process] = False - settings.processes['particle advection'] += 1 - settings.processes['fluid advection'] += 1 - simulation = Simulation(settings, Storage(), None) - simulation.reinit() - simulation.run(controller=FlowFieldAsserts(simulation)) diff --git a/tests/unit_tests/attributes/test_terminal_velocity.py b/tests/unit_tests/attributes/test_terminal_velocity.py index 8e056e75b..9446b3e16 100644 --- a/tests/unit_tests/attributes/test_terminal_velocity.py +++ b/tests/unit_tests/attributes/test_terminal_velocity.py @@ -2,7 +2,8 @@ import matplotlib.pyplot as plt import numpy as np from PySDM.physics import constants as const -from PySDM.physics.terminal_velocity.gunn_and_kinzer import RogersYau, Interpolation +from PySDM.physics.terminal_velocity.gunn_and_kinzer import Interpolation +from PySDM.physics.terminal_velocity.rogers_and_yau import RogersYau from tests.backends_fixture import backend_class from tests.unit_tests.dummy_particulator import DummyParticulator diff --git a/tests/unit_tests/test_products.py b/tests/unit_tests/test_products.py index 62f82240b..357837491 100644 --- a/tests/unit_tests/test_products.py +++ b/tests/unit_tests/test_products.py @@ -12,7 +12,9 @@ from PySDM.products import (AqueousMassSpectrum, AqueousMoleFraction, TotalDryMassMixingRatio, ParticleSizeSpectrumPerMass, GaseousMoleFraction, FreezableSpecificConcentration, DynamicWallTime, - ParticleSizeSpectrumPerVolume, ParticlesVolumeSpectrum) + ParticleSizeSpectrumPerVolume, + ParticleVolumeVersusRadiusLogarithmSpectrum, + RadiusBinnedNumberAveragedTerminalVelocity) _ARGUMENTS = { AqueousMassSpectrum: {'key': 'S_VI', 'dry_radius_bins_edges': (0, np.inf)}, @@ -23,7 +25,8 @@ FreezableSpecificConcentration: {'temperature_bins_edges': (0, 300)}, DynamicWallTime: {'dynamic': 'Condensation'}, ParticleSizeSpectrumPerVolume: {'radius_bins_edges': (0, np.inf)}, - ParticlesVolumeSpectrum: {'radius_bins_edges': (0, np.inf)} + ParticleVolumeVersusRadiusLogarithmSpectrum: {'radius_bins_edges': (0, np.inf)}, + RadiusBinnedNumberAveragedTerminalVelocity: {'radius_bin_edges': (0, np.inf)} }