Nektar++
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  break;
186 
187  case 2:
188  {
189  const int nq = inarray[0].num_elements();
190  const int vx = (int)vecLocs[i][0];
191  const int vy = (int)vecLocs[i][1];
192 
193  Vmath::Vmul (nq, inarray [vx], 1, normals [0], 1,
194  outarray[vx], 1);
195  Vmath::Vvtvp(nq, inarray [vy], 1, normals [1], 1,
196  outarray[vx], 1, outarray[vx], 1);
197  Vmath::Vmul (nq, inarray [vx], 1, normals [1], 1,
198  outarray[vy], 1);
199  Vmath::Vvtvm(nq, inarray [vy], 1, normals [0], 1,
200  outarray[vy], 1, outarray[vy], 1);
201  break;
202  }
203 
204  case 3:
205  {
206  const int nq = inarray[0].num_elements();
207  const int vx = (int)vecLocs[i][0];
208  const int vy = (int)vecLocs[i][1];
209  const int vz = (int)vecLocs[i][2];
210 
211  // Generate matrices if they don't already exist.
212  if (m_rotMat.num_elements() == 0)
213  {
214  GenerateRotationMatrices(normals);
215  }
216 
217  // Apply rotation matrices.
218  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[0], 1,
219  inarray [vy], 1, m_rotMat[1], 1,
220  outarray[vx], 1);
221  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[2], 1,
222  outarray[vx], 1, outarray[vx], 1);
223  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[3], 1,
224  inarray [vy], 1, m_rotMat[4], 1,
225  outarray[vy], 1);
226  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[5], 1,
227  outarray[vy], 1, outarray[vy], 1);
228  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[6], 1,
229  inarray [vy], 1, m_rotMat[7], 1,
230  outarray[vz], 1);
231  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[8], 1,
232  outarray[vz], 1, outarray[vz], 1);
233  break;
234  }
235 
236  default:
237  ASSERTL1(false, "Invalid space dimension.");
238  break;
239  }
240  }
241  }
242 
243  /**
244  * @brief Rotate a vector field from trace normal.
245  *
246  * This function performs a rotation of the triad of vector components
247  * provided in inarray so that the first component aligns with the
248  * Cartesian components; it performs the inverse operation of
249  * RiemannSolver::rotateToNormal.
250  */
252  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
253  const Array<OneD, const Array<OneD, NekDouble> > &normals,
254  const Array<OneD, const Array<OneD, NekDouble> > &vecLocs,
255  Array<OneD, Array<OneD, NekDouble> > &outarray)
256  {
257  for (int i = 0; i < inarray.num_elements(); ++i)
258  {
259  Vmath::Vcopy(inarray[i].num_elements(), inarray[i], 1,
260  outarray[i], 1);
261  }
262 
263  for (int i = 0; i < vecLocs.num_elements(); i++)
264  {
265  ASSERTL1(vecLocs[i].num_elements() == normals.num_elements(),
266  "vecLocs[i] element count mismatch");
267 
268  switch (normals.num_elements())
269  {
270  case 1:
271  // do nothing
272  break;
273 
274  case 2:
275  {
276  const int nq = normals[0].num_elements();
277  const int vx = (int)vecLocs[i][0];
278  const int vy = (int)vecLocs[i][1];
279 
280  Vmath::Vmul (nq, inarray [vy], 1, normals [1], 1,
281  outarray[vx], 1);
282  Vmath::Vvtvm(nq, inarray [vx], 1, normals [0], 1,
283  outarray[vx], 1, outarray[vx], 1);
284  Vmath::Vmul (nq, inarray [vx], 1, normals [1], 1,
285  outarray[vy], 1);
286  Vmath::Vvtvp(nq, inarray [vy], 1, normals [0], 1,
287  outarray[vy], 1, outarray[vy], 1);
288  break;
289  }
290 
291  case 3:
292  {
293  const int nq = normals[0].num_elements();
294  const int vx = (int)vecLocs[i][0];
295  const int vy = (int)vecLocs[i][1];
296  const int vz = (int)vecLocs[i][2];
297 
298  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[0], 1,
299  inarray [vy], 1, m_rotMat[3], 1,
300  outarray[vx], 1);
301  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[6], 1,
302  outarray[vx], 1, outarray[vx], 1);
303  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[1], 1,
304  inarray [vy], 1, m_rotMat[4], 1,
305  outarray[vy], 1);
306  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[7], 1,
307  outarray[vy], 1, outarray[vy], 1);
308  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[2], 1,
309  inarray [vy], 1, m_rotMat[5], 1,
310  outarray[vz], 1);
311  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[8], 1,
312  outarray[vz], 1, outarray[vz], 1);
313  break;
314  }
315 
316  default:
317  ASSERTL1(false, "Invalid space dimension.");
318  break;
319  }
320  }
321  }
322 
323  /**
324  * @brief Determine whether a scalar has been defined in #m_scalars.
325  *
326  * @param name Scalar name.
327  */
328  bool RiemannSolver::CheckScalars(std::string name)
329  {
331  m_scalars.find(name);
332 
333  return it != m_scalars.end();
334  }
335 
336  /**
337  * @brief Determine whether a vector has been defined in #m_vectors.
338  *
339  * @param name Vector name.
340  */
341  bool RiemannSolver::CheckVectors(std::string name)
342  {
344  m_vectors.find(name);
345 
346  return it != m_vectors.end();
347  }
348 
349  /**
350  * @brief Determine whether a parameter has been defined in #m_params.
351  *
352  * @param name Parameter name.
353  */
354  bool RiemannSolver::CheckParams(std::string name)
355  {
357  m_params.find(name);
358 
359  return it != m_params.end();
360  }
361 
362  /**
363  * @brief Determine whether a scalar has been defined in #m_auxScal.
364  *
365  * @param name Scalar name.
366  */
367  bool RiemannSolver::CheckAuxScal(std::string name)
368  {
370  m_auxScal.find(name);
371 
372  return it != m_auxScal.end();
373  }
374 
375  /**
376  * @brief Determine whether a vector has been defined in #m_auxVec.
377  *
378  * @param name Vector name.
379  */
380  bool RiemannSolver::CheckAuxVec(std::string name)
381  {
383  m_auxVec.find(name);
384 
385  return it != m_auxVec.end();
386  }
387 
388  /**
389  * @brief Generate rotation matrices for 3D expansions.
390  */
392  const Array<OneD, const Array<OneD, NekDouble> > &normals)
393  {
394  Array<OneD, NekDouble> xdir(3,0.0);
395  Array<OneD, NekDouble> tn (3);
396  NekDouble tmp[9];
397  const int nq = normals[0].num_elements();
398  int i, j;
399  xdir[0] = 1.0;
400 
401  // Allocate storage for rotation matrices.
403 
404  for (i = 0; i < 9; ++i)
405  {
407  }
408  for (i = 0; i < normals[0].num_elements(); ++i)
409  {
410  // Generate matrix which takes us from (1,0,0) vector to trace
411  // normal.
412  tn[0] = normals[0][i];
413  tn[1] = normals[1][i];
414  tn[2] = normals[2][i];
415  FromToRotation(tn, xdir, tmp);
416 
417  for (j = 0; j < 9; ++j)
418  {
419  m_rotMat[j][i] = tmp[j];
420  }
421  }
422  }
423 
424  /**
425  * @brief A function for creating a rotation matrix that rotates a
426  * vector @a from into another vector @a to.
427  *
428  * Authors: Tomas Möller, John Hughes
429  * "Efficiently Building a Matrix to Rotate One Vector to
430  * Another" Journal of Graphics Tools, 4(4):1-4, 1999
431  *
432  * @param from Normalised 3-vector to rotate from.
433  * @param to Normalised 3-vector to rotate to.
434  * @param out Resulting 3x3 rotation matrix (row-major order).
435  */
439  NekDouble *mat)
440  {
441  NekDouble v[3];
442  NekDouble e, h, f;
443 
444  CROSS(v, from, to);
445  e = DOT(from, to);
446  f = (e < 0)? -e:e;
447  if (f > 1.0 - EPSILON)
448  {
449  NekDouble u[3], v[3];
450  NekDouble x[3];
451  NekDouble c1, c2, c3;
452  int i, j;
453 
454  x[0] = (from[0] > 0.0)? from[0] : -from[0];
455  x[1] = (from[1] > 0.0)? from[1] : -from[1];
456  x[2] = (from[2] > 0.0)? from[2] : -from[2];
457 
458  if (x[0] < x[1])
459  {
460  if (x[0] < x[2])
461  {
462  x[0] = 1.0; x[1] = x[2] = 0.0;
463  }
464  else
465  {
466  x[2] = 1.0; x[0] = x[1] = 0.0;
467  }
468  }
469  else
470  {
471  if (x[1] < x[2])
472  {
473  x[1] = 1.0; x[0] = x[2] = 0.0;
474  }
475  else
476  {
477  x[2] = 1.0; 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  for (j = 0; j < 3; j++) {
494  mat[3*i+j] = - c1 * u[i] * u[j]
495  - c2 * v[i] * v[j]
496  + c3 * v[i] * u[j];
497  }
498  mat[i+3*i] += 1.0;
499  }
500  }
501  else
502  {
503  NekDouble hvx, hvz, hvxy, hvxz, hvyz;
504  h = 1.0/(1.0 + e);
505  hvx = h * v[0];
506  hvz = h * v[2];
507  hvxy = hvx * v[1];
508  hvxz = hvx * v[2];
509  hvyz = hvz * v[1];
510  mat[0] = e + hvx * v[0];
511  mat[1] = hvxy - v[2];
512  mat[2] = hvxz + v[1];
513  mat[3] = hvxy + v[2];
514  mat[4] = e + h * v[1] * v[1];
515  mat[5] = hvyz - v[0];
516  mat[6] = hvxz - v[1];
517  mat[7] = hvyz + v[0];
518  mat[8] = e + hvz * v[2];
519  }
520  }
521  }
522 }
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:428
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:523
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:451
#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:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
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:169
Provides a generic Factory class.
Definition: NekFactory.hpp:116