This Repository contains the code for the N-Body project from the UZH course on Computational Astrophysics.
You'll need Eigen for the code to compile. If needed change the EIGEN variable in the makefile to the path of your eigen source files. You'll need MathGL 2.0 or higher to plot the data of the first task and Python 3 or higher to plot the data of the second task.
Simply compile the code with make and run it with ./main.out. The plots for the first task are generated directly. The second task produces output files with data. To plot that data use python3 plot.py. The code will perform the necessary tasks for the project but the solvers are fully encapsulated and can be used on their own.
Computes the n-body forces from a set of initial conditions with the direct summation algorithm in O(n²) where n is the number of particles. The initial conditions can be read in two ways
Reads a set of initial conditions in the format according to the data provided by the project. The data file is made of a header and sequential arrays, as follows: Number of particles(N), Number of gas particles(NG), Number of star particles (NS).
N NG NS
Masses[i]
x[i]
y[i]
z[i]
Vx[i]
Vy[i]
Vz[i]
Note that the data first runns through all N masses then through all N x coordinates, all N y coordinates, all N z coordinates etc. The function returns true when successful and false otherwise. In the case of an unsuccessful read the console should give some output.
Reads a header in the format N eps where N is the number of particles and eps is the softening. It then reads the data particle by particle so mass, x, y, z, vx, vy, vz of particle i then mass, x, y, z, vx, vy, vz of particle i+1 and so on. The function returns true when successful and false otherwise. In the case of an unsuccessful read the console should give some output.
Writes the particles to a file in the same format as being used by the readData() function above. The function returns true when successfully writing the file and false otherwise. In the case of an unsuccessful write the console should give some output.
Returns the currently used softening for the solver.
Sets the softening value used by the solver to the function's input.
Computes the n-body forces acting on each particle with the direct summation algorithm.
Returns the current state (mass, position, velocity, force) of all particles in the form of a 10xN Eigen Matrix Matrix<float, 10, Dynamic> where the i-th column holds the state of the i-th particle in the form of [mass, x, y, z, vx, vz, vz, fx, fy, fz]. Note that until either readDataOld() or readData() have been called the output of the function remains uninitialized.
- Note this solver has not been tested for correctness and should therefore not be used! Right now it uses the unsupported Eigen Tensor for the mesh which is not nearly as efficient as the Eigen matrices.
Computes the n-body forces from a set of initial conditions by constructing a particle mesh whose size in one dimension is specified in the constructor. It computes the forces by solving the poisson equation with a fast fourier transform. It scales as O(n³log(n)) where n is the number of mesh cells along one dimension. The interface is identical to the direct solver class except that the methods to get and set the softening are missing.
Reads a set of initial conditions, evolves the system for a specified amount of time and writes the state of the system and the relative energy error over the entire time (for each timestep) to ouput files. One file for the particle positions and one file for the energy error.
Reads a set of initial conditions from a file with the same format as specified in Gravitysolver::Direct::readDataOld. The very last value of the file can hold the softening if you whish to do so.
Sets the softening for the evolution to the function's input value
Only writes the energy error to a file. Can be useful if the file for the positions becomes extremely large and one want's to avoid writing it.
Inverse of enableLean(). If calling the simulation function now the particles positions are written to a file again.
Only keeps track of every k-th iteration in the simulation where k is the function's input. This allows to save a lot of memory when usign very small step sizes.
Calling integrate(dt, numSteps) evolves the system with timestep dt for numSteps steps and (depending on the lean mode) writes the result of the simulation to a file.
Returns the state of the system over the entire simulation time. Each column contains the [x, y, z] position of the particles. The columns 0 to n-1 contain the state of n particles at the first time step, columns n to 2n-1 contain the state of n particles at the next (or a multiple of next depending on the block size) time step and so on.
#include <gravitysolvers.hpp>
#include <hermite.hpp>
int main()
{
Gravitysolver::Direct *direct = new Gravitysolver::Direct();
direct->readData("data.txt");
direct->setSoftening(0.001);
direct->solve();
direct->writeData("forces.txt");
Hermite *hermite = new Hermite();
hermite->readData("data.txt");
hermite->setSoftening(0.001);
hermite->setBlockSize(20);
hermite->enableLean();
hermite->integrate(0.01, 1000);
return 0;
}