Sitemap
Open in app

Sign in

Medium Logo
Write
Search

Sign in

Develop an SPH Thermal Simulation Program in C++

13 min readAug 19, 2024

--

Smoothed Particle Hydrodynamics (SPH) is a mesh-free computational method representing continuous fields through kernel functions centered on discrete particles carrying the system’s physical properties. It utilizes a kernel function W to interpolate field equations by considering neighboring particles within a smoothing length h, see the figure below. This kernel interpolation technique enables SPH to accurately approximate field variables and their derivatives, resulting in smooth and continuous solutions without needing a predefined mesh structure.

SPH kernel function (for 2D space)

Today, we will explore building an SPH thermal simulation program from scratch. We’ll begin with a simple problem and gradually introduce new features and optimizations throughout the upcoming series.

Prerequisites: 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_thermal_CPU.git

Get Ahmed Elbossily’s stories in your inbox

Join Medium for free to get updates from this writer.

This article is structured as follows: It begins by defining a thermal problem that will be simulated using the C++ program we will develop today. Next, we will discuss the governing equations and their discretization using the SPH method. After that, we will dive into the code, walking through it step by step to build the project. Additionally, we will cover how to build and compile the code. Finally, the article concludes by highlighting the current limitations of the code and potential optimizations, which will be explored in our next article.

Problem Definition

A 0.1 x 0.1 x 0.1-meter aluminum cube has one face with an initial temperature of 500°C. What is the temperature distribution within the cube after 10 seconds? assuming there is no convection to the surroundings.

Problem initial state

Problem Understanding

The energy equation governing the current problem can be simplified to

e˙ is the change in energy over time. ρ is the density. k is the thermal conductivity and T is the temperature. By employing the SPH method, the above equation can be approximated to be applied to discrete particles as follows

cp_i represents the heat capacity and rho represents the density, while N denotes the total number of neighboring particles for particle i. The distance between particle i and its neighboring particle j is expressed as x_ij. The smoothing function, W_ij, is employed in the calculation. A “Cubic B-spline” function is utilized for this purpose and is defined as follows:

The gradient of the smoothing function is given as

The simulation time of 10 seconds will be divided into small intervals of 0.01 seconds each. At every interval, the temperature time derivative for each particle will be calculated and used to update its current temperature. The following chart illustrates this process.

Press enter or click to view image in full size
SPH thermal simulation diagram

Project Structure

Before diving into the code, let’s have an overview of the file structure we will develop.

```
.
├── readme.md
├── CMakeLists.txt
├── include
│   ├── glm
│   ├── interaction.h
│   ├── kernel.h
│   ├── particle.h
│   ├── type.h
│   └── vtk_writer.h
└── src
    ├── interaction.cpp
    ├── kernel.cpp
    ├── main.cpp
    ├── particle.cpp
    └── vtk_writer.cpp
```

Dive into The Code

Where do we start? We’ll begin by selecting the data type for our simulation. Should we use double-precision accuracy or single-precision? We’ll opt for a generic data type that can be easily switched between single and double precision. Additionally, we’ll use the GLM library to provide vector and matrix data types. This will be implemented as follows:

// type.h file
#ifndef TYPES_H_
#define TYPES_H_

#include <stdio.h>
#include <glm/glm.hpp>

#define USE_DOUBLE
#ifdef USE_DOUBLE
// built-in vector types
#define float_t double
#define float2_t glm::dvec2
#define float3_t glm::dvec3
#define float4_t glm::dvec4

#define make_float2_t glm::dvec2
#define make_float3_t glm::dvec3
#define make_float4_t glm::dvec4

// mathematical functions
#define sqrt_t sqrt
#define exp_t exp
#define log_t log

#else

// built-in vector types
#define float_t float
#define float2_t glm::vec2
#define float3_t glm::vec3
#define float4_

#define make_float2_t glm::vec2
#define make_float3_t glm::vec3
#define make_float4_t glm::vec4

// mathematical functions
#define sqrt_t sqrtf
#define exp_t expf
#define log_t logf
#endif

#endif

Declaring an Object for Particles

We will create a struct to hold all particle data, including positions, thermal conductivities, temperatures, …etc. The struct will include a constructor method to allocate memory for these member variables on the heap and a destructor to free the memory. Additionally, a setter method will be used to initialize the struct’s member variables.

// particle.h file
#ifndef PARTICLE
#define PARTICLE

#include "type.h"
#include <vector>

struct Particles
{
 // position of particles in 3D x,y,z
 float3_t *pos = 0;
 float_t *T = 0;    // current temperature
 float_t *k;     // thermal conductivity
 float_t *rho = 0;  // density of particles
 float_t *mass = 0; // mass of particles
 float_t *cp;    // specific heat capacity
 float_t *h = 0;    // smoothing length
 float_t *T_t = 0;  // time derivative of temperature
 unsigned int N;    // count number of patricles

