342{
  343    
  344    
  345    
  346    
  347    
  348    
  350    
  352    riemannSolver.SetParam("gamma",
  353                           [&gamma]() -> 
NekDouble & { 
return gamma; });
 
  354 
  355    size_t spaceDim = 3;
  356    size_t npts     = 11; 
  357 
  358    
  361    for (size_t i = 0; i < spaceDim; ++i)
  362    {
  363        vecLocs[0][i] = 1 + i;
  364    }
  365    riemannSolver.SetAuxVec(
  366        "vecLocs",
  367        [&vecLocs]() -> const Array<OneD, const Array<OneD, NekDouble>> & {
  368            return vecLocs;
  369        });
  370 
  371    
  372    Array<OneD, Array<OneD, NekDouble>> normals(spaceDim);
  373    for (size_t i = 0; i < spaceDim; ++i)
  374    {
  375        normals[i] = Array<OneD, NekDouble>(npts);
  376    }
  377    riemannSolver.SetVector(
  378        "N", [&normals]() -> const Array<OneD, const Array<OneD, NekDouble>> & {
  379            return normals;
  380        });
  381 
  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)
  386    {
  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);
  391    }
  392 
  393    
  394    
  395 
  396    for (size_t i = 0; i < npts; ++i)
  397    {
  398        
  401        fwd[0][i]      = rhoL;
  402        bwd[0][i]      = rhoR;
  403        
  405        fwd[1][i]      = rhou;
  406        bwd[1][i]      = rhou;
  407        
  409        fwd[2][i]      = rhov;
  410        bwd[2][i]      = rhov;
  411        
  413        fwd[3][i]      = rhow;
  414        bwd[3][i]      = rhow;
  415        
  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;
  424        
  425        normals[0][i] = 1.0;
  426        
  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;
  432    }
  433 
  434    riemannSolver.Solve(spaceDim, fwd, bwd, flx);
  435 
  436    
  437    for (size_t i = 0; i < npts; ++i)
  438    {
  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);
  444    }
  445}