Nektar++
Public Member Functions | Public Attributes | Protected Member Functions | Protected Attributes | List of all members
Nektar::SolverUtils::RiemannSolver Class Referenceabstract

The RiemannSolver class provides an abstract interface under which solvers for various Riemann problems can be implemented. More...

#include <RiemannSolver.h>

Inheritance diagram for Nektar::SolverUtils::RiemannSolver:
[legend]

Public Member Functions

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. More...
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetScalar (std::string name, FuncPointerT func, ObjectPointerT obj)
 
void SetScalar (std::string name, RSScalarFuncType fp)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetVector (std::string name, FuncPointerT func, ObjectPointerT obj)
 
void SetVector (std::string name, RSVecFuncType fp)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetParam (std::string name, FuncPointerT func, ObjectPointerT obj)
 
void SetALEFlag (bool &ALE)
 
void SetParam (std::string name, RSParamFuncType fp)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetAuxScal (std::string name, FuncPointerT func, ObjectPointerT obj)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetAuxVec (std::string name, FuncPointerT func, ObjectPointerT obj)
 
void SetAuxVec (std::string name, RSVecFuncType fp)
 
std::map< std::string, RSScalarFuncType > & GetScalars ()
 
std::map< std::string, RSVecFuncType > & GetVectors ()
 
std::map< std::string, RSParamFuncType > & GetParams ()
 
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. More...
 

Public Attributes

int m_spacedim
 

Protected Member Functions

SOLVER_UTILS_EXPORT RiemannSolver ()
 
SOLVER_UTILS_EXPORT RiemannSolver (const LibUtilities::SessionReaderSharedPtr &pSession)
 
virtual SOLVER_UTILS_EXPORT ~RiemannSolver ()
 
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 GenerateRotationMatrices (const Array< OneD, const Array< OneD, NekDouble > > &normals)
 Generate rotation matrices for 3D expansions. More...
 
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. More...
 
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. More...
 
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. More...
 
SOLVER_UTILS_EXPORT bool CheckScalars (std::string name)
 Determine whether a scalar has been defined in m_scalars. More...
 
SOLVER_UTILS_EXPORT bool CheckVectors (std::string name)
 Determine whether a vector has been defined in m_vectors. More...
 
SOLVER_UTILS_EXPORT bool CheckParams (std::string name)
 Determine whether a parameter has been defined in m_params. More...
 
SOLVER_UTILS_EXPORT bool CheckAuxScal (std::string name)
 Determine whether a scalar has been defined in m_auxScal. More...
 
SOLVER_UTILS_EXPORT bool CheckAuxVec (std::string name)
 Determine whether a vector has been defined in m_auxVec. More...
 
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)
 

Protected Attributes

bool m_requiresRotation
 Indicates whether the Riemann solver requires a rotation to be applied to the velocity fields. More...
 
std::map< std::string, RSScalarFuncTypem_scalars
 Map of scalar function types. More...
 
std::map< std::string, RSVecFuncTypem_vectors
 Map of vector function types. More...
 
std::map< std::string, RSParamFuncTypem_params
 Map of parameter function types. More...
 
std::map< std::string, RSScalarFuncTypem_auxScal
 Map of auxiliary scalar function types. More...
 
std::map< std::string, RSVecFuncTypem_auxVec
 Map of auxiliary vector function types. More...
 
Array< OneD, Array< OneD, NekDouble > > m_rotMat
 Rotation matrices for each trace quadrature point. More...
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_rotStorage
 Rotation storage. More...
 
bool m_ALESolver = false
 Flag if using the ALE formulation. More...
 

Detailed Description

The RiemannSolver class provides an abstract interface under which solvers for various Riemann problems can be implemented.

Definition at line 57 of file RiemannSolver.h.

Constructor & Destructor Documentation

◆ RiemannSolver() [1/2]

Nektar::SolverUtils::RiemannSolver::RiemannSolver ( )
protected

Definition at line 76 of file RiemannSolver.cpp.

77{
78}
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.

◆ RiemannSolver() [2/2]

