Commit 5302ee67 authored by Ivo Roghair's avatar Ivo Roghair
Browse files

Initial commit

parent f33f3833
Pipeline #3576 canceled with stages
all:
g++ -o summation -O3 -fopenmp -std=c++14 main.cpp
# Summations
# Neighborlists
## Intro
It has been found that a parallel implementation of the DPM did not provide
the same (exact) results as its serial counterpart. We suspected that the
addition of finite precision floating point numbers would be the root cause.
The (serial) DPM model loops over all the Lagrangian particles in sequence
to account for their contribution to various field variables, e.g. porosity,
pressure, momentum exchange etc. The parallel code adds these contributions
in a different order, and this accounts for a numerical error that perturbs
the simulation so that the results are not exactly identical compared to the
serial code. As a matter of fact, multiple runs of the same simulation may
yield slightly different results due to different scheduling.
## Components
* Python script `generateBinaryPositions.py`: Generates a binary file that
contains the number of particles (1 integer) and consequently their (x,y,z)
positions.
* A main file `main.cpp` that loads the binary file into an N x 3 array, and
proceeds to sum all the positions to a single summation:
```math
s_1 = \sum_n=0^n=\mathrm{npart} x_n + y_n + z_n
```
It accounts for a summation using OpenMP parallelization, a serial summation
(counting from 0 up to npart), another serial (counting down from npart to 0)
and finally a sum that accounts for minimizing the numerical error using the
[Kahan summation algorithm](https://en.wikipedia.org/wiki/Kahan_summation_algorithm)
## Sample output
The program loads an input file given as argument to the program:
```
$ ./summation example.bin
Trying to open file example.bin ...
Trying to read 10000000 positions from file...
Four different summations:
OpenMP: 1.4999288710279200e+07 -- Serial 0->nparts: 1.4999288710278464e+07 -- Serial nparts-->0: 1.4999288710281001e+07 -- Kahan: 1.4999288710279370e+07
Exiting normally...
```
## To do
File added
#!/usr/bin/env python
import numpy as np
from array import array
npart = np.intc(10000000)
fname = 'example.bin'
print ("Generating {} particle positions...".format(npart))
a = np.random.rand(npart, 3)
np.set_printoptions(precision=16)
print(a)
print ("Writing to file ....")
output_file = open(fname, 'wb')
output_file.write(npart.tobytes())
output_file.write(a.tobytes())
output_file.close()
print ("Done!")
#include <iostream>
#include <fstream>
#include <vector>
#include <limits>
#include <chrono>
#include <omp.h>
typedef std::chrono::high_resolution_clock Clock;
using namespace std;
int main(int argc, char* argv[])
{
int npart;
double** pos;
cout.precision(numeric_limits<double>::max_digits10);
cout << "Trying to open file " << argv[1] << " ..." << endl;
ifstream file(argv[1], ios::binary);
if (file.is_open())
{
// Read number of particles
file.read((char*) &npart, sizeof(int));
cout << "Trying to read " << npart << " positions from file..." << endl;
// Allocate memory
pos = new double*[npart];
for (int i = 0; i < npart; i++)
pos[i] = new double[3];
// Read all particle positions into array
for (int i = 0; i < npart; i++) {
file.read((char*) pos[i], sizeof(double[3]));
}
file.close();
double c = 0.0;
double sum1 = 0.0, sum2 = 0.0;
double sum3 = 0.0, sum4 = 0.0;
double y, t;
#pragma omp parallel for schedule(dynamic) reduction(+:sum1)
for (int i=0;i<npart;i++) {
sum1 +=pos[i][0]+pos[i][1]+pos[i][2];
}
for (int i=0;i<npart;i++) {
sum2 +=pos[i][0]+pos[i][1]+pos[i][2];
}
for (int i=npart-1;i>=0;i--) {
sum3 +=pos[i][0]+pos[i][1]+pos[i][2];
}
for (int i=0;i<npart;i++) {
y = pos[i][0]+pos[i][1]+pos[i][2] - c;
t = sum4+y;
c = (t - sum4) - y;
sum4 = t;
}
printf("Four different summations:\n");
printf("OpenMP: %1.16e -- Serial 0->nparts: %1.16e -- Serial nparts-->0: %1.16e -- Kahan: %1.16e\n", sum1,sum2,sum3,sum4);
}
else cout << "Error opening file " << argv[1] << ".\n\n" << endl;
cout << "Exiting normally..." << endl;
}
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment