Page MenuHomec4science
Files F63180420

simulation.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, May 18, 07:16

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_, 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::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;
}
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();
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;
}
//if(i == 186 && j == 186){
// std::cout << Current << std::endl;
//}
}
}
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;
//std::cout << C << " " << G << std::endl;
for(std::size_t i(1);i<m()-1;i++){
for(std::size_t j(1);j<n()-1;j++) {
/*if(i == 486 && j == 486){
std::cout << oH(i,j) << " ------------- " << oHU(i,j) << " ----------------------- " << oHV(i,j) << std::endl;
}*/
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 ) );
/*if(i == 486 && j == 486){
std::cout << "HU CONTRIB : ----------------------------- "
<< + 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 ) ) << std::endl;
std::cout << dt << " ----------------------------- " << G << " ----------------------------- " << cH( i , j ) << " ----------------------------- " << Zdx( i , j ) << std::endl;
}*/
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

c4science · Help