How to update the neighbours, which the positions are being updated?

[dineshadepu]

Hey everyone,

I'm currently working on implementing an SPH-DEM solver in Cabana to analyze particle dispersion in fluid tanks. I've managed to write the code for both the SPH and DEM components, where the DEM particles follow rigid body dynamics and interact with SPH fluid particles through discretized dummy particles.

This is not a bug but a question I am dealing with. I'm facing a common challenge in both components during the time loop. I need to update the neighbors for each particle, but I'm struggling with the implementation. Here is part of my code:

  ListType verlet_list( aosoa_position, 0,
            aosoa_position.size(), neighborhood_radius,
            cell_ratio, grid_min, grid_max );

  // Main timestep loop
  for ( int step = 0; step < steps; step++ )
    {
      rigid_body_stage_1(rb, dt);
      rigid_body_particles_stage_1(aosoa, rb, dt, rigid_limits);

      rigid_body_stage_2(rb, dt);
      rigid_body_particles_stage_2(aosoa, rb, dt, rigid_limits);

      // Compute the neighbours
      ListType verlet_list( aosoa_position, 0,
                aosoa_position.size(), neighborhood_radius,
                cell_ratio, grid_min, grid_max );

      compute_force_on_rigid_body_particles(aosoa, dt,
                        &verlet_list,
                        rigid_limits);

      body_force_rigid_body(rb, gravity, dt);
      compute_effective_force_and_torque_on_rigid_body(rb,
                               aosoa);

      rigid_body_stage_3(rb, dt);
      rigid_body_particles_stage_3(aosoa, rb, dt, rigid_limits);

      if ( step % print_freq == 0 )
    {
      std::cout << "Time is:" << time << std::endl;
      output_data(aosoa, num_particles, step, time);
      output_rb_data(rb, no_of_bodies, step, time);
    }

      time += dt;

    }

Here, I am creating the verlet_list variable every time in the loop, and I am not sure if it is hampering the performance. Is there any way to just update it with the updated positions, provided it working on all the architectures independently?

I will share the fully code in the next comment. Thank you so much for the help :) .

[dineshadepu]

https://gist.github.com/dineshadepu/a5df26833844fea204b406dae3d578f2

A small demonstration of the design:

I have created two types AoSoA's, where the one AoSoA is for the particles each rigid body composed of, and the other AoSoA (RBAoSoA) is to handle the dynamics of the center of rigid bodies. The rigid bodies interact using the discretized points on the rigid body. We only use functions to handle the update, interaction or any part of the algorithm. Easy way to look at the code is to see the time loop (https://gist.github.com/dineshadepu/a5df26833844fea204b406dae3d578f2#file-bodies_settling_in_tank_under_gravity_cabana-cpp-L1690). I need to seperate out the functions into respective files and make it a package. The current version is not commented, sorry.

One can run this example using the following command:

./01RBFreelyRotatingBodies 0.1 1.0 1.0 6 3.0 100

By changing number (6) corresponding to the no of bodies, will increase the no of bodies. This will output, two kinds of files, one is particles_*.xmf and rigid_bodies_*.xmf, which can be visualised using paraview. Save the following script in a python file (visualise.py) and run paraview visualise.py, and visualise :)

# trace generated using paraview version 5.12.0
#import paraview
#paraview.compatibility.major = 5
#paraview.compatibility.minor = 12

#### import the simple module from the paraview
from paraview.simple import *
#### disable automatic camera reset on 'Show'
paraview.simple._DisableFirstRenderCameraReset()

import os

# =====================================
# start: get the files and sort
# =====================================
files = [filename for filename in os.listdir('.') if filename.startswith("particles") and filename.endswith("xmf") ]
files.sort()
files_num = []
for f in files:
    f_last = f[10:]
    files_num.append(int(f_last[:-4]))
files_num.sort()

sorted_files = []
for num in files_num:
    sorted_files.append("particles_" + str(num) + ".xmf")
files = sorted_files
# =====================================
# end: get the files and sort
# =====================================

# =====================================
# start: get the files and sort
# =====================================
files_rb = [filename for filename in os.listdir('.') if filename.startswith("rigid_bodies") and filename.endswith("xmf") ]
files_rb.sort()
# print(files_rb)
files_rb_num = []
for f in files_rb:
    f_last = f[13:]
    files_rb_num.append(int(f_last[:-4]))
files_rb_num.sort()

sorted_files_rb = []
for num in files_rb_num:
    sorted_files_rb.append("rigid_bodies_" + str(num) + ".xmf")
# print(sorted_files_rb)
files_rb = sorted_files_rb
# =====================================
# end: get the files and sort
# =====================================
# print(files)
# print(files_rb)

# create a new 'XDMF Reader'
rigid_bodies_0xmf = XDMFReader(registrationName='rigid_bodies_0.xmf*', FileNames=files_rb)

# get animation scene
animationScene1 = GetAnimationScene()

# update animation scene based on data timesteps
animationScene1.UpdateAnimationUsingDataTimeSteps()

# set active source
SetActiveSource(rigid_bodies_0xmf)

# get active view
renderView1 = GetActiveViewOrCreate('RenderView')

# show data in view
rigid_bodies_0xmfDisplay = Show(rigid_bodies_0xmf, renderView1, 'UnstructuredGridRepresentation')

# trace defaults for the display properties.
rigid_bodies_0xmfDisplay.Representation = 'Surface'

# show color bar/color legend
rigid_bodies_0xmfDisplay.SetScalarBarVisibility(renderView1, True)

#changing interaction mode based on data extents
renderView1.InteractionMode = '2D'
renderView1.CameraPosition = [1.15, 4.855, 0.0]
renderView1.CameraFocalPoint = [1.15, 0.5, 0.0]
renderView1.CameraViewUp = [1.0, 0.0, 0.0]

