346 auto riemannSolver = RoeSolverSIMD();
349 riemannSolver.SetParam(
"gamma",
350 [&gamma]() ->
NekDouble & {
return gamma; });
356 Array<OneD, Array<OneD, NekDouble>> vecLocs(1);
357 vecLocs[0] = Array<OneD, NekDouble>(spaceDim);
358 for (
size_t i = 0; i < spaceDim; ++i)
360 vecLocs[0][i] = 1 + i;
362 riemannSolver.SetAuxVec(
364 [&vecLocs]() ->
const Array<OneD,
const Array<OneD, NekDouble>> & {
369 Array<OneD, Array<OneD, NekDouble>> normals(spaceDim);
370 for (
size_t i = 0; i < spaceDim; ++i)
372 normals[i] = Array<OneD, NekDouble>(npts);
374 riemannSolver.SetVector(
375 "N", [&normals]() ->
const Array<OneD,
const Array<OneD, NekDouble>> & {
379 size_t nFields = spaceDim + 2;
380 Array<OneD, Array<OneD, NekDouble>> fwd(nFields), bwd(nFields),
381 flx(nFields), flxRef(nFields);
382 for (
size_t i = 0; i < nFields; ++i)
384 fwd[i] = Array<OneD, NekDouble>(npts);
385 bwd[i] = Array<OneD, NekDouble>(npts);
386 flx[i] = Array<OneD, NekDouble>(npts);
387 flxRef[i] = Array<OneD, NekDouble>(npts);
393 for (
size_t i = 0; i < npts; ++i)
416 rhoe + 0.5 * (rhou * rhou + rhov * rhov + rhow * rhow) / rhoL;
418 rhoe + 0.5 * (rhou * rhou + rhov * rhov + rhow * rhow) / rhoR;
419 fwd[nFields - 1][i] = EL;
420 bwd[nFields - 1][i] = ER;
424 flxRef[0][i] = 0.87858599768171342;
425 flxRef[1][i] = 2.0449028304431223;
426 flxRef[2][i] = 1.8282946712594808;
427 flxRef[3][i] = 2.7424420068892208;
428 flxRef[nFields - 1][i] = 9.8154698039903128;
431 riemannSolver.Solve(spaceDim, fwd, bwd, flx);
434 for (
size_t i = 0; i < npts; ++i)
436 BOOST_CHECK_CLOSE(flxRef[0][i], flx[0][i], 1e-10);
437 BOOST_CHECK_CLOSE(flxRef[1][i], flx[1][i], 1e-10);
438 BOOST_CHECK_CLOSE(flxRef[2][i], flx[2][i], 1e-10);
439 BOOST_CHECK_CLOSE(flxRef[3][i], flx[3][i], 1e-10);
440 BOOST_CHECK_CLOSE(flxRef[4][i], flx[4][i], 1e-10);