Page MenuHomec4science
Files F64173892

Simulation.cc
No OneTemporary
Actions

File Metadata

Created
Sat, May 25, 02:02

Simulation.cc
View Options

//
// Created by Joachim Koerfer on 10.05.2019.
//
#include <cmath>
#include <iostream>
#include <iomanip>
#include "Simulation.h"
Simulation::Simulation(std::string path, int nx, int size, double Tend) : nx_(nx), size_(size), Tend_(Tend), delta_nx_((double)size/nx) {
//Begin to measure time
start_time_ = Clock::now();
std::string file_name = path + "Data_nx" + std::to_string(nx) + "_" + std::to_string(size) + "km_T0.2";
H_ = std::make_unique<Matrix>(file_name + "_h.bin", nx, nx);
Hu_ = std::make_unique<Matrix>(file_name + "_hu.bin", nx, nx);
Hv_ = std::make_unique<Matrix>(file_name + "_hv.bin", nx, nx);
Zdx_ = std::make_unique<Matrix>(file_name + "_Zdx.bin", nx, nx);
Zdy_ = std::make_unique<Matrix>(file_name + "_Zdy.bin", nx, nx);
//Temp storage
Ht_ = std::make_unique<Matrix>(*H_);
Hut_ = std::make_unique<Matrix>(*Hu_);
Hvt_ = std::make_unique<Matrix>(*Hv_);
//Initialisation of the variables necessary to the simulation
std::cout << std::setprecision(16);
T_ = 0.;
niter_ = 0;
delta_t_ = compute_delta_t();
cst_ = delta_t_/(2*delta_nx_);
}
void Simulation::compute() {
while(compute_step());
}
bool Simulation::compute_step() {
if (T_ > Tend_) {
//End the measurement of time
end_time_ = Clock::now();
measured_time = std::chrono::duration_cast<std::chrono::milliseconds>(end_time_ - start_time_).count();
std::cout << niter_ << std::endl;
nflops_ = (15 + 2 + 11 + 30 + 30 + 1 )*nx_*nx_*(double)niter_/measured_time*1e3;
return false;
}
Clock::time_point startTime = Clock::now();
enforce_boundary_conditions();
//Avoid dereferencing the pointers for each x,y
Matrix &Zdx = *Zdx_;
Matrix &Zdy = *Zdy_;
Matrix &H = *H_;
Matrix &Hu = *Hu_;
Matrix &Hv = *Hv_;
Matrix &Ht = *Ht_;
Matrix &Hut = *Hut_;
Matrix &Hvt = *Hvt_;
//*
//Use Lax-Friedrichs method to discretize the problem
for(int y = 1;y < nx_ - 1;y++) {
for (int x = 1; x < nx_ - 1; x++) {
Ht(x, y) = 0.25 * (H(x, y - 1) + H(x, y + 1) + H(x - 1, y) + H(x + 1, y))
+ cst_*(Hu(x, y - 1) - Hu(x, y + 1) + Hv(x - 1, y) - Hv(x + 1, y));
}
}
for(int y = 1;y < nx_ - 1;y++) {
for (int x = 1; x < nx_ - 1; x++) {
Hut(x, y) = 0.25 * (Hu(x, y - 1) + Hu(x, y + 1) + Hu(x - 1, y) + Hu(x + 1, y))
- delta_t_ * g_ * Ht(x, y) * Zdx(x, y)
+ cst_ * (Hu(x, y - 1) * Hu(x, y - 1) / H(x, y - 1) + .5 * g_ * H(x, y - 1) * H(x, y - 1)
- Hu(x, y + 1) * Hu(x, y + 1) / H(x, y + 1) - .5 * g_ * H(x, y + 1) * H(x, y + 1))
+ cst_ * (Hu(x - 1, y) * Hv(x - 1, y) / H(x - 1, y) -
Hu(x + 1, y) * Hv(x + 1, y) / H(x + 1, y));
}
}
for(int y = 1;y < nx_ - 1;y++) {
for (int x = 1; x < nx_ - 1; x++) {
Hvt(x,y) = 0.25*(Hv(x,y-1) + Hv(x,y+1) + Hv(x-1,y) + Hv(x+1,y))
- delta_t_*g_*Ht(x,y)*Zdy(x,y)
+ cst_*(Hv(x-1,y)*Hv(x-1,y)/H(x-1,y)+ .5*g_*H(x-1,y)*H(x-1,y)
- Hv(x+1,y)*Hv(x+1,y)/H(x+1,y) - .5*g_*H(x+1,y)*H(x+1,y)
+ Hu(x,y-1)*Hv(x,y-1)/H(x,y-1) - Hu(x,y+1)*Hv(x,y+1)/H(x,y+1));
}
}
copy_old_borders();
H_.swap(Ht_);
Hu_.swap(Hut_);
Hv_.swap(Hvt_);//*/
impose_tolerances();
niter_++;
if (T_ + delta_t_ > Tend_ && Tend_ != T_)
delta_t_ = Tend_ - T_;
T_ += delta_t_;
delta_t_ = compute_delta_t();
cst_ = delta_t_/(2*delta_nx_);
auto endTime = std::chrono::high_resolution_clock::now();
double speed = std::chrono::duration_cast<std::chrono::milliseconds>(endTime - startTime).count();
double ops = (15 + 2 + 11 + 30 + 30 + 1 )*nx_*nx_;
std::cout << "Current T = " << T_ << " , perf. = " << ops/speed/1e6 << "Gflops" << " , time = " << speed << " , " << T_/Tend_*100 << "%" << std::endl;
return true;
}
double Simulation::compute_delta_t() {
double nu;
double max_nu = 0;
Matrix &H = *H_;
Matrix &Hu = *Hu_;
Matrix &Hv = *Hv_;
for(int y = 0;y < nx_;y++) {
for (int x = 0; x < nx_; x++) {
double frac_hu = Hu(x,y)/H(x,y);
double frac_hv = Hv(x,y)/H(x,y);
double gh = std::sqrt(H(x,y)*g_);
nu = std::sqrt(pow(std::max(std::abs(frac_hu - gh), std::abs(frac_hu + gh)),2) + pow(std::max(std::abs(frac_hv - gh), std::abs(frac_hv + gh)),2));
if (nu > max_nu) {
max_nu = nu;
}
}
}
return delta_nx_/(std::sqrt(2)*max_nu);
}
void Simulation::enforce_boundary_conditions() {
for(int y = 0; y < nx_; y++){
//Boundary conditions
(*H_)(0,y) = (*H_)(1,y);
(*H_)(nx_ - 1, y) = (*H_)(nx_ - 2, y);
(*Hu_)(0,y) = (*Hu_)(1,y);
(*Hu_)(nx_ - 1, y) = (*Hu_)(nx_ - 2, y);
(*Hv_)(0,y) = (*Hv_)(1,y);
(*Hv_)(nx_ - 1, y) = (*Hv_)(nx_ - 2, y);
}
for(int x = 0; x < nx_; x++){
//Boundary conditions
(*H_)(x,0) = (*H_)(x,1);
(*H_)(x, nx_ - 1) = (*H_)(x, nx_ - 2);
(*Hu_)(x,0) = (*Hu_)(x,1);
(*Hu_)(x, nx_ - 1) = (*Hu_)(x, nx_ - 2);
(*Hv_)(x,0) = (*Hv_)(x,1);
(*Hv_)(x, nx_ - 1) = (*Hv_)(x, nx_ - 2);
}
}
void Simulation::copy_old_borders() {
for(int y = 0; y < nx_; y++){
//Old borders have to be copied from temp. storage before swapping pointers
(*H_)(0,y) = (*Ht_)(0, y);
(*Hu_)(0,y) = (*Hut_)(0, y);
(*Hv_)(0,y) = (*Hvt_)(0, y);
(*H_)(nx_ - 1,y) = (*Ht_)(nx_ - 1, y);
(*Hu_)(nx_ - 1,y) = (*Hut_)(nx_ - 1, y);
(*Hv_)(nx_ - 1,y) = (*Hvt_)(nx_ - 1, y);
}
for(int x = 0; x < nx_; x++){
//Old borders have to be copied from temp. storage before swapping pointers
(*H_)(x,0) = (*Ht_)(x, 0);
(*Hu_)(x,0) = (*Hut_)(x, 0);
(*Hv_)(x,0) = (*Hvt_)(x, 0);
(*H_)(x, nx_ - 1) = (*Ht_)(x, nx_ - 1);
(*Hu_)(x, nx_ - 1) = (*Hut_)(x, nx_ - 1);
(*Hv_)(x, nx_ - 1) = (*Hvt_)(x, nx_ - 1);
}
}
void Simulation::save(const std::string &file_name) {
(*H_).save(file_name);
}
void Simulation::impose_tolerances() {
Matrix &H = *H_;
Matrix &Hu = *Hu_;
Matrix &Hv = *Hv_;
for(int y = 0;y < nx_;y++) {
for (int x = 0; x < nx_; x++) {
if(H(x,y) <= 5*1e-4) {
Hu(x,y) = 0;
Hv(x,y) = 0;
if(H(x,y) < 0)
H(x,y) = 1e-5;
}
}
}
}

Event Timeline

c4science · Help