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 (1 d, 21 h)
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