349 auto riemannSolver = RoeSolverSIMD();
352 riemannSolver.SetParam(
"gamma",
353 [&gamma]() ->
NekDouble & {
return gamma; });
359 Array<OneD, Array<OneD, NekDouble>> vecLocs(1);
360 vecLocs[0] = Array<OneD, NekDouble>(spaceDim);
361 for (
size_t i = 0; i < spaceDim; ++i)
363 vecLocs[0][i] = 1 + i;
365 riemannSolver.SetAuxVec(
367 [&vecLocs]() ->
const Array<OneD,
const Array<OneD, NekDouble>> & {
372 Array<OneD, Array<OneD, NekDouble>> normals(spaceDim);
373 for (
size_t i = 0; i < spaceDim; ++i)
375 normals[i] = Array<OneD, NekDouble>(npts);
377 riemannSolver.SetVector(
378 "N", [&normals]() ->
const Array<OneD,
const Array<OneD, NekDouble>> & {
382 size_t nFields = spaceDim + 2;
383 Array<OneD, Array<OneD, NekDouble>> fwd(nFields), bwd(nFields),
384 flx(nFields), flxRef(nFields);
385 for (
size_t i = 0; i < nFields; ++i)
387 fwd[i] = Array<OneD, NekDouble>(npts);
388 bwd[i] = Array<OneD, NekDouble>(npts);
389 flx[i] = Array<OneD, NekDouble>(npts);
390 flxRef[i] = Array<OneD, NekDouble>(npts);
396 for (
size_t i = 0; i < npts; ++i)
419 rhoe + 0.5 * (rhou * rhou + rhov * rhov + rhow * rhow) / rhoL;
421 rhoe + 0.5 * (rhou * rhou + rhov * rhov + rhow * rhow) / rhoR;
422 fwd[nFields - 1][i] = EL;
423 bwd[nFields - 1][i] = ER;
427 flxRef[0][i] = 0.87858599768171342;
428 flxRef[1][i] = 2.0449028304431223;
429 flxRef[2][i] = 1.8282946712594808;
430 flxRef[3][i] = 2.7424420068892208;
431 flxRef[nFields - 1][i] = 9.8154698039903128;
434 riemannSolver.Solve(spaceDim, fwd, bwd, flx);
437 for (
size_t i = 0; i < npts; ++i)
439 BOOST_CHECK_CLOSE(flxRef[0][i], flx[0][i], 1e-10);
440 BOOST_CHECK_CLOSE(flxRef[1][i], flx[1][i], 1e-10);
441 BOOST_CHECK_CLOSE(flxRef[2][i], flx[2][i], 1e-10);
442 BOOST_CHECK_CLOSE(flxRef[3][i], flx[3][i], 1e-10);
443 BOOST_CHECK_CLOSE(flxRef[4][i], flx[4][i], 1e-10);