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]; \
80 : m_requiresRotation(false), m_rotStorage(3)
114 int nFields = Fwd.size();
115 int nPts = Fwd[0].size();
120 for (
int i = 0; i < 3; ++i)
123 for (
int j = 0; j < nFields; ++j)
141 for (
int i = 0; i < flux.size(); ++i)
143 for (
int j = 0; j < flux[i].size(); ++j)
146 for (
int k = 0; k < nDim; ++k)
148 tmp += N[k][j] * vgt[k][j];
150 flux[i][j] -= 0.5 * (Fwd[i][j] + Bwd[i][j]) * tmp;
186 for (
int i = 0; i < inarray.size(); ++i)
188 Vmath::Vcopy(inarray[i].size(), inarray[i], 1, outarray[i], 1);
191 for (
int i = 0; i < vecLocs.size(); i++)
193 ASSERTL1(vecLocs[i].size() == normals.size(),
194 "vecLocs[i] element count mismatch");
196 switch (normals.size())
200 const int nq = inarray[0].size();
201 const int vx = (int)vecLocs[i][0];
202 Vmath::Vmul(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1);
207 const int nq = inarray[0].size();
208 const int vx = (int)vecLocs[i][0];
209 const int vy = (int)vecLocs[i][1];
211 Vmath::Vmul(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1);
212 Vmath::Vvtvp(nq, inarray[vy], 1, normals[1], 1, outarray[vx], 1,
214 Vmath::Vmul(nq, inarray[vx], 1, normals[1], 1, outarray[vy], 1);
215 Vmath::Vvtvm(nq, inarray[vy], 1, normals[0], 1, outarray[vy], 1,
222 const int nq = inarray[0].size();
223 const int vx = (int)vecLocs[i][0];
224 const int vy = (int)vecLocs[i][1];
225 const int vz = (int)vecLocs[i][2];
235 1,
m_rotMat[1], 1, outarray[vx], 1);
239 1,
m_rotMat[4], 1, outarray[vy], 1);
243 1,
m_rotMat[7], 1, outarray[vz], 1);
250 ASSERTL1(
false,
"Invalid space dimension.");
270 for (
int i = 0; i < inarray.size(); ++i)
272 Vmath::Vcopy(inarray[i].size(), inarray[i], 1, outarray[i], 1);
275 for (
int i = 0; i < vecLocs.size(); i++)
277 ASSERTL1(vecLocs[i].size() == normals.size(),
278 "vecLocs[i] element count mismatch");
280 switch (normals.size())
284 const int nq = normals[0].size();
285 const int vx = (int)vecLocs[i][0];
286 Vmath::Vmul(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1);
291 const int nq = normals[0].size();
292 const int vx = (int)vecLocs[i][0];
293 const int vy = (int)vecLocs[i][1];
295 Vmath::Vmul(nq, inarray[vy], 1, normals[1], 1, outarray[vx], 1);
296 Vmath::Vvtvm(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1,
298 Vmath::Vmul(nq, inarray[vx], 1, normals[1], 1, outarray[vy], 1);
299 Vmath::Vvtvp(nq, inarray[vy], 1, normals[0], 1, outarray[vy], 1,
306 const int nq = normals[0].size();
307 const int vx = (int)vecLocs[i][0];
308 const int vy = (int)vecLocs[i][1];
309 const int vz = (int)vecLocs[i][2];
312 1,
m_rotMat[3], 1, outarray[vx], 1);
316 1,
m_rotMat[4], 1, outarray[vy], 1);
320 1,
m_rotMat[5], 1, outarray[vz], 1);
327 ASSERTL1(
false,
"Invalid space dimension.");
392 const int nq = normals[0].size();
399 for (i = 0; i < 9; ++i)
403 for (i = 0; i < normals[0].size(); ++i)
407 tn[0] = normals[0][i];
408 tn[1] = normals[1][i];
409 tn[2] = normals[2][i];
412 for (j = 0; j < 9; ++j)
440 f = (e < 0) ? -e : e;
448 x[0] = (from[0] > 0.0) ? from[0] : -from[0];
449 x[1] = (from[1] > 0.0) ? from[1] : -from[1];
450 x[2] = (from[2] > 0.0) ? from[2] : -from[2];
479 u[0] = x[0] - from[0];
480 u[1] = x[1] - from[1];
481 u[2] = x[2] - from[2];
486 c1 = 2.0 /
DOT(u, u);
487 c2 = 2.0 /
DOT(v, v);
488 c3 = c1 * c2 *
DOT(u, v);
490 for (i = 0; i < 3; i++)
492 for (j = 0; j < 3; j++)
495 -c1 * u[i] * u[j] - c2 * v[i] * v[j] + c3 * v[i] * u[j];
497 mat[i + 3 * i] += 1.0;
509 mat[0] = e + hvx * v[0];
510 mat[1] = hvxy - v[2];
511 mat[2] = hvxz + v[1];
512 mat[3] = hvxy + v[2];
513 mat[4] = e + h * v[1] * v[1];
514 mat[5] = hvyz - v[0];
515 mat[6] = hvxz - v[1];
516 mat[7] = hvyz + v[0];
517 mat[8] = e + hvz * v[2];
533 int nPts = Fwd[0].size();
549 for (
int i = 0; i < nDim; i++)
560 [[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()
bool m_ALESolver
Flag if using the ALE formulation.
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)