Smooth Particle Hydrodynamics (SPH) and Hashing
Today, we’ll enhance the SPH simulation program developed in the previous article by speeding it up 24x using hashing. Previously, we identified that the neighbor interaction method was inefficient because each particle iterates through all the other particles in the domain. By using hashing, we can reduce the number of iterations and avoid this inefficiency. So, what exactly is hashing?
Prerequisites: Read the previous article, Basic C++ programming skills, and familiarity with the SPH method.
Required tools: Linux operating system (Ubuntu), g++ compiler, VS Code, and Paraview (all open-source).
Project Repository: https://github.com/AhmedElbossily/SPH_Hashing.git
Introduction
Hashing is a big topic in computer science and is used for many applications including cryptography, password Storage, data storage and retrieval, and more. What we are concerned about today is data storage and retrieval or spatial hashing. The core idea of hashing is to transform input data of arbitrary size into a fixed size using what is known as a hash value. The function that performs this transformation is called a hash function. But how does this definition apply to our problem? To simplify, we’ll illustrate the concept in 2D space, which can then be easily extended to 3D.
Spatial Hashing
Let’s break it down step by step. As shown in the figure below, we begin with 36 particles and then draw a box around them. This small adjustment allows us to maintain a fixed data size. What do you mean by fixed data size? It’s simple: regardless of the number of particles inside the rectangle, the rectangle, itself, remains unchanged. This is the core idea of transforming input data of varying sizes into a fixed size.
Hash Function and Hash Values
After defining the fixed-size data — represented by the rectangle — a hash function is created based on the rectangle and its grid, which is why a grid is displayed in the figure below. This hash function is then used to assign a hash value (hv) to each particle within the rectangle as follows:
where (ix,iy) represent the grid cell indices in which particle i is located. nx and ny denote the number of grid cells in the x and y directions, respectively. px and py correspond to the current coordinates of the particle’s position.
Magic of Hashing
As shown above, particles within the same cell are represented with the same color, indicating that each particle is now aware of its neighbors. But how can this be useful in SPH? If we set the grid cell size equal to the smoothing length (h), we can limit our iteration to the cells surrounding the cell containing our particle of interest i, as illustrated below.
But how can we identify the neighboring cells marked in yellow in the figure above? The magic of hashing allows us to do this. Starting with the point of interest i, we can determine its grid cell coordinates (ix, iy) from its hash value by “unhashing” it using the following equations:
With grid cell indices (ix, iy), we can move one step backward, one step forward, one step up, and one step down to define the domain where we’ll search for neighbors. Why search in all four directions? Because we don’t know the exact location of particle i within its grid cell. However, this approach can be further optimized. Anyway, you can notice with hashing, the search domain is reduced to only 9 cells.
Neighbor Retrieving
But how can we retrieve the particles inside these cells? The answer lies in the hash value. By iterating through the 9 neighboring cells and calculating the hash value for each cell based on its indices (ix, iy), we can identify all the particles inside that cell.
Get Ahmed Elbossily’s stories in your inbox
Join Medium for free to get updates from this writer.
To achieve this efficiently, we will sort the particle data by their hash values. Then, we will store the start and end indexes for each hash value in two additional arrays. For example, in the figure below, if we want to retrieve all the particles within the cell with a hash value of (1), the cell_start array will indicate that we should begin retrieving data from the sorted particles at index (3), while the cell_end array will show where to stop, which is at index (5).
Dive into Code
We will start from the code developed in the previous article, by declaring an object for the Grid that contains all the necessary information about it.
class Grid
{
public:
// Constructor and destructor
Grid(float3_t minCorner, float3_t maxCornerInit, float_t h);
~Grid() {};
// cell start and end indices arrays
int *cellStart = 0;
int *cellEnd = 0;
// helper methods
void compute_hashes(Particles *particles); // compute hash values for all particles
int compute_hash(int i, int j, int k); // compute hash values for a cell
void compute_unhashes(int hv, int &i, int &j, int &k); // unhash a hash value to get the cell indices
void hash_sort(Particles *particles); // sort particles based on hash values
void compute_cellStartEnd(Particles *particles); // compute cell start and end indices for particles
// Getter methods
long int get_nx() const;
long int get_ny() const;
long int get_nz() const;
long int get_numCell() const;
float3_t get_minCorner() const;
float3_t get_maxCorner() const;
private:
float_t dx; // cell size
float3_t l; // grid size
long int nx = 0, ny = 0, nz = 0; // number of cells in each direction
long int numCell = 0; // total number of cells
float3_t minCorner, maxCornerInit; // min and max corners of the grid
};Grid Class Constructor
The constructor method of the Grid class initiates the grid size and calculates the number of grid cells in each dimension, and the total number of cells. Additionally, it allocates memory on the heap for the two arrays: cellStart and cellEnd.
Grid::Grid(float3_t minCorner, float3_t maxCorner, float_t h)
{
// set min and max corners of the grid
// add a small value to max and subtract a small value
// from min to avoid particles on the boundary
this->minCorner = minCorner - 1.e-6;
this->maxCorner = maxCorner + 1.e-6;
dx = h; // set cell size equal to the smoothin length
l = maxCorner - minCorner; // compute grid size
// compute number of cells in each direction
nx = ceil(l.x / dx);
ny = ceil(l.y / dx);
nz = ceil(l.z / dx);
// compute total number of cells
numCell = nx * ny * nz;
// allocate memory for cellStart and cellEnd
cellStart = new int[numCell];
cellEnd = new int[numCell];
// initialize cellStart with max value to indicate empty cells
std::fill(cellStart, cellStart + numCell, std::numeric_limits<int>::max());
}Compute Hash Values
The compute_hashes method generates hash values for each particle based on their coordinates, grid dimensions, and grid cell indices in which the particle is located.
void Grid::compute_hashes(Particles *particles)
{
for (size_t i = 0; i < particles->N; i++)
{
// get the position of the particle
float_t _px = particles->pos[i].x;
float_t _py = particles->pos[i].y;
float_t _pz = particles->pos[i].z;
// compute the grid cell indices
int ix = (_px - minCorner.x) / dx;
int iy = (_py - minCorner.y) / dx;
int iz = (_pz - minCorner.z) / dx;
// compute the hash value
particles->hash[i] = iz * ny * nx + iy * nx + ix;
}
}Unhashing
compute_unhashes method determines a grid cell indices from a given hash value.
void Grid::compute_unhashes(int hv, int &ix, int &iy, int &iz)
{
ix = hv % nx; // compute grid cell x-indice
iy = (hv / nx) % ny; // compute grid cell y-indices
iz = hv / (nx * ny);// compute grid cell z-indices
}Particles sorting
The hash_sort method sorts the properties of the particles based on their hash values. It uses a separate vector of indexes to keep track of the sorting order and then rearranges the properties accordingly.
void Grid::hash_sort(Particles *particles)
{
// create a separate vector of indices
std::vector<size_t> indices(particles->N);
// fill the vector with increasing values starting from 0
std::iota(indices.begin(), indices.end(), 0);
// sort the indices based on the hash values
std::sort(indices.begin(), indices.end(), [&](size_t i, size_t j)
{ return particles->hash[i] < particles->hash[j]; });
// buffers to store sorted properties
std::vector<float3_t> sortedPos(particles->N);
std::vector<float_t> sortedT(particles->N);
std::vector<float_t> sortedCp(particles->N);
std::vector<float_t> sortedK(particles->N);
std::vector<float_t> sortedRho(particles->N);
std::vector<float_t> sortedH(particles->N);
std::vector<float_t> sortedMass(particles->N);
std::vector<int> sortedHash(particles->N);
// sort the properties based on the sorted indices and save them in buffers
for (size_t i = 0; i < particles->N; i++)
{
size_t index = indices[i];
sortedPos[i] = particles->pos[index];
sortedT[i] = particles->T[index];
sortedCp[i] = particles->cp[index];
sortedK[i] = particles->k[index];
sortedRho[i] = particles->rho[index];
sortedH[i] = particles->h[index];
sortedMass[i] = particles->mass[index];
sortedHash[i] = particles->hash[index];
}
// copy the sorted properties back to the original arrays
for (size_t i = 0; i < particles->N; i++)
{
particles->pos[i] = sortedPos[i];
particles->T[i] = sortedT[i];
particles->cp[i] = sortedCp[i];
particles->k[i] = sortedK[i];
particles->rho[i] = sortedRho[i];
particles->h[i] = sortedH[i];
particles->mass[i] = sortedMass[i];
particles->hash[i] = sortedHash[i];
}
}Filling cellStart and cellEnd arrays
The compute_cellStartEnd method finds the start and end indexes of every hash value and stored in cellStart and cellEnd arrays, respectively.
void Grid::compute_cellStartEnd(Particles *particles)
{
// get the first hash value (first cell)
int currentHash = particles->hash[0];
// set the start index of the first cell
cellStart[currentHash] = 0;
// iterate over the particles to compute the start and end indices of each cell
for (size_t i = 1; i < particles->N; i++)
{
// if the hash value of the current particle is different from the previous one
if (particles->hash[i] != currentHash)
{
// set the end index of the previous cell
cellEnd[currentHash] = i - 1;
// update the current hash value
currentHash = particles->hash[i];
// set the start index of the current cell
cellStart[currentHash] = i;
}
}
// set the end index of the last cell
cellEnd[currentHash] = particles->N - 1;
}Integrating Grid with SPH Simulation Program
Now, it is time to utilize the Grid object by creating an instance of it using the minimum and maximum particle positions within the grid, along with the smoothing length.
- After creating the grid, we can verify that the grid box was created correctly by visualizing it with the vtk_writer_write_box method.
- Then, the hash values for all particles are calculated using the compute_hashes method.
- Next, the particles’ properties are rearranged (sorted) in ascending order based on their hash values using the hash_sort method.
- Moreover, the cellStart and cellEnds arrays are populated using the compute_cellStartEnd method.
- Finally, the grid is passed to run_simulation to be used within the interaction method.
// main function
Particles *particles = new Particles(N);
initialize_particles(particles, nx, ny, nz, dx, dy, dz, hf);
// create grid with min and max corners and cell size
Grid* grid = new Grid(particles->pos[0], particles->pos[N-1], hf * dx);
vtk_writer_write_box(grid);
grid->compute_hashes(particles);
grid->hash_sort(particles);
grid->compute_cellStartEnd(particles);
run_simulation(particles, grid,total_time, dt); // run simulationGrid and Particles Interactions
The interaction method iterates through all particles. For each particle i,
- It retrieves its hash value.
- Using this hash value, the particle’s grid cell indices (ix, iy, iz) are determined with the compute_unhashes method.
- Next, the method identifies the neighboring grid cells by the lowest and highest corners of the grid neighboring cells. while accounting for boundary particles.
- By determining the lowest and highest corners of these neighboring grid cells, the method can iterate through them.
- For each cell: The hash value is computed using compute_hash. This hash value is then used to retrieve the cellStart and cellEnd indexes, which serve as the start and end indexes for searching for neighbors within the entire particles.
void interactions_heat(Particles *particles, Grid *grid)
{
unsigned int N = particles->N;
for (size_t i = 0; i < N; i++)
{
float3_t pi = particles->pos[i]; // position of particle i
float_t Ti = particles->T[i]; // temperature of particle i
float_t cpi = particles->cp[i]; // specific heat capacity of particle i
float_t ki = particles->k[i]; // thermal conductivity of particle i
float_t rhoi = particles->rho[i]; // density of particle i
float_t hi = particles->h[i]; // smoothing length of particle i
float_t mass = particles->mass[i]; // mass of particle
float_t T_ti = 0.; // time derivative of temperature for particle i
float_t hash = particles->hash[i]; // hash value of particle i
// compute cell indices for particle i
int ix, iy, iz;
grid->compute_unhashes(hash, ix, iy, iz);
// lowest indice of the neighboring cells
// -1 to go one step back
// adding condition to avoid negative indices
// and indices greater than the grid size for boundary particles
int low_i = ix - 1 < 0 ? 0 : ix - 1;
int low_j = iy - 1 < 0 ? 0 : iy - 1;
int low_k = iz - 1 < 0 ? 0 : iz - 1;
// highest indice of the neighboring cells
// +2 to go two steps forward
// Why +2? Because we are not considering the last cell in the loop below
int high_i = ix + 2 > grid->get_nx() ? grid->get_nx() : ix + 2;
int high_j = iy + 2 > grid->get_ny() ? grid->get_ny() : iy + 2;
int high_k = iz + 2 > grid->get_nz() ? grid->get_nz() : iz + 2;
for (size_t ii = low_i; ii < high_i; ii++)
{
for (size_t jj = low_j; jj < high_j; jj++)
{
for (size_t kk = low_k; kk < high_k; kk++)
{
// compute hash value for the neighboring cell
int h_idx = grid->compute_hash(ii, jj, kk);
// get the start and end indices of the particles in the neighboring cell
int start = grid->cellStart[h_idx];
int end = grid->cellEnd[h_idx];
// if the cell is empty, continue to the next cell
if (start == std::numeric_limits<int>::max())
continue;
for (size_t j = start; j <= end; j++)
{
float3_t pj = particles->pos[j]; // position of particle j
float_t Tj = particles->T[j]; // temperature of particle j
float_t cpi = particles->cp[j]; // specific heat capacity of particle j
float_t kj = particles->k[j]; // thermal conductivity of particle j
float_t rhoj = particles->rho[j]; // density of particle j
// distance between particle i and j
float_t xij = pi.x - pj.x;
float_t yij = pi.y - pj.y;
float_t zij = pi.z - pj.z;
float_t rij_2 = xij * xij + yij * yij + zij * zij;
float_t rij = sqrt(rij_2);
float4_t ww = cubic_spline(pi, pj, hi); // cubic spline kernel
float_t w = ww.x; // weight
float_t w_x = ww.y; // weight gradient
float_t w_y = ww.z; // weight gradient
float_t w_z = ww.w; // weight gradient
// avoid division by zero
if (rij > 1e-8)
{
float_t eijx = xij / rij;
float_t eijy = yij / rij;
float_t eijz = zij / rij;
float_t rij1 = 1. / rij;
T_ti += ((4 * ki * kj) / (ki + kj)) * (mass / rhoj) * (Ti - Tj) * rij1 *
(eijx * w_x + eijy * w_y + eijz * w_z);
}
}
}
}
}
float_t alpha = 1.0 / (rhoi * cpi); // thermal diffusivity of particle i
particles->T_t[i] = alpha * T_ti; // time derivative of temperature for particle i
}
}Running Time
Running the current code on an 11th Gen Intel® Core™ i7–11800H @ 2.30GHz with 16 GB of RAM takes ~ 1 minute, 24 times faster than the initial version of the code developed in the previous article.
Let’s Recap
We started by exploring the hashing concept specifically spatial hashing. Spatial hashing efficiently maps spatial data (sparticles) to a discrete grid. The main idea is to divide the continuous space into a grid of cells, where each cell is associated with a unique hash value. This allows for fast lookup and storage of spatial data, making it easier to perform neighbor searches.
What is next?!
Can we make it even faster? While 24 times faster is good, there’s always better. We’ll explore this further in our next topic:
“SPH and CPU Multi-threading.”
Acknowledgment
- This project makes use of the excellent iwf_mfree_gpu_3d open-source project developed at ETHZ.
Useful Links
- Project Repository: https://github.com/AhmedElbossily/SPH_Hashing.git
- Paper: “Contemporary Meshfree Methods for Three Dimensional Heat Conduction Problems”
- Paper: “Metal cutting simulations using smoothed particle hydrodynamics on the GPU”
- Paraview: https://www.paraview.org/download/
- VTK file format documentation: https://docs.vtk.org/en/latest/design_documents/VTKFileFormats.html
