Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F86287385
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
Sat, Oct 5, 14:01
Size
6 KB
Mime Type
text/x-c
Expires
Mon, Oct 7, 14:01 (2 d)
Engine
blob
Format
Raw Data
Handle
21392517
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 <omp.h>
#include <sstream>
#include <limits>
#include <chrono>
Simulation::Simulation (std::size_t NX_, double G_, int SIZE_, double TEND_, double DX_, std::string dir_)
: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)){};
void Simulation::readInitialConditionsFromFile () {
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
H.current () = Reader::readGridFromFile(filenameH.str(), NX,NX);
HU.current() = Reader::readGridFromFile(filenameHU.str(), NX,NX);
HV.current() = Reader::readGridFromFile(filenameHV.str(), NX,NX);
// Load topography slopes from data files
Zdx = Reader::readGridFromFile(filenameZdx.str(), NX,NX);
Zdy = Reader::readGridFromFile(filenameZdy.str(), NX,NX);
//std::cout << "ZDX (486,486) " << Zdx(486,486) << std::endl;
//std::cout << "ZDy (486,486) " << Zdy(486,486) << std::endl;
// Allocate memory for computing
H.old () = Grid(NX,NX);
HU.old() = Grid(NX,NX);
HV.old() = Grid(NX,NX);
std::cout << "Data Loaded - Start of simulation" << 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 << "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 << "Time for a step : " << ms << std::endl;
std::cout << "\r" << " dt : " << dt << " Computing T: " << T << ". " << (100*T/TEND) <<"%" << " --------- Time for a step : " << ms << " " << std::flush;
}
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 () {
double Max = 0.0;
double Current;
Grid& cH = H.current();
Grid& cHU = HU.current();
Grid& cHV = HV.current();
#pragma omp parallel for reduction(max: Max)
for (std::size_t i (0); i < m (); i++) {
for (std::size_t j (0); j < n (); j++) {
Current = sqrt (+pow (std::max (abs (cHU (i, j) / cH (i, j) - sqrt (cH (i, j) * G)),
abs (cHU (i, j) / cH (i, j) + sqrt (cH (i, j) * G))), 2)
+ pow (std::max (abs (cHV (i, j) / cH (i, j) - sqrt (cH (i, j) * G)),
abs (cHV (i, j) / cH (i, j) + sqrt (cH (i, j) * G))), 2));
if (Current > Max) {
Max = Current;
}
}
}
dt = DX/(sqrt(2)*Max);
//std::cout << "DT : " << dt << " --------------------------" << H.current()(486,486) << std::endl;
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();
Grid& cH = H.current();
double C = 0.5*dt/DX;
int id;
#pragma omp parallel for
for(std::size_t i(1);i<m()-1;i++){
for(std::size_t j(1);j<n()-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 ) );
}
}
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);
}
Event Timeline
Log In to Comment