Nektar::SolverUtils::RiemannSolver::RiemannSolver ( const LibUtilities::SessionReaderSharedPtr pSession)
protected

Definition at line 80 of file RiemannSolver.cpp.

83{
84}

◆ ~RiemannSolver()

virtual SOLVER_UTILS_EXPORT Nektar::SolverUtils::RiemannSolver::~RiemannSolver ( )
inlineprotectedvirtual

Definition at line 167 of file RiemannSolver.h.

167{};

Member Function Documentation

◆ CalcFluxJacobian()

void Nektar::SolverUtils::RiemannSolver::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.

Parameters
FwdForwards trace space.
BwdBackwards trace space.
fluxResultant flux along trace space.

Definition at line 530 of file RiemannSolver.cpp.

534{
535 int nPts = Fwd[0].size();
536
538 {
539 ASSERTL1(CheckVectors("N"), "N not defined.");
540 ASSERTL1(CheckAuxVec("vecLocs"), "vecLocs not defined.");
541 const Array<OneD, const Array<OneD, NekDouble>> normals =
542 m_vectors["N"]();
543 const Array<OneD, const Array<OneD, NekDouble>> vecLocs =
544 m_auxVec["vecLocs"]();
545
546 v_CalcFluxJacobian(nDim, Fwd, Bwd, normals, FJac, BJac);
547 }
548 else
549 {
550 Array<OneD, Array<OneD, NekDouble>> normals(nDim);
551 for (int i = 0; i < nDim; i++)
552 {
553 normals[i] = Array<OneD, NekDouble>(nPts, 0.0);
554 }
555 Vmath::Fill(nPts, 1.0, normals[0], 1);
556
557 v_CalcFluxJacobian(nDim, Fwd, Bwd, normals, FJac, BJac);
558 }
559}
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
SOLVER_UTILS_EXPORT bool CheckAuxVec(std::string name)
Determine whether a vector has been defined in m_auxVec.
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)
std::map< std::string, RSVecFuncType > m_vectors
Map of vector function types.
SOLVER_UTILS_EXPORT bool CheckVectors(std::string name)
Determine whether a vector has been defined in m_vectors.
std::map< std::string, RSVecFuncType > m_auxVec
Map of auxiliary vector function types.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54

References ASSERTL1, CheckAuxVec(), CheckVectors(), Vmath::Fill(), m_auxVec, m_requiresRotation, m_vectors, and v_CalcFluxJacobian().

◆ CheckAuxScal()

bool Nektar::SolverUtils::RiemannSolver::CheckAuxScal ( std::string  name)
protected

Determine whether a scalar has been defined in m_auxScal.

Parameters
nameScalar name.

Definition at line 370 of file RiemannSolver.cpp.

371{
372 return m_auxScal.find(name) != m_auxScal.end();
373}
std::map< std::string, RSScalarFuncType > m_auxScal
Map of auxiliary scalar function types.

References m_auxScal, and CellMLToNektar.pycml::name.

◆ CheckAuxVec()

bool Nektar::SolverUtils::RiemannSolver::CheckAuxVec ( std::string  name)
protected

Determine whether a vector has been defined in m_auxVec.

Parameters
nameVector name.

Definition at line 380 of file RiemannSolver.cpp.

381{
382 return m_auxVec.find(name) != m_auxVec.end();
383}

References m_auxVec, and CellMLToNektar.pycml::name.

Referenced by CalcFluxJacobian(), and Solve().

◆ CheckParams()

bool Nektar::SolverUtils::RiemannSolver::CheckParams ( std::string  name)
protected

Determine whether a parameter has been defined in m_params.

Parameters
nameParameter name.

Definition at line 360 of file RiemannSolver.cpp.

361{
362 return m_params.find(name) != m_params.end();
363}
std::map< std::string, RSParamFuncType > m_params
Map of parameter function types.

References m_params, and CellMLToNektar.pycml::name.

Referenced by Nektar::UpwindPulseSolver::RiemannSolverUpwind().

◆ CheckScalars()