 // class constructor and destructor
 Particles(unsigned int N);
 ~Particles();

 // set properties of particles
 void set_properties(std::vector<float3_t> pos_vec, std::vector<float_t> T_vec,
      std::vector<float_t> cp_vec, std::vector<float_t> k_vec,
      std::vector<float_t> rho_vec, std::vector<float_t> h_vec,
      std::vector<float_t> mass_vec);
};

#endif
// particle.cpp file
#include "particle.h"
// Constructor method
Particles::Particles(unsigned int N)
{
    this->N = N; 
    pos = new float3_t[N]; h = new float_t[N]; rho = new float_t[N];
    T = new float_t[N]; T_t = new float_t[N]; cp = new float_t[N]; 
    k = new float_t[N]; mass = new float_t[N];
}
// Destructor method
Particles::~Particles()
{
    delete[] pos; delete[] h; delete[] rho; delete[] T;
    delete[] T_t; delete[] cp; delete[] k; delete[] mass;
}
// setter method
void Particles::set_properties(std::vector<float3_t> pos_vec, std::vector<float_t> T_vec,
                               std::vector<float_t> cp_vec, std::vector<float_t> k_vec,
                               std::vector<float_t> rho_vec, std::vector<float_t> h_vec,
                               std::vector<float_t> mass_vec)
{
    for (size_t i = 0; i < N; i++)
    {
        pos[i] = pos_vec[i]; cp[i] = cp_vec[i];
        k[i] = k_vec[i]; rho[i] = rho_vec[i];
        T[i] = T_vec[i]; T_t[i] = 0.;
        h[i] = h_vec[i]; mass[i] = mass_vec[i];
    }
}

Initializing The Problem

After creating the particle struct, we’ll start setting up the problem in the main function. The first step will be to create a results directory and clear it if any files already exist inside.

// main.cpp file
// inside the main function
// make results directoy if not present
struct stat st = {0};
if (stat("../results", &st) == -1) mkdir("../results", 0777);
// clear files from result directory
int ret;
ret = system("rm ../results/*");

Next, we will define the cube dimensions and discretization size, which will allow us to calculate the number of particles in the simulation. Additionally, we will specify the smoothing length factor.

// main.cpp file
// cube size
float_t lx = 0.1; float_t ly = 0.1; float_t lz = 0.1;
// Space between particles
float_t dx = .005; float_t dy = dx; float_t dz = dx;
// number of particles in each direction
int nx = lx / dx; int ny = ly / dy; int nz = lz / dz;
// total number of particles
int N = nx * ny * nz;
// smoothing length factor
float_t hf = 1.2;

Following that, we will allocate an instance of particles on the heap and call the initialize_particles function. This function will generate particles based on the specified dimensions and properties, storing them in the provided particles object. Finally, the run_simulation method will be called

// main.cpp
// allocating an instance of particles and initializing them
Particles *particles = new Particles(N);
initialize_particles(particles, nx, ny, nz, dx, dy, dz, hf);
run_simulation(particles, total_time, dt); // run simulation
// main.cpp
void initialize_particles(Particles *particles, int nx, int ny, int nz,
                          float_t dx, float_t dy, float_t dz, float_t hf)
{
    float_t thermal_cp = 850.;           // units: J/kgK
    float_t thermal_k = 157.;            // units: W/mK
    float_t rho = 2830.;                 // units: kg/m^3
    float_t mass = 2830. * dx * dy * dz; // units: kg
    float_t h = hf * dx;                 // units: m
    float_t T = 20.;                     // units: C (initial temperature)
    float_t T_boundary = 500.;           // units: C  (boundary temperature)

    // particles containers
    std::vector<float3_t> pos_vec;
    std::vector<float_t> T_vec;
    std::vector<float_t> cp_vec;
    std::vector<float_t> k_vec;
    std::vector<float_t> rho_vec;
    std::vector<float_t> h_vec;
    std::vector<float_t> mass_vec;

    // create particles
    for (size_t i = 0; i < nx; i++)
    {
        for (size_t j = 0; j < ny; j++)
        {
            for (size_t k = 0; k < nz; k++)
            {
                pos_vec.push_back(make_float3_t(i * dx, j * dx, k * dx)); // set position of particles
                cp_vec.push_back(thermal_cp);
                k_vec.push_back(thermal_k);
                rho_vec.push_back(rho);
                mass_vec.push_back(mass);
                h_vec.push_back(h);
                T_vec.push_back(i == 0 ? T_boundary : T); // set boundary temperature and initial temperature
            }
        }
    }

    particles->set_properties(pos_vec, T_vec, cp_vec, k_vec, rho_vec, h_vec, mass_vec); // set properties of particles
}

