Monte Carlo Method for Displacement Dynamic #1647

Bodzio-2 · 2025-06-10T13:17:49Z

No description provided.

…-Carlo/PySDM_END into monte-carlo-impl

…nt function

Updated the Explicit in space displacement function to be on one line

Updated displacement methods so that the new displacemnt functions compile

Added two numpy functions so the displacment functions uses

add test_single_cell

add the rest of the tests

Eniterusx

minor fixes

Eniterusx · 2025-06-13T10:11:48Z

PySDM/physics/particle_advection/implicit_in_space.py

+    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))


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

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

Eniterusx · 2025-06-13T10:14:44Z

PySDM/physics/particle_advection/explicit_in_space.py

+    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)


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

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

pepe5p · 2025-06-13T18:41:32Z

tests/unit_tests/dynamics/displacement/test_monte_carlo.py

+        sut, particulator = settings.get_displacement(
+            backend_class, scheme="ImplicitInSpace"
+        )


What do you think about parametrizing tests in this file with ["ImplicitInSpace", "ExplicitInSpace"]? You use them selectively, but probably both cases should be tested

Added parametrization in this test!

422d166

pepe5p · 2025-06-13T18:44:18Z

tests/unit_tests/dynamics/displacement/test_monte_carlo.py

+        value_a = 0.1
+        value_b = 0.2
+        weight = 0.125
+        settings.courant_field_data = (
+            np.array([[0, 0]]).T,
+            np.array([[value_a, value_b]]),
+        )


Could you test behaviour with different Courant number e.g 0, 1, 2, 3, 1000?

I think that for basic checking thats enough but it can be added of course in a future

kacpermajchrzak · 2025-06-14T18:19:56Z

tests/unit_tests/dynamics/displacement/test_advection.py

@@ -127,15 +128,16 @@ def test_calculate_displacement_dim1(backend_class):
        )


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 ):

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, )

kacpermajchrzak · 2025-06-14T18:44:18Z

tests/unit_tests/dynamics/displacement/test_monte_carlo.py

+        )
+
+    @staticmethod
+    def test_advection(backend_class):


I see the DisplacementSettings setup is repeated in most tests. Could we extract this into a shared pytest fixture to keep the tests DRY (Don't Repeat Yourself)?

stabor705 · 2025-06-16T09:40:11Z

PySDM/dynamics/displacement.py

    ):
+        if self.enable_sedimentation and self.enable_monte_carlo:
+            dt = self.particulator.dt / self._n_substeps
+            dt_over_dz = dt / self.particulator.mesh.dz


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?

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.

stabor705 · 2025-06-16T09:43:59Z

tests/unit_tests/dynamics/displacement/test_monte_carlo.py

+        if backend_class.__name__ == "Numba":
+            np.testing.assert_equal(
+                sut.displacement[0, slice(0, 1)].to_ndarray(),
+                0.125


So for different backends we get different results? Isn't that a little bit concerning?

This is because our solution didn't implement Monte Carlo in ThrustRTC (GPU), hence why it tests against values that you would get without Monte Carlo, this is surely something to add later down the line, just went out of scope of this PR.

stabor705 · 2025-06-16T09:49:01Z

PySDM/particulator.py

                position_in_cell=position_in_cell,
                n_substeps=n_substeps,
+                enable_monte_carlo=enable_monte_carlo,
+                rng=self.Random(1,1).generator,


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.

jakub-ml · 2025-06-16T17:05:06Z

PySDM/particulator.py

                position_in_cell=position_in_cell,
                n_substeps=n_substeps,
+                enable_monte_carlo=enable_monte_carlo,
+                rng=self.Random(1,1).generator,


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.

jakub-ml · 2025-06-16T17:06:55Z

PySDM/backends/impl_thrust_rtc/methods/displacement_methods.py

                position_in_cell.data,
                trtc.DVInt64(n_substeps),
+                trtc.DVBool(enable_monte_carlo),
+                trtc.DVDouble(1.0) # TODO


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.

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

pepe5p · 2025-06-16T18:08:12Z

PySDM/backends/impl_numba/methods/displacement_methods.py

@@ -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


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)

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

