Skip to content

Commit

Permalink
Merge pull request #712 from slayoo/freezing
Browse files Browse the repository at this point in the history
adding test_freezing for ICMW (with both singular & time-dependent runs and an assert for non-zero ice-water content)
  • Loading branch information
slayoo authored Dec 27, 2021
2 parents e4d89f7 + 6901491 commit 86a41d7
Show file tree
Hide file tree
Showing 39 changed files with 404 additions and 259 deletions.
4 changes: 2 additions & 2 deletions .github/workflows/julia.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions .github/workflows/matlab.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
4 changes: 2 additions & 2 deletions .github/workflows/python.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 2 additions & 0 deletions PySDM/attributes/physics/volume.py
Original file line number Diff line number Diff line change
@@ -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

Expand Down
51 changes: 3 additions & 48 deletions PySDM/backends/impl_numba/methods/collisions_methods.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
5 changes: 3 additions & 2 deletions PySDM/backends/impl_numba/methods/freezing_methods.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
56 changes: 56 additions & 0 deletions PySDM/backends/impl_numba/methods/terminal_velocity_methods.py
Original file line number Diff line number Diff line change
@@ -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
)
3 changes: 2 additions & 1 deletion PySDM/backends/impl_numba/storage_impl.py
Original file line number Diff line number Diff line change
Expand Up @@ -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}})
Expand Down
74 changes: 0 additions & 74 deletions PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py
Original file line number Diff line number Diff line change
Expand Up @@ -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'),
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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
Expand Down
Original file line number Diff line number Diff line change
@@ -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)
)
Loading

0 comments on commit 86a41d7

Please sign in to comment.