Skip to content
Open
Show file tree
Hide file tree
Changes from 13 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 26 additions & 5 deletions PySDM/backends/impl_numba/methods/displacement_methods.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,12 +14,15 @@
@numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}})
# pylint: disable=too-many-arguments
def calculate_displacement_body_common(
dim, droplet, scheme, _l, _r, displacement, courant, position_in_cell, n_substeps
dim, droplet, scheme, _l, _r, displacement, courant, position_in_cell, cell_id, n_substeps, enable_monte_carlo, rng
Copy link

@pepe5p pepe5p Jun 16, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see you pass rng object which is numpy Random Generator for this jit function. I think it would be better to pass already created random uniform distribution as array (see supported types for numba)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have added your changes and the tests that I run locally run correctly :)

):
displacement[dim, droplet] = scheme(
position_in_cell[dim, droplet],
cell_id[droplet],
courant[_l] / n_substeps,
courant[_r] / n_substeps,
enable_monte_carlo,
rng.uniform(0.,1.)
)


Expand All @@ -28,7 +31,7 @@ class DisplacementMethods(BackendMethods):
@numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False, "cache": False}})
# pylint: disable=too-many-arguments
def calculate_displacement_body_1d(
dim, scheme, displacement, courant, cell_origin, position_in_cell, n_substeps
dim, scheme, displacement, courant, cell_origin, position_in_cell, cell_id, n_substeps, enable_monte_carlo, rng
):
length = displacement.shape[1]
for droplet in numba.prange(length): # pylint: disable=not-an-iterable
Expand All @@ -44,14 +47,17 @@ def calculate_displacement_body_1d(
displacement,
courant,
position_in_cell,
cell_id,
n_substeps,
enable_monte_carlo,
rng,
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here it would be something like this:

Suggested change
rng,
rng[droplet],

The same for 2d and 3d

)

@staticmethod
@numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False, "cache": False}})
# pylint: disable=too-many-arguments
def calculate_displacement_body_2d(
dim, scheme, displacement, courant, cell_origin, position_in_cell, n_substeps
dim, scheme, displacement, courant, cell_origin, position_in_cell, cell_id, n_substeps, enable_monte_carlo, rng
):
length = displacement.shape[1]
for droplet in numba.prange(length): # pylint: disable=not-an-iterable
Expand All @@ -73,14 +79,17 @@ def calculate_displacement_body_2d(
displacement,
courant,
position_in_cell,
cell_id,
n_substeps,
enable_monte_carlo,
rng,
)

@staticmethod
@numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False, "cache": False}})
# pylint: disable=too-many-arguments
def calculate_displacement_body_3d(
dim, scheme, displacement, courant, cell_origin, position_in_cell, n_substeps
dim, scheme, displacement, courant, cell_origin, position_in_cell, cell_id, n_substeps, enable_monte_carlo, rng
):
n_sd = displacement.shape[1]
for droplet in numba.prange(n_sd): # pylint: disable=not-an-iterable
Expand All @@ -104,11 +113,14 @@ def calculate_displacement_body_3d(
displacement,
courant,
position_in_cell,
cell_id,
n_substeps,
enable_monte_carlo,
rng,
)

