Code overview
The code in this repository contains a serial C++ implementation of the SPH methodology described earlier. The variables associated with the particles' positions, velocities and forces, are stored as members of an object called Fluid, while the methods of another class called SphSolver perform the steps of the algorithm.
Files
The src directory comprises five *.cpp files and their corresponding header (*.h) files, as well as the files to build the code.
src/SPH-main.cpp: Contains the source code of the main program which is responsible for instrumenting the execution of the SPH simulation.src/initial_conditions.{h, cpp}: Contains the functions which are required in order to create the initial conditions.src/main_prog_func.h: Header file for the declarations of the main program functions.src/particles.{h, cpp}: Contains a class describing the positions and velocities of particles.src/fluid.{h, cpp}: Contains a class which extends the functionality described in thesrc/particles.{h, cpp}to represent a cluster of particles which form a fluid.src/sph_solver.{h, cpp}: Contains a class which manifests the implementation of the SPH methodology described in this project.src/CMakeLists.txt: Contains the instructions to compile the code.
Reading Inputs
The program expects the parameters which are specified in the exec/input directory. When the function initialise() is invoked by the main program to read the .txt files, the <boost/program_options.hpp> library is utilised. This library facilitates the mapping of these parameters to their corresponding values, which are then stored in their respective variable maps in the retrieveInputsFromFile() function. This approach enhances the flexibility and robustness of the input reading process. Users can specify input parameters in the .txt files in any order, provided they are presented as key = value pairs.
/* **************************** SPH_main.cpp **************************** */
void initialise(std::unique_ptr<Fluid>& fluidPtr, SphSolver& sphSolver) {
// Process to obtain the inputs provided by the user
po::options_description desc("Allowed options");
desc.add_options()("init_condition", po::value<std::string>(),
"take an initial condition")("T", po::value<double>(),
"take integration time")(
"dt", po::value<double>(), "take time-step")(
"coeffCfl1", po::value<double>(), "take lamda nu")(
"coeffCfl2", po::value<double>(), "take lamda f")(
"adaptive_timestep", po::value<bool>(),
"take flag for adaptive time-step")("h", po::value<double>(),
"take radius of influence")(
"gas_constant", po::value<double>(), "take gas constant")(
"density_resting", po::value<double>(), "take resting density")(
"viscosity", po::value<double>(), "take viscosity")(
"acceleration_gravity", po::value<double>(), "take acc due to gravity")(
"coeff_restitution", po::value<double>(), "take coeff of restitution")(
"left_wall", po::value<double>(), "take left wall position")(
"right_wall", po::value<double>(), "take right wall position")(
"bottom_wall", po::value<double>(), "take bottom wall position")(
"top_wall", po::value<double>(), "take top wall position")(
"length", po::value<double>(), "take length of the block")(
"width", po::value<double>(), "take width of the block")(
"radius", po::value<double>(), "take radius of the droplet")(
"n", po::value<int>(), "take number of particles")(
"center_x", po::value<double>(), "take center of the particle mass in x")(
"center_y", po::value<double>(), "take center of the particle mass in y")(
"init_x_1", po::value<double>(), "take x_1")(
"init_y_1", po::value<double>(), "take y_2")(
"init_x_2", po::value<double>(), "take x_2")(
"init_y_2", po::value<double>(), "take y_2")(
"init_x_3", po::value<double>(), "take x_3")(
"init_y_3", po::value<double>(), "take y_3")(
"init_x_4", po::value<double>(), "take x_4")(
"init_y_4", po::value<double>(), "take y_4")(
"output_frequency", po::value<int>(),
"take frequency that output will be written to file");
...
icCase = caseVm["init_condition"].as<std::string>();
std::string fileName = icCase + ".txt";
retrieveInputsFromFile(fileName, icCase, desc, icVm);
handleInputErrors(caseVm, domainVm, constantsVm, icVm);
...
}
void retrieveInputsFromFile(const std::string& fileName,
const std::string& icCase,
const po::options_description& desc,
po::variables_map& vm) {
std::ifstream caseFile;
std::string errorMessage = "Error opening file: " + fileName;
if (fileName == icCase + ".txt") {
errorMessage +=
" Make sure that the value of the init_condition in the case.txt "
"file is one of the following: ic-one-particle, ic-two-particles, "
"ic-three-particles, ic-four-particles, ic-droplet, ic-block-drop.";
}
// Try to open the file
try {
caseFile.open("../input/" + fileName);
// Throw an exception if the file cannot be opened
if (!caseFile.is_open()) {
throw std::runtime_error(errorMessage);
}
po::store(po::parse_config_file(caseFile, desc), vm);
} catch (std::runtime_error& e) {
// Handle the exception by printing the error message and exiting the
// program
std::cerr << e.what() << std::endl;
exit(1);
}
po::notify(vm);
}
Additionally, using the handleInputErrors() function, error handling is integrated into the input file reading process to guarantee that the provided values conform to the constraints imposed by the underlying mathematical models and the physical meaning of each variable. For example, if the user attempts to set a negative value for the timestep, or a value that is greater than the integration time, the program will throw an error and instruct the user to choose a more suitable value.
/* **************************** SPH_main.cpp **************************** */
void handleInputErrors(const po::variables_map& caseVm,
const po::variables_map& domainVm,
const po::variables_map& constantsVm,
const po::variables_map& icVm) {
try {
// Error handling for the total integration time
if (caseVm["T"].as<double>() <= 0) {
throw std::runtime_error(
"Error: Total integration time must be positive!");
// Error handling for the time step
} else if (caseVm["dt"].as<double>() <= 0 or
caseVm["dt"].as<double>() > caseVm["T"].as<double>()) {
throw std::runtime_error(
"Error: Time step must be positive and lower than the total "
"integration time!");
// Error handling for the output frequency
} else if (caseVm["output_frequency"].as<int>() <= 0 or
caseVm["output_frequency"].as<int>() >
ceil(caseVm["T"].as<double>() / caseVm["dt"].as<double>())) {
throw std::runtime_error(
"Error: Output frequency must be positive and lower than the total "
"number of iterations!");
// Error handling for the CFL coefficients
} else if (caseVm["coeffCfl1"].as<double>() <= 0 or
caseVm["coeffCfl1"].as<double>() >= 1 or
caseVm["coeffCfl2"].as<double>() <= 0 or
caseVm["coeffCfl2"].as<double>() >= 1) {
throw std::runtime_error(
"Error: The CFL coefficients must be positive and less than 1");
// Error handling for the domain boundaries
} else if (domainVm["left_wall"].as<double>() >=
domainVm["right_wall"].as<double>() ||
domainVm["bottom_wall"].as<double>() >=
domainVm["top_wall"].as<double>()) {
throw std::runtime_error(
"Error: Please adjust your domain boundaries so that left_wall < "
"right wall and bottom_wall < top_wall.");
// Error handling for the number of particles
} else if (icVm["n"].as<int>() <= 0) {
throw std::runtime_error("Error: Number of particles must be positive!");
}
} catch (std::runtime_error& e) {
// Handle the exception by printing the error message and exiting the
// program
std::cerr << e.what() << std::endl;
exit(1);
}
}
Class Initialisation
The code makes use of three different classes which are purposed to represent the fluid and the SPH algorithm deployed in this project.
Firstly, one SphSolver object and one pointer to a Fluid object are declared in the main program. The pointer declaration is used for the Fluid, because to initialise the object properly the number of particles is required in the user defined constructor and this information is not yet available since the input files have not been read. These objects are passed as a reference to the initialise() function.
Once the input values are read and stored, the provided initial condition (IC) is used to determine the number of particles. This means that although the user has already provided a number of particles, this is just an indication, since the IC (droplet and block drop) require specific formation and the particles to be distributed uniformly. These two conditions cannot be satisfied simultaneously by any number of particles and therefore several adjustments need to be made. The functions closestIntegerSqrt() and rectangleN() from initial_conditions.h are functions suitable for this purpose.
The IC functions are called within the setInitialConditions() function and a reference to the pointer of the fluid object is passed as an argument*, as well as the updated number of particles. Inside these functions the user defined constructor of the Fluid class is called and the memory allocation process for the object's containers is invoked. In this, the containers are declared as new raw pointers to arrays, dynamically allocating memory proportional to the number of particles. The function used to initialise the Fluid class for the simple cases of 1,2,3 and 4 particles is demonstrated below.
* The rationale behind passing the pointer to the fluid object as a reference is identical to the reason an object is passed as a reference to a function. In these functions we are allocating new memory that the pointer should point to. If the pointer was passed by value (Fluid *fluidPtr) then a copy of the pointer would be pointing to the new memory and our original fluid pointer will still be a nullptr. Of course, since the allocation is happening inside the function, the caller is responsible to delete the object manually.
/* **************************** initial_conditions.cpp **************************** */
void icBasic(std::unique_ptr<Fluid> &fluidPtr, int nbParticles,
std::vector<double> &positionX, std::vector<double> &positionY) {
fluidPtr = std::make_unique<Fluid>(nbParticles);
for (int i = 0; i < nbParticles; i++) {
fluid.setPositionX(i, positionX[i]);
fluid.setPositionY(i, positionY[i]);
fluid.setVelocityX(i, 0.0);
fluid.setVelocityY(i, 0.0);
}
}
Inside the setInitialConditions() function, an std::map object is used to map the initial condition case with a fixed number of particles to their corresponding number of particles, to avoid the use of multiple if statements and make the code cleaner. Then, the initial condition case is set based on the user's input, including additional error handling for the input parameters specific to the selected IC. The workflow for the droplet case is demonstrated below.
/* **************************** SPH_main.cpp **************************** */
void setInitialConditions(const std::string& icCase,
std::unique_ptr<Fluid>& fluidPtr,
const po::variables_map& icVm,
const po::variables_map& domainVm) {
// Fixed nbParticles ic cases map
std::map<std::string, int> initConditionToParticlesMap = {
{"ic-one-particle", 1},
{"ic-two-particles", 2},
{"ic-three-particles", 3},
{"ic-four-particles", 4}};
// Get the number of particles based on the ic case (for the more complex ic)
int nbParticles;
if (icCase == "ic-droplet" || icCase == "ic-block-drop") {
nbParticles = icVm["n"].as<int>();
} else {
nbParticles = initConditionToParticlesMap[icCase];
}
// Block of code ...
} else if (icCase == "ic-droplet") {
// Get the droplet radius and center coordinates from the ic file
double radius = icVm["radius"].as<double>();
// Error handling for the droplet radius
try {
if (radius <= 0) {
throw std::runtime_error("Error: Radius must be positive!");
}
} catch (std::runtime_error& e) {
// Handle the exception by printing the error message and exiting the
// program
std::cerr << e.what() << std::endl;
exit(1);
}
double centerX = icVm["center_x"].as<double>();
double centerY = icVm["center_y"].as<double>();
// Error handling for the droplet initial position (centerX, centerY)
try {
if (centerX - radius < domainVm["left_wall"].as<double>() ||
centerX + radius > domainVm["right_wall"].as<double>() ||
centerY - radius < domainVm["bottom_wall"].as<double>() ||
centerY + radius > domainVm["top_wall"].as<double>()) {
throw std::runtime_error(
"Error: The droplet must be within the domain boundaries! Please "
"adjust the center coordinates.");
}
} catch (std::runtime_error& e) {
// Handle the exception by printing the error message and exiting the
// program
std::cerr << e.what() << std::endl;
exit(1);
}
icDroplet(fluidPtr, nbParticles, radius, centerX, centerY);
} else {
std::cerr << "Error: Initial condition function not found! Make sure "
<< "that the value of the init_condition in the case.txt file is "
<< "one of the following: ic-one-particle, ic-two-particles, "
<< "ic-three-particles, ic-four-particles, ic-droplet, "
<< "ic-block-drop." << std::endl;
exit(1);
}
}
Finally, after the object initialisation, the rest of the parameters which are required by the SphSolver and the Fluid objects are being set with the use of setter functions, the Neighbour Searching grid is created, the initialised particles are placed in their corresponding grid cells, and their neighbours are identified.
/* **************************** SPH_main.cpp **************************** */
// void initialise(fluid** fluidPtr, SphSolver& sphSolver) { ...
// Set the parameters of the solver for the specific simulation
sphSolver.setTimestep(caseVm["dt"].as<double>());
sphSolver.setAdaptiveTimestep(caseVm["adaptive_timestep"].as<bool>());
sphSolver.setCflCoefficients(caseVm["coeffCfl1"].as<double>(),
caseVm["coeffCfl2"].as<double>());
sphSolver.setTotalTime(caseVm["T"].as<double>());
sphSolver.setOutputFrequency(caseVm["output_frequency"].as<int>());
sphSolver.setCoeffRestitution(constantsVm["coeff_restitution"].as<double>());
sphSolver.setLeftWall(domainVm["left_wall"].as<double>());
sphSolver.setRightWall(domainVm["right_wall"].as<double>());
sphSolver.setTopWall(domainVm["top_wall"].as<double>());
sphSolver.setBottomWall(domainVm["bottom_wall"].as<double>());
sphSolver.setPrecalculatedValues(constantsVm["h"].as<double>());
// Define the fluid based on the inputs
fluidPtr->setRadInfl(constantsVm["h"].as<double>());
fluidPtr->setGasConstant(constantsVm["gas_constant"].as<double>());
fluidPtr->setDensityResting(constantsVm["density_resting"].as<double>());
fluidPtr->setViscosity(constantsVm["viscosity"].as<double>());
fluidPtr->setAccelerationGravity(
constantsVm["acceleration_gravity"].as<double>());
// Calculate the mass of the particles
sphSolver.createGrid(*fluidPtr);
sphSolver.neighbourParticlesSearch(*fluidPtr);
std::vector<std::vector<std::pair<int, double>>> neighboursPerParticle =
sphSolver.getNeighbourParticles();
fluidPtr->calculateMass(neighboursPerParticle);
Output files
The output files are initialised by the function initOutputFiles(), after the Fluid class initialisation and before the time integration procedure is invoked. This is done because during the time integration, the code will produce runtime outputs and therefore the corresponding files need to have been created beforehand. The outputs are exported in .csv format which displays good readability and facilitates data manipulation compared to .txt files. They are stored in a centralised location, specifically within the /exec/output/ directory. This centralisation simplifies data organisation and retrieval, making it easier for users to access and analyse output data.
Upon successful execution, the program generates two types of files:
-
Energies File: This file, containing Total, Kinetic, and Potential energies, is updated at each timestep. The results can be visualized by using the script
post/plot_energies.py. -
Particle Positions File: This file captures the positions of the particles at a given timestep, and it can be visualized using the scripts
post/visualize_particles.py.
/* **************************** SPH_main.cpp **************************** */
void storeToFile(Fluid &fluid, std::string type, std::ofstream &targetFile,
double currentIntegrationTime) {
if (type == "energy") {
// Write energies in the Energy-File
targetFile << currentIntegrationTime << "," << fluid.getKineticEnergy()
<< "," << fluid.getPotentialEnergy() << ","
<< fluid.getPotentialEnergy() + fluid.getKineticEnergy() << "\n";
} else if (type == "position") {
// Write positions in the position file
for (size_t k = 0; k < fluid.getNumberOfParticles(); k++) {
targetFile << currentIntegrationTime << "," << fluid.getPositionX(k)
<< "," << fluid.getPositionY(k) << "\n";
}
}
}
 |