pepe5p · 2025-06-16T18:29:44Z

PySDM/backends/impl_numba/methods/displacement_methods.py

                n_substeps,
+                enable_monte_carlo,
+                rng,


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]),

pepe5p · 2025-06-16T18:32:21Z

PySDM/backends/impl_numba/methods/displacement_methods.py

                n_substeps,
+                enable_monte_carlo,
+                rng,


Here it would be something like this:

Suggested change

rng,

rng[droplet],

The same for 2d and 3d

slayoo · 2025-06-18T06:59:26Z

CI fails (https://github.com/open-atmos/PySDM/actions/runs/15562732280/job/43825166589?pr=1647) due to pylint-reported issues:

************* Module PySDM.backends.impl_numba.methods.displacement_methods
PySDM/backends/impl_numba/methods/displacement_methods.py:17:0: C0301: Line too long (119/100) (line-too-long)
PySDM/backends/impl_numba/methods/displacement_methods.py:34:0: C0301: Line too long (119/100) (line-too-long)
PySDM/backends/impl_numba/methods/displacement_methods.py:60:0: C0301: Line too long (119/100) (line-too-long)
PySDM/backends/impl_numba/methods/displacement_methods.py:92:0: C0301: Line too long (119/100) (line-too-long)
PySDM/backends/impl_numba/methods/displacement_methods.py:123:0: C0301: Line too long (120/100) (line-too-long)
************* Module PySDM.backends.impl_thrust_rtc.methods.displacement_methods
PySDM/backends/impl_thrust_rtc/methods/displacement_methods.py:105:0: C0301: Line too long (120/100) (line-too-long)
PySDM/backends/impl_thrust_rtc/methods/displacement_methods.py:105:117: W0613: Unused argument 'rng' (unused-argument)
************* Module PySDM.dynamics.displacement
PySDM/dynamics/displacement.py:144:69: C0303: Trailing whitespace (trailing-whitespace)
PySDM/dynamics/displacement.py:127:4: R0913: Too many arguments (6/5) (too-many-arguments)
************* Module PySDM.particulator
PySDM/particulator.py:429:0: C0301: Line too long (110/100) (line-too-long)
************* Module PySDM.physics.particle_advection.explicit_in_space
PySDM/physics/particle_advection/explicit_in_space.py:12:0: C0301: Line too long (282/100) (line-too-long)
PySDM/physics/particle_advection/explicit_in_space.py:11:4: R0913: Too many arguments (7/5) (too-many-arguments)
PySDM/physics/particle_advection/explicit_in_space.py:11:42: W0613: Unused argument 'cell_id' (unused-argument)
************* Module PySDM.physics.particle_advection.implicit_in_space
PySDM/physics/particle_advection/implicit_in_space.py:13:0: C0301: Line too long (302/100) (line-too-long)
PySDM/physics/particle_advection/implicit_in_space.py:11:4: R0913: Too many arguments (7/5) (too-many-arguments)
PySDM/physics/particle_advection/implicit_in_space.py:11:42: W0613: Unused argument 'cell_id' (unused-argument)
************* Module tests.unit_tests.dynamics.displacement.displacement_settings
tests/unit_tests/dynamics/displacement/displacement_settings.py:42:0: C0301: Line too long (125/100) (line-too-long)
************* Module tests.unit_tests.dynamics.displacement.test_monte_carlo
tests/unit_tests/dynamics/displacement/test_monte_carlo.py:161:0: C0303: Trailing whitespace (trailing-whitespace)

-----------------------------------
Your code has been rated at 9.99/10

and due to pre-commit not being run or installed (https://github.com/open-atmos/PySDM/actions/runs/15562732326/job/43825167611?pr=1647):

black....................................................................Failed
- hook id: black
- files were modified by this hook

reformatted PySDM/backends/impl_numba/methods/displacement_methods.py
reformatted PySDM/backends/impl_thrust_rtc/methods/displacement_methods.py
reformatted PySDM/dynamics/displacement.py
reformatted PySDM/physics/particle_advection/explicit_in_space.py
reformatted PySDM/physics/particle_advection/implicit_in_space.py
reformatted PySDM/particulator.py
reformatted tests/unit_tests/dynamics/displacement/displacement_settings.py
reformatted tests/unit_tests/dynamics/displacement/test_advection.py
reformatted tests/unit_tests/dynamics/displacement/test_monte_carlo.py

All done! ✨ 🍰 ✨
9 files reformatted, 819 files left unchanged.

For pre-commit help see e.g. https://github.com/open-atmos/python-dev-hints/wiki/Pre-commit-FAQ

…ySDM_END into monte-carlo-impl

Bodzio-2 and others added 13 commits June 10, 2025 15:10

add monte carlo flags for numba displacement_methods.py

8f7f437

update get_displacement for Monte Carlo solution

6cb90ef

monte carlo method for displacement.py

ecc44b2

Merge branch 'monte-carlo-impl' of https://github.com/MoPS-2025-Monte…

d9cb6f9

…-Carlo/PySDM_END into monte-carlo-impl

Update current tests to work with Monte Carlo flag

0490971

changes for monte carlo solution in particulator calculate_displaceme…

76c7d45

…nt function

added calculate displacement tests

3582bd8

Update explicit_in_space.py

2e9bacf

Updated the Explicit in space displacement function to be on one line

Update displacement_methods.py

91c2705

Updated displacement methods so that the new displacemnt functions compile

Update cpp2python.py

3b5cb31

Added two numpy functions so the displacment functions uses

Update test_monte_carlo.py

6a1bd13

add test_single_cell

Update test_monte_carlo.py

6106384

add the rest of the tests

Update implicit_in_space.py

43d1078

Bodzio-2 mentioned this pull request Jun 10, 2025

Monte-Carlo particle transport scheme for PySDM slayoo/AGH-MPS-2025#3

Open

Eniterusx reviewed Jun 13, 2025

View reviewed changes

pepe5p reviewed Jun 13, 2025

View reviewed changes

kacpermajchrzak reviewed Jun 14, 2025

View reviewed changes

stabor705 reviewed Jun 16, 2025

View reviewed changes

jakub-ml reviewed Jun 16, 2025

View reviewed changes

pepe5p reviewed Jun 16, 2025

View reviewed changes

ForNeus57 and others added 5 commits June 19, 2025 23:21

Changes to make reviewers happy

4c78097

Merge branch 'main' into monte-carlo-impl

4f10f47

Fix to make the tests pass :)))

715f02d

Merge branch 'monte-carlo-impl' of github.com:MoPS-2025-Monte-Carlo/P…

34b666e

…ySDM_END into monte-carlo-impl

parametrize displacement scheme in test_single_cell

422d166

		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):
-        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

		@@ -127,15 +128,16 @@ def test_calculate_displacement_dim1(backend_class):
		)

Monte Carlo Method for Displacement Dynamic #1647

Are you sure you want to change the base?

Monte Carlo Method for Displacement Dynamic #1647

Uh oh!

Conversation

Bodzio-2 commented Jun 10, 2025

Uh oh!

Eniterusx left a comment

Choose a reason for hiding this comment

Uh oh!

Eniterusx Jun 13, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Eniterusx Jun 13, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

kubijaku Jun 23, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Variable naming

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pepe5p Jun 16, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

slayoo commented Jun 18, 2025

Uh oh!

Uh oh!

Eniterusx Jun 13, 2025 •

edited

Loading

Eniterusx Jun 13, 2025 •

edited

Loading

kubijaku Jun 23, 2025 •

edited

Loading

pepe5p Jun 16, 2025 •

edited

Loading