def calculate_displacement(
self, *, dim, displacement, courant, cell_origin, position_in_cell, n_substeps
self, *, dim, displacement, courant, cell_origin, position_in_cell, cell_id, n_substeps, enable_monte_carlo, rng
):
n_dims = len(courant.shape)
scheme = self.formulae.particle_advection.displacement
Expand All @@ -120,7 +132,10 @@ def calculate_displacement(
courant.data,
cell_origin.data,
position_in_cell.data,
cell_id.data,
n_substeps,
enable_monte_carlo,
rng,
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

According to previous command it is how the code would look like here. It is the same for 2d and 3d

Suggested change
rng,
rng.uniform(0., 1., displacement.data.shape[1]),

)
elif n_dims == 2:
DisplacementMethods.calculate_displacement_body_2d(
Expand All @@ -130,7 +145,10 @@ def calculate_displacement(
courant.data,
cell_origin.data,
position_in_cell.data,
cell_id.data,
n_substeps,
enable_monte_carlo,
rng,
)
elif n_dims == 3:
DisplacementMethods.calculate_displacement_body_3d(
Expand All @@ -140,7 +158,10 @@ def calculate_displacement(
courant.data,
cell_origin.data,
position_in_cell.data,
cell_id.data,
n_substeps,
enable_monte_carlo,
rng,
)
else:
raise NotImplementedError()
Expand Down
13 changes: 11 additions & 2 deletions PySDM/backends/impl_thrust_rtc/methods/displacement_methods.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,11 @@ def __calculate_displacement_body(self):
"courant_shape_0",
"courant_shape_1",
"cell_origin",
"cell_id",
"position_in_cell",
"n_substeps",
"enable_monte_carlo",
"u01",
),
name_iter="i",
body=f"""
Expand All @@ -52,8 +55,11 @@ def __calculate_displacement_body(self):
displacement[i + n_sd * dim] = {
self.formulae.particle_advection.displacement.c_inline(
position_in_cell="position_in_cell[i + n_sd * dim]",
cell_id="cell_id[i + n_sd * dim]",
c_l="courant[_l] / n_substeps",
c_r="courant[_r] / n_substeps"
c_r="courant[_r] / n_substeps",
enable_monte_carlo="false",
u01="u01"
)
};
""".replace(
Expand Down Expand Up @@ -96,7 +102,7 @@ def __flag_precipitated_body(self):

@nice_thrust(**NICE_THRUST_FLAGS)
def calculate_displacement(
self, *, dim, displacement, courant, cell_origin, position_in_cell, n_substeps
self, *, dim, displacement, courant, cell_origin, position_in_cell, cell_id, n_substeps, enable_monte_carlo, rng
):
n_dim = len(courant.shape)
n_sd = position_in_cell.shape[1]
Expand All @@ -110,8 +116,11 @@ def calculate_displacement(
trtc.DVInt64(courant.shape[0]),
trtc.DVInt64(courant.shape[1] if n_dim > 2 else -1),
cell_origin.data,
cell_id.data,
position_in_cell.data,
trtc.DVInt64(n_substeps),
trtc.DVBool(enable_monte_carlo),
trtc.DVDouble(1.0) # TODO

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In line 123, the value 1.0 is hardcoded:
trtc.DVDouble(1.0) # TODO

This leads to deterministic behavior in Monte Carlo simulations, which should instead use a random value from a uniform distribution. The propossed solution:
trtc.DVDouble(rng.uniform(0.0, 1.0))

This is the same issue as in the CPU backend, where the RNG is deterministically seeded – here the randomness is replaced by a constant.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for mentioning this, when I was debugging the GPU backend I must have missed this. I think now the problem is solved.

),
)

Expand Down
3 changes: 2 additions & 1 deletion PySDM/backends/impl_thrust_rtc/test_helpers/cpp2python.py
Original file line number Diff line number Diff line change
Expand Up @@ -199,7 +199,8 @@ def make(self):
arcsinh as asinh,
sinh,
maximum, minimum,
where, # TODO #1295
where, # TODO #1295
absolute, sign
)
import numba

Expand Down
16 changes: 13 additions & 3 deletions PySDM/dynamics/displacement.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,12 +21,14 @@ class Displacement: # pylint: disable=too-many-instance-attributes
def __init__(
self,
enable_sedimentation=False,
enable_monte_carlo=False,
precipitation_counting_level_index: int = 0,
adaptive=DEFAULTS.adaptive,
rtol=DEFAULTS.rtol,
): # pylint: disable=too-many-arguments
self.particulator = None
self.enable_sedimentation = enable_sedimentation
self.enable_monte_carlo = enable_monte_carlo
self.dimension = None
self.grid = None
self.courant = None
Expand Down Expand Up @@ -101,11 +103,12 @@ def __call__(self):
# TIP: not need all array only [idx[:sd_num]]
cell_origin = self.particulator.attributes["cell origin"]
position_in_cell = self.particulator.attributes["position in cell"]
cell_id = self.particulator.attributes["cell id"]

