Skip to content

Commit 88033b8

Browse files
pgretepgrete
committed
Fix un-inline Riemann (didn't help). Fix bounds.
1 parent 70896cd commit 88033b8

1 file changed

Lines changed: 7 additions & 355 deletions

File tree

‎src/hydro/hydro.cpp‎

Lines changed: 7 additions & 355 deletions
Original file line numberDiff line numberDiff line change
@@ -1142,357 +1142,11 @@ TaskStatus CalculateFluxes(MeshData<Real> *u0_data, MeshData<Real> *u1_data,
11421142
});
11431143
}
11441144
pmb->par_for(
1145-
"x1 Riemann", 0, u0_cons_pack.GetDim(5) - 1, 0, u0_cons_pack.GetDim(4) - 1, kb.s,
1146-
kb.e, jb.s, jb.e, ib.s, ib.e + 1,
1147-
KOKKOS_LAMBDA(const int b, const int n, const int k, const int j, const int i) {
1145+
"x1 Riemann", 0, u0_cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e + 1,
1146+
KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
11481147
auto &wl = u0_wl_pack(b);
11491148
const auto &wr = u0_wr_pack(b);
1150-
1151-
const int ivx = IV1;
1152-
const int ivy = IV1 + ((ivx - IV1) + 1) % 3;
1153-
const int ivz = IV1 + ((ivx - IV1) + 2) % 3;
1154-
const int iBx = ivx - 1 + NHYDRO;
1155-
const int iBy = ivy - 1 + NHYDRO;
1156-
const int iBz = ivz - 1 + NHYDRO;
1157-
1158-
const auto gamma = eos.GetGamma();
1159-
const auto gm1 = gamma - 1.0;
1160-
const auto igm1 = 1.0 / gm1;
1161-
1162-
// TODO(pgrete) move to a more central center and add logic
1163-
constexpr int NGLMMHD = 9;
1164-
1165-
Real wli[NGLMMHD], wri[NGLMMHD];
1166-
Real spd[5]; // signal speeds, left to right
1167-
Cons1D ul, ur; // L/R states, conserved variables (computed)
1168-
Cons1D ulst, uldst, urdst, urst; // Conserved variable for all states
1169-
Cons1D fl, fr; // Fluxes for left & right states
1170-
1171-
//--- Step 1. Load L/R states into local variables
1172-
1173-
wli[IDN] = wl(IDN, k, j, i);
1174-
wli[IV1] = wl(ivx, k, j, i);
1175-
wli[IV2] = wl(ivy, k, j, i);
1176-
wli[IV3] = wl(ivz, k, j, i);
1177-
wli[IPR] = wl(IPR, k, j, i);
1178-
wli[IB1] = wl(iBx, k, j, i);
1179-
wli[IB2] = wl(iBy, k, j, i);
1180-
wli[IB3] = wl(iBz, k, j, i);
1181-
wli[IPS] = wl(IPS, k, j, i);
1182-
1183-
wri[IDN] = wr(IDN, k, j, i);
1184-
wri[IV1] = wr(ivx, k, j, i);
1185-
wri[IV2] = wr(ivy, k, j, i);
1186-
wri[IV3] = wr(ivz, k, j, i);
1187-
wri[IPR] = wr(IPR, k, j, i);
1188-
wri[IB1] = wr(iBx, k, j, i);
1189-
wri[IB2] = wr(iBy, k, j, i);
1190-
wri[IB3] = wr(iBz, k, j, i);
1191-
wri[IPS] = wr(IPS, k, j, i);
1192-
1193-
// first solve the decoupled state, see eq (24) in Mignone & Tzeferacos (2010)
1194-
Real bxi = 0.5 * (wli[IB1] + wri[IB1]) - 0.5 / c_h * (wri[IPS] - wli[IPS]);
1195-
Real psii = 0.5 * (wli[IPS] + wri[IPS]) - 0.5 * c_h * (wri[IB1] - wli[IB1]);
1196-
// and store flux
1197-
wl(iBx, k, j, i) = psii;
1198-
wl(IPS, k, j, i) = SQR(c_h) * bxi;
1199-
1200-
// Compute L/R states for selected conserved variables
1201-
Real bxsq = bxi * bxi;
1202-
// (KGF): group transverse vector components for floating-point associativity
1203-
// symmetry
1204-
Real pbl =
1205-
0.5 * (bxsq + (SQR(wli[IB2]) + SQR(wli[IB3]))); // magnetic pressure (l/r)
1206-
Real pbr = 0.5 * (bxsq + (SQR(wri[IB2]) + SQR(wri[IB3])));
1207-
Real kel = 0.5 * wli[IDN] * (SQR(wli[IV1]) + (SQR(wli[IV2]) + SQR(wli[IV3])));
1208-
Real ker = 0.5 * wri[IDN] * (SQR(wri[IV1]) + (SQR(wri[IV2]) + SQR(wri[IV3])));
1209-
1210-
ul.d = wli[IDN];
1211-
ul.mx = wli[IV1] * ul.d;
1212-
ul.my = wli[IV2] * ul.d;
1213-
ul.mz = wli[IV3] * ul.d;
1214-
ul.e = wli[IPR] * igm1 + kel + pbl;
1215-
ul.by = wli[IB2];
1216-
ul.bz = wli[IB3];
1217-
1218-
ur.d = wri[IDN];
1219-
ur.mx = wri[IV1] * ur.d;
1220-
ur.my = wri[IV2] * ur.d;
1221-
ur.mz = wri[IV3] * ur.d;
1222-
ur.e = wri[IPR] * igm1 + ker + pbr;
1223-
ur.by = wri[IB2];
1224-
ur.bz = wri[IB3];
1225-
1226-
//--- Step 2. Compute L & R wave speeds according to Miyoshi & Kusano, eqn. (67)
1227-
1228-
const auto cfl =
1229-
eos.FastMagnetosonicSpeed(wli[IDN], wli[IPR], wli[IB1], wli[IB2], wli[IB3]);
1230-
const auto cfr =
1231-
eos.FastMagnetosonicSpeed(wri[IDN], wri[IPR], wri[IB1], wri[IB2], wri[IB3]);
1232-
1233-
spd[0] = std::min(wli[IV1] - cfl, wri[IV1] - cfr);
1234-
spd[4] = std::max(wli[IV1] + cfl, wri[IV1] + cfr);
1235-
1236-
// Real cfmax = std::max(cfl,cfr);
1237-
// if (wli[IV1] <= wri[IV1]) {
1238-
// spd[0] = wli[IV1] - cfmax;
1239-
// spd[4] = wri[IV1] + cfmax;
1240-
// } else {
1241-
// spd[0] = wri[IV1] - cfmax;
1242-
// spd[4] = wli[IV1] + cfmax;
1243-
// }
1244-
1245-
//--- Step 3. Compute L/R fluxes
1246-
1247-
Real ptl = wli[IPR] + pbl; // total pressures L,R
1248-
Real ptr = wri[IPR] + pbr;
1249-
1250-
fl.d = ul.mx;
1251-
fl.mx = ul.mx * wli[IV1] + ptl - bxsq;
1252-
fl.my = ul.my * wli[IV1] - bxi * ul.by;
1253-
fl.mz = ul.mz * wli[IV1] - bxi * ul.bz;
1254-
fl.e =
1255-
wli[IV1] * (ul.e + ptl - bxsq) - bxi * (wli[IV2] * ul.by + wli[IV3] * ul.bz);
1256-
fl.by = ul.by * wli[IV1] - bxi * wli[IV2];
1257-
fl.bz = ul.bz * wli[IV1] - bxi * wli[IV3];
1258-
1259-
fr.d = ur.mx;
1260-
fr.mx = ur.mx * wri[IV1] + ptr - bxsq;
1261-
fr.my = ur.my * wri[IV1] - bxi * ur.by;
1262-
fr.mz = ur.mz * wri[IV1] - bxi * ur.bz;
1263-
fr.e =
1264-
wri[IV1] * (ur.e + ptr - bxsq) - bxi * (wri[IV2] * ur.by + wri[IV3] * ur.bz);
1265-
fr.by = ur.by * wri[IV1] - bxi * wri[IV2];
1266-
fr.bz = ur.bz * wri[IV1] - bxi * wri[IV3];
1267-
1268-
//--- Step 4. Compute middle and Alfven wave speeds
1269-
1270-
Real sdl = spd[0] - wli[IV1]; // S_i-u_i (i=L or R)
1271-
Real sdr = spd[4] - wri[IV1];
1272-
1273-
// S_M: eqn (38) of Miyoshi & Kusano
1274-
// (KGF): group ptl, ptr terms for floating-point associativity symmetry
1275-
spd[2] = (sdr * ur.mx - sdl * ul.mx + (ptl - ptr)) / (sdr * ur.d - sdl * ul.d);
1276-
1277-
Real sdml = spd[0] - spd[2]; // S_i-S_M (i=L or R)
1278-
Real sdmr = spd[4] - spd[2];
1279-
Real sdml_inv = 1.0 / sdml;
1280-
Real sdmr_inv = 1.0 / sdmr;
1281-
// eqn (43) of Miyoshi & Kusano
1282-
ulst.d = ul.d * sdl * sdml_inv;
1283-
urst.d = ur.d * sdr * sdmr_inv;
1284-
Real ulst_d_inv = 1.0 / ulst.d;
1285-
Real urst_d_inv = 1.0 / urst.d;
1286-
Real sqrtdl = std::sqrt(ulst.d);
1287-
Real sqrtdr = std::sqrt(urst.d);
1288-
1289-
// eqn (51) of Miyoshi & Kusano
1290-
spd[1] = spd[2] - std::abs(bxi) / sqrtdl;
1291-
spd[3] = spd[2] + std::abs(bxi) / sqrtdr;
1292-
1293-
//--- Step 5. Compute intermediate states
1294-
// eqn (23) explicitly becomes eq (41) of Miyoshi & Kusano
1295-
// TODO(felker): place an assertion that ptstl==ptstr
1296-
Real ptstl = ptl + ul.d * sdl * (spd[2] - wli[IV1]);
1297-
Real ptstr = ptr + ur.d * sdr * (spd[2] - wri[IV1]);
1298-
// Real ptstl = ptl + ul.d*sdl*(sdl-sdml); // these equations had issues when
1299-
// averaged Real ptstr = ptr + ur.d*sdr*(sdr-sdmr);
1300-
Real ptst = 0.5 * (ptstr + ptstl); // total pressure (star state)
1301-
1302-
// ul* - eqn (39) of M&K
1303-
ulst.mx = ulst.d * spd[2];
1304-
if (std::abs(ul.d * sdl * sdml - bxsq) < (SMALL_NUMBER)*ptst) {
1305-
// Degenerate case
1306-
ulst.my = ulst.d * wli[IV2];
1307-
ulst.mz = ulst.d * wli[IV3];
1308-
1309-
ulst.by = ul.by;
1310-
ulst.bz = ul.bz;
1311-
} else {
1312-
// eqns (44) and (46) of M&K
1313-
Real tmp = bxi * (sdl - sdml) / (ul.d * sdl * sdml - bxsq);
1314-
ulst.my = ulst.d * (wli[IV2] - ul.by * tmp);
1315-
ulst.mz = ulst.d * (wli[IV3] - ul.bz * tmp);
1316-
1317-
// eqns (45) and (47) of M&K
1318-
tmp = (ul.d * SQR(sdl) - bxsq) / (ul.d * sdl * sdml - bxsq);
1319-
ulst.by = ul.by * tmp;
1320-
ulst.bz = ul.bz * tmp;
1321-
}
1322-
// v_i* dot B_i*
1323-
// (KGF): group transverse momenta terms for floating-point associativity symmetry
1324-
Real vbstl =
1325-
(ulst.mx * bxi + (ulst.my * ulst.by + ulst.mz * ulst.bz)) * ulst_d_inv;
1326-
// eqn (48) of M&K
1327-
// (KGF): group transverse by, bz terms for floating-point associativity symmetry
1328-
ulst.e =
1329-
(sdl * ul.e - ptl * wli[IV1] + ptst * spd[2] +
1330-
bxi * (wli[IV1] * bxi + (wli[IV2] * ul.by + wli[IV3] * ul.bz) - vbstl)) *
1331-
sdml_inv;
1332-
1333-
// ur* - eqn (39) of M&K
1334-
urst.mx = urst.d * spd[2];
1335-
if (std::abs(ur.d * sdr * sdmr - bxsq) < (SMALL_NUMBER)*ptst) {
1336-
// Degenerate case
1337-
urst.my = urst.d * wri[IV2];
1338-
urst.mz = urst.d * wri[IV3];
1339-
1340-
urst.by = ur.by;
1341-
urst.bz = ur.bz;
1342-
} else {
1343-
// eqns (44) and (46) of M&K
1344-
Real tmp = bxi * (sdr - sdmr) / (ur.d * sdr * sdmr - bxsq);
1345-
urst.my = urst.d * (wri[IV2] - ur.by * tmp);
1346-
urst.mz = urst.d * (wri[IV3] - ur.bz * tmp);
1347-
1348-
// eqns (45) and (47) of M&K
1349-
tmp = (ur.d * SQR(sdr) - bxsq) / (ur.d * sdr * sdmr - bxsq);
1350-
urst.by = ur.by * tmp;
1351-
urst.bz = ur.bz * tmp;
1352-
}
1353-
// v_i* dot B_i*
1354-
// (KGF): group transverse momenta terms for floating-point associativity symmetry
1355-
Real vbstr =
1356-
(urst.mx * bxi + (urst.my * urst.by + urst.mz * urst.bz)) * urst_d_inv;
1357-
// eqn (48) of M&K
1358-
// (KGF): group transverse by, bz terms for floating-point associativity symmetry
1359-
urst.e =
1360-
(sdr * ur.e - ptr * wri[IV1] + ptst * spd[2] +
1361-
bxi * (wri[IV1] * bxi + (wri[IV2] * ur.by + wri[IV3] * ur.bz) - vbstr)) *
1362-
sdmr_inv;
1363-
// ul** and ur** - if Bx is near zero, same as *-states
1364-
if (0.5 * bxsq < (SMALL_NUMBER)*ptst) {
1365-
uldst = ulst;
1366-
urdst = urst;
1367-
} else {
1368-
Real invsumd = 1.0 / (sqrtdl + sqrtdr);
1369-
Real bxsig = (bxi > 0.0 ? 1.0 : -1.0);
1370-
1371-
uldst.d = ulst.d;
1372-
urdst.d = urst.d;
1373-
1374-
uldst.mx = ulst.mx;
1375-
urdst.mx = urst.mx;
1376-
1377-
// eqn (59) of M&K
1378-
Real tmp =
1379-
invsumd * (sqrtdl * (ulst.my * ulst_d_inv) +
1380-
sqrtdr * (urst.my * urst_d_inv) + bxsig * (urst.by - ulst.by));
1381-
uldst.my = uldst.d * tmp;
1382-
urdst.my = urdst.d * tmp;
1383-
1384-
// eqn (60) of M&K
1385-
tmp = invsumd * (sqrtdl * (ulst.mz * ulst_d_inv) +
1386-
sqrtdr * (urst.mz * urst_d_inv) + bxsig * (urst.bz - ulst.bz));
1387-
uldst.mz = uldst.d * tmp;
1388-
urdst.mz = urdst.d * tmp;
1389-
1390-
// eqn (61) of M&K
1391-
tmp = invsumd * (sqrtdl * urst.by + sqrtdr * ulst.by +
1392-
bxsig * sqrtdl * sqrtdr *
1393-
((urst.my * urst_d_inv) - (ulst.my * ulst_d_inv)));
1394-
uldst.by = urdst.by = tmp;
1395-
1396-
// eqn (62) of M&K
1397-
tmp = invsumd * (sqrtdl * urst.bz + sqrtdr * ulst.bz +
1398-
bxsig * sqrtdl * sqrtdr *
1399-
((urst.mz * urst_d_inv) - (ulst.mz * ulst_d_inv)));
1400-
uldst.bz = urdst.bz = tmp;
1401-
1402-
// eqn (63) of M&K
1403-
tmp = spd[2] * bxi + (uldst.my * uldst.by + uldst.mz * uldst.bz) / uldst.d;
1404-
uldst.e = ulst.e - sqrtdl * bxsig * (vbstl - tmp);
1405-
urdst.e = urst.e + sqrtdr * bxsig * (vbstr - tmp);
1406-
}
1407-
1408-
//--- Step 6. Compute flux
1409-
uldst.d = spd[1] * (uldst.d - ulst.d);
1410-
uldst.mx = spd[1] * (uldst.mx - ulst.mx);
1411-
uldst.my = spd[1] * (uldst.my - ulst.my);
1412-
uldst.mz = spd[1] * (uldst.mz - ulst.mz);
1413-
uldst.e = spd[1] * (uldst.e - ulst.e);
1414-
uldst.by = spd[1] * (uldst.by - ulst.by);
1415-
uldst.bz = spd[1] * (uldst.bz - ulst.bz);
1416-
1417-
ulst.d = spd[0] * (ulst.d - ul.d);
1418-
ulst.mx = spd[0] * (ulst.mx - ul.mx);
1419-
ulst.my = spd[0] * (ulst.my - ul.my);
1420-
ulst.mz = spd[0] * (ulst.mz - ul.mz);
1421-
ulst.e = spd[0] * (ulst.e - ul.e);
1422-
ulst.by = spd[0] * (ulst.by - ul.by);
1423-
ulst.bz = spd[0] * (ulst.bz - ul.bz);
1424-
1425-
urdst.d = spd[3] * (urdst.d - urst.d);
1426-
urdst.mx = spd[3] * (urdst.mx - urst.mx);
1427-
urdst.my = spd[3] * (urdst.my - urst.my);
1428-
urdst.mz = spd[3] * (urdst.mz - urst.mz);
1429-
urdst.e = spd[3] * (urdst.e - urst.e);
1430-
urdst.by = spd[3] * (urdst.by - urst.by);
1431-
urdst.bz = spd[3] * (urdst.bz - urst.bz);
1432-
1433-
urst.d = spd[4] * (urst.d - ur.d);
1434-
urst.mx = spd[4] * (urst.mx - ur.mx);
1435-
urst.my = spd[4] * (urst.my - ur.my);
1436-
urst.mz = spd[4] * (urst.mz - ur.mz);
1437-
urst.e = spd[4] * (urst.e - ur.e);
1438-
urst.by = spd[4] * (urst.by - ur.by);
1439-
urst.bz = spd[4] * (urst.bz - ur.bz);
1440-
1441-
if (spd[0] >= 0.0) {
1442-
// return Fl if flow is supersonic
1443-
wl(IDN, k, j, i) = fl.d;
1444-
wl(ivx, k, j, i) = fl.mx;
1445-
wl(ivy, k, j, i) = fl.my;
1446-
wl(ivz, k, j, i) = fl.mz;
1447-
wl(IEN, k, j, i) = fl.e;
1448-
wl(iBy, k, j, i) = fl.by;
1449-
wl(iBz, k, j, i) = fl.bz;
1450-
} else if (spd[4] <= 0.0) {
1451-
// return Fr if flow is supersonic
1452-
wl(IDN, k, j, i) = fr.d;
1453-
wl(ivx, k, j, i) = fr.mx;
1454-
wl(ivy, k, j, i) = fr.my;
1455-
wl(ivz, k, j, i) = fr.mz;
1456-
wl(IEN, k, j, i) = fr.e;
1457-
wl(iBy, k, j, i) = fr.by;
1458-
wl(iBz, k, j, i) = fr.bz;
1459-
} else if (spd[1] >= 0.0) {
1460-
// return Fl*
1461-
wl(IDN, k, j, i) = fl.d + ulst.d;
1462-
wl(ivx, k, j, i) = fl.mx + ulst.mx;
1463-
wl(ivy, k, j, i) = fl.my + ulst.my;
1464-
wl(ivz, k, j, i) = fl.mz + ulst.mz;
1465-
wl(IEN, k, j, i) = fl.e + ulst.e;
1466-
wl(iBy, k, j, i) = fl.by + ulst.by;
1467-
wl(iBz, k, j, i) = fl.bz + ulst.bz;
1468-
} else if (spd[2] >= 0.0) {
1469-
// return Fl**
1470-
wl(IDN, k, j, i) = fl.d + ulst.d + uldst.d;
1471-
wl(ivx, k, j, i) = fl.mx + ulst.mx + uldst.mx;
1472-
wl(ivy, k, j, i) = fl.my + ulst.my + uldst.my;
1473-
wl(ivz, k, j, i) = fl.mz + ulst.mz + uldst.mz;
1474-
wl(IEN, k, j, i) = fl.e + ulst.e + uldst.e;
1475-
wl(iBy, k, j, i) = fl.by + ulst.by + uldst.by;
1476-
wl(iBz, k, j, i) = fl.bz + ulst.bz + uldst.bz;
1477-
} else if (spd[3] > 0.0) {
1478-
// return Fr**
1479-
wl(IDN, k, j, i) = fr.d + urst.d + urdst.d;
1480-
wl(ivx, k, j, i) = fr.mx + urst.mx + urdst.mx;
1481-
wl(ivy, k, j, i) = fr.my + urst.my + urdst.my;
1482-
wl(ivz, k, j, i) = fr.mz + urst.mz + urdst.mz;
1483-
wl(IEN, k, j, i) = fr.e + urst.e + urdst.e;
1484-
wl(iBy, k, j, i) = fr.by + urst.by + urdst.by;
1485-
wl(iBz, k, j, i) = fr.bz + urst.bz + urdst.bz;
1486-
} else {
1487-
// return Fr*
1488-
wl(IDN, k, j, i) = fr.d + urst.d;
1489-
wl(ivx, k, j, i) = fr.mx + urst.mx;
1490-
wl(ivy, k, j, i) = fr.my + urst.my;
1491-
wl(ivz, k, j, i) = fr.mz + urst.mz;
1492-
wl(IEN, k, j, i) = fr.e + urst.e;
1493-
wl(iBy, k, j, i) = fr.by + urst.by;
1494-
wl(iBz, k, j, i) = fr.bz + urst.bz;
1495-
}
1149+
riemann.Solve(k, j, i, IV1, wl, wr, eos, c_h);
14961150
});
14971151
pmb->par_for(
14981152
"UpdateFluxDiv x1", 0, u0_cons_pack.GetDim(5) - 1, 0, u0_cons_pack.GetDim(4) - 1,
@@ -1527,9 +1181,8 @@ TaskStatus CalculateFluxes(MeshData<Real> *u0_data, MeshData<Real> *u1_data,
15271181
});
15281182
}
15291183
pmb->par_for(
1530-
"x2 Riemann", 0, u0_cons_pack.GetDim(5) - 1, 0, u0_cons_pack.GetDim(4) - 1, kb.s,
1531-
kb.e, jb.s, jb.e + 1, ib.s, ib.e,
1532-
KOKKOS_LAMBDA(const int b, const int n, const int k, const int j, const int i) {
1184+
"x2 Riemann", 0, u0_cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e + 1, ib.s,
1185+
ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
15331186
auto &wl = u0_wl_pack(b);
15341187
const auto &wr = u0_wr_pack(b);
15351188
riemann.Solve(k, j, i, IV2, wl, wr, eos, c_h);
@@ -1569,9 +1222,8 @@ TaskStatus CalculateFluxes(MeshData<Real> *u0_data, MeshData<Real> *u1_data,
15691222
});
15701223
}
15711224
pmb->par_for(
1572-
"x3 Riemann", 0, u0_cons_pack.GetDim(5) - 1, 0, u0_cons_pack.GetDim(4) - 1, kb.s,
1573-
kb.e + 1, jb.s, jb.e, ib.s, ib.e,
1574-
KOKKOS_LAMBDA(const int b, const int n, const int k, const int j, const int i) {
1225+
"x3 Riemann", 0, u0_cons_pack.GetDim(5) - 1, kb.s, kb.e + 1, jb.s, jb.e, ib.s,
1226+
ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
15751227
auto &wl = u0_wl_pack(b);
15761228
const auto &wr = u0_wr_pack(b);
15771229
riemann.Solve(k, j, i, IV3, wl, wr, eos, c_h);

0 commit comments

Comments
 (0)