Hello!

I am working on a tool in Julia which allows one to redistribute an initial set of particles inside a profile based on SPH (smoothed particle hydrodynamics), to get a more uniform distribution. The algorithm is based on this paper, https://www.researchgate.net/publication/256688486_Particle_packing_algorithm_for_SPH_schemes, for the interested. In pictures this means that one goes from left to right, so that an initial circular distribution inside a square profile becomes distributed homogeneously:

I have uploaded my Julia file and one relevant data file here: https://github.com/AhmedSalih3d/ParticlePacking.jl , but what I specifically want help with is the “PackStep” function I made:

```
#----------------------------------- PackStep -----------------------------------#
# Put one time step of the algorithm in a function
# Inputs are:
# 1. Initial position of particles, pg
# 2. Initial velocity of particles, updated
# 3. Total number of particles, nCoords
function PackStep(pg,u,nCoords)
# Focus on one particle at at a time, iter
for iter = 1:nCoords
# Calculate difference between position of particle iter minus all other particles
rij = map(x -> pg[iter] .- x, pg)
# Using LinearAlgebra, calculate the Euclidean distance between iter and each point
RIJ = norm.(rij)
# Transform it to a normalized distance
q = RIJ ./ H
# For some reason it becomes a 2D array, I prefer it as 1D..
q = dropdims(q,dims=2)
# Find all particles in near vicinity
indices = q.<=2;
# Set the influence of the current particle in focus, iter, to zero
indices[iter] = 0
# Extract only the normalized distance from the closest particles to iter
qc = q[indices]
# Calculate main body of the Wendland Kernel derivative, dW/dq
qq3 = qc .* (qc .- 2).^3;
# Calculate the actual dW/dq
Wq = AD*FAC*qq3
# Extract the distance, between,
# (x_i[1],x_i[2],x_i[3]) - (x_j[1],x_j[2],x_j[3]),
# which was calculated earlier in "rij"
x_ij = first.(rij[indices]);
z_ij = last.(rij[indices]);
# From C++ I have learned to divide before hand.. maybe it doesn't make sense to do so in Julia
RIJ1 = 1 ./ RIJ[indices];
# Calculate actual gradient values for particle iter
Wgx = sum(Wq .* (x_ij .* RIJ1) * H1);
Wgz = sum(Wq .* (z_ij .* RIJ1) * H1);
# Extract ux and uz components and make it 1D array
ux = dropdims(first.(u),dims=2);
uz = dropdims(last.(u),dims=2);
# Calculate change of velocity for particle iter
dux = (-BETA .* Wgx * V - ZETA .* ux[iter])*DT;
duz = (-BETA .* Wgz * V - ZETA .* uz[iter])*DT;
# Calculate displacement of particle iter
dx = dux*DT;
dz = duz*DT;
# Set new velocity and displacement
# NOTE: I am aware I "cheat" since I base calculations
# on updated values over time, so I halfway through
# would have 50% updated and 50% non-updated values.
# Since displacements are so small in this method, it
# really does not do a big difference, if one sticks
# to the original method
u[iter] = u[iter] .+ tuple(dux,0,duz);
pg[iter] = pg[iter] .+ tuple(dx,0,dz);
end
end
```

“pg” is a array of tuple’s which holds the initial circular distribution and fixed square - u is an initial velocity field to initialize the algorithm and nCoords, is the length of number of particles in the circular distribution to ensure that one only calculates for the particles inside of the square, since the square is fixed. Performance wise it is not too bad, since for 100 time steps, it outputs:

```
# Time and run the main body
@time RunAlgorithm()
14.160436 seconds (14.46 M allocations: 50.042 GiB, 6.42% gc time)
```

But of course way too many allocations and memory usage for my liking. I have made my initial implementation using “tuple” since I read that they should be more effective than arrays, but this makes it a bit harder to alter “pg” based on in-place operations.

The code in the Github repository and the datafile should run seamlessly since no dependence on niche packages is used.

I have read the performance tips section, but am still struggling with optimizing this further, without starting to test @threads and parallel computing etc.

Looking forward to hear any feedback/suggestions.

Kind regards