self.precipitation_mass_in_last_step = 0.0
for _ in range(self._n_substeps):
self.calculate_displacement(
self.displacement, self.courant, cell_origin, position_in_cell
self.displacement, self.courant, cell_origin, position_in_cell, cell_id
)
self.update_position(position_in_cell, self.displacement)
if self.enable_sedimentation:
Expand All @@ -122,16 +125,23 @@ def __call__(self):
self.particulator.attributes.mark_updated(key)

def calculate_displacement(
self, displacement, courant, cell_origin, position_in_cell
self, displacement, courant, cell_origin, position_in_cell, cell_id
):
if self.enable_sedimentation and self.enable_monte_carlo:
dt = self.particulator.dt / self._n_substeps
dt_over_dz = dt / self.particulator.mesh.dz

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

From above to variables names it is not obvious what they represent. Maybe delta_time and delta_altitude (if z is altitude) would be better?

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Variable naming

Thank you for the proposed update concerning the naming of the certain variables!

I agree that clarity in variable naming is important, for both maintainability and readability of the code. In this case, however, the names dt and dt_over_dz are intentionally kept concise and consistent with standard nomenclature widely used in numerical methods and literature.

With that being said, we can consider adding a short inline comment to clarify their meaning if that would help in understanding the functionality.

courant += self.particulator.attributes["relative fall velocity"] / dt_over_dz

self.particulator.calculate_displacement(
displacement=displacement,
courant=courant,
cell_origin=cell_origin,
cell_id=cell_id,
position_in_cell=position_in_cell,
n_substeps=self._n_substeps,
enable_monte_carlo=self.enable_monte_carlo,
)
if self.enable_sedimentation:
if self.enable_sedimentation and not self.enable_monte_carlo:
displacement_z = displacement[self.dimension - 1, :]
dt = self.particulator.dt / self._n_substeps
dt_over_dz = dt / self.particulator.mesh.dz
Expand Down
5 changes: 4 additions & 1 deletion PySDM/particulator.py
Original file line number Diff line number Diff line change
Expand Up @@ -426,16 +426,19 @@ def flag_out_of_column(self):
self.attributes.sanitize()

def calculate_displacement(
self, *, displacement, courant, cell_origin, position_in_cell, n_substeps
self, *, displacement, courant, cell_origin, cell_id, position_in_cell, n_substeps, enable_monte_carlo
):
for dim in range(len(self.environment.mesh.grid)):
self.backend.calculate_displacement(
dim=dim,
displacement=displacement,
courant=courant[dim],
cell_origin=cell_origin,
cell_id=cell_id,
position_in_cell=position_in_cell,
n_substeps=n_substeps,
enable_monte_carlo=enable_monte_carlo,
rng=self.Random(1,1).generator,

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function slowly starts being too complex. There are at least two flags used in a little bit confusing way (if flag1 and flag2 and after that if flag1 and not flag2).
The fact that we are passing rng which I presume we only need when enable_monte_carlo == True already sounds like reason enough to create a new function such as calculate_displacement_monte_carlo.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the displacement calculation, the Random generator is instantiated with a hardcoded seed (seed=1) for each backend call
rng=self.Random(1,1).generator
This makes the backend behavior fully deterministic, which is incorrect for Monte Carlo simulations.

)

def isotopic_fractionation(self, heavy_isotopes: tuple):
Expand Down
6 changes: 3 additions & 3 deletions PySDM/physics/particle_advection/explicit_in_space.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
"""
basic explicit-in-space Euler scheme
"""

import numpy as np

class ExplicitInSpace: # pylint: disable=too-few-public-methods
def __init__(self, _):
pass

@staticmethod
def displacement(_, position_in_cell, c_l, c_r):
return c_l * (1 - position_in_cell) + c_r * position_in_cell
def displacement(_, position_in_cell, cell_id, c_l, c_r, enable_monte_carlo, u01):
return (position_in_cell + (np.floor(np.abs(max(c_l, c_r))) * np.sign(np.abs(max(c_l, c_r)) / max(c_l, c_r))) + (np.abs(max(c_l, c_r)) > u01) * np.sign(np.abs(max(c_l, c_r)) / max(c_l, c_r))) if enable_monte_carlo else (c_l * (1 - position_in_cell) + c_r * position_in_cell)
Copy link

@Eniterusx Eniterusx Jun 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above, unreadable.

