To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit 3bca17d1 authored by ursular's avatar ursular
Browse files

solution ex10

introduction ex11
parent 1bec2905
#ifndef ARRAYOFPARTICLE_H
#define ARRAYOFPARTICLE_H
#include "common.h"
class ArrayOfParticles
{
public:
int Np;
double * x;
double * y;
double * u;
double * v;
double * gamma;
ArrayOfParticles(int Np) : Np(Np)
{
x = new double[Np];
y = new double[Np];
u = new double[Np];
v = new double[Np];
gamma = new double[Np];
for (int i=0; i<Np; i++)
{
x[i] = 0.;
y[i] = 0.;
u[i] = 0.;
v[i] = 0.;
gamma[i] = 0.;
}
}
~ArrayOfParticles()
{
delete [] x;
delete [] y;
delete [] u;
delete [] v;
delete [] gamma;
}
ArrayOfParticles& operator=(ArrayOfParticles &p)
{
Np = p.Np;
for (int i=0; i<Np; i++)
{
x[i] = p.x[i];
y[i] = p.y[i];
u[i] = p.u[i];
v[i] = p.v[i];
gamma[i] = p.gamma[i];
}
return *this;
}
void AdvectEuler(double dt)
{
for (int i=0; i<Np; i++)
{
x[i] += u[i]*dt;
y[i] += v[i]*dt;
}
}
void ClearVelocities()
{
for (int i=0; i<Np; i++)
{
u[i] = 0.;
v[i] = 0.;
}
}
};
#endif /* ARRAYOFPARTICLE_H */
CC = mpic++
.PHONY: all
all: sheet
sheet: mainSheet.o VelocitySolverNSquared.o
$(CC) -o $@ $^
%.o: %.cpp
$(CC) -c $^ -o $@
.PHONY: clean
clean:
rm -f *.o
rm -f sheet
#include "VelocitySolverNSquared.h"
VelocitySolverNSquared::VelocitySolverNSquared(ArrayOfParticles & dstParticles, ArrayOfParticles & srcParticles, const int rank, const int size)
: dstParticles(dstParticles), srcParticles(srcParticles), rank(rank), size(size), timeC(0), timeT(0)
{
}
void VelocitySolverNSquared::ComputeVelocity()
{
const double i2pi = 0.5/M_PI;
const int N = dstParticles.Np;
for (int pass=0; pass<size; pass++)
{
// 1. exchange
Timer timerT;
timerT.start();
if (pass!=0)
{
// no need to send velocities
MPI_Request reqs[6];
MPI_Irecv(srcParticles.x , N, MPI_DOUBLE, (rank-pass+size)%size, 0, MPI_COMM_WORLD, &reqs[0]);
MPI_Isend(dstParticles.x , N, MPI_DOUBLE, (rank+pass )%size, 0, MPI_COMM_WORLD, &reqs[1]);
MPI_Irecv(srcParticles.y , N, MPI_DOUBLE, (rank-pass+size)%size, 1, MPI_COMM_WORLD, &reqs[2]);
MPI_Isend(dstParticles.y , N, MPI_DOUBLE, (rank+pass )%size, 1, MPI_COMM_WORLD, &reqs[3]);
MPI_Irecv(srcParticles.gamma, N, MPI_DOUBLE, (rank-pass+size)%size, 2, MPI_COMM_WORLD, &reqs[4]);
MPI_Isend(dstParticles.gamma, N, MPI_DOUBLE, (rank+pass )%size, 2, MPI_COMM_WORLD, &reqs[5]);
MPI_Status status[6];
MPI_Waitall(6, reqs, status);
}
timerT.stop();
timeT += timerT.get_timing();
// 2. compute local
Timer timerC;
timerC.start();
for (int i=0; i<N; i++) // loop over dst (particles of this rank)
{
double u=0;
double v=0;
for (int j=0; j<N; j++) // loop over src (received particles)
{
if(i+rank*N==j+((rank+pass)%size)*N) continue;
const double dx = dstParticles.x[i] - srcParticles.x[j];
const double dy = dstParticles.y[i] - srcParticles.y[j];
const double rsq = dx*dx + dy*dy;
const double irsq = 1./rsq;
u += -srcParticles.gamma[j] * dy * irsq;
v += srcParticles.gamma[j] * dx * irsq;
}
dstParticles.u[i] += i2pi * u;
dstParticles.v[i] += i2pi * v;
}
timerC.stop();
timeC += timerC.get_timing();
}
}
#ifndef VELOCITYSOLVERNSQUARED_H
#define VELOCITYSOLVERNSQUARED_H
#include "common.h"
#include "ArrayOfParticles.h"
class VelocitySolverNSquared
{
private:
ArrayOfParticles & dstParticles;
ArrayOfParticles & srcParticles;
const int rank, size;
public:
VelocitySolverNSquared(ArrayOfParticles & dstParticles, ArrayOfParticles & srcParticles,const int rank,const int size);
~VelocitySolverNSquared(){}
void ComputeVelocity();
double timeC, timeT;
};
#endif /* VELOCITYSOLVERNSQUARED_H */
#ifndef COMMON_H
#define COMMON_H
#include <iostream>
#include <vector>
#include <cmath>
#include <mpi.h>
#include "timer.hpp"
#endif /* COMMON_H */
#include <fstream>
#include <sstream>
#include <iomanip>
#include "common.h"
#include "ArrayOfParticles.h"
#include "VelocitySolverNSquared.h"
//#define WEAK_SCALING
using namespace std;
void WriteTimeToFile(const int Nranks, const double time, const char * sFileName)
{
ofstream outfile;
outfile.open(sFileName, ios::out | ios::app);
outfile << Nranks << " " << scientific << setprecision(6) << time << "\n";
outfile.close();
}
void Dump(ArrayOfParticles& dstParticles, ArrayOfParticles& allParticles, const int Np, const int NpProcess, const int step, const int rank)
{
// need to gather all particles here
MPI_Gather(dstParticles.x , NpProcess, MPI_DOUBLE, &allParticles.x[rank*NpProcess] , NpProcess, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(dstParticles.y , NpProcess, MPI_DOUBLE, &allParticles.y[rank*NpProcess] , NpProcess, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(dstParticles.u , NpProcess, MPI_DOUBLE, &allParticles.u[rank*NpProcess] , NpProcess, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(dstParticles.v , NpProcess, MPI_DOUBLE, &allParticles.v[rank*NpProcess] , NpProcess, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(dstParticles.gamma, NpProcess, MPI_DOUBLE, &allParticles.gamma[rank*NpProcess], NpProcess, MPI_DOUBLE, 0, MPI_COMM_WORLD);
if (rank==0)
{
char filename[500];
sprintf(filename, "vortexline_%5.5d.csv", step);
const double zero=0.0;
ofstream outfile;
outfile.open(filename);
outfile << "x,y,z,vx,vy,vz,g\n";
outfile << scientific << setprecision(6);
for (int i=0; i<Np;i++)
{
outfile << allParticles.x[i] << ",";
outfile << allParticles.y[i] << ",";
outfile << zero << ",";
outfile << allParticles.u[i] << ",";
outfile << allParticles.v[i] << ",";
outfile << zero << ",";
outfile << allParticles.gamma[i];
outfile << "\n";
}
outfile.close();
}
}
void DumpMPI(ArrayOfParticles& dstParticles, const int NpProcess, const int step, const int rank)
{
char filename[500];
sprintf(filename, "mpi_vortexline_%5.5d.csv", step);
MPI_File f;
MPI_File_open(MPI_COMM_WORLD, filename , MPI_MODE_WRONLY | MPI_MODE_CREATE, MPI_INFO_NULL, &f);
MPI_File_set_size (f, 0);
MPI_Offset base;
MPI_File_get_position(f, &base);
std::stringstream ss;
// only rank 0 writes the header
if (rank==0)
ss << "x,y,z,vx,vy,vz,g\n";
// all ranks write data
const double zero=0.0;
ss << scientific << setprecision(6);
for (int i = 0; i < NpProcess; ++i)
{
ss << dstParticles.x[i] << ",";
ss << dstParticles.y[i] << ",";
ss << zero << ",";
ss << dstParticles.u[i] << ",";
ss << dstParticles.v[i] << ",";
ss << zero << ",";
ss << dstParticles.gamma[i];
ss << "\n";
}
string content = ss.str();
MPI_Offset len = content.size();
MPI_Offset offset = 0;
MPI_Exscan(&len, &offset, 1, MPI_OFFSET, MPI_SUM, MPI_COMM_WORLD);
MPI_Status status;
MPI_File_write_at_all(f, base + offset, const_cast<char *>(content.c_str()), len, MPI_CHAR, &status);
MPI_File_close(&f);
}
int main (int argc, char ** argv)
{
MPI_Init(&argc,&argv);
int rank, size;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
// number of processes
if (rank==0)
std::cout << "Running with " << size << " MPI processes\n";
// timer setup
Timer timerIC, timerSim;
// time integration setup
const double dt = 0.001;
const double tfinal = 2.5;
const int ndump = 10;
// output
const bool bAnimation = true;
const bool bVerbose = true;
const bool dump_mpi = true;
// number of particles
const int N = 10000;
#ifndef WEAK_SCALING
size_t Np = (N/size)*size; // ensures that the number of particles is divisible by the number of workers
#else
size_t Np = N*sqrt(size); // keep number of work per rank fixed
#endif
timerIC.start();
// initialization
// each worker gets a part of the whole array
size_t NpProcess = Np/size;
// particle vectors
// dstParticles: particles owned by rank
// srcParticles: particles with which the interaction has to be computed
ArrayOfParticles dstParticles(NpProcess), srcParticles(NpProcess);
ArrayOfParticles allParticles(rank==0 ? Np : 0); // list of all particles for output
const double dx = 1./Np;
double totGamma=0.; // total circulation is sum over gamma of all particles
const double Gamma_s = 1.;
// initialize particles: position and circulation
for (size_t i=0; i<NpProcess; i++)
{
const double x = -0.5 + ((double)i+(double)rank*(double)NpProcess+.5)*dx;
dstParticles.x[i] = x;
dstParticles.y[i] = 0.0;
const double g = 4.*Gamma_s*x/sqrt(1.-4.*x*x);
dstParticles.gamma[i] = g*dx;
totGamma += dstParticles.gamma[i];
// each rank has also to compute the interactions with its own particles
srcParticles.x[i] = dstParticles.x[i];
srcParticles.y[i] = dstParticles.y[i];
srcParticles.gamma[i] = dstParticles.gamma[i];
}
MPI_Reduce(rank==0 ? MPI_IN_PLACE : &totGamma, &totGamma, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
timerIC.stop();
if (rank==0)
{
std::cout << "Number of particles: " << Np << std::endl;
std::cout << "Number of particles per process: " << NpProcess << std::endl;
std::cout << "Initial circulation: " << totGamma << std::endl;
std::cout << "IC time: " << timerIC.get_timing() << std::endl;
}
// initialize velocity solver
VelocitySolverNSquared VelocitySolver(dstParticles, srcParticles, rank, size);
timerSim.start();
double t=0;
int it=0;
for (it=1; it<=std::ceil(tfinal/dt); it++)
{
// reset particle velocities
dstParticles.ClearVelocities();
// compute velocities corresponding to time n
VelocitySolver.ComputeVelocity();
// dump the particles
if ((it-1)%ndump==0 && bAnimation)
{
if (!dump_mpi)
Dump(dstParticles,allParticles,Np,NpProcess,it-1,rank);
else
DumpMPI(dstParticles,NpProcess,it-1,rank);
}
// update time
if (rank==0)
t += dt;
MPI_Bcast(&t, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
if (it%ndump==0 && rank==0 && bVerbose)
std::cout << "Iteration " << it << " time " << t << std::endl;
// advance particles in time to n+1 using forward Euler
dstParticles.AdvectEuler(dt);
// update "source" particles
srcParticles = dstParticles;
}
// dump final state
if ((it-1)%ndump==0 && bAnimation)
{
dstParticles.ClearVelocities();
// compute velocities corresponding to time n
VelocitySolver.ComputeVelocity();
if (!dump_mpi)
Dump(dstParticles,allParticles,Np,NpProcess,it-1,rank);
else
DumpMPI(dstParticles,NpProcess,it-1,rank);
}
if (bVerbose)
std::cout << "Bye from rank " << rank << std::endl;
timerSim.stop();
if (rank==0)
{
char buf[500];
sprintf(buf, "timing.dat");
WriteTimeToFile(size,timerSim.get_timing(),buf);
std::cout << "#Ranks, Time - " << size << "\t" << timerSim.get_timing() << "\t( " << VelocitySolver.timeT << "\t" << VelocitySolver.timeC << " )\n";
}
MPI_Finalize();
return 0;
}
// Timer for the HPCSE I course
//
#ifndef HPCSE_TIMER_HPP
#define HPCSE_TIMER_HPP
#include <sys/time.h>
class Timer {
public:
Timer() {
start_time.tv_sec = 0;
start_time.tv_usec = 0;
stop_time.tv_sec = 0;
stop_time.tv_usec = 0;
}
inline void start() {
gettimeofday(&start_time, NULL);
}
inline void stop() {
gettimeofday(&stop_time, NULL);
}
double get_timing() const {
return (stop_time.tv_sec - start_time.tv_sec) + (stop_time.tv_usec - start_time.tv_usec)*1e-6;
}
private:
struct timeval start_time, stop_time;
};
#endif //HPCSE_TIMER_HPP
# TODO: complete the compile commands
all: nbody_serial_nocells nbody_omp_cells
nbody_serial_nocells: nbody_serial_nocells.cpp
#FILL IN COMPILE COMMAND HERE
nbody_omp_cells: nbody_omp_cells.cpp
#FILL IN COMPILE COMMAND HERE
clean:
rm -f nbody_serial_nocells nbody_omp_cells
// Example codes for HPC course
// (c) 2016 ETH Zurich
#include <array>
#include <vector>
#include <algorithm>
#include <functional>
#include <random>
#include <iostream>
#include <iterator>
#include <cassert>
#include <tuple>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <timer.hpp>
const unsigned DIMENSIONS = 2;
typedef std::size_t size_type;
typedef double scalar_type;
typedef std::array<scalar_type,DIMENSIONS> position;
struct potential
{
potential(scalar_type rm, scalar_type epsilon):
rm2_(rm*rm),
eps_(epsilon),
rc2_(0.0),
shift_(0)
{
rc2_ = 2.5*rm/std::pow(2,1/6.);
position x = {{}};
position y = {{rc2_}};
assert( x[0] == 0 && x[1] == 0 );
assert( y[1] == 0 );
shift_ = -(*this)(x,y,position());
rc2_ = rc2_*rc2_;
std::cout << "# Potential shift -V(rc=" << rc2_ << ")=" << shift_ << std::endl;
}
/// compute potential V(x,y)
scalar_type operator()(const position& x, const position& y, const position& extent) const
{
// TODO 1a
}
/// compute the Lennard-Jones force F(x,y) which particle y exerts on x and add it to f
void add_force(position& f, const position& x, const position& y, const position& extent) const
{
// TODO 1a
}
scalar_type cutoff_radius() const { return std::sqrt(rc2_); }
private:
/// compute distance between particle x and y
/// consider periodic boundaries
scalar_type dist(scalar_type x, scalar_type y, scalar_type extent) const
{
// TODO 1a
}
scalar_type rm2_; // r_m^2
scalar_type eps_; // \epsilon
scalar_type rc2_; // cut-off radius r_c^2
scalar_type shift_; // potential shift -V(r_c)
};
class simulation
{
public:
/// Initialize simulation in rectangular box with corners (0,0) and extent.
/// Initial positions and velocities are given as x, v.
simulation(const position& extent, const potential& pot,
const std::vector<position>& x, const std::vector<position>& v ):
extent_(extent),
potential_(pot),
x_(x),
v_(v),
a_(x.size())
{
assert( x.size() == v.size() );
calculate_forces(a_,x_);
}
/// evolve the system for [steps] time steps of size [dt]
void evolve(scalar_type dt, size_type steps)
{
using std::swap;
configuration aold(a_);
for( size_type s = 0; s < steps; ++s )
{
update_positions(x_,v_,a_,dt);
swap(a_,aold);