344 auto riemannSolver = RoeSolverSIMD();
347 riemannSolver.SetParam(
"gamma", [&gamma]()
354 Array<OneD, Array<OneD, NekDouble>> vecLocs(1);
355 vecLocs[0] = Array<OneD, NekDouble>(spaceDim);
356 for (
size_t i = 0; i < spaceDim; ++i)
358 vecLocs[0][i] = 1 + i;
360 riemannSolver.SetAuxVec(
"vecLocs", [&vecLocs]()
361 ->
const Array<OneD,
const Array<OneD, NekDouble>>&
365 Array<OneD, Array<OneD, NekDouble>> normals(spaceDim);
366 for (
size_t i = 0; i < spaceDim; ++i)
368 normals[i] = Array<OneD, NekDouble>(npts);
370 riemannSolver.SetVector(
"N", [&normals]()
371 ->
const Array<OneD,
const Array<OneD, NekDouble>>&
374 size_t nFields = spaceDim + 2;
375 Array<OneD, Array<OneD, NekDouble>> fwd(nFields), bwd(nFields),
376 flx(nFields), flxRef(nFields);
377 for (
size_t i = 0; i < nFields; ++i)
379 fwd[i] = Array<OneD, NekDouble>(npts);
380 bwd[i] = Array<OneD, NekDouble>(npts);
381 flx[i] = Array<OneD, NekDouble>(npts);
382 flxRef[i] = Array<OneD, NekDouble>(npts);
388 for (
size_t i = 0; i < npts; ++i)
410 NekDouble EL = rhoe + 0.5*(rhou*rhou + rhov*rhov + rhow*rhow) / rhoL;
411 NekDouble ER = rhoe + 0.5*(rhou*rhou + rhov*rhov + rhow*rhow) / rhoR;
412 fwd[nFields-1][i] = EL;
413 bwd[nFields-1][i] = ER;
417 flxRef[0][i] = 0.87858599768171342;
418 flxRef[1][i] = 2.0449028304431223;
419 flxRef[2][i] = 1.8282946712594808;
420 flxRef[3][i] = 2.7424420068892208;
421 flxRef[nFields-1][i] = 9.8154698039903128;
424 riemannSolver.Solve(spaceDim, fwd, bwd, flx);
427 for (
size_t i = 0; i < npts; ++i)
429 BOOST_CHECK_CLOSE(flxRef[0][i], flx[0][i], 1e-10);
430 BOOST_CHECK_CLOSE(flxRef[1][i], flx[1][i], 1e-10);
431 BOOST_CHECK_CLOSE(flxRef[2][i], flx[2][i], 1e-10);
432 BOOST_CHECK_CLOSE(flxRef[3][i], flx[3][i], 1e-10);
433 BOOST_CHECK_CLOSE(flxRef[4][i], flx[4][i], 1e-10);