Suggested change
def displacement(_, position_in_cell, cell_id, c_l, c_r, enable_monte_carlo, u01):
return (position_in_cell + (np.floor(np.abs(max(c_l, c_r))) * np.sign(np.abs(max(c_l, c_r)) / max(c_l, c_r))) + (np.abs(max(c_l, c_r)) > u01) * np.sign(np.abs(max(c_l, c_r)) / max(c_l, c_r))) if enable_monte_carlo else (c_l * (1 - position_in_cell) + c_r * position_in_cell)
def displacement(_, position_in_cell, cell_id, c_l, c_r, enable_monte_carlo, u01):
if enable_monte_carlo:
max_c = max(c_l, c_r)
abs_max_c = np.abs(max_c)
sign_max_c = np.sign(abs_max_c / max_c)
floor_term = np.floor(abs_max_c) * sign_max_c
jump_term = (abs_max_c > u01) * sign_max_c
return position_in_cell + floor_term + jump_term
else:
numerator = c_l * (1 - position_in_cell) + c_r * position_in_cell
return numerator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have changed the file to make it more readable :)

7 changes: 4 additions & 3 deletions PySDM/physics/particle_advection/implicit_in_space.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,13 @@
"""
eqs. 14-16 in [Arabas et al. 2015](https://doi.org/10.5194/gmd-8-1677-2015)
"""

import numpy as np

class ImplicitInSpace: # pylint: disable=too-few-public-methods
def __init__(self, _):
pass

@staticmethod
def displacement(_, position_in_cell, c_l, c_r):
return (c_l * (1 - position_in_cell) + c_r * position_in_cell) / (1 - c_r + c_l)
def displacement(_, position_in_cell, cell_id, c_l, c_r, enable_monte_carlo, u01):

return (position_in_cell + (np.floor(np.abs(max(c_l, c_r))) * np.sign(np.abs(max(c_l, c_r)) / max(c_l, c_r))) + (np.abs(max(c_l, c_r)) > u01) * np.sign(np.abs(max(c_l, c_r)) / max(c_l, c_r))) if enable_monte_carlo else ((c_l * (1 - position_in_cell) + c_r * position_in_cell) / (1 - c_r + c_l))
Copy link

@Eniterusx Eniterusx Jun 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line is absolutely unreadable. Please consider splitting it into a few lines for the sake of clarity.

Suggested change
def displacement(_, position_in_cell, cell_id, c_l, c_r, enable_monte_carlo, u01):
return (position_in_cell + (np.floor(np.abs(max(c_l, c_r))) * np.sign(np.abs(max(c_l, c_r)) / max(c_l, c_r))) + (np.abs(max(c_l, c_r)) > u01) * np.sign(np.abs(max(c_l, c_r)) / max(c_l, c_r))) if enable_monte_carlo else ((c_l * (1 - position_in_cell) + c_r * position_in_cell) / (1 - c_r + c_l))
def displacement(_, position_in_cell, cell_id, c_l, c_r, enable_monte_carlo, u01):
if enable_monte_carlo:
max_c = max(c_l, c_r)
abs_max_c = np.abs(max_c)
sign_max_c = np.sign(abs_max_c / max_c)
floor_term = np.floor(abs_max_c) * sign_max_c
jump_term = (abs_max_c > u01) * sign_max_c
return position_in_cell + floor_term + jump_term
else:
numerator = c_l * (1 - position_in_cell) + c_r * position_in_cell
denominator = 1 - c_r + c_l
normalized_position = numerator / denominator
return normalized_position

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have changed the file to make it more readable :)

Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ def __init__(self, n_sd=1, grid=None, positions=None, courant_field_data=None):
self.sedimentation = False
self.dt = None

def get_displacement(self, backend, scheme, adaptive=True):
def get_displacement(self, backend, scheme, adaptive=True, enable_monte_carlo=False):
formulae = Formulae(particle_advection=scheme)
particulator = DummyParticulator(backend, n_sd=len(self.n), formulae=formulae)
particulator.environment = DummyEnvironment(
Expand All @@ -39,7 +39,7 @@ def get_displacement(self, backend, scheme, adaptive=True):
"position in cell": position_in_cell,
}
particulator.build(attributes)
sut = Displacement(enable_sedimentation=self.sedimentation, adaptive=adaptive)
sut = Displacement(enable_sedimentation=self.sedimentation, adaptive=adaptive, enable_monte_carlo=enable_monte_carlo)
sut.register(particulator)
sut.upload_courant_field(self.courant_field_data)