Running The Simulation

The run_simulation function will be called from the main function after the initialize_particles function.

  • This function starts by saving the initial state of the problem in VTK file format, enabling visualization in the open-source program ParaView.
  • It then iterates through the entire simulation duration, calculating the time derivative of the particles’ temperatures and updating their temperatures in each iteration.
  • Additionally, about every 1% of completion, it will save the current state in a VTK file.
// main.cpp file
void run_simulation(Particles *particles, float_t total_time, float_t dt)
{
    // for tracking the progress of the simulation
    float_t percent_increment = 0.01;
    float_t next_percent = percent_increment;
    unsigned int step = 0;

    // write initial state
    vtk_writer_write(particles, step);

    for (float_t t = 0; t < total_time; t += dt)
    {
        // calculate heat transfer between particles T_t
        interactions_heat(particles);
        // update temperature T of particles
        update(particles, dt);

        // tracking progress
        float_t current_percent = t / total_time;
        if (current_percent >= next_percent)
        {
            // printing on terminal the percentage of Completion
            std::cout << "Completion: " << int(current_percent * 100) << "%" << std::endl;
            next_percent += percent_increment;

            // write vtk file
            vtk_writer_write(particles, step);
        }
        step++;
    }
}

VTK Files and Visualization

To gain a better understanding of the vtk_writer_write function, you should review the VTK file format documentation. However, it’s relatively straightforward, and you can use the following function:

// vtk_writer.cpp file
void vtk_writer_write(Particles *particles, unsigned int step)
{
    // number of particles in the simulation
    unsigned int N_particles = particles->N;

    char buf[256];
    // output file name
    sprintf(buf, "../results/output_%012d.vtk", step);
    FILE *fp = fopen(buf, "w+");

    // write vtk file header
    fprintf(fp, "# vtk DataFile Version 2.0\n");
    fprintf(fp, "mfree iwf\n");
    fprintf(fp, "ASCII\n");
    fprintf(fp, "\n");

    fprintf(fp, "DATASET UNSTRUCTURED_GRID\n");
    fprintf(fp, "POINTS %d float\n", N_particles);
    for (unsigned int i = 0; i < N_particles; i++)
    {
        fprintf(fp, "%f %f %f\n", particles->pos[i].x,
                particles->pos[i].y, particles->pos[i].z);
    }
    fprintf(fp, "\n");

    fprintf(fp, "CELLS %d %d\n", N_particles, 2 * N_particles);
    for (unsigned int i = 0; i < N_particles; i++)
    {

        fprintf(fp, "%d %d\n", 1, i);
    }
    fprintf(fp, "\n");

    fprintf(fp, "CELL_TYPES %d\n", N_particles);
    for (unsigned int i = 0; i < N_particles; i++)
    {

        fprintf(fp, "%d\n", 1);
    }
    fprintf(fp, "\n");

    fprintf(fp, "POINT_DATA %d\n", N_particles);

    /************** Here temperature is printed****************/
    // Current particle temperature
    fprintf(fp, "SCALARS Temperature float 1\n");
    fprintf(fp, "LOOKUP_TABLE default\n");
    for (unsigned int i = 0; i < N_particles; i++)
    {
        fprintf(fp, "%f\n", particles->T[i]);
    }
    fprintf(fp, "\n");
    /****Use the above code block to save any other property*****/

    fclose(fp);
}

The function above will save the simulation state, i.e. particles positions and temperatures in VTK file format. You can then open these files using the ParaView program for visualization.

Press enter or click to view image in full size

Interaction between particles

Next is the interactions_heat function, which will be called during each iteration of the run_simulation function to compute the time derivative of the temperature for every particle in the simulation.

// interaction.cpp
void interactions_heat(Particles *particles)
{
    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

        // iterate over neighboring particles for particle i
        for (size_t j = 0; j < N; 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;                       // Kernel value
            float_t w_x = ww.y;                     // Kernel gradient
            float_t w_y = ww.z;                     // Kernel gradient
            float_t w_z = ww.w;                     // Kernel 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 *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
    }
}

The function above iterates through all particles to calculate the time derivative of their temperatures.

  • It starts by retrieving each particle’s position, temperature, thermal conductivity, and other properties.
  • For each particle, an additional iteration is performed to find its neighboring particles and gather their position and property values.
  • The interaction effect of the neighboring particles is then calculated. This effect is multiplied by a weight, which determines the influence strength of each neighbor.
  • The weight increases as the neighbor particle gets closer and decreases as it moves farther away. This weight is computed in the code using a function called cubic_spline, as shown below.
Press enter or click to view image in full size
SPH thermal interaction illustration

The cubic_spline function is defined as follows


