Skip to content

Commit af03ba8

Browse files
pgretepgrete
committed
Try reflect in kernel
1 parent c5c6bb3 commit af03ba8

1 file changed

Lines changed: 115 additions & 81 deletions

File tree

‎src/hydro/hydro.cpp‎

Lines changed: 115 additions & 81 deletions
Original file line numberDiff line numberDiff line change
@@ -1066,40 +1066,14 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshData<Real>> &md) {
10661066
parthenon::ScratchPad2D<Real>::shmem_size(num_scratch_vars, nx1) * 2;
10671067

10681068
auto riemann = Riemann<fluid, rsolver>();
1069-
1070-
const auto num_blocks = md->NumBlocks();
10711069
auto const &outside_pack = md->PackVariables(std::vector<std::string>{"outside"});
1072-
pmb->par_for(
1073-
"Set x1 reflect", 0, num_blocks - 1, kl, ku, jl, ju, ib.s - 1, ib.e + 1,
1074-
KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
1075-
auto &cons = cons_in(b);
1076-
auto &outside = outside_pack(b);
1077-
1078-
// check for left and right sided reflection
1079-
for (int dir = -1; dir <= 1; dir += 2) {
1080-
if (outside(0, k, j, i) > 0 && outside(0, k, j, i + dir) == 0) {
1081-
PARTHENON_REQUIRE(outside(0, k, j, i + 2 * dir) == 0,
1082-
"Corner case, not sure how to handle corner case");
1083-
PARTHENON_REQUIRE(outside(0, k, j, i - dir) > 0,
1084-
"Corner case, not sure how to handle corner case");
1085-
// mirror all components
1086-
for (int n = 0; n < prim_in.GetDim(4); n++) {
1087-
const bool reflect = n == IV1;
1088-
prim_in(b, n, k, j, i) =
1089-
(reflect ? -1.0 : 1.0) * prim_in(b, n, k, j, i + dir);
1090-
prim_in(b, n, k, j, i - dir) =
1091-
(reflect ? -1.0 : 1.0) * prim_in(b, n, k, j, i + 2 * dir);
1092-
}
1093-
}
1094-
}
1095-
});
1096-
10971070
parthenon::par_for_outer(
10981071
DEFAULT_OUTER_LOOP_PATTERN, "x1 flux", DevExecSpace(), scratch_size_in_bytes,
10991072
scratch_level, 0, cons_in.GetDim(5) - 1, kl, ku, jl, ju,
11001073
KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k, const int j) {
11011074
const auto &prim = prim_in(b);
11021075
auto &cons = cons_in(b);
1076+
auto &outside = outside_pack(b);
11031077
parthenon::ScratchPad2D<Real> wl(member.team_scratch(scratch_level),
11041078
num_scratch_vars, nx1);
11051079
parthenon::ScratchPad2D<Real> wr(member.team_scratch(scratch_level),
@@ -1108,10 +1082,44 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshData<Real>> &md) {
11081082
Reconstruct<recon, X1DIR>(member, k, j, ib.s - 1, ib.e + 1, prim, wl, wr);
11091083
// Sync all threads in the team so that scratch memory is consistent
11101084
member.team_barrier();
1111-
1085+
#if 1
1086+
parthenon::par_for_inner(member, ib.s, ib.e + 1, [&](const int i) {
1087+
// if the cell left of the current reconstructed interface is outside and the
1088+
// right (current)cell inside, replace reconstructed variable on both sides with
1089+
// (reflected) first order versions
1090+
if (outside(0, k, j, i - 1) > 0 && outside(0, k, j, i) == 0) {
1091+
for (auto n = 0; n < prim.GetDim(4); ++n) {
1092+
const bool reflect = n == IV1;
1093+
wl(n, i) = (reflect ? -1.0 : 1.0) * prim(n, k, j, i);
1094+
wr(n, i) = prim(n, k, j, i);
1095+
}
1096+
}
1097+
// if the cell left of the current reconstructed interface is inside and the
1098+
// right (current) cell outside, replace reconstructed variable on both sides
1099+
// with (reflected) first order versions
1100+
if (outside(0, k, j, i) > 0 && outside(0, k, j, i - 1) == 0) {
1101+
for (auto n = 0; n < prim.GetDim(4); ++n) {
1102+
const bool reflect = n == IV1;
1103+
wl(n, i) = prim(n, k, j, i - 1);
1104+
wr(n, i) = (reflect ? -1.0 : 1.0) * prim(n, k, j, i - 1);
1105+
}
1106+
}
1107+
});
1108+
member.team_barrier();
1109+
#endif
11121110
riemann.Solve(member, k, j, ib.s, ib.e + 1, IV1, wl, wr, cons, eos, c_h);
11131111
member.team_barrier();
1114-
1112+
#if 0
1113+
parthenon::par_for_inner(member, ib.s, ib.e + 1, [&](const int i) {
1114+
if (outside(0, k, j, i) > 0 ||
1115+
(outside(0, k, j, i) == 0 && outside(0, k, j, i - 1) > 0)) {
1116+
for (auto n = 0; n < prim.GetDim(4); ++n) {
1117+
cons.flux(IV1, n, k, j, i) = 0;
1118+
}
1119+
}
1120+
});
1121+
member.team_barrier();
1122+
#endif
11151123
// Passive scalar fluxes
11161124
for (auto n = nhydro; n < nhydro + nscalars; ++n) {
11171125
parthenon::par_for_inner(member, ib.s, ib.e + 1, [&](const int i) {
@@ -1136,36 +1144,13 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshData<Real>> &md) {
11361144
else // 3D
11371145
kl = kb.s - 1, ku = kb.e + 1;
11381146

1139-
pmb->par_for(
1140-
"Set x2 reflect", 0, num_blocks - 1, kl, ku, jb.s - 1, jb.e + 1, il, iu,
1141-
KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
1142-
auto &cons = cons_in(b);
1143-
auto &outside = outside_pack(b);
1144-
1145-
// check for left and right sided reflection
1146-
for (int dir = -1; dir <= 1; dir += 2) {
1147-
if (outside(0, k, j, i) > 0 && outside(0, k, j + dir, i) == 0) {
1148-
PARTHENON_REQUIRE(outside(0, k, j + 2 * dir, i) == 0,
1149-
"Corner case, not sure how to handle corner case");
1150-
PARTHENON_REQUIRE(outside(0, k, j - dir, i) > 0,
1151-
"Corner case, not sure how to handle corner case");
1152-
// mirror all components
1153-
for (int n = 0; n < prim_in.GetDim(4); n++) {
1154-
const bool reflect = n == IV2;
1155-
prim_in(b, n, k, j, i) =
1156-
(reflect ? -1.0 : 1.0) * prim_in(b, n, k, j + dir, i);
1157-
prim_in(b, n, k, j - dir, i) =
1158-
(reflect ? -1.0 : 1.0) * prim_in(b, n, k, j + 2 * dir, i);
1159-
}
1160-
}
1161-
}
1162-
});
11631147
parthenon::par_for_outer(
11641148
DEFAULT_OUTER_LOOP_PATTERN, "x2 flux", DevExecSpace(), scratch_size_in_bytes,
11651149
scratch_level, 0, cons_in.GetDim(5) - 1, kl, ku,
11661150
KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k) {
11671151
const auto &prim = prim_in(b);
11681152
auto &cons = cons_in(b);
1153+
auto &outside = outside_pack(b);
11691154
parthenon::ScratchPad2D<Real> wl(member.team_scratch(scratch_level),
11701155
num_scratch_vars, nx1);
11711156
parthenon::ScratchPad2D<Real> wr(member.team_scratch(scratch_level),
@@ -1179,8 +1164,44 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshData<Real>> &md) {
11791164
member.team_barrier();
11801165

11811166
if (j > jb.s - 1) {
1167+
#if 1
1168+
parthenon::par_for_inner(member, il, iu, [&](const int i) {
1169+
// if the cell left of the current reconstructed interface is outside and
1170+
// the right (current)cell inside, replace reconstructed variable on both
1171+
// sides with (reflected) first order versions
1172+
if (outside(0, k, j - 1, i) > 0 && outside(0, k, j, i) == 0) {
1173+
for (auto n = 0; n < prim.GetDim(4); ++n) {
1174+
const bool reflect = n == IV2;
1175+
wl(n, i) = (reflect ? -1.0 : 1.0) * prim(n, k, j, i);
1176+
wr(n, i) = prim(n, k, j, i);
1177+
}
1178+
}
1179+
// if the cell left of the current reconstructed interface is inside and
1180+
// the right (current) cell outside, replace reconstructed variable on
1181+
// both sides with (reflected) first order versions
1182+
if (outside(0, k, j, i) > 0 && outside(0, k, j - 1, i) == 0) {
1183+
for (auto n = 0; n < prim.GetDim(4); ++n) {
1184+
const bool reflect = n == IV2;
1185+
wl(n, i) = prim(n, k, j - 1, i);
1186+
wr(n, i) = (reflect ? -1.0 : 1.0) * prim(n, k, j - 1, i);
1187+
}
1188+
}
1189+
});
1190+
member.team_barrier();
1191+
#endif
11821192
riemann.Solve(member, k, j, il, iu, IV2, wl, wr, cons, eos, c_h);
11831193
member.team_barrier();
1194+
#if 0
1195+
parthenon::par_for_inner(member, il, iu, [&](const int i) {
1196+
if (outside(0, k, j, i) > 0 ||
1197+
(outside(0, k, j, i) == 0 && outside(0, k, j - 1, i) > 0)) {
1198+
for (auto n = 0; n < prim.GetDim(4); ++n) {
1199+
cons.flux(IV2, n, k, j, i) = 0;
1200+
}
1201+
}
1202+
});
1203+
member.team_barrier();
1204+
#endif
11841205

11851206
// Passive scalar fluxes
11861207
for (auto n = nhydro; n < nhydro + nscalars; ++n) {
@@ -1207,37 +1228,14 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshData<Real>> &md) {
12071228
if (pmb->pmy_mesh->ndim >= 3) {
12081229
// set the loop limits
12091230
il = ib.s - 1, iu = ib.e + 1, jl = jb.s - 1, ju = jb.e + 1;
1210-
pmb->par_for(
1211-
"Set x3 reflect", 0, num_blocks - 1, kb.s - 1, kb.e + 1, jl, ju, il, iu,
1212-
KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
1213-
auto &cons = cons_in(b);
1214-
auto &outside = outside_pack(b);
1215-
1216-
// check for left and right sided reflection
1217-
for (int dir = -1; dir <= 1; dir += 2) {
1218-
if (outside(0, k, j, i) > 0 && outside(0, k + dir, j, i) == 0) {
1219-
PARTHENON_REQUIRE(outside(0, k + 2 * dir, j, i) == 0,
1220-
"Corner case, not sure how to handle corner case");
1221-
PARTHENON_REQUIRE(outside(0, k - dir, j, i) > 0,
1222-
"Corner case, not sure how to handle corner case");
1223-
// mirror all components
1224-
for (int n = 0; n < prim_in.GetDim(4); n++) {
1225-
const bool reflect = n == IV3;
1226-
prim_in(b, n, k, j, i) =
1227-
(reflect ? -1.0 : 1.0) * prim_in(b, n, k + dir, j, i);
1228-
prim_in(b, n, k - dir, j, i) =
1229-
(reflect ? -1.0 : 1.0) * prim_in(b, n, k + 2 * dir, j, i);
1230-
}
1231-
}
1232-
}
1233-
});
12341231

12351232
parthenon::par_for_outer(
12361233
DEFAULT_OUTER_LOOP_PATTERN, "x3 flux", DevExecSpace(), scratch_size_in_bytes,
12371234
scratch_level, 0, cons_in.GetDim(5) - 1, jl, ju,
12381235
KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int j) {
12391236
const auto &prim = prim_in(b);
12401237
auto &cons = cons_in(b);
1238+
auto &outside = outside_pack(b);
12411239
parthenon::ScratchPad2D<Real> wl(member.team_scratch(scratch_level),
12421240
num_scratch_vars, nx1);
12431241
parthenon::ScratchPad2D<Real> wr(member.team_scratch(scratch_level),
@@ -1249,11 +1247,47 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshData<Real>> &md) {
12491247
Reconstruct<recon, X3DIR>(member, k, j, il, iu, prim, wlb, wr);
12501248
// Sync all threads in the team so that scratch memory is consistent
12511249
member.team_barrier();
1252-
12531250
if (k > kb.s - 1) {
1251+
#if 1
1252+
parthenon::par_for_inner(member, il, iu, [&](const int i) {
1253+
// if the cell left of the current reconstructed interface is outside and
1254+
// the right (current)cell inside, replace reconstructed variable on both
1255+
// sides with (reflected) first order versions
1256+
if (outside(0, k - 1, j, i) > 0 && outside(0, k, j, i) == 0) {
1257+
for (auto n = 0; n < prim.GetDim(4); ++n) {
1258+
const bool reflect = n == IV3;
1259+
wl(n, i) = (reflect ? -1.0 : 1.0) * prim(n, k, j, i);
1260+
wr(n, i) = prim(n, k, j, i);
1261+
}
1262+
}
1263+
// if the cell left of the current reconstructed interface is inside and
1264+
// the right (current) cell outside, replace reconstructed variable on
1265+
// both sides with (reflected) first order versions
1266+
if (outside(0, k, j, i) > 0 && outside(0, k - 1, j, i) == 0) {
1267+
for (auto n = 0; n < prim.GetDim(4); ++n) {
1268+
const bool reflect = n == IV3;
1269+
wl(n, i) = prim(n, k - 1, j, i);
1270+
wr(n, i) = (reflect ? -1.0 : 1.0) * prim(n, k - 1, j, i);
1271+
}
1272+
}
1273+
});
1274+
member.team_barrier();
1275+
#endif
12541276
riemann.Solve(member, k, j, il, iu, IV3, wl, wr, cons, eos, c_h);
12551277
member.team_barrier();
12561278

1279+
#if 0
1280+
parthenon::par_for_inner(member, il, iu, [&](const int i) {
1281+
if (outside(0, k, j, i) > 0 ||
1282+
(outside(0, k, j, i) == 0 && outside(0, k - 1, j, i) > 0)) {
1283+
for (auto n = 0; n < prim.GetDim(4); ++n) {
1284+
cons.flux(IV3, n, k, j, i) = 0;
1285+
}
1286+
}
1287+
});
1288+
member.team_barrier();
1289+
#endif
1290+
12571291
// Passive scalar fluxes
12581292
for (auto n = nhydro; n < nhydro + nscalars; ++n) {
12591293
parthenon::par_for_inner(member, il, iu, [&](const int i) {
@@ -1333,9 +1367,9 @@ TaskStatus FirstOrderFluxCorrect(MeshData<Real> *u0_data, MeshData<Real> *u1_dat
13331367

13341368
std::int64_t num_corrected, num_need_floor;
13351369
// Potentially need multiple attempts as flux correction corrects 6 (in 3D) fluxes
1336-
// of a single cell at the same time. So the neighboring cells need to be rechecked with
1337-
// the corrected fluxes as the corrected fluxes in one cell may result in the need to
1338-
// correct all the fluxes of an originally "good" neighboring cell.
1370+
// of a single cell at the same time. So the neighboring cells need to be rechecked
1371+
// with the corrected fluxes as the corrected fluxes in one cell may result in the
1372+
// need to correct all the fluxes of an originally "good" neighboring cell.
13391373
size_t num_attempts = 0;
13401374
do {
13411375
num_corrected = 0;

0 commit comments

Comments
 (0)