Monte-Carlo particle transport scheme for PySDM

[slayoo]

Displacement of particles in PySDM is handled using:

• https://github.com/open-atmos/PySDM/blob/main/PySDM/dynamics/displacement.py
• https://github.com/open-atmos/PySDM/blob/main/PySDM/backends/impl_numba/methods/displacement_methods.py
• https://github.com/open-atmos/PySDM/blob/main/PySDM/backends/impl_thrust_rtc/methods/displacement_methods.py

Tests are here: https://github.com/open-atmos/PySDM/tree/main/tests/unit_tests/dynamics/displacement

The new Monte-Carlo formulation could either be an option of the existing "dynamic" or a new "dynamic" class.

[kubijaku]

So we should just add one method for MonteCarlo solutions in displacement.py? or what exactly are we supposed to do?

[slayoo]

So we should just add one method for MonteCarlo solutions in displacement.py? or what exactly are we supposed to do?

I would suggest:

• adding a flag to the existing Displacement constructor which would trigger switching from the existing deterministic transport to a probabilistic one
• adding pytest.mark.parametrize decorator (with appropriate True and False values for this flag) to the displacement tests thus making them run for both types of displacement handling (original deterministic and the new probabilistic)
• implementing the Monte-Carlo displacement for CPU and/or GPU backend
• introducing logic which would drop tracking of x, y and z co-ordinates for the particles if probabilistic transport is used - cell id is then enough
• adding a new notebook to the examples which would depict a comparison between the original and new behaviour for a 1D or 2D simulation
• adding a smoke test that verifies if the comparison presented in that notebook produces an expected result

[MarcinZ20]

Hi,
I seem to have a problem with assigning myself to the issue.
Could someone please help me with this?

[Bodzio-2]

@slayoo Questions about implementation:

• If i understand correctly, we're ignoring displacement within cell in this scenario? Because monte carlo method will only move droplets between cells if the probability check is met?
• Assuming yes to above, where is the movement across cells implemented in PySDM, or how would we output the displacement to force a movement to another cell?

For any of these pointing to specific parts of code for explanation would be of great help.

I'll ask more questions if anything comes to mind

[slayoo]

* If i understand correctly, we're ignoring displacement within cell in this scenario? Because monte carlo method will only move droplets between cells if the probability check is met?

Droplet positions (x,y,z) are no longer needed with Monte-Carlo transport, hence one cannot even define a displacement within a cell. Instead of using (x,y,z) as droplet attributes, and deriving cell id from it, it should be enough to use the cell id as the only location attribute.

For simplicity, though, an initial MC transport implementation could use the existing setups with (x,y,z) attributes, and shift positions by dx/dy/dz.

* Assuming yes to above, where is the movement across cells implemented in PySDM, or how would we output the displacement to force a movement to another cell?

First, Particulator.step() calls the __call__() method of the Displacement dynamic here:
https://github.com/open-atmos/PySDM/blob/a955f80b15b195caa551e4afb975426c7a502b3a/PySDM/particulator.py#L53

Then, the Displacement dynamic calculates the displacement and updates particle positions:
https://github.com/open-atmos/PySDM/blob/a955f80b15b195caa551e4afb975426c7a502b3a/PySDM/dynamics/displacement.py#L107-L110

Then, in the Displacement.calculate_displacement() method, the Particulator.calculate_displacement() is called:
https://github.com/open-atmos/PySDM/blob/a955f80b15b195caa551e4afb975426c7a502b3a/PySDM/dynamics/displacement.py#L127

Then, the backend's calculate_displacement is called here:
https://github.com/open-atmos/PySDM/blob/a955f80b15b195caa551e4afb975426c7a502b3a/PySDM/particulator.py#L432

In the case of CPU backend:
https://github.com/open-atmos/PySDM/blob/a955f80b15b195caa551e4afb975426c7a502b3a/PySDM/backends/impl_numba/methods/displacement_methods.py#L110

In the case of GPU backend:
https://github.com/open-atmos/PySDM/blob/a955f80b15b195caa551e4afb975426c7a502b3a/PySDM/backends/impl_thrust_rtc/methods/displacement_methods.py#L98

Unit tests are here:
https://github.com/open-atmos/PySDM/tree/main/tests/unit_tests/dynamics/displacement