# get the material library
materialLibrary1 = GetMaterialLibrary()

# reset view to fit data
renderView1.ResetCamera(False, 0.9)

# get color transfer function/color map for 'rb_ids'
rb_idsLUT = GetColorTransferFunction('rb_ids')

# get opacity transfer function/opacity map for 'rb_ids'
rb_idsPWF = GetOpacityTransferFunction('rb_ids')

# get 2D transfer function for 'rb_ids'
rb_idsTF2D = GetTransferFunction2D('rb_ids')

# change representation type
rigid_bodies_0xmfDisplay.SetRepresentationType('Point Gaussian')

# set scalar coloring
ColorBy(rigid_bodies_0xmfDisplay, ('POINTS', 'rb_force', 'Magnitude'))

# Hide the scalar bar for this color map if no visible data is colored by it.
HideScalarBarIfNotNeeded(rb_idsLUT, renderView1)

# rescale color and/or opacity maps used to include current data range
rigid_bodies_0xmfDisplay.RescaleTransferFunctionToDataRange(True, False)

# show color bar/color legend
rigid_bodies_0xmfDisplay.SetScalarBarVisibility(renderView1, True)

# get color transfer function/color map for 'rb_force'
rb_forceLUT = GetColorTransferFunction('rb_force')

# get opacity transfer function/opacity map for 'rb_force'
rb_forcePWF = GetOpacityTransferFunction('rb_force')

# get 2D transfer function for 'rb_force'
rb_forceTF2D = GetTransferFunction2D('rb_force')

# create a new 'XDMF Reader'
particles_0xmf = XDMFReader(registrationName='particles_0.xmf*', FileNames=files)

# set active source
SetActiveSource(particles_0xmf)

# show data in view
particles_0xmfDisplay = Show(particles_0xmf, renderView1, 'UnstructuredGridRepresentation')

# trace defaults for the display properties.
particles_0xmfDisplay.Representation = 'Surface'

# show color bar/color legend
particles_0xmfDisplay.SetScalarBarVisibility(renderView1, True)

# get color transfer function/color map for 'ids'
idsLUT = GetColorTransferFunction('ids')

# get opacity transfer function/opacity map for 'ids'
idsPWF = GetOpacityTransferFunction('ids')

# get 2D transfer function for 'ids'
idsTF2D = GetTransferFunction2D('ids')

# change representation type
particles_0xmfDisplay.SetRepresentationType('Point Gaussian')

# set scalar coloring
ColorBy(particles_0xmfDisplay, ('POINTS', 'frc_dem', 'Magnitude'))

# Hide the scalar bar for this color map if no visible data is colored by it.
HideScalarBarIfNotNeeded(idsLUT, renderView1)

# rescale color and/or opacity maps used to include current data range
particles_0xmfDisplay.RescaleTransferFunctionToDataRange(True, False)

# show color bar/color legend
particles_0xmfDisplay.SetScalarBarVisibility(renderView1, True)

# get color transfer function/color map for 'frc_dem'
frc_demLUT = GetColorTransferFunction('frc_dem')

# get opacity transfer function/opacity map for 'frc_dem'
frc_demPWF = GetOpacityTransferFunction('frc_dem')

# get 2D transfer function for 'frc_dem'
frc_demTF2D = GetTransferFunction2D('frc_dem')

renderView1.ResetActiveCameraToNegativeZ()

# reset view to fit data
renderView1.ResetCamera(False, 0.9)

#================================================================
# addendum: following script captures some of the application
# state to faithfully reproduce the visualization during playback
#================================================================

# get layout
layout1 = GetLayout()

#--------------------------------
# saving layout sizes for layouts

# layout/tab size in pixels
layout1.SetSize(1542, 751)

#-----------------------------------
# saving camera placements for views

# current camera placement for renderView1
renderView1.InteractionMode = '2D'
renderView1.CameraPosition = [1.15, 0.5, 4.845059295778354]
renderView1.CameraFocalPoint = [1.15, 0.5, 0.0]
renderView1.CameraParallelScale = 1.2539936203984452

##--------------------------------------------
## You may need to add some code at the end of this python script depending on your usage, eg:
#
## Render all views to see them appears
# RenderAllViews()
#
## Interact with the view, usefull when running from pvpython
# Interact()
#
## Save a screenshot of the active view
# SaveScreenshot("path/to/screenshot.png")
#
## Save a screenshot of a layout (multiple splitted view)
# SaveScreenshot("path/to/screenshot.png", GetLayout())
#
## Save all "Extractors" from the pipeline browser
# SaveExtracts()
#
## Save a animation of the current active view
# SaveAnimation()
#
## Please refer to the documentation of paraview.simple
## https://kitware.github.io/paraview-docs/latest/python/paraview.simple.html
##--------------------------------------------

.

[streeve]

This will absolutely affect performance in two main ways: first constructing every step will do a lot of initialization memory allocation that doesn't need to be re-done, second presuming a relatively consistent neighborhood size re-iteration during search can be minimized.

After construction, the loop can call: verlet_list.build( aosoa_position, 0, aosoa_position.size(), neighborhood_radius, cell_ratio, grid_min, grid_max );

You can also call:

max_neigh_guess = 100; // could also depend on previous max_neighbor count
verlet_list.build( aosoa_position, 0, aosoa_position.size(), neighborhood_radius, cell_ratio, grid_min, grid_max, max_neigh_guess );

for the additional benefits if you have some way to estimate how many will be found in the next neighbor search. The more you overallocate the more unused memory you will have, but you will avoid reallocation more often.

[dineshadepu]

Thank you Sam. I will close the issue as it addresses the problem. However, I have a performance related problem, which I will raise in a different issue.