38#define EPSILON 0.000001
40#define CROSS(dest, v1, v2) \
42 dest[0] = v1[1] * v2[2] - v1[2] * v2[1]; \
43 dest[1] = v1[2] * v2[0] - v1[0] * v2[2]; \
44 dest[2] = v1[0] * v2[1] - v1[1] * v2[0]; \
47#define DOT(v1, v2) (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2])
49#define SUB(dest, v1, v2) \
51 dest[0] = v1[0] - v2[0]; \
52 dest[1] = v1[1] - v2[1]; \
53 dest[2] = v1[2] - v2[2]; \
82 : m_requiresRotation(false), m_rotStorage(3)
116 int nFields = Fwd.size();
117 int nPts = Fwd[0].size();
122 for (
int i = 0; i < 3; ++i)
125 for (
int j = 0; j < nFields; ++j)
168 for (
int i = 0; i < inarray.size(); ++i)
170 Vmath::Vcopy(inarray[i].size(), inarray[i], 1, outarray[i], 1);
173 for (
int i = 0; i < vecLocs.size(); i++)
175 ASSERTL1(vecLocs[i].size() == normals.size(),
176 "vecLocs[i] element count mismatch");
178 switch (normals.size())
182 const int nq = inarray[0].size();
183 const int vx = (int)vecLocs[i][0];
184 Vmath::Vmul(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1);
189 const int nq = inarray[0].size();
190 const int vx = (int)vecLocs[i][0];
191 const int vy = (int)vecLocs[i][1];
193 Vmath::Vmul(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1);
194 Vmath::Vvtvp(nq, inarray[vy], 1, normals[1], 1, outarray[vx], 1,
196 Vmath::Vmul(nq, inarray[vx], 1, normals[1], 1, outarray[vy], 1);
197 Vmath::Vvtvm(nq, inarray[vy], 1, normals[0], 1, outarray[vy], 1,
204 const int nq = inarray[0].size();
205 const int vx = (int)vecLocs[i][0];
206 const int vy = (int)vecLocs[i][1];
207 const int vz = (int)vecLocs[i][2];
217 1,
m_rotMat[1], 1, outarray[vx], 1);
221 1,
m_rotMat[4], 1, outarray[vy], 1);
225 1,
m_rotMat[7], 1, outarray[vz], 1);
232 ASSERTL1(
false,
"Invalid space dimension.");
252 for (
int i = 0; i < inarray.size(); ++i)
254 Vmath::Vcopy(inarray[i].size(), inarray[i], 1, outarray[i], 1);
257 for (
int i = 0; i < vecLocs.size(); i++)
259 ASSERTL1(vecLocs[i].size() == normals.size(),
260 "vecLocs[i] element count mismatch");
262 switch (normals.size())
266 const int nq = normals[0].size();
267 const int vx = (int)vecLocs[i][0];
268 Vmath::Vmul(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1);
273 const int nq = normals[0].size();
274 const int vx = (int)vecLocs[i][0];
275 const int vy = (int)vecLocs[i][1];
277 Vmath::Vmul(nq, inarray[vy], 1, normals[1], 1, outarray[vx], 1);
278 Vmath::Vvtvm(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1,
280 Vmath::Vmul(nq, inarray[vx], 1, normals[1], 1, outarray[vy], 1);
281 Vmath::Vvtvp(nq, inarray[vy], 1, normals[0], 1, outarray[vy], 1,
288 const int nq = normals[0].size();
289 const int vx = (int)vecLocs[i][0];
290 const int vy = (int)vecLocs[i][1];
291 const int vz = (int)vecLocs[i][2];
294 1,
m_rotMat[3], 1, outarray[vx], 1);
298 1,
m_rotMat[4], 1, outarray[vy], 1);
302 1,
m_rotMat[5], 1, outarray[vz], 1);
309 ASSERTL1(
false,
"Invalid space dimension.");
374 const int nq = normals[0].size();
381 for (i = 0; i < 9; ++i)
385 for (i = 0; i < normals[0].size(); ++i)
389 tn[0] = normals[0][i];
390 tn[1] = normals[1][i];
391 tn[2] = normals[2][i];
394 for (j = 0; j < 9; ++j)
422 f = (e < 0) ? -e : e;
430 x[0] = (from[0] > 0.0) ? from[0] : -from[0];
431 x[1] = (from[1] > 0.0) ? from[1] : -from[1];
432 x[2] = (from[2] > 0.0) ? from[2] : -from[2];
461 u[0] = x[0] - from[0];
462 u[1] = x[1] - from[1];
463 u[2] = x[2] - from[2];
468 c1 = 2.0 /
DOT(u, u);
469 c2 = 2.0 /
DOT(v, v);
470 c3 = c1 * c2 *
DOT(u, v);
472 for (i = 0; i < 3; i++)
474 for (j = 0; j < 3; j++)
477 -c1 * u[i] * u[j] - c2 * v[i] * v[j] + c3 * v[i] * u[j];
479 mat[i + 3 * i] += 1.0;
491 mat[0] = e + hvx * v[0];
492 mat[1] = hvxy - v[2];
493 mat[2] = hvxz + v[1];
494 mat[3] = hvxy + v[2];
495 mat[4] = e + h * v[1] * v[1];
496 mat[5] = hvyz - v[0];
497 mat[6] = hvxz - v[1];
498 mat[7] = hvyz + v[0];
499 mat[8] = e + hvz * v[2];
515 int nPts = Fwd[0].size();
531 for (
int i = 0; i < nDim; i++)
542 [[maybe_unused]]
const int nDim,
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
#define CROSS(dest, v1, v2)
Provides a generic Factory class.
SOLVER_UTILS_EXPORT bool CheckAuxVec(std::string name)
Determine whether a vector has been defined in m_auxVec.
SOLVER_UTILS_EXPORT void Solve(const int nDim, const Array< OneD, const Array< OneD, NekDouble > > &Fwd, const Array< OneD, const Array< OneD, NekDouble > > &Bwd, Array< OneD, Array< OneD, NekDouble > > &flux)
Perform the Riemann solve given the forwards and backwards spaces.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_rotStorage
Rotation storage.
bool m_requiresRotation
Indicates whether the Riemann solver requires a rotation to be applied to the velocity fields.
virtual SOLVER_UTILS_EXPORT void v_CalcFluxJacobian(const int nDim, const Array< OneD, const Array< OneD, NekDouble > > &Fwd, const Array< OneD, const Array< OneD, NekDouble > > &Bwd, const Array< OneD, const Array< OneD, NekDouble > > &normals, DNekBlkMatSharedPtr &FJac, DNekBlkMatSharedPtr &BJac)
SOLVER_UTILS_EXPORT bool CheckParams(std::string name)
Determine whether a parameter has been defined in m_params.
std::map< std::string, RSScalarFuncType > m_auxScal
Map of auxiliary scalar function types.
SOLVER_UTILS_EXPORT void rotateToNormal(const Array< OneD, const Array< OneD, NekDouble > > &inarray, const Array< OneD, const Array< OneD, NekDouble > > &normals, const Array< OneD, const Array< OneD, NekDouble > > &vecLocs, Array< OneD, Array< OneD, NekDouble > > &outarray)
Rotate a vector field to trace normal.
virtual void v_Solve(const int nDim, const Array< OneD, const Array< OneD, NekDouble > > &Fwd, const Array< OneD, const Array< OneD, NekDouble > > &Bwd, Array< OneD, Array< OneD, NekDouble > > &flux)=0
SOLVER_UTILS_EXPORT void rotateFromNormal(const Array< OneD, const Array< OneD, NekDouble > > &inarray, const Array< OneD, const Array< OneD, NekDouble > > &normals, const Array< OneD, const Array< OneD, NekDouble > > &vecLocs, Array< OneD, Array< OneD, NekDouble > > &outarray)
Rotate a vector field from trace normal.
Array< OneD, Array< OneD, NekDouble > > m_rotMat
Rotation matrices for each trace quadrature point.
SOLVER_UTILS_EXPORT bool CheckScalars(std::string name)
Determine whether a scalar has been defined in m_scalars.
SOLVER_UTILS_EXPORT void CalcFluxJacobian(const int nDim, const Array< OneD, const Array< OneD, NekDouble > > &Fwd, const Array< OneD, const Array< OneD, NekDouble > > &Bwd, DNekBlkMatSharedPtr &FJac, DNekBlkMatSharedPtr &BJac)
Calculate the flux jacobian of Fwd and Bwd.
std::map< std::string, RSVecFuncType > m_vectors
Map of vector function types.
std::map< std::string, RSScalarFuncType > m_scalars
Map of scalar function types.
SOLVER_UTILS_EXPORT bool CheckVectors(std::string name)
Determine whether a vector has been defined in m_vectors.
void FromToRotation(Array< OneD, const NekDouble > &from, Array< OneD, const NekDouble > &to, NekDouble *mat)
A function for creating a rotation matrix that rotates a vector from into another vector to.
SOLVER_UTILS_EXPORT void GenerateRotationMatrices(const Array< OneD, const Array< OneD, NekDouble > > &normals)
Generate rotation matrices for 3D expansions.
SOLVER_UTILS_EXPORT bool CheckAuxScal(std::string name)
Determine whether a scalar has been defined in m_auxScal.
std::map< std::string, RSParamFuncType > m_params
Map of parameter function types.
SOLVER_UTILS_EXPORT RiemannSolver()
std::map< std::string, RSVecFuncType > m_auxVec
Map of auxiliary vector function types.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
RiemannSolverFactory & GetRiemannSolverFactory()
std::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
void Vvtvm(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvm (vector times vector minus vector): z = w*x - y
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)