Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F87586017
simulation.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Oct 13, 13:35
Size
9 KB
Mime Type
text/x-c
Expires
Tue, Oct 15, 13:35 (2 d)
Engine
blob
Format
Raw Data
Handle
21622125
Attached To
rTSUNAMIPROJECT PHPC-TSUNAMI-PROJECT
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>
Simulation::Simulation (std::size_t NX_, double G_, std::size_t 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 << " " << global_m << " " << global_n << " "
// << local_m << " " << local_n << " " << offset_m << " "
// << offset_n << " " << north_prank << " " << south_prank
// << std::endl;
};
void Simulation::readInitialConditionsFromFile () {
std::size_t x_start = (prank == 0 ? 0 : 1);
std::size_t 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 (std::size_t x = x_start; x < x_end; x++) {
for (std::size_t 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 << " 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
std::size_t x_start = (prank == 0 ? 0 : 1);
std::size_t x_end = (prank == psize - 1 ? local_x : local_x - 1);
double Max = 0.0;
double Current;
Grid& cH = H.current();
Grid& cHU = HU.current();
Grid& cHV = HV.current();
for (std::size_t x = x_start; x < x_end; x++) {
for (std::size_t 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;
}
}
}
// 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 () {
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();
// Taking care of communications going up (so receiving from bottom)
// oH
MPI_Sendrecv(&oH(1, 0), y(), MPI_DOUBLE, north_prank, 0,
&oH(x() - 1, 0), y(), MPI_DOUBLE, south_prank, 0,
communicator, MPI_STATUS_IGNORE);
// oHU
MPI_Sendrecv(&oHV(1, 0), y(), MPI_DOUBLE, north_prank, 0,
&oHV(x() - 1, 0), y(), MPI_DOUBLE, south_prank, 0,
communicator, MPI_STATUS_IGNORE);
//oHV
MPI_Sendrecv(&oHU(1, 0), y(), MPI_DOUBLE, north_prank, 0,
&oHU(x() - 1, 0), y(), MPI_DOUBLE, south_prank, 0,
communicator, MPI_STATUS_IGNORE);
// Taking care of communications going down (so receiving from top)
MPI_Sendrecv(&oH(x() - 2, 0), y(), MPI_DOUBLE, south_prank, 0,
&oH(0, 0), y(), MPI_DOUBLE, north_prank, 0,
communicator, MPI_STATUS_IGNORE);
MPI_Sendrecv(&oHU(x() - 2, 0), y(), MPI_DOUBLE, south_prank, 0,
&oHU(0, 0), y(), MPI_DOUBLE, north_prank, 0,
communicator, MPI_STATUS_IGNORE);
MPI_Sendrecv(&oHV(x() - 2, 0), y(), MPI_DOUBLE, south_prank, 0,
&oHV(0, 0), y(), MPI_DOUBLE, north_prank, 0,
communicator, MPI_STATUS_IGNORE);
for(std::size_t i(1);i<x()-1;i++){
compute_row (i);
}
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 (size_t nx, size_t 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 (size_t i) {
Grid& oH = H.old();
Grid& oHU = HU.old();
Grid& oHV = HV.old();
Grid& cH = H.current();
for(std::size_t 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
Log In to Comment