// kernel.cpp file
// This function calculates the cubic spline weight and its derivatives between two particles.
// It takes the position of particle i (posi), the position of particle j (posj), and the smoothing length (hi) as inputs.
// It returns a float4_t variable (w) that contains the weight and its derivatives.
float4_t cubic_spline(float3_t posi, float3_t posj, float_t hi) {
    float4_t w = {0.0, 0.0, 0.0, 0.0};

    // Calculate the difference between the positions of particle i and particle j
    float3_t diff = {posi.x - posj.x, posi.y - posj.y, posi.z - posj.z};
    float_t rr2 = diff.x * diff.x + diff.y * diff.y + diff.z * diff.z;
    float_t h1 = 1 / hi;
    float_t fourh2 = 4 * hi * hi;

    // Check if the particles are too far apart or too close
    if (rr2 >= fourh2 || rr2 < 1e-8) {
        return w;
    }

    float_t rad = sqrtf(rr2);
    float_t q = rad * h1;
    float_t fac = (M_1_PI) * h1 * h1 * h1;

    // Calculate the weight and its derivatives based on the value of q
    if (q > 1) {
        float_t _2mq = 2 - q;
        float_t _2mq2 = _2mq * _2mq;
        float_t val = 0.25 * _2mq2 * _2mq;
        float_t der = -0.75 * _2mq2 * h1 / rad;

        // Assign the calculated values to the components of the float4_t variable w
        w.x = val * fac;         // weight
        w.y = der * diff.x * fac;    // derivative with respect to x
        w.z = der * diff.y * fac;    // derivative with respect to y
        w.w = der * diff.z * fac;    // derivative with respect to z
    } else {
        float_t val = 1 - 1.5 * q * q * (1 - 0.5 * q);
        float_t der = -3.0 * q * (1 - 0.75 * q) * h1 / rad;

        // Assign the calculated values to the components of the float4_t variable w
        w.x = val * fac;         // weight
        w.y = der * diff.x * fac;    // derivative with respect to x
        w.z = der * diff.y * fac;    // derivative with respect to y
        w.w = der * diff.z * fac;    // derivative with respect to z
    }

    return w;
}

Updating problem state

Finally, in each iteration of the run_simulation function, the temperatures of the particles are updated as follows:

// main.cpp file
void update(Particles *particles, float_t dt)
{
    for (size_t i = 0; i < particles->N; i++)
    {
        // update temperature of particles
        particles->T[i] += particles->T_t[i] * dt;
    }
}

Project Building

let’s create first a CMakeList.txt file

# CMakeLists file
cmake_minimum_required(VERSION 3.12)
project(SPH_thermal_CPU)

# Set the C++ standard
set(CMAKE_CXX_STANDARD 17)

# Set the include directories
include_directories(include)

# Add your source files here
file(GLOB SOURCES "src/*.cpp")

# Create the executable
add_executable(${PROJECT_NAME} ${SOURCES})

To build and compile the project

## Building the project

1. Open the terminal and navigate to the project directory:
```
cd SPH_thermal_CPU
```
2. Create a `build` directory, Change to the `build` directory:
```
mkdir build && cd build
```
4. Build and compile the project:
```
cmake .. && make
```
5. Run it 

```
./SPH_thermal_CPU
```

Expected output

  • You can track the completion percentage on the terminal.
  • The output files will be saved in the results folder within the project directory and can be dragged and dropped into the ParaView program for visualization.

Let’s Recap

  • We start by defining our own generic data type, which allows us to switch between single and double precision.
  • Next, we initialize our problem by creating particles with their properties.
  • We then develop the `run_simulation` function, which prints the initial state of the problem and iterates through the simulation time. In each iteration, the function:
    - Calculates the interactions between particles
    - Updates the particles’ temperatures
    - Saves the current state of the problem every 1% of completion

What is next

Running the current code on an 11th Gen Intel® Core™ i7–11800H @ 2.30GHz with 16 GB of RAM takes ~ 24 minutes. However, there is potential to optimize the code for improved performance. Let’s examine the interactions_heat method to explore possible enhancements.

Press enter or click to view image in full size
SPH interaction function
  • Interaction function time complexity is O(N²).
  • We can introduce hashing to speed it up
  • What is hashing?!
  • This will be discussed in the next article: SPH and Hashing.

Acknowledgment

Useful Links

Sph
Programming
Simulation Software
Thermal Conductivity
C Plus Plus Language

--

--

Ahmed Elbossily

Written by Ahmed Elbossily

3 followers
1 following

Research Assistant & Ph.D. Candidate | GPU Computing & Computational Mechanics | ahmedelbossily.github.io

Help

Status

About

Careers

Press

Blog

Store

Privacy

Rules

Terms

Text to speech