bool Nektar::SolverUtils::RiemannSolver::CheckScalars ( std::string  name)
protected

Determine whether a scalar has been defined in m_scalars.

Parameters
nameScalar name.

Definition at line 340 of file RiemannSolver.cpp.

341{
342 return m_scalars.find(name) != m_scalars.end();
343}
std::map< std::string, RSScalarFuncType > m_scalars
Map of scalar function types.

References m_scalars, and CellMLToNektar.pycml::name.

Referenced by Nektar::SolverUtils::UpwindSolver::v_Solve(), and Nektar::UpwindPulseSolver::v_Solve().

◆ CheckVectors()

bool Nektar::SolverUtils::RiemannSolver::CheckVectors ( std::string  name)
protected

Determine whether a vector has been defined in m_vectors.

Parameters
nameVector name.

Definition at line 350 of file RiemannSolver.cpp.

351{
352 return m_vectors.find(name) != m_vectors.end();
353}

References m_vectors, and CellMLToNektar.pycml::name.

Referenced by CalcFluxJacobian(), Nektar::AcousticSolver::GetRotBasefield(), Solve(), and Nektar::RoeSolverSIMD::v_Solve().

◆ FromToRotation()

void Nektar::SolverUtils::RiemannSolver::FromToRotation ( Array< OneD, const NekDouble > &  from,
Array< OneD, const NekDouble > &  to,
NekDouble mat 
)
protected

A function for creating a rotation matrix that rotates a vector from into another vector to.

Authors: Tomas Möller, John Hughes "Efficiently Building a Matrix to Rotate One Vector to Another" Journal of Graphics Tools, 4(4):1-4, 1999

Parameters
fromNormalised 3-vector to rotate from.
toNormalised 3-vector to rotate to.
outResulting 3x3 rotation matrix (row-major order).

Definition at line 433 of file RiemannSolver.cpp.

436{
437 NekDouble v[3];
438 NekDouble e, h, f;
439
440 CROSS(v, from, to);
441 e = DOT(from, to);
442 f = (e < 0) ? -e : e;
443 if (f > 1.0 - EPSILON)
444 {
445 NekDouble u[3], v[3];
446 NekDouble x[3];
447 NekDouble c1, c2, c3;
448 int i, j;
449
450 x[0] = (from[0] > 0.0) ? from[0] : -from[0];
451 x[1] = (from[1] > 0.0) ? from[1] : -from[1];
452 x[2] = (from[2] > 0.0) ? from[2] : -from[2];
453
454 if (x[0] < x[1])
455 {
456 if (x[0] < x[2])
457 {
458 x[0] = 1.0;
459 x[1] = x[2] = 0.0;
460 }
461 else
462 {
463 x[2] = 1.0;
464 x[0] = x[1] = 0.0;
465 }
466 }
467 else
468 {
469 if (x[1] < x[2])
470 {
471 x[1] = 1.0;
472 x[0] = x[2] = 0.0;
473 }
474 else
475 {
476 x[2] = 1.0;
477 x[0] = x[1] = 0.0;
478 }
479 }
480
481 u[0] = x[0] - from[0];
482 u[1] = x[1] - from[1];
483 u[2] = x[2] - from[2];
484 v[0] = x[0] - to[0];
485 v[1] = x[1] - to[1];
486 v[2] = x[2] - to[2];
487
488 c1 = 2.0 / DOT(u, u);
489 c2 = 2.0 / DOT(v, v);
490 c3 = c1 * c2 * DOT(u, v);
491
492 for (i = 0; i < 3; i++)
493 {
494 for (j = 0; j < 3; j++)
495 {
496 mat[3 * i + j] =
497 -c1 * u[i] * u[j] - c2 * v[i] * v[j] + c3 * v[i] * u[j];
498 }
499 mat[i + 3 * i] += 1.0;
500 }
501 }
502 else
503 {
504 NekDouble hvx, hvz, hvxy, hvxz, hvyz;
505 h = 1.0 / (1.0 + e);
506 hvx = h * v[0];
507 hvz = h * v[2];
508 hvxy = hvx * v[1];
509 hvxz = hvx * v[2];
510 hvyz = hvz * v[1];
511 mat[0] = e + hvx * v[0];
512 mat[1] = hvxy - v[2];
513 mat[2] = hvxz + v[1];
514 mat[3] = hvxy + v[2];
515 mat[4] = e + h * v[1] * v[1];
516 mat[5] = hvyz - v[0];
517 mat[6] = hvxz - v[1];
518 mat[7] = hvyz + v[0];
519 mat[8] = e + hvz * v[2];
520 }
521}
#define EPSILON
#define CROSS(dest, v1, v2)
#define DOT(v1, v2)
double NekDouble

References CROSS, DOT, and EPSILON.

Referenced by GenerateRotationMatrices().

◆ GenerateRotationMatrices()

void Nektar::SolverUtils::RiemannSolver::GenerateRotationMatrices ( const Array< OneD, const Array< OneD, NekDouble > > &  normals)
protected

Generate rotation matrices for 3D expansions.

Definition at line 388 of file RiemannSolver.cpp.

390{
391 Array<OneD, NekDouble> xdir(3, 0.0);
392 Array<OneD, NekDouble> tn(3);
393 NekDouble tmp[9];
394 const int nq = normals[0].size();
395 int i, j;
396 xdir[0] = 1.0;
397
398 // Allocate storage for rotation matrices.
399 m_rotMat = Array<OneD, Array<OneD, NekDouble>>(9);
400
401 for (i = 0; i < 9; ++i)
402 {
403 m_rotMat[i] = Array<OneD, NekDouble>(nq);
404 }
405 for (i = 0; i < normals[0].size(); ++i)
406 {
407 // Generate matrix which takes us from (1,0,0) vector to trace
408 // normal.
409 tn[0] = normals[0][i];
410 tn[1] = normals[1][i];
411 tn[2] = normals[2][i];
412 FromToRotation(tn, xdir, tmp);
413
414 for (j = 0; j < 9; ++j)
415 {
416 m_rotMat[j][i] = tmp[j];
417 }
418 }
419}
Array< OneD, Array< OneD, NekDouble > > m_rotMat
Rotation matrices for each trace quadrature point.
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.

References FromToRotation(), and m_rotMat.

Referenced by rotateToNormal(), and Nektar::RoeSolverSIMD::v_Solve().

◆ GetParams()

std::map< std::string, RSParamFuncType > & Nektar::SolverUtils::RiemannSolver::GetParams ( )
inline

Definition at line 130 of file RiemannSolver.h.

131 {
132 return m_params;
133 }

References m_params.

◆ GetScalars()

std::map< std::string, RSScalarFuncType > & Nektar::SolverUtils::RiemannSolver::GetScalars ( )
inline

Definition at line 120 of file RiemannSolver.h.

121 {
122 return m_scalars;
123 }

References m_scalars.

◆ GetVectors()

std::map< std::string, RSVecFuncType > & Nektar::SolverUtils::RiemannSolver::GetVectors ( )
inline

Definition at line 125 of file RiemannSolver.h.

126 {
127 return m_vectors;
128 }

References m_vectors.

◆ rotateFromNormal()

void Nektar::SolverUtils::RiemannSolver::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 
)
protected

Rotate a vector field from trace normal.

This function performs a rotation of the triad of vector components provided in inarray so that the first component aligns with the Cartesian components; it performs the inverse operation of RiemannSolver::rotateToNormal.

Definition at line 266 of file RiemannSolver.cpp.

271{
272 for (int i = 0; i < inarray.size(); ++i)
273 {
274 Vmath::Vcopy(inarray[i].size(), inarray[i], 1, outarray[i], 1);
275 }
276
277 for (int i = 0; i < vecLocs.size(); i++)
278 {
279 ASSERTL1(vecLocs[i].size() == normals.size(),
280 "vecLocs[i] element count mismatch");
281
282 switch (normals.size())
283 {
284 case 1:
285 { // do nothing
286 const int nq = normals[0].size();
287 const int vx = (int)vecLocs[i][0];
288 Vmath::Vmul(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1);
289 break;
290 }
291 case 2:
292 {
293 const int nq = normals[0].size();
294 const int vx = (int)vecLocs[i][0];
295 const int vy = (int)vecLocs[i][1];
296
297 Vmath::Vmul(nq, inarray[vy], 1, normals[1], 1, outarray[vx], 1);
298 Vmath::Vvtvm(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1,
299 outarray[vx], 1);
300 Vmath::Vmul(nq, inarray[vx], 1, normals[1], 1, outarray[vy], 1);
301 Vmath::Vvtvp(nq, inarray[vy], 1, normals[0], 1, outarray[vy], 1,
302 outarray[vy], 1);
303 break;
304 }
305
306 case 3:
307 {
308 const int nq = normals[0].size();
309 const int vx = (int)vecLocs[i][0];
310 const int vy = (int)vecLocs[i][1];
311 const int vz = (int)vecLocs[i][2];
312
313 Vmath::Vvtvvtp(nq, inarray[vx], 1, m_rotMat[0], 1, inarray[vy],
314 1, m_rotMat[3], 1, outarray[vx], 1);
315 Vmath::Vvtvp(nq, inarray[vz], 1, m_rotMat[6], 1, outarray[vx],
316 1, outarray[vx], 1);
317 Vmath::Vvtvvtp(nq, inarray[vx], 1, m_rotMat[1], 1, inarray[vy],
318 1, m_rotMat[4], 1, outarray[vy], 1);
319 Vmath::Vvtvp(nq, inarray[vz], 1, m_rotMat[7], 1, outarray[vy],
320 1, outarray[vy], 1);
321 Vmath::Vvtvvtp(nq, inarray[vx], 1, m_rotMat[2], 1, inarray[vy],
322 1, m_rotMat[5], 1, outarray[vz], 1);
323 Vmath::Vvtvp(nq, inarray[vz], 1, m_rotMat[8], 1, outarray[vz],
324 1, outarray[vz], 1);
325 break;
326 }
327
328 default:
329 ASSERTL1(false, "Invalid space dimension.");
330 break;
331 }
332 }
333}
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.
Definition: Vmath.hpp:72
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
Definition: Vmath.hpp:366
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
Definition: Vmath.hpp:381
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):
Definition: Vmath.hpp:439
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825

References ASSERTL1, m_rotMat, Vmath::Vcopy(), Vmath::Vmul(), Vmath::Vvtvm(), Vmath::Vvtvp(), and Vmath::Vvtvvtp().

Referenced by Solve().

◆ rotateToNormal()

void Nektar::SolverUtils::RiemannSolver::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 
)
protected

Rotate a vector field to trace normal.

This function performs a rotation of a vector so that the first component aligns with the trace normal direction.

The vectors components are stored in inarray. Their locations must be specified in the "vecLocs" array. vecLocs[0] contains the locations of the first vectors components, vecLocs[1] those of the second and so on.

In 2D, this is accomplished through the transform:

\[ (u_x, u_y) = (n_x u_x + n_y u_y, -n_x v_x + n_y v_y) \]

In 3D, we generate a (non-unique) transformation using RiemannSolver::fromToRotation.

Definition at line 182 of file RiemannSolver.cpp.

187{
188 for (int i = 0; i < inarray.size(); ++i)
189 {
190 Vmath::Vcopy(inarray[i].size(), inarray[i], 1, outarray[i], 1);
191 }
192
193 for (int i = 0; i < vecLocs.size(); i++)
194 {
195 ASSERTL1(vecLocs[i].size() == normals.size(),
196 "vecLocs[i] element count mismatch");
197
198 switch (normals.size())
199 {
200 case 1:
201 { // do nothing
202 const int nq = inarray[0].size();
203 const int vx = (int)vecLocs[i][0];
204 Vmath::Vmul(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1);
205 break;
206 }
207 case 2:
208 {
209 const int nq = inarray[0].size();
210 const int vx = (int)vecLocs[i][0];
211 const int vy = (int)vecLocs[i][1];
212
213 Vmath::Vmul(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1);
214 Vmath::Vvtvp(nq, inarray[vy], 1, normals[1], 1, outarray[vx], 1,
215 outarray[vx], 1);
216 Vmath::Vmul(nq, inarray[vx], 1, normals[1], 1, outarray[vy], 1);
217 Vmath::Vvtvm(nq, inarray[vy], 1, normals[0], 1, outarray[vy], 1,
218 outarray[vy], 1);
219 break;
220 }
221
222 case 3:
223 {
224 const int nq = inarray[0].size();
225 const int vx = (int)vecLocs[i][0];
226 const int vy = (int)vecLocs[i][1];
227 const int vz = (int)vecLocs[i][2];
228
229 // Generate matrices if they don't already exist.
230 if (m_rotMat.size() == 0)
231 {
233 }
234
235 // Apply rotation matrices.
236 Vmath::Vvtvvtp(nq, inarray[vx], 1, m_rotMat[0], 1, inarray[vy],
237 1, m_rotMat[1], 1, outarray[vx], 1);
238 Vmath::Vvtvp(nq, inarray[vz], 1, m_rotMat[2], 1, outarray[vx],
239 1, outarray[vx], 1);
240 Vmath::Vvtvvtp(nq, inarray[vx], 1, m_rotMat[3], 1, inarray[vy],
241 1, m_rotMat[4], 1, outarray[vy], 1);
242 Vmath::Vvtvp(nq, inarray[vz], 1, m_rotMat[5], 1, outarray[vy],
243 1, outarray[vy], 1);
244 Vmath::Vvtvvtp(nq, inarray[vx], 1, m_rotMat[6], 1, inarray[vy],
245 1, m_rotMat[7], 1, outarray[vz], 1);
246 Vmath::Vvtvp(nq, inarray[vz], 1, m_rotMat[8], 1, outarray[vz],
247 1, outarray[vz], 1);
248 break;
249 }
250
251 default:
252 ASSERTL1(false, "Invalid space dimension.");
253 break;
254 }
255 }
256}
SOLVER_UTILS_EXPORT void GenerateRotationMatrices(const Array< OneD, const Array< OneD, NekDouble > > &normals)
Generate rotation matrices for 3D expansions.

References ASSERTL1, GenerateRotationMatrices(), m_rotMat, Vmath::Vcopy(), Vmath::Vmul(), Vmath::Vvtvm(), Vmath::Vvtvp(), and Vmath::Vvtvvtp().

Referenced by Nektar::AcousticSolver::GetRotBasefield(), and Solve().

◆ SetALEFlag()

void Nektar::SolverUtils::RiemannSolver::SetALEFlag ( bool &  ALE)
inline

Definition at line 93 of file RiemannSolver.h.

94 {
95 m_ALESolver = ALE;
96 }
bool m_ALESolver
Flag if using the ALE formulation.

References m_ALESolver.

◆ SetAuxScal()

template<typename FuncPointerT , typename ObjectPointerT >
void Nektar::SolverUtils::RiemannSolver::SetAuxScal ( std::string  name,
FuncPointerT  func,
ObjectPointerT  obj 
)
inline

Definition at line 104 of file RiemannSolver.h.

105 {
106 m_auxScal[name] = std::bind(func, obj);
107 }

References m_auxScal, and CellMLToNektar.pycml::name.

◆ SetAuxVec() [1/2]

template<typename FuncPointerT , typename ObjectPointerT >
void Nektar::SolverUtils::RiemannSolver::SetAuxVec ( std::string  name,
FuncPointerT  func,
ObjectPointerT  obj 
)
inline

Definition at line 110 of file RiemannSolver.h.

111 {
112 m_auxVec[name] = std::bind(func, obj);
113 }

References m_auxVec, and CellMLToNektar.pycml::name.

◆ SetAuxVec() [2/2]

void Nektar::SolverUtils::RiemannSolver::SetAuxVec ( std::string  name,
RSVecFuncType  fp 
)
inline

Definition at line 115 of file RiemannSolver.h.

116 {
117 m_auxVec[name] = fp;
118 }

References m_auxVec, and CellMLToNektar.pycml::name.

◆ SetParam() [1/2]

template<typename FuncPointerT , typename ObjectPointerT >
void Nektar::SolverUtils::RiemannSolver::SetParam ( std::string  name,
FuncPointerT  func,
ObjectPointerT  obj 
)
inline

Definition at line 88 of file RiemannSolver.h.

89 {
90 m_params[name] = std::bind(func, obj);
91 }

References m_params, and CellMLToNektar.pycml::name.

◆ SetParam() [2/2]

void Nektar::SolverUtils::RiemannSolver::SetParam ( std::string  name,
RSParamFuncType  fp 
)
inline

Definition at line 98 of file RiemannSolver.h.

99 {
100 m_params[name] = fp;
101 }

References m_params, and CellMLToNektar.pycml::name.

◆ SetScalar() [1/2]

template<typename FuncPointerT , typename ObjectPointerT >
void Nektar::SolverUtils::RiemannSolver::SetScalar ( std::string  name,
FuncPointerT  func,
ObjectPointerT  obj 
)
inline

Definition at line 66 of file RiemannSolver.h.

67 {
68 m_scalars[name] = std::bind(func, obj);
69 }

References m_scalars, and CellMLToNektar.pycml::name.

◆ SetScalar() [2/2]

void Nektar::SolverUtils::RiemannSolver::SetScalar ( std::string  name,
RSScalarFuncType  fp 
)
inline

Definition at line 71 of file RiemannSolver.h.

72 {
73 m_scalars[name] = fp;
74 }

References m_scalars, and CellMLToNektar.pycml::name.

◆ SetVector() [1/2]

template<typename FuncPointerT , typename ObjectPointerT >
void Nektar::SolverUtils::RiemannSolver::SetVector ( std::string  name,
FuncPointerT  func,
ObjectPointerT  obj 
)
inline

Definition at line 77 of file RiemannSolver.h.

78 {
79 m_vectors[name] = std::bind(func, obj);
80 }

References m_vectors, and CellMLToNektar.pycml::name.

◆ SetVector() [2/2]

void Nektar::SolverUtils::RiemannSolver::SetVector ( std::string  name,
RSVecFuncType  fp 
)
inline

Definition at line 82 of file RiemannSolver.h.

83 {
84 m_vectors[name] = fp;
85 }

References m_vectors, and CellMLToNektar.pycml::name.

◆ Solve()

void Nektar::SolverUtils::RiemannSolver::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.

This routine calls the virtual function v_Solve to perform the Riemann solve. If the flag m_requiresRotation is set, then the velocity field is rotated to the normal direction to perform dimensional splitting, and the resulting fluxes are rotated back to the Cartesian directions before being returned. For the Rotation to work, the normal vectors "N" and the location of the vector components in Fwd "vecLocs"must be set via the SetAuxVec() method.

Parameters
FwdForwards trace space.
BwdBackwards trace space.
fluxResultant flux along trace space.

Definition at line 102 of file RiemannSolver.cpp.

106{
108 {
109 ASSERTL1(CheckVectors("N"), "N not defined.");
110 ASSERTL1(CheckAuxVec("vecLocs"), "vecLocs not defined.");
111 const Array<OneD, const Array<OneD, NekDouble>> normals =
112 m_vectors["N"]();
113 const Array<OneD, const Array<OneD, NekDouble>> vecLocs =
114 m_auxVec["vecLocs"]();
115
116 int nFields = Fwd.size();
117 int nPts = Fwd[0].size();
118
119 if (m_rotStorage[0].size() != nFields ||
120 m_rotStorage[0][0].size() != nPts)
121 {
122 for (int i = 0; i < 3; ++i)
123 {
124 m_rotStorage[i] = Array<OneD, Array<OneD, NekDouble>>(nFields);
125 for (int j = 0; j < nFields; ++j)
126 {
127 m_rotStorage[i][j] = Array<OneD, NekDouble>(nPts);
128 }
129 }
130 }
131
132 rotateToNormal(Fwd, normals, vecLocs, m_rotStorage[0]);
133 rotateToNormal(Bwd, normals, vecLocs, m_rotStorage[1]);
134
136 rotateFromNormal(m_rotStorage[2], normals, vecLocs, flux);
137
138 // Attempt to subtract (\vec{U}\vec{vg})\dot n for ALE
139 if (m_ALESolver)
140 {
141 auto vgt = m_vectors["vgt"]();
142 auto N = m_vectors["N"]();
143 for (int i = 0; i < flux.size(); ++i)
144 {
145 for (int j = 0; j < flux[i].size(); ++j)
146 {
147 NekDouble tmp = 0;
148 for (int k = 0; k < nDim; ++k)
149 {
150 tmp += N[k][j] * vgt[k][j];
151 }
152 flux[i][j] -= 0.5 * (Fwd[i][j] + Bwd[i][j]) * tmp;
153 }
154 }
155 }
156 }
157 else
158 {
159 v_Solve(nDim, Fwd, Bwd, flux);
160 }
161}
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.

References ASSERTL1, CheckAuxVec(), CheckVectors(), m_ALESolver, m_auxVec, m_requiresRotation, m_rotStorage, m_vectors, rotateFromNormal(), rotateToNormal(), and v_Solve().

◆ v_CalcFluxJacobian()

void Nektar::SolverUtils::RiemannSolver::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 
)
protectedvirtual

Definition at line 561 of file RiemannSolver.cpp.

568{
569 NEKERROR(ErrorUtil::efatal, "v_CalcFluxJacobian not specified.");
570}
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202

References Nektar::ErrorUtil::efatal, and NEKERROR.

Referenced by CalcFluxJacobian().

◆ v_Solve()

virtual void Nektar::SolverUtils::RiemannSolver::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 
)
protectedpure virtual

Member Data Documentation

◆ m_ALESolver

bool Nektar::SolverUtils::RiemannSolver::m_ALESolver = false
protected

Flag if using the ALE formulation.

Definition at line 161 of file RiemannSolver.h.

Referenced by SetALEFlag(), and Solve().

◆ m_auxScal

std::map<std::string, RSScalarFuncType> Nektar::SolverUtils::RiemannSolver::m_auxScal
protected

Map of auxiliary scalar function types.

Definition at line 153 of file RiemannSolver.h.

Referenced by CheckAuxScal(), and SetAuxScal().

◆ m_auxVec

std::map<std::string, RSVecFuncType> Nektar::SolverUtils::RiemannSolver::m_auxVec
protected

Map of auxiliary vector function types.

Definition at line 155 of file RiemannSolver.h.

Referenced by CalcFluxJacobian(), CheckAuxVec(), SetAuxVec(), and Solve().

◆ m_params

std::map<std::string, RSParamFuncType> Nektar::SolverUtils::RiemannSolver::m_params
protected

◆ m_requiresRotation

bool Nektar::SolverUtils::RiemannSolver::m_requiresRotation
protected

◆ m_rotMat

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::RiemannSolver::m_rotMat
protected

Rotation matrices for each trace quadrature point.

Definition at line 157 of file RiemannSolver.h.

Referenced by GenerateRotationMatrices(), rotateFromNormal(), rotateToNormal(), and Nektar::RoeSolverSIMD::v_Solve().

◆ m_rotStorage

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::RiemannSolver::m_rotStorage
protected

Rotation storage.

Definition at line 159 of file RiemannSolver.h.

Referenced by Solve().

◆ m_scalars

std::map<std::string, RSScalarFuncType> Nektar::SolverUtils::RiemannSolver::m_scalars
protected

◆ m_spacedim

int Nektar::SolverUtils::RiemannSolver::m_spacedim

Definition at line 135 of file RiemannSolver.h.

◆ m_vectors

std::map<std::string, RSVecFuncType> Nektar::SolverUtils::RiemannSolver::m_vectors
protected