For a "hello-world" example case with visualisation, I'd suggest the single-dimensional setup, like:
https://github.com/open-atmos/PySDM/blob/main/examples/PySDM_examples/Shipway_and_Hill_2012/fig_1.ipynb

HTH

[Bodzio-2]

The advection implementation in the deterministic solution uses c_l and c_r to calculate the displacement. During the labs, the calculations were done using just one provided Courant number. Should we adapt the lab method to work with c_l and c_r (if so then how should we go about doing it?) or would we acquire the Courant number in a different way?

• Monte Carlo Advection from lab

def advect(*, params_p, params_t, particles, rng):
    """ performs Monte-Carlo advection of the particles """
    probability_of_shift = abs(params_t.courant_number)
    assert probability_of_shift < 1

    sign = int(abs(params_t.courant_number) / params_t.courant_number)
    for _ in range(params_t.n_step):
        u01 = rng.uniform(0, 1, params_p.n_part)
        for part in particles.values():
            part['cell'] += (probability_of_shift > u01) * sign

• Displacement calculation formula from PySDM

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)

[Bodzio-2]

Also, should we worry about the cell_origin attribute or can we ignore it? I'm not entirely sure what it represents exactly.

[slayoo]

@Bodzio-2, in the implicit_in_space.py file, just above the code you quoted, there's a link to a paper with equation numbers:

So, c_l and c_r are the Courant numbers at left-hand and right-hand edges of a cell around point i, respectively.
During the lab, we assumed c_l=c_r, and of course this is a valid initial test case.

Also, should we worry about the cell_origin attribute or can we ignore it? I'm not entirely sure what it represents exactly.

"Cell origin" are the coordinates of the corner of a grid cell, so adding "cell origin" to "position within cell" we get the position of a particle in the domain. Note that all positions and cell origins are expressed in units of cell sizes (uniform grid is assumed), hence cell origin is integer, and "position within cell" is a float between 0 and 1. Recomputing "cell origin" after displacement is done here - should not be needed to worry about it:
https://github.com/open-atmos/PySDM/blob/44a9b05d4da0fd04ae683cd2e1f9eea008ebfcad/PySDM/dynamics/displacement.py#L146-L150

Also, cell IDs are defined using "cell origin:
https://github.com/open-atmos/PySDM/blob/44a9b05d4da0fd04ae683cd2e1f9eea008ebfcad/PySDM/impl/mesh.py#L47-L49

HTH

[ForNeus57]

Prerequisites:

• We need to provide randomness when calculating the displacement of particles. The code from the laboratory creates a u01 vector at each time step simulating if a given super particle should move. There is a need to provide in some form the random generator when evaluating particle movement.

  # Laboratory code
  def advect(*, params_p, params_t, particles, rng):
      probability_of_shift = abs(params_t.courant_number)
      assert probability_of_shift < 1
  
      sign = int(abs(params_t.courant_number) / params_t.courant_number)
  
      for _ in range(params_t.n_step):
          u01 = rng.uniform(0., 1., params_p.n_part)
          for part in particles.values():
              part['cell'] += (probability_of_shift > u01) * sign
  
  advect(particles=PARTICLES, params_t=PARAMS_T, params_p=PARAMS_P, rng=RNG)

• It appears that the ExplicitInSpace.displacement and ImplicitInSpace.displacement static methods and corresponding classes should be stateless. Meaning that they should not have any reference to backend randomness generators. If we create an alternative class (thank you @Bodzio-2) for the Monte Carlo displacement, than propagating a u01 vector directly to displacement static method arguments would not comply with ExplicitInSpace.displacement and ImplicitInSpace.displacement interfaces.

  class MonteCarlo:
      def __init__(self, _):
          pass
  
      @staticmethod
      def displacement(_, position_in_cell, cell_id, c_l, c_r):
          c = max(c_l, c_r)
          probability_of_shift = abs(c)
          assert probability_of_shift < 1, ""
  
          sign = int(abs(c) / c)
          u01 = 0 # rng call here
          return cell_id + (probability_of_shift > u01) * sign

• Our idea was to modify the interface so that some state (backend randomness generators) can be bounded to DisplacementMethods.formulae.particle_advection object, by utilizing the __init__ of the MonteCarlo?

Question:

• Is there other special way the randomness should be supplied to the displacement staticmethod in MonteCarlo?