Implementing a neighbor list like data structure
Stephan
Dear developers,
I am trying to implement something which requires me to use a variable per particle pair at the previous timestep (we will refer to this variable as n_ij from here). The reason that I need n_ij on the previous time step is because I can compute n_ij on the current timestep only by making use of the previous n_ij value. You can look at it as a variable similar to RIJ. RIJ is always computed from the current particle positions but it is never stored. If you would like to know RIJ on the previous time step one way to do this would be to store additional variables x_old, y_old and z_old in the particle arrays and make use of the cached neighbor list information to compute the RIJ on the previous time step. Unfortunately, for our variable of interest this doesn’t work because n_ij at the current time is computed by n_ij at the previous time. The most efficient way to do this would be to include a neighborlist like structure to store the n_ij variable. That would mean that we do not only store an array of neighbor indices for all our particles but also an array containing the corresponding n_ij values. This can be explained as follows:
Say the neighborlist of particle i looks like [2, 5, 32], then we would like to store the n_ij values in an array of same size in which the values correspond to the pair, i.e. [n_i2, n_i5, n_i32]. These values then need to be cached, so it would look something very similar to creating a neighborlist and setting the caching option to True.
Are there any suggestions on how to do this most efficiently? I can imagine that it would be possible by adjusting the NNPS methods of the code. There are plenty of examples and documentations about how to write your own equations, integrators and kernels and use them via create_solver() and create_equations(), but I couldn’t find anything on how to write your own NNPS methods. I did see that users can overload the create_NNPS() method but I have no idea on how to use this method and if this is the most efficient or preferable way to go (not sure if it is even possible to do what we want with the create_NNPS() method.
Thanks in advance!
Best,
Stephan
A Dinesh
class UsingOldValuesSampleEquation(Equation):
def loop(self, d_idx, s_id, d_x_old, d_y_old, d_z_old, s_x_old, s_y_old,
s_z_old, RIJ, XIJ):
# ========================================== #
# Compute the old 'nij' vector i.e. 'nij_old'
# ========================================== #
dx_old = d_x_old[d_idx] - s_x_old[s_idx]
dy_old = d_y_old[d_idx] - s_y_old[s_idx]
dz_old = d_z_old[d_idx] - s_z_old[s_idx]
rij_old = (dx_old**2. + dy_old**2. + dz_old**2.)**0.5
nij_x_old = 0.
nij_y_old = 0.
nij_z_old = 0.
if rij_old > 1e-12:
nij_x_old = dx_old / rij_old
nij_y_old = dy_old / rij_old
nij_z_old = dz_old / rij_old
# ========================================== #
# Compute the current 'nij' vector i.e. 'nij'
# ========================================== #
nij_x = 0.
nij_y = 0.
nij_z = 0.
if RIJ > 1e-12:
nij_x = XIJ[0] / RIJ
nij_y = XIJ[1] / RIJ
nij_z = XIJ[2] / RIJ
- There is no guarantee that the previous neighbor is a current neighbor.
- The number of neighbors could change at the next time step. (We can create the length of the array equal to maximum neighbours (len_of_particles * 30 (possibly maximum number of neighbours))
class TrackingTemplateEquation(Equation):
def __init__(self, dest, sources):
super(TrackingTemplateEquation, self).__init__(dest, sources)
def loop(self, d_idx, d_m, d_u, d_v, d_w,
d_cnt_idx, d_nij_x_old, d_nij_y_old, d_nij_z_old,
d_total_contacts, d_max_contacts_limit,
XIJ, RIJ, s_idx, s_m, s_u, s_v, s_w, dt, t):
p, q1, tot_ctcs, j, found_at, found = declare('int', 6)
# total number of contacts of particle i in destination
tot_ctcs = d_total_contacts[d_idx]
# d_idx has a range of tracking indices with sources
# starting index is p
p = d_idx * d_max_contacts_limit[0]
# ending index is q -1
q1 = p + tot_ctcs
# ===================================================
# Here we are dynamically adding the new particles
# if they are previously not a neighbour
# ===================================================
# check if the particle is in the tracking list
# if so, then save the location at found_at
found = 0
for j in range(p, q1):
if s_idx == d_cnt_idx[j]:
found_at = j
found = 1
break
# if the particle is not been tracked then assign an index in
# tracking history.
if found == 0:
found_at = q1
d_cnt_idx[found_at] = s_idx
d_total_contacts[d_idx] += 1
# implies we are tracking the particle
else:
# get the previous `nij` values
nij_x_old = d_nij_x_old[found_at]
nij_y_old = d_nij_y_old[found_at]
nij_z_old = d_nij_z_old[found_at]
# Compute current 'nij' values
nij_x = 0.
nij_y = 0.
nij_z = 0.
if RIJ > 1e-12:
nij_x = XIJ[0] / RIJ
nij_y = XIJ[1] / RIJ
nij_z = XIJ[2] / RIJ
# ================================
# Do whatever computation you want here
# using nij, and nij_old
# ================================
# xxxxxxxxxxxxxxxxxxxxxxxxxx
# ================================
# ================================
# now save the current n_ij values for next time step computation
d_nij_x_old[found_at] = nij_x
d_nij_y_old[found_at] = nij_y
d_nij_z_old[found_at] = nij_z
--
You received this message because you are subscribed to the Google Groups "pysph-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/68c59e0b-d8b5-42f3-9750-c03092bd035an%40googlegroups.com.
Stephan
Dear Dinesh,
Thanks for your reply. This really got us on the right track. I have several remaining questions. The property named cnt_idx should be defined as an integer property with stride 30 I suppose. Does that also mean that the maximum amount of neighbors used in PySPH is 30 in general? What would be the best way to incorporate cnt_idx into get_particle_arrays? I currently use fluid.add_property(name='cnt_idx', type='int', stride=30).
For the n_ij variable I use fluid.add_property(name='n_ij', type='double', stride=30) but after the simulation I get an error message:
double free or corruption (!prev)
Aborted (core dumped)
Could you tell if I'm doing something wrong?
Cheers,
Stephan
A Dinesh
To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/ca085284-910b-4b0e-97a2-f4386de7fde7n%40googlegroups.com.
Prabhu Ramachandran
Dear Dinesh,
Thanks for your reply. This really got us on the right track. I have several remaining questions. The property named cnt_idx should be defined as an integer property with stride 30 I suppose. Does that also mean that the maximum amount of neighbors used in PySPH is 30 in general? What would be the best way to incorporate cnt_idx into get_particle_arrays? I currently use fluid.add_property(name='cnt_idx', type='int', stride=30).
cheers,
Prabhu