Page MenuHomec4science
Files F64560781

simulation.cpp
No OneTemporary
Actions

File Metadata

Created
Mon, May 27, 19:28

simulation.cpp
View Options

//
// Created by Arnaud Pannatier on 06.05.18.
// Based on the code of :
// - Nicolas Richart <nicolas.richart@epfl.ch>
// - Vincent Keller <vincent.keller@epfl.ch>
// - Vittoria Rezzonico <vittoria.rezzonico@epfl.ch>
// See the files AUTHORS and COPYRIGHT for the concerning information
//
/* -------------------------------------------------------------------------- */
#include "simulation.h"
#include "reader.h"
/* -------------------------------------------------------------------------- */
#include <math.h>
#include <sstream>
#include <chrono>
Simulation::Simulation (int NX_, double G_, int SIZE_, double TEND_, double DX_, std::string dir_, MPI_Comm communicator_)
:NX(NX_), NY(NX_),dir(dir_), G(G_),TEND(TEND_),DX(DX_),dt(0),T(0),nt(0),
SIZE(SIZE_), H(DoubleBuffer(NX,NX)), HU(DoubleBuffer(NX,NX)), HV(DoubleBuffer(NX,NX)),
Zdx(Grid(NX,NX)), Zdy(Grid(NX,NX)){
communicator = communicator_;
// retrieving the number of proc and the rank in the proc pool
MPI_Comm_rank(communicator, &prank);
MPI_Comm_size(communicator, &psize);
// computation of the local size of the grid the remainder is spread equally
// on the first processors
local_x = NX / psize + (prank < NX % psize ? 1 : 0);
local_y = NY;
// adding the ghosts lines if needed
if (psize > 1)
local_x += (prank == 0 || prank == psize - 1) ? 1 : 2;
// computing the offsets of the local grid in the global one
offset_x = (NX / psize) * prank + (prank < NX % psize ? prank : NX % psize);
offset_y = 0;
// resizing the different grids
resizeGrids(local_x, local_y);
// determining the rank of the neighbors
north_prank = (prank == 0 ? MPI_PROC_NULL : prank - 1);
south_prank = (prank == (psize - 1) ? MPI_PROC_NULL : prank + 1);
// Some info if needed to debug
std::cout << prank << " " << NX << " " << NY << " "
<< local_x << " " << local_y << " " << offset_x << " "
<< offset_y << " " << north_prank << " " << south_prank
<< std::endl;
};
void Simulation::readInitialConditionsFromFile () {
int x_start = (prank == 0 ? 0 : 1);
int x_end = (prank == psize - 1 ? local_x : local_x - 1);
std::stringstream filename, filenameH, filenameHU,filenameHV, filenameZdx, filenameZdy;
filename << dir << "Data_nx" << NX << "_" << SIZE << "km_T" <<"0.2";
filenameH << filename.str() << "_h.bin";
filenameHU << filename.str() << "_hu.bin";
filenameHV << filename.str() << "_hv.bin";
filenameZdx << filename.str() << "_Zdx.bin";
filenameZdy << filename.str() << "_Zdy.bin";
// Load initial condition from data files
Grid HRead = Reader::readGridFromFile(filenameH.str(), NX,NX);
Grid HURead = Reader::readGridFromFile(filenameHU.str(), NX,NX);
Grid HVRead = Reader::readGridFromFile(filenameHV.str(), NX,NX);
// Load topography slopes from data files
Grid ZdxRead = Reader::readGridFromFile(filenameZdx.str(), NX,NX);
Grid ZdyRead = Reader::readGridFromFile(filenameZdy.str(), NX,NX);
for (int x = x_start; x < x_end; x++) {
for (int y = 0; y < local_y; y++) {
H.current()(x,y) = HRead(offset_x + x - x_start,offset_y + y);
HU.current()(x,y) = HURead(offset_x + x - x_start,offset_y + y);
HV.current()(x,y) = HVRead(offset_x + x - x_start,offset_y + y);
Zdx(x,y) = ZdxRead(offset_x + x - x_start,offset_y + y);
Zdy(x,y) = ZdyRead(offset_x + x - x_start,offset_y + y);
}
}
std::cout << " Process " << prank+1 << " Initialized ! " << std::endl;
}
/* -------------------------------------------------------------------------- */
std::tuple<double, int> Simulation::compute() {
int s = 0;
double totalTime = 0;
auto begin = std::chrono::high_resolution_clock::now ();
readInitialConditionsFromFile();
while(T < TEND){
s++;
auto start = std::chrono::high_resolution_clock::now ();
compute_mu_and_set_dt ();
std::cout << " Process " << prank+1 << "dt : " << dt << " Computing T: " << (T+dt) << ". " << (100*(T+dt)/TEND) <<"%" << std::endl;
compute_step();
T = T + dt;
nt = nt + 1;
auto end = std::chrono::high_resolution_clock::now();
auto duration = end-start;
auto ms = std::chrono::duration_cast<std::chrono::milliseconds>(duration).count();
std::cout <<" Process " << prank+1 << " Time for a step : " << ms << std::endl;
}
writeResultsToFile ();
auto end = std::chrono::high_resolution_clock::now();
auto duration = end-begin;
totalTime = std::chrono::duration_cast<std::chrono::milliseconds>(duration).count();
return std::make_tuple(totalTime, s);
}
/* -------------------------------------------------------------------------- */
void Simulation::compute_mu_and_set_dt () {
// No need to take account of the ghost cell for computing the timesteps
int x_start = (prank == 0 ? 0 : 1);
int x_end = (prank == psize - 1 ? local_x : local_x - 1);
double Max = 0.0;
double Current = 0.0;
Grid& cH = H.current();
Grid& cHU = HU.current();
Grid& cHV = HV.current();
for (int x = x_start; x < x_end; x++) {
for (int y = 0; y < local_y; y++) {
Current = sqrt( +pow(std::max( abs(cHU(x,y)/cH(x,y) - sqrt(cH(x,y)*G)) ,
abs(cHU(x,y)/cH(x,y) + sqrt(cH(x,y)*G)) ) ,2)
+pow(std::max( abs(cHV(x,y)/cH(x,y) - sqrt(cH(x,y)*G)) ,
abs(cHV(x,y)/cH(x,y) + sqrt(cH(x,y)*G)) ),2));
if(Current > Max){
Max = Current;
}
if(cH(x,y) < 5e-6){
std::cout << " Process " << prank+1 << " Low value : " << cH(x,y) << std::endl;
}
}
}
std::cout << " Process " << prank+1 << " Max : " << Max << std::endl;
// Getting the max the value of all the processors together
MPI_Allreduce(MPI_IN_PLACE, & Max, 1, MPI_DOUBLE, MPI_MAX, communicator);
dt = DX/(sqrt(2)*Max);
C = 0.5*dt/DX;
if(T+dt > TEND){
dt = TEND-T;
}
}
void Simulation::compute_step () {
std::cout << " Compute Step " << std::endl;
H.swap();
HU.swap();
HV.swap();
H.old().applyBoundaryConditions();
HU.old().applyBoundaryConditions();
HV.old().applyBoundaryConditions();
Grid & oH = H.old();
Grid & oHU = HU.old();
Grid & oHV = HV.old();
std::cout<< " Process " << prank+1 << "Adress of oH(1,0) : " << &oH(1,0) << std::endl;
// Taking care of communications going up (so receiving from bottom)
// oH
MPI_Sendrecv(&(oH(1, 0)), local_y, MPI_DOUBLE, north_prank, 0,
&(oH(local_x - 1, 0)), local_y, MPI_DOUBLE, south_prank, 0,
communicator, MPI_STATUS_IGNORE);
// oHU
MPI_Sendrecv(&oHV(1, 0), local_y, MPI_DOUBLE, north_prank, 0,
&oHV(local_x - 1, 0), local_y, MPI_DOUBLE, south_prank, 0,
communicator, MPI_STATUS_IGNORE);
//oHV
MPI_Sendrecv(&(oHU(1, 0)), local_y, MPI_DOUBLE, north_prank, 0,
&(oHU(local_x - 1, 0)), local_y, MPI_DOUBLE, south_prank, 0,
communicator, MPI_STATUS_IGNORE);
// Taking care of communications going down (so receiving from top)
MPI_Sendrecv(&oH(local_x - 2, 0), local_y, MPI_DOUBLE, south_prank, 0,
&oH(0, 0), local_y, MPI_DOUBLE, north_prank, 0,
communicator, MPI_STATUS_IGNORE);
MPI_Sendrecv(&oHU(local_x - 2, 0), local_y, MPI_DOUBLE, south_prank, 0,
&oHU(0, 0), local_y, MPI_DOUBLE, north_prank, 0,
communicator, MPI_STATUS_IGNORE);
MPI_Sendrecv(&oHV(local_x - 2, 0), local_y, MPI_DOUBLE, south_prank, 0,
&oHV(0, 0), local_y, MPI_DOUBLE, north_prank, 0,
communicator, MPI_STATUS_IGNORE);
std::cout<< " AFTER EXCHANGE : Process " << prank+1 << "Adress of oH(1,0) : " << &oH(1,0) << std::endl;
int x_start = (prank == 0 ? 0 : 1);
int x_end = (prank == psize - 1 ? local_x : local_x - 1);
for (int x = x_start; x < x_end; x++) {
compute_row (x);
}
H.current().imposeTolerances (0.0, 1e-5, H.current());
HU.current().imposeTolerances (1e-4, 0.0, H.current());
HV.current().imposeTolerances (1e-4, 0.0, H.current());
}
void Simulation::writeResultsToFile () {
std::stringstream file;
file << dir << "Solution_nx" << NX << "_" << SIZE << "km_T"<<TEND<<"_h.bin";
Reader::writeGridInFile (H.current (), file.str(), NX,NX);
}
void Simulation::resizeGrids (int nx, int ny) {
H = DoubleBuffer(nx, ny);
HU = DoubleBuffer(nx, ny);
HV = DoubleBuffer(nx, ny);
Zdx = Grid(nx,ny);
Zdy = Grid(nx,ny);
}
void Simulation::compute_row (int i) {
Grid& oH = H.old();
Grid& oHU = HU.old();
Grid& oHV = HV.old();
Grid& cH = H.current();
for(int j(1);j<y()-1;j++) {
H.current()(i,j) = + 0.25 * ( + oH ( i , j-1 ) + oH ( i , j+1 ) + oH ( i-1 , j ) + oH ( i+1 , j))
+ C * ( + oHU( i , j-1 ) - oHU( i , j+1 )
+ oHV( i-1 , j ) - oHV( i+1 , j ));
HU.current()(i,j) = + 0.25 * ( + oHU( i , j-1 ) + oHU( i , j+1 ) + oHU( i-1 , j ) + oHU( i+1 , j )) - dt * G * cH( i , j ) * Zdx( i , j )
+ C * ( + pow( oHU( i , j-1 ) ,2) / oH( i , j-1 ) + 0.5 * G * pow( oH( i , j-1 ) ,2)
- pow( oHU( i , j+1 ) ,2) / oH( i , j+1 ) - 0.5 * G * pow( oH( i , j+1 ) ,2) )
+ C * ( + oHU( i-1 , j ) * oHV( i-1 , j ) / oH( i-1 , j )
- oHU( i+1 , j ) * oHV( i+1 , j ) / oH( i+1 , j ) );
HV.current ()(i,j) = + 0.25 * ( + oHV( i , j-1 ) + oHV( i , j+1 ) + oHV( i-1 , j ) + oHV( i+1 , j )) - dt * G * cH( i , j ) * Zdy( i , j )
+ C * ( + pow( oHV( i-1 , j ),2) / oH( i-1 , j ) + 0.5 * G * pow( oH( i-1 , j ),2)
- pow( oHV( i+1 , j ),2) / oH( i+1 , j ) - 0.5 * G * pow( oH( i+1 , j ),2) )
+ C * ( + oHU( i , j-1 ) * oHV( i , j-1 ) / oH( i , j-1 )
- oHU( i , j+1 ) * oHV( i , j+1 ) / oH( i , j+1 ) );
}
}

Event Timeline

c4science · Help