This is just a rework of the SPH Dam break simulation optimized to get better performance we will focus on the optimization and differences with the previous example
Simulation video 1The first optimization that we operate is the usage of verlet list. The verlet are reconstructed only when the maximum displacement is bigger than the half skin. Because we have to calculate the maximum displacement the verlet and euler integration has been modified to do this. The function accept a reference to max_disp that is filled with the maximum displacement calculated from these functions.
The variable is reset inside verlet and euler time integration function
while iteration across particle the maximum displacement is saved inside the variable max_disp
We also have to be careful that if we are removing particles we have to reconstruct the verlet list, so we set it to a really big number
Because the maximum displacement has to be calculated across processors, we use the function max in Vcluster to calculate the maximum displacement across processors.
We also calculate the skin part and ghost plus skin. Consider also that the ghost must be extended to ghost + skin so r_gskin
As we explained before, we update the verlet only if particles move more than the skin. In case they move more than the skin we do first a map to redistribute the particles and in the meanwhile we check if it is a good moment to rebalance. We decided to combine these two steps because in case we rebalance we have anyway to reconstruct the Verler-list. Than we calculate the pressure for all the particles, refresh the ghost, update the Verlet-list and reset the total displacement. In case the the total displacement does not overshoot the skin we just calculate the pressure for all the particles and refresh the ghost. We must use the option SKIP_LABELLING
We pass the max_displacement variable to verlet_int and euler_int function. We also add the maximum displacement per iteration to the total maximum displacement
Symmetric interaction give the possibility to reduce the computation by half and speed-up your simulation. To do this we have to do some changes into the function calc forces. Symmetric interaction require to write on the ghost area. So at the beginning of the function we reset the ghost part. In the meanwhile because we have the external force gravity that operate per particles, we set this force at the beginning.
With this code we set the per particle external force to gravity and reset the derivative of the density for the domain particles
With this code we reset the force and derivative of the density of the particles on the ghost part
Small changes must be done to iterate over the neighborhood particles
skip the self interaction
This is instead an important change (and honestly it took quite some hour of debuging to discover the problem). In case we are on boundary particle (p = boundary particle) and calculating an interaction with a particle q = fluid particle we have to remeber that we have also to calculate the force for q (not only drho)
for a fluid particle instead we calculate p-q interaction and we add the contribution to p and q. Because we do not integrate over the boundary particles we can also avoid to check that q is a boundary particle
After the calculation cycle we have to merge the forces and delta density calculated on the ghost with the real particles.
It is important when we construct our vector of particles to pass the option BIND_DEC_TO_GHOST. To use symmetric calculation in parallel environment the decomposition must be consistent with the cell decomposition of the space.
To construct a Verlet-list using the CRS scheme we use the following function
while to update the verlet list we use the following
Using redecompose instead of decompose produce less jumping decomposition during the simulation