|---|
| Energy plots (left) and initial position (right) for a droplet of 60 particles. |
Time integration
Following the initialisation of the class and the output files, the function SphSolver::timeIntegration() is invoked. Within this function, the steps of the SPH algorithm are executed, and the outputs are exported.
/* ***************************** SPH-main.cpp ****************************** */
// Time integration loop
sphSolver.timeIntegration(*sphFluid, simulationPositionsFile,
finalPositionsFile, energiesFile);
/* **************************** sph_solver.cpp **************************** */
void SphSolver::timeIntegration(Fluid &data,
std::ofstream &simulationPositionsFile,
std::ofstream &finalPositionsFile,
std::ofstream &energiesFile) {
std ::cout << "Time integration started -- OK"
<< "\n";
while (currentIntegrationTime < totalTime) {
if (adaptiveTimestepBool) {
// Reset the adaptive timestep related variables
maxVelocity = 0.0;
maxAcceleration = 0.0;
}
// In each iteration the distances between the particles are recalculated,
// as well as their density and pressure
neighbourParticlesSearch(data);
data.calculateDensity(neighbourParticles);
data.calculatePressure();
particleIterations(data);
currentIntegrationTime += dt;
if (t % outputFrequency == 0) {
storeToFile(data, "energy", energiesFile, currentIntegrationTime);
storeToFile(data, "position", simulationPositionsFile, currentIntegrationTime);
}
t++;
if (adaptiveTimestepBool) {
adaptiveTimestep(data);
}
}
// Store particles' positions after integration is completed
storeToFile(data, "position", finalPositionsFile, currentIntegrationTime);
std ::cout << "Time integration finished -- OK"
<< "\n";
}
Implementation of the SPH algorithm
The steps of the SPH algorithm are executed by the function SphSolver::timeIntegration() and are implemented as follows.
Density
The density of the fluid associated with each particle is approximated as:
/* **************************** fluid.cpp **************************** */
void Fluid::calculateDensity(
const std::vector<std::vector<std::pair<int, double>>>& neighbours) {
double phi, normalisedDistance, normalisedDistanceSqr;
for (size_t i = 0; i < nbParticles; i++) {
density[i] = mass * fourPih2;
for (size_t j = 0; j < neighbours[i].size(); j++) {
normalisedDistance = neighbours[i][j].second * hInverse;
normalisedDistanceSqr = (1.0 - normalisedDistance * normalisedDistance);
phi = fourPih2 * normalisedDistanceSqr * normalisedDistanceSqr *
normalisedDistanceSqr;
density[i] += mass * phi;
}
}
}
where:
massis the mass of a particle ().phiis the kernel density function for density ().normalisedDistanceis the distance between particle and particle , normalised by the interaction radius ().hInverseis the inverse of the interaction radius () which is constant throughout the simulation and is calculated once and stored in a variable.
Pressure
The pressure is calculated based on the ideal gas law
/* **************************** fluid.cpp **************************** */
void Fluid::calculatePressure() {
for (size_t i = 0; i < nbParticles; i++) {
pressure[i] = gasConstant * (density[i] - densityResting);
}
}
where:
densityRestingis the fluid's resting density ().gasConstantis the gas constant ().
Pressure force
The force exerted on the particle due to pressure from neighboring fluid particles is calculated as
/* **************************** sph_solver.cpp **************************** */
double SphSolver::calculatePressureForce(Fluid &data,
std::function<double(int)> getPosition,
int particleIndex) {
double sum = 0.0; // Initializing the summation
double normalisedDistance;
double position = getPosition(particleIndex);
double pressure = data.getPressure(particleIndex);
double mass = data.getMass();
double radiusOfInfluence = data.getRadInfl();
size_t neighbourIndex;
for (size_t j = 0; j < neighbourParticles[particleIndex].size(); j++) {
neighbourIndex = neighbourParticles[particleIndex][j].first;
if (particleIndex != neighbourIndex) {
try {
normalisedDistance =
neighbourParticles[particleIndex][j].second / radiusOfInfluence;
// Throw an exception if a singularity is about to appear
if (normalisedDistance == 0.0) {
throw std::runtime_error(
"A singularity appeared. Consider reducing the number of "
"particles or increasing the initial distance between the "
"particles.");
}
} catch (std::runtime_error &e) {
// Handle the exception by printing the error message and exiting the
// program
std::cerr << e.what() << std::endl;
exit(1);
}
sum += (mass / data.getDensity(neighbourIndex)) *
((pressure + data.getPressure(neighbourIndex)) / 2.0) *
(thirtyPih3 * (position - getPosition(neighbourIndex))) *
(((1.0 - normalisedDistance) * (1.0 - normalisedDistance)) /
normalisedDistance);
}
}
return -sum;
}
Where:
thirtyPih3is a constant precalculated value () which assists in the calculation of the nabla of the kernel density function for pressure .pressureis the pressure of the particle.data.getPressure(neighbourIndex)is used to retrieve the pressure of the neighbouring particle.data.getDensity(neighbourIndex)is used to retrieve the density of the neighbouring particle.
It is important to note that the function admits as an argument the std::function<double(int)> getPosition which is a pointer to a function of the Fluid class. This is done because the SphSolver::calculatePressureForce() function needs to be called twice for each particle - once for the x component and once for the y component of the pressure force.
forcePressureX = calculatePressureForce(data, ptrGetPositionX, i);
forcePressureY = calculatePressureForce(data, ptrGetPositionY, i);
In this way, the function for the calculation of pressure forces remains generic and direction agnostic which means that we do not have to duplicate any piece of code.
Also, in this function, an error handling occurs, to check whether a singularity appears if two particles come very close to each other. This would make the normalized distance nearly zero and NaN values would appear in the code. Herein we choose to inform the users about that incident and advise them to adjust the case parameters in order to avoid such behaviour.
Viscous force
The force acting on each particle due to viscous effects is calculated as
/* **************************** sph_solver.cpp **************************** */
double SphSolver::calcViscousForce(Fluid &data,
std::function<double(int)> getVelocity,
int particleIndex) {
double sum = 0.0; // Initializing the summation
double normalisedDistance;
double velocity = getVelocity(particleIndex);
double mass = data.getMass();
double radiusOfInfluence = data.getRadInfl();
size_t neighbourIndex;
for (size_t j = 0; j < neighbourParticles[particleIndex].size(); j++) {
neighbourIndex = neighbourParticles[particleIndex][j].first;
if (particleIndex != neighbourIndex) {
normalisedDistance =
neighbourParticles[particleIndex][j].second / radiusOfInfluence;
sum += (mass / data.getDensity(neighbourIndex)) *
(velocity - getVelocity(neighbourIndex)) *
(fourtyPih4 * (1.0 - normalisedDistance));
}
}
return -data.getViscosity() * sum;
}
Where:
fourtyPih4is a constant precalculated value () which assists in the calculation of the nabla of the kernel density function for the visocsity ().velocityis the velocity of the particle.getVelocity(neighbourIndex)is used to retrieve the velocity of the neighbouring particle.
For the reasons explained for SphSolver::calculatePressureForce(), a pointer to a function is required as an argument.
Gravity force:
Finally, the force due to gravity is calculated as:
/* **************************** sph_solver.cpp **************************** */
double SphSolver::calcGravityForce(Fluid &data, int particleIndex) {
return -data.getDensity(particleIndex) * data.getAccelerationGravity();
}
Where:
data.getAccelerationGravity()is used to retrieve the value of the acceleration of gravity.
Acceleration
The acceleration of each particle is calculated as:
/* **************************** sph_solver.cpp **************************** */
double acceleration = (forcePressure + forceViscous + forceGravity) /
data.getDensity(particleIndex);
Time marching
We solve the equation as a function of time by finding the velocity and position of each particle at each of a number of time steps. We denote a property of particle at time step as . The state of the property half way between time steps and is denoted as .
We begin with the initial conditions of the system, which are the positions and velocities of the particles at time . We iteratively use the state of the system at time step to find the state of the system at time step using a leap-frog scheme, which provides improved stability characteristics. For the x-directions we do the following:
/* **************************** sph_solver.cpp **************************** */
// x-direction
newVelocity = data.getVelocityX(particleIndex) +
integrationCoeff *
velocityIntegration(data, particleIndex, forcePressureX,
forceViscousX, forceGravityX);
data.setVelocityX(particleIndex, newVelocity);
newPosition = data.getPositionX(particleIndex) + newVelocity * dt;
data.setPositionX(particleIndex, newPosition);
We normally use integrationCoeff=1.0, but because the velocity is calculated at half-steps, we need to initialise the scheme on the first time step using integrationCoeff=0.5:
/* **************************** sph_solver.cpp **************************** */
// First step to initialise the scheme
if (t == 0) {
integrationCoeff = 0.5;
}