Expand Down
22 changes: 12 additions & 10 deletions tests/unit_tests/dynamics/displacement/test_advection.py
Original file line number Diff line number Diff line change
Expand Up @@ -97,15 +97,16 @@ def test_calculate_displacement(backend_class):
)
settings.positions = [[weight], [0]]
sut, particulator = settings.get_displacement(
backend_class, scheme="ExplicitInSpace", adaptive=False
backend_class, scheme="ExplicitInSpace", adaptive=False, enable_monte_carlo=False
)

# Act
sut.calculate_displacement(
sut.displacement,
sut.courant,
particulator.attributes["cell origin"],
particulator.attributes["position in cell"],
displacement=sut.displacement,
courant=sut.courant,
cell_origin=particulator.attributes["cell origin"],
position_in_cell=particulator.attributes["position in cell"],
cell_id=particulator.attributes["cell id"],
)

# Assert
Expand All @@ -127,15 +128,16 @@ def test_calculate_displacement_dim1(backend_class):
)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The tests test_calculate_displacement and test_calculate_displacement_dim1 are very similar. They both test the calculate_displacement method, but with different input data to check behavior along different dimensions. This structural duplication can be reduced.

@staticmethod @pytest.mark.parametrize( "courant_field_data, positions, asserted_dimension, asserted_slice", [ ( (np.array([[0.1, 0.2]]).T, np.array([[0, 0]])), [[0.25], [0]], 0, slice(0, 1) ), ( (np.array([[0, 0]]).T, np.array([[0.1, 0.2]])), [[0], [0.25]], 1, slice(0, 1) ), ], ) def test_calculate_displacement_by_dimension( backend_class, courant_field_data, positions, asserted_dimension, asserted_slice ):

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are absolutely right! Those two test methods are nearly identical and can definitely be refactored into a single parametrized test. Here is how I would implement your suggestion:

@staticmethod
@pytest.mark.parametrize(
    "courant_field_data, positions, asserted_dimension, asserted_slice",
    [
        (
            (np.array([[0.1, 0.2]]).T, np.array([[0, 0]])),
            [[0.25], [0]],
            0,
            slice(0, 1)
        ),
        (
            (np.array([[0, 0]]).T, np.array([[0.1, 0.2]])),
            [[0], [0.25]],
            1,
            slice(0, 1)
        ),
    ],
)
def test_calculate_displacement_by_dimension(
    backend_class, courant_field_data, positions, asserted_dimension, asserted_slice
):
    # Arrange
    settings = DisplacementSettings()
    value_a = 0.1
    value_b = 0.2
    weight = 0.25
    settings.courant_field_data = courant_field_data
    settings.positions = positions
    sut, particulator = settings.get_displacement(
        backend_class, scheme="ExplicitInSpace", adaptive=False, enable_monte_carlo=False
    )

    # Act
    sut.calculate_displacement(
        displacement=sut.displacement,
        courant=sut.courant,
        cell_origin=particulator.attributes["cell origin"],
        position_in_cell=particulator.attributes["position in cell"],
        cell_id=particulator.attributes["cell id"],
    )

    # Assert
    np.testing.assert_equal(
        sut.displacement[asserted_dimension, asserted_slice].to_ndarray(),
        (1 - weight) * value_a + weight * value_b,
    )

settings.positions = [[0], [weight]]
sut, particulator = settings.get_displacement(
backend_class, scheme="ExplicitInSpace", adaptive=False
backend_class, scheme="ExplicitInSpace", adaptive=False, enable_monte_carlo=False
)

# Act
sut.calculate_displacement(
sut.displacement,
sut.courant,
particulator.attributes["cell origin"],
particulator.attributes["position in cell"],
displacement=sut.displacement,
courant=sut.courant,
cell_origin=particulator.attributes["cell origin"],
position_in_cell=particulator.attributes["position in cell"],
cell_id=particulator.attributes["cell id"],
)

# Assert
Expand Down
Loading
Loading