Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
RiemannSolver.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: RiemannSolver.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Abstract base class for Riemann solvers with factory.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #define LOKI_CLASS_LEVEL_THREADING
39 
40 #define EPSILON 0.000001
41 
42 #define CROSS(dest, v1, v2){ \
43  dest[0] = v1[1] * v2[2] - v1[2] * v2[1]; \
44  dest[1] = v1[2] * v2[0] - v1[0] * v2[2]; \
45  dest[2] = v1[0] * v2[1] - v1[1] * v2[0];}
46 
47 #define DOT(v1, v2) (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2])
48 
49 #define SUB(dest, v1, v2){ \
50  dest[0] = v1[0] - v2[0]; \
51  dest[1] = v1[1] - v2[1]; \
52  dest[2] = v1[2] - v2[2];}
53 
54 using namespace std;
55 
56 namespace Nektar
57 {
58  namespace SolverUtils
59  {
60  /**
61  * Retrieves the singleton instance of the Riemann solver factory.
62  */
64  {
65  typedef Loki::SingletonHolder<RiemannSolverFactory,
66  Loki::CreateUsingNew,
67  Loki::NoDestroy,
68  Loki::SingleThreaded> Type;
69  return Type::Instance();
70  }
71 
72  /**
73  * @class RiemannSolver
74  *
75  * @brief The RiemannSolver class provides an abstract interface under
76  * which solvers for various Riemann problems can be implemented.
77  */
78 
79  RiemannSolver::RiemannSolver() : m_requiresRotation(false),
80  m_rotStorage (3)
81  {
82 
83  }
84 
85  /**
86  * @brief Perform the Riemann solve given the forwards and backwards
87  * spaces.
88  *
89  * This routine calls the virtual function #v_Solve to perform the
90  * Riemann solve. If the flag #m_requiresRotation is set, then the
91  * velocity field is rotated to the normal direction to perform
92  * dimensional splitting, and the resulting fluxes are rotated back to
93  * the Cartesian directions before being returned. For the Rotation to
94  * work, the normal vectors "N" and the location of the vector
95  * components in Fwd "vecLocs"must be set via the SetAuxVec() method.
96  *
97  * @param Fwd Forwards trace space.
98  * @param Bwd Backwards trace space.
99  * @param flux Resultant flux along trace space.
100  */
102  const int nDim,
103  const Array<OneD, const Array<OneD, NekDouble> > &Fwd,
104  const Array<OneD, const Array<OneD, NekDouble> > &Bwd,
106  {
107  if (m_requiresRotation)
108  {
109  ASSERTL1(CheckVectors("N"), "N not defined.");
110  ASSERTL1(CheckAuxVec("vecLocs"), "vecLocs not defined.");
112  m_vectors["N"]();
114  m_auxVec["vecLocs"]();
115 
116  int nFields = Fwd .num_elements();
117  int nPts = Fwd[0].num_elements();
118 
119  if (m_rotStorage[0].num_elements() != nFields ||
120  m_rotStorage[0][0].num_elements() != nPts)
121  {
122  for (int i = 0; i < 3; ++i)
123  {
124  m_rotStorage[i] =
126  for (int j = 0; j < nFields; ++j)
127  {
128  m_rotStorage[i][j] = Array<OneD, NekDouble>(nPts);
129  }
130  }
131  }
132 
133  rotateToNormal (Fwd, normals, vecLocs, m_rotStorage[0]);
134  rotateToNormal (Bwd, normals, vecLocs, m_rotStorage[1]);
135  v_Solve (nDim, m_rotStorage[0], m_rotStorage[1],
136  m_rotStorage[2]);
137  rotateFromNormal(m_rotStorage[2], normals, vecLocs, flux);
138  }
139  else
140  {
141  v_Solve(nDim, Fwd, Bwd, flux);
142  }
143  }
144 
145  /**
146  * @brief Rotate a vector field to trace normal.
147  *
148  * This function performs a rotation of a vector so that the first
149  * component aligns with the trace normal direction.
150  *
151  * The vectors components are stored in inarray. Their locations must
152  * be specified in the "vecLocs" array. vecLocs[0] contains the locations
153  * of the first vectors components, vecLocs[1] those of the second and
154  * so on.
155  *
156  * In 2D, this is accomplished through the transform:
157  *
158  * \f[ (u_x, u_y) = (n_x u_x + n_y u_y, -n_x v_x + n_y v_y) \f]
159  *
160  * In 3D, we generate a (non-unique) transformation using
161  * RiemannSolver::fromToRotation.
162  *
163  */
165  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
166  const Array<OneD, const Array<OneD, NekDouble> > &normals,
167  const Array<OneD, const Array<OneD, NekDouble> > &vecLocs,
168  Array<OneD, Array<OneD, NekDouble> > &outarray)
169  {
170  for (int i = 0; i < inarray.num_elements(); ++i)
171  {
172  Vmath::Vcopy(inarray[i].num_elements(), inarray[i], 1,
173  outarray[i], 1);
174  }
175 
176  for (int i = 0; i < vecLocs.num_elements(); i++)
177  {
178  ASSERTL1(vecLocs[i].num_elements() == normals.num_elements(),
179  "vecLocs[i] element count mismatch");
180 
181  switch (normals.num_elements())
182  {
183  case 1:
184  { // do nothing
185  const int nq = inarray[0].num_elements();
186  const int vx = (int)vecLocs[i][0];
187  Vmath::Vmul (nq, inarray [vx], 1, normals [0], 1,
188  outarray[vx], 1);
189  break;
190  }
191  case 2:
192  {
193  const int nq = inarray[0].num_elements();
194  const int vx = (int)vecLocs[i][0];
195  const int vy = (int)vecLocs[i][1];
196 
197  Vmath::Vmul (nq, inarray [vx], 1, normals [0], 1,
198  outarray[vx], 1);
199  Vmath::Vvtvp(nq, inarray [vy], 1, normals [1], 1,
200  outarray[vx], 1, outarray[vx], 1);
201  Vmath::Vmul (nq, inarray [vx], 1, normals [1], 1,
202  outarray[vy], 1);
203  Vmath::Vvtvm(nq, inarray [vy], 1, normals [0], 1,
204  outarray[vy], 1, outarray[vy], 1);
205  break;
206  }
207 
208  case 3:
209  {
210  const int nq = inarray[0].num_elements();
211  const int vx = (int)vecLocs[i][0];
212  const int vy = (int)vecLocs[i][1];
213  const int vz = (int)vecLocs[i][2];
214 
215  // Generate matrices if they don't already exist.
216  if (m_rotMat.num_elements() == 0)
217  {
218  GenerateRotationMatrices(normals);
219  }
220 
221  // Apply rotation matrices.
222  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[0], 1,
223  inarray [vy], 1, m_rotMat[1], 1,
224  outarray[vx], 1);
225  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[2], 1,
226  outarray[vx], 1, outarray[vx], 1);
227  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[3], 1,
228  inarray [vy], 1, m_rotMat[4], 1,
229  outarray[vy], 1);
230  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[5], 1,
231  outarray[vy], 1, outarray[vy], 1);
232  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[6], 1,
233  inarray [vy], 1, m_rotMat[7], 1,
234  outarray[vz], 1);
235  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[8], 1,
236  outarray[vz], 1, outarray[vz], 1);
237  break;
238  }
239 
240  default:
241  ASSERTL1(false, "Invalid space dimension.");
242  break;
243  }
244  }
245  }
246 
247  /**
248  * @brief Rotate a vector field from trace normal.
249  *
250  * This function performs a rotation of the triad of vector components
251  * provided in inarray so that the first component aligns with the
252  * Cartesian components; it performs the inverse operation of
253  * RiemannSolver::rotateToNormal.
254  */
256  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
257  const Array<OneD, const Array<OneD, NekDouble> > &normals,
258  const Array<OneD, const Array<OneD, NekDouble> > &vecLocs,
259  Array<OneD, Array<OneD, NekDouble> > &outarray)
260  {
261  for (int i = 0; i < inarray.num_elements(); ++i)
262  {
263  Vmath::Vcopy(inarray[i].num_elements(), inarray[i], 1,
264  outarray[i], 1);
265  }
266 
267  for (int i = 0; i < vecLocs.num_elements(); i++)
268  {
269  ASSERTL1(vecLocs[i].num_elements() == normals.num_elements(),
270  "vecLocs[i] element count mismatch");
271 
272  switch (normals.num_elements())
273  {
274  case 1:
275  { // do nothing
276  const int nq = normals[0].num_elements();
277  const int vx = (int)vecLocs[i][0];
278  Vmath::Vmul (nq, inarray [vx], 1, normals [0], 1,
279  outarray[vx], 1);
280  break;
281  }
282  case 2:
283  {
284  const int nq = normals[0].num_elements();
285  const int vx = (int)vecLocs[i][0];
286  const int vy = (int)vecLocs[i][1];
287 
288  Vmath::Vmul (nq, inarray [vy], 1, normals [1], 1,
289  outarray[vx], 1);
290  Vmath::Vvtvm(nq, inarray [vx], 1, normals [0], 1,
291  outarray[vx], 1, outarray[vx], 1);
292  Vmath::Vmul (nq, inarray [vx], 1, normals [1], 1,
293  outarray[vy], 1);
294  Vmath::Vvtvp(nq, inarray [vy], 1, normals [0], 1,
295  outarray[vy], 1, outarray[vy], 1);
296  break;
297  }
298 
299  case 3:
300  {
301  const int nq = normals[0].num_elements();
302  const int vx = (int)vecLocs[i][0];
303  const int vy = (int)vecLocs[i][1];
304  const int vz = (int)vecLocs[i][2];
305 
306  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[0], 1,
307  inarray [vy], 1, m_rotMat[3], 1,
308  outarray[vx], 1);
309  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[6], 1,
310  outarray[vx], 1, outarray[vx], 1);
311  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[1], 1,
312  inarray [vy], 1, m_rotMat[4], 1,
313  outarray[vy], 1);
314  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[7], 1,
315  outarray[vy], 1, outarray[vy], 1);
316  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[2], 1,
317  inarray [vy], 1, m_rotMat[5], 1,
318  outarray[vz], 1);
319  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[8], 1,
320  outarray[vz], 1, outarray[vz], 1);
321  break;
322  }
323 
324  default:
325  ASSERTL1(false, "Invalid space dimension.");
326  break;
327  }
328  }
329  }
330 
331  /**
332  * @brief Determine whether a scalar has been defined in #m_scalars.
333  *
334  * @param name Scalar name.
335  */
336  bool RiemannSolver::CheckScalars(std::string name)
337  {
339  m_scalars.find(name);
340 
341  return it != m_scalars.end();
342  }
343 
344  /**
345  * @brief Determine whether a vector has been defined in #m_vectors.
346  *
347  * @param name Vector name.
348  */
349  bool RiemannSolver::CheckVectors(std::string name)
350  {
352  m_vectors.find(name);
353 
354  return it != m_vectors.end();
355  }
356 
357  /**
358  * @brief Determine whether a parameter has been defined in #m_params.
359  *
360  * @param name Parameter name.
361  */
362  bool RiemannSolver::CheckParams(std::string name)
363  {
365  m_params.find(name);
366 
367  return it != m_params.end();
368  }
369 
370  /**
371  * @brief Determine whether a scalar has been defined in #m_auxScal.
372  *
373  * @param name Scalar name.
374  */
375  bool RiemannSolver::CheckAuxScal(std::string name)
376  {
378  m_auxScal.find(name);
379 
380  return it != m_auxScal.end();
381  }
382 
383  /**
384  * @brief Determine whether a vector has been defined in #m_auxVec.
385  *
386  * @param name Vector name.
387  */
388  bool RiemannSolver::CheckAuxVec(std::string name)
389  {
391  m_auxVec.find(name);
392 
393  return it != m_auxVec.end();
394  }
395 
396  /**
397  * @brief Generate rotation matrices for 3D expansions.
398  */
400  const Array<OneD, const Array<OneD, NekDouble> > &normals)
401  {
402  Array<OneD, NekDouble> xdir(3,0.0);
403  Array<OneD, NekDouble> tn (3);
404  NekDouble tmp[9];
405  const int nq = normals[0].num_elements();
406  int i, j;
407  xdir[0] = 1.0;
408 
409  // Allocate storage for rotation matrices.
411 
412  for (i = 0; i < 9; ++i)
413  {
415  }
416  for (i = 0; i < normals[0].num_elements(); ++i)
417  {
418  // Generate matrix which takes us from (1,0,0) vector to trace
419  // normal.
420  tn[0] = normals[0][i];
421  tn[1] = normals[1][i];
422  tn[2] = normals[2][i];
423  FromToRotation(tn, xdir, tmp);
424 
425  for (j = 0; j < 9; ++j)
426  {
427  m_rotMat[j][i] = tmp[j];
428  }
429  }
430  }
431 
432  /**
433  * @brief A function for creating a rotation matrix that rotates a
434  * vector @a from into another vector @a to.
435  *
436  * Authors: Tomas Möller, John Hughes
437  * "Efficiently Building a Matrix to Rotate One Vector to
438  * Another" Journal of Graphics Tools, 4(4):1-4, 1999
439  *
440  * @param from Normalised 3-vector to rotate from.
441  * @param to Normalised 3-vector to rotate to.
442  * @param out Resulting 3x3 rotation matrix (row-major order).
443  */
447  NekDouble *mat)
448  {
449  NekDouble v[3];
450  NekDouble e, h, f;
451 
452  CROSS(v, from, to);
453  e = DOT(from, to);
454  f = (e < 0)? -e:e;
455  if (f > 1.0 - EPSILON)
456  {
457  NekDouble u[3], v[3];
458  NekDouble x[3];
459  NekDouble c1, c2, c3;
460  int i, j;
461 
462  x[0] = (from[0] > 0.0)? from[0] : -from[0];
463  x[1] = (from[1] > 0.0)? from[1] : -from[1];
464  x[2] = (from[2] > 0.0)? from[2] : -from[2];
465 
466  if (x[0] < x[1])
467  {
468  if (x[0] < x[2])
469  {
470  x[0] = 1.0; x[1] = x[2] = 0.0;
471  }
472  else
473  {
474  x[2] = 1.0; x[0] = x[1] = 0.0;
475  }
476  }
477  else
478  {
479  if (x[1] < x[2])
480  {
481  x[1] = 1.0; x[0] = x[2] = 0.0;
482  }
483  else
484  {
485  x[2] = 1.0; x[0] = x[1] = 0.0;
486  }
487  }
488 
489  u[0] = x[0] - from[0];
490  u[1] = x[1] - from[1];
491  u[2] = x[2] - from[2];
492  v[0] = x[0] - to [0];
493  v[1] = x[1] - to [1];
494  v[2] = x[2] - to [2];
495 
496  c1 = 2.0 / DOT(u, u);
497  c2 = 2.0 / DOT(v, v);
498  c3 = c1 * c2 * DOT(u, v);
499 
500  for (i = 0; i < 3; i++) {
501  for (j = 0; j < 3; j++) {
502  mat[3*i+j] = - c1 * u[i] * u[j]
503  - c2 * v[i] * v[j]
504  + c3 * v[i] * u[j];
505  }
506  mat[i+3*i] += 1.0;
507  }
508  }
509  else
510  {
511  NekDouble hvx, hvz, hvxy, hvxz, hvyz;
512  h = 1.0/(1.0 + e);
513  hvx = h * v[0];
514  hvz = h * v[2];
515  hvxy = hvx * v[1];
516  hvxz = hvx * v[2];
517  hvyz = hvz * v[1];
518  mat[0] = e + hvx * v[0];
519  mat[1] = hvxy - v[2];
520  mat[2] = hvxz + v[1];
521  mat[3] = hvxy + v[2];
522  mat[4] = e + h * v[1] * v[1];
523  mat[5] = hvyz - v[0];
524  mat[6] = hvxz - v[1];
525  mat[7] = hvyz + v[0];
526  mat[8] = e + hvz * v[2];
527  }
528  }
529  }
530 }
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 bool CheckScalars(std::string name)
Determine whether a scalar has been defined in m_scalars.
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.cpp:442
STL namespace.
SOLVER_UTILS_EXPORT bool CheckVectors(std::string name)
Determine whether a vector has been defined in m_vectors.
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.
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
LibUtilities::NekFactory< std::string, RiemannSolver > RiemannSolverFactory
Datatype of the NekFactory used to instantiate classes derived from the RiemannSolver class...
#define CROSS(dest, v1, v2)
std::map< std::string, RSScalarFuncType > m_scalars
Map of scalar function types.
void GenerateRotationMatrices(const Array< OneD, const Array< OneD, NekDouble > > &normals)
Generate rotation matrices for 3D expansions.
RiemannSolverFactory & GetRiemannSolverFactory()
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.
std::map< std::string, RSScalarFuncType > m_auxScal
Map of auxiliary scalar function types.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_rotStorage
Rotation storage.
double NekDouble
SOLVER_UTILS_EXPORT bool CheckAuxScal(std::string name)
Determine whether a scalar has been defined in m_auxScal.
#define EPSILON
Array< OneD, Array< OneD, NekDouble > > m_rotMat
Rotation matrices for each trace quadrature point.
SOLVER_UTILS_EXPORT bool CheckAuxVec(std::string name)
Determine whether a vector has been defined in m_auxVec.
std::map< std::string, RSVecFuncType > m_auxVec
Map of auxiliary vector function types.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
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.cpp:537
bool m_requiresRotation
Indicates whether the Riemann solver requires a rotation to be applied to the velocity fields...
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 plus vector): z = w*x - y
Definition: Vmath.cpp:465
#define DOT(v1, v2)
std::map< std::string, RSParamFuncType > m_params
Map of parameter function types.
SOLVER_UTILS_EXPORT bool CheckParams(std::string name)
Determine whether a parameter has been defined in m_params.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
std::map< std::string, RSVecFuncType > m_vectors
Map of vector 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.
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.cpp:183
Provides a generic Factory class.
Definition: NekFactory.hpp:116