Skip to content

Commit 0462723

Browse files
mdelormemdelorme
committed
First commit
0 parents  commit 0462723

26 files changed

+1630
-0
lines changed

‎.gitignore

+2
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,2 @@
1+
.vscode
2+
build

‎.gitmodules

+3
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
[submodule "external/inih"]
2+
path = external/inih
3+
url = [email protected]:jtilly/inih.git

‎BoundaryConditions.h

+107
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
#pragma once
2+
3+
#include <map>
4+
#include <cassert>
5+
6+
#include "SimInfo.h"
7+
8+
namespace fv1d {
9+
10+
State fill_absorbing(Array &Q, int i, int j, int iref, int jref, real_t dt, IDir dir) {
11+
return Q[jref][iref];
12+
}
13+
14+
State fill_reflecting(Array &Q, int i, int j, int iref, int jref, real_t dt, IDir dir) {
15+
int isym, jsym;
16+
if (dir == IX) {
17+
int ipiv = (i < iref ? ibeg : iend);
18+
isym = 2*ipiv - i - 1;
19+
jsym = j;
20+
}
21+
else {
22+
int jpiv = (j < jref ? jbeg : jend);
23+
isym = i;
24+
jsym = 2*jpiv - j - 1;
25+
}
26+
27+
State q = Q[jsym][isym];
28+
29+
if (dir == IX)
30+
q[IU] *= -1.0;
31+
else
32+
q[IV] *= -1.0;
33+
34+
return q;
35+
}
36+
37+
State fill_periodic(Array &Q, int i, int j, int iref, int jref, real_t dt, IDir dir) {
38+
if (dir == IX) {
39+
if (i < ibeg)
40+
i += Nx;
41+
else
42+
i -= Nx;
43+
}
44+
else {
45+
if (j < jbeg)
46+
j += Ny;
47+
else
48+
j -= Ny;
49+
}
50+
51+
return Q[j][i];
52+
}
53+
54+
/**
55+
* One time initialization for the boundary conditions
56+
**/
57+
void init_boundaries() {
58+
assert(Ng == 2); // Most functions have been coded with Ng = 2 !!!!
59+
60+
switch (boundary_x) {
61+
case BC_ABSORBING: bc_function_x = fill_absorbing; break;
62+
case BC_PERIODIC: bc_function_x = fill_periodic; break;
63+
default: bc_function_x = fill_reflecting; break;
64+
}
65+
66+
switch (boundary_y) {
67+
case BC_ABSORBING: bc_function_y = fill_absorbing; break;
68+
case BC_PERIODIC: bc_function_y = fill_periodic; break;
69+
default: bc_function_y = fill_reflecting; break;
70+
}}
71+
72+
/**
73+
* Fills the boundary with a value. The value is returned from a user-defed method.
74+
* Note that this function is overly complicated for this type of code.
75+
* There is a real reason to that: to replicate the way it is called in dyablo
76+
*/
77+
void fill_boundaries(Array &Q, real_t dt) {
78+
// Filling boundaries along X
79+
for (int j=jbeg; j < jend; ++j) {
80+
for (int i=0; i < Ng; ++i) {
81+
int ileft = i;
82+
int iright = iend+i;
83+
int iref_left = ibeg;
84+
int iref_right = iend-1;
85+
86+
// And apply
87+
Q[j][ileft] = bc_function_x(Q, ileft, j, iref_left, j, dt, IX);
88+
Q[j][iright] = bc_function_x(Q, iright, j, iref_right, j, dt, IX);
89+
}
90+
}
91+
92+
for (int j=0; j < Ng; ++j) {
93+
for (int i=0; i < Ntx; ++i) {
94+
int jtop = j;
95+
int jbot = jend+j;
96+
int jref_top = jbeg;
97+
int jref_bot = jend-1;
98+
99+
// And apply
100+
Q[jtop][i] = bc_function_y(Q, i, jtop, i, jref_top, dt, IY);
101+
Q[jbot][i] = bc_function_y(Q, i, jbot, i, jref_bot, dt, IY);
102+
}
103+
}
104+
}
105+
106+
107+
}

‎CMakeLists.txt

+22
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
cmake_minimum_required(VERSION 3.16)
2+
3+
project(fv2d LANGUAGES C CXX)
4+
5+
find_package(OpenMP)
6+
7+
add_subdirectory(external/HighFive)
8+
9+
include_directories(external/inih)
10+
11+
# Various sets
12+
set(CMAKE_CXX_STANDARD 17)
13+
14+
add_executable(fv2d
15+
main.cpp)
16+
17+
if (OpenMP_CXX_FOUND)
18+
target_link_libraries(fv2d PRIVATE HighFive OpenMP::OpenMP_CXX)
19+
else()
20+
message(WARNING "OpenM not found, compiling in serial")
21+
target_link_libraries(fv2d PRIVATE HighFive)
22+
endif()

‎ComputeDt.h

+49
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,49 @@
1+
#pragma once
2+
3+
#include "SimInfo.h"
4+
#include "ThermalConduction.h"
5+
#include "Viscosity.h"
6+
7+
namespace fv1d {
8+
9+
real_t compute_dt(Array &Q, real_t max_dt, real_t t, bool diag) {
10+
11+
real_t inv_dt_hyp = 0.0;
12+
real_t inv_dt_tc = 0.0;
13+
real_t inv_dt_visc = 0.0;
14+
15+
#pragma omp parallel reduction(max:inv_dt_hyp) reduction(max:inv_dt_tc) reduction(max:inv_dt_visc)
16+
for (int j=jbeg; j<jend; ++j) {
17+
for (int i=ibeg; i<iend; ++i) {
18+
Pos pos = get_pos(i, j);
19+
20+
// Hyperbolic CFL
21+
real_t cs = speed_of_sound(Q[j][i]);
22+
inv_dt_hyp = std::max(inv_dt_hyp, (cs + std::fabs(Q[j][i][IU]))/dx);
23+
inv_dt_hyp = std::max(inv_dt_hyp, (cs + std::fabs(Q[j][i][IV]))/dy);
24+
25+
// Parabolic
26+
if (thermal_conductivity_active) {
27+
inv_dt_tc = std::max(inv_dt_tc, 2.0 * compute_kappa(pos[IX], pos[IY]) / (dx*dx));
28+
inv_dt_tc = std::max(inv_dt_tc, 2.0 * compute_kappa(pos[IX], pos[IY]) / (dy*dy));
29+
}
30+
31+
if (viscosity_active) {
32+
inv_dt_visc = std::max(inv_dt_visc, 2.0 * compute_mu(pos[IX], pos[IY]) / (dx*dx));
33+
inv_dt_visc = std::max(inv_dt_visc, 2.0 * compute_mu(pos[IY], pos[IY]) / (dy*dy));
34+
}
35+
}
36+
}
37+
38+
if (diag) {
39+
std::cout << "Computing dts at (t=" << t << ") : dt_hyp=" << 1.0/inv_dt_hyp
40+
<< "; dt_TC=" << 1.0/inv_dt_tc
41+
<< "; dt_visc=" << 1.0/inv_dt_visc << std::endl;
42+
}
43+
44+
real_t dt = CFL / (std::max({inv_dt_hyp, inv_dt_tc, inv_dt_visc}));
45+
46+
return std::min(dt, max_dt);
47+
}
48+
49+
}

‎Init.h

+187
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,187 @@
1+
#pragma once
2+
3+
#include <fstream>
4+
5+
#include "SimInfo.h"
6+
#include "BoundaryConditions.h"
7+
8+
9+
10+
namespace fv1d {
11+
12+
void init_sod_x(Array &Q) {
13+
std::cout << "Initializing SOD (x)" << std::endl;
14+
15+
for (int j=jbeg; j < jend; ++j) {
16+
for (int i=ibeg; i < iend; ++i) {
17+
if (get_pos(i, j)[IX] <= 0.5) {
18+
Q[j][i][IR] = 1.0;
19+
Q[j][i][IP] = 1.0;
20+
Q[j][i][IU] = 0.0;
21+
}
22+
else {
23+
Q[j][i][IR] = 0.125;
24+
Q[j][i][IP] = 0.1;
25+
Q[j][i][IU] = 0.0;
26+
}
27+
}
28+
}
29+
}
30+
31+
void init_sod_y(Array &Q) {
32+
std::cout << "Initializing SOD (y)" << std::endl;
33+
34+
for (int j=jbeg; j < jend; ++j) {
35+
for (int i=ibeg; i < iend; ++i) {
36+
if (get_pos(i, j)[IY] <= 0.5) {
37+
Q[j][i][IR] = 1.0;
38+
Q[j][i][IP] = 1.0;
39+
Q[j][i][IU] = 0.0;
40+
}
41+
else {
42+
Q[j][i][IR] = 0.125;
43+
Q[j][i][IP] = 0.1;
44+
Q[j][i][IU] = 0.0;
45+
}
46+
}
47+
}
48+
}
49+
50+
void init_blast(Array &Q) {
51+
std::cout << "Initializing blast" << std::endl;
52+
53+
real_t xmid = 0.5 * (xmin+xmax);
54+
real_t ymid = 0.5 * (ymin+ymax);
55+
56+
for (int j=jbeg; j < jend; ++j) {
57+
58+
for (int i=ibeg; i < iend; ++i) {
59+
Pos pos = get_pos(i, j);
60+
real_t x = pos[IX];
61+
real_t y = pos[IY];
62+
63+
real_t xr = xmid - x;
64+
real_t yr = ymid - y;
65+
real_t r = sqrt(xr*xr+yr*yr);
66+
67+
if (r < 0.2) {
68+
Q[j][i][IR] = 1.0;
69+
Q[j][i][IU] = 0.0;
70+
Q[j][i][IP] = 10.0;
71+
}
72+
else {
73+
Q[j][i][IR] = 1.2;
74+
Q[j][i][IU] = 0.0;
75+
Q[j][i][IP] = 0.1;
76+
}
77+
}
78+
}
79+
}
80+
81+
void init_C91(Array &Q) {
82+
std::cout << "Initializing C91" << std::endl;
83+
84+
srand(seed);
85+
86+
for (int j=jbeg; j < jend; ++j) {
87+
for (int i=ibeg; i < iend; ++i) {
88+
Pos pos = get_pos(i, j);
89+
real_t x = pos[IX];
90+
real_t y = pos[IY];
91+
92+
real_t T = (1.0 + theta1*y);
93+
real_t rho = std::pow(T, m1);
94+
real_t prs = std::pow(T, m1+1.0);
95+
96+
real_t pert = c91_pert * ((float)rand() / RAND_MAX * 2.0 - 1.0);
97+
98+
prs = prs * (1.0 + pert);
99+
100+
Q[j][i][IR] = rho;
101+
Q[j][i][IU] = 0.0;
102+
Q[j][i][IV] = 0.0;
103+
Q[j][i][IP] = prs;
104+
}
105+
}
106+
}
107+
108+
void init_diff(Array &Q) {
109+
std::cout << "Initializing Diffusion test" << std::endl;
110+
111+
real_t xmid = 0.5 * (xmin + xmax);
112+
real_t ymid = 0.5 * (ymin + ymax);
113+
114+
for (int j=jbeg; j < jend; ++j) {
115+
for (int i=ibeg; i < iend; ++i) {
116+
Pos pos = get_pos(i, j);
117+
118+
real_t x0 = (pos[IX]-xmid);
119+
real_t y0 = (pos[IY]-ymid);
120+
121+
real_t r = sqrt(x0*x0+y0*y0);
122+
123+
if (r < 0.2)
124+
Q[j][i][IR] = 1.0;
125+
else
126+
Q[j][i][IR] = 0.1;
127+
128+
Q[j][i][IP] = 1.0;
129+
Q[j][i][IU] = 1.0;
130+
Q[j][i][IV] = 0.0;
131+
}
132+
}
133+
}
134+
135+
void init_rayleigh_taylor(Array &Q) {
136+
std::cout << "Initializing Rayleigh-Taylor" << std::endl;
137+
138+
srand(seed);
139+
140+
real_t ymid = 0.5*(ymin + ymax);
141+
142+
for (int j=ibeg; j < jend; ++j) {
143+
for (int i=ibeg; i < iend; ++i) {
144+
Pos pos = get_pos(i, j);
145+
real_t x = pos[IX];
146+
real_t y = pos[IY];
147+
148+
const real_t P0 = 2.5;
149+
150+
if (y < ymid) {
151+
Q[j][i][IR] = 1.0;
152+
Q[j][i][IU] = 0.0;
153+
Q[j][i][IP] = P0 + 0.1 * g * y;
154+
}
155+
else {
156+
Q[j][i][IR] = 2.0;
157+
Q[j][i][IU] = 0.0;
158+
Q[j][i][IP] = P0 + 0.1 * g * y;
159+
}
160+
161+
if (y > -1.0/3.0 && y < 1.0/3.0)
162+
Q[j][i][IV] = 0.01 * (1.0 + cos(4*M_PI*x)) * (1 + cos(3.0*M_PI*y))/4.0;
163+
}
164+
}
165+
}
166+
167+
168+
void init(Array &Q) {
169+
if (problem == "sod_x")
170+
init_sod_x(Q);
171+
else if (problem == "sod_y")
172+
init_sod_y(Q);
173+
else if (problem == "blast")
174+
init_blast(Q);
175+
else if (problem == "C91")
176+
init_C91(Q);
177+
else if (problem == "rayleigh-taylor")
178+
init_rayleigh_taylor(Q);
179+
else if (problem == "diffusion")
180+
init_diff(Q);
181+
182+
// One time init stuff
183+
init_boundaries();
184+
}
185+
186+
187+
}

0 commit comments

Comments
 (0)