Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 #include <LibUtilities/BasicUtils/VmathArray.hpp>
37 #include <SolverUtils/RiemannSolvers/RiemannSolver.h>
38 
39 #define EPSILON 0.000001
40 
41 #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];}
45 
46 #define DOT(v1, v2) (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2])
47 
48 #define SUB(dest, v1, v2){ \
49  dest[0] = v1[0] - v2[0]; \
50  dest[1] = v1[1] - v2[1]; \
51  dest[2] = v1[2] - v2[2];}
52 
53 using namespace std;
54 
55 namespace Nektar
56 {
57  namespace SolverUtils
58  {
59  /**
60  * Retrieves the singleton instance of the Riemann solver factory.
61  */
62  RiemannSolverFactory& GetRiemannSolverFactory()
63  {
64  typedef Loki::SingletonHolder<RiemannSolverFactory,
65  Loki::CreateUsingNew,
66  Loki::NoDestroy > Type;
67  return Type::Instance();
68  }
69 
70  /**
71  * @class RiemannSolver
72  *
73  * @brief The RiemannSolver class provides an abstract interface under
74  * which solvers for various Riemann problems can be implemented.
75  */
76 
77  RiemannSolver::RiemannSolver() : m_requiresRotation(false),
78  m_rotStorage (3)
79  {
80 
81  }
82 
83  /**
84  * @brief Perform the Riemann solve given the forwards and backwards
85  * spaces.
86  *
87  * This routine calls the virtual function #v_Solve to perform the
88  * Riemann solve. If the flag #m_requiresRotation is set, then the
89  * velocity field is rotated to the normal direction to perform
90  * dimensional splitting, and the resulting fluxes are rotated back to
91  * the Cartesian directions before being returned. For the Rotation to
92  * work, the normal vectors "N" and the location of the vector
93  * components in Fwd "vecLocs"must be set via the SetAuxVec() method.
94  *
95  * @param Fwd Forwards trace space.
96  * @param Bwd Backwards trace space.
97  * @param flux Resultant flux along trace space.
98  */
99  void RiemannSolver::Solve(
100  const int nDim,
101  const Array<OneD, const Array<OneD, NekDouble> > &Fwd,
102  const Array<OneD, const Array<OneD, NekDouble> > &Bwd,
103  Array<OneD, Array<OneD, NekDouble> > &flux)
104  {
105  if (m_requiresRotation)
106  {
107  ASSERTL1(CheckVectors("N"), "N not defined.");
108  ASSERTL1(CheckAuxVec("vecLocs"), "vecLocs not defined.");
109  const Array<OneD, const Array<OneD, NekDouble> > normals =
110  m_vectors["N"]();
111  const Array<OneD, const Array<OneD, NekDouble> > vecLocs =
112  m_auxVec["vecLocs"]();
113 
114  int nFields = Fwd .num_elements();
115  int nPts = Fwd[0].num_elements();
116 
117  if (m_rotStorage[0].num_elements() != nFields ||
118  m_rotStorage[0][0].num_elements() != nPts)
119  {
120  for (int i = 0; i < 3; ++i)
121  {
122  m_rotStorage[i] =
123  Array<OneD, Array<OneD, NekDouble> >(nFields);
124  for (int j = 0; j < nFields; ++j)
125  {
126  m_rotStorage[i][j] = Array<OneD, NekDouble>(nPts);
127  }
128  }
129  }
130 
131  rotateToNormal (Fwd, normals, vecLocs, m_rotStorage[0]);
132  rotateToNormal (Bwd, normals, vecLocs, m_rotStorage[1]);
133  v_Solve (nDim, m_rotStorage[0], m_rotStorage[1],
134  m_rotStorage[2]);
135  rotateFromNormal(m_rotStorage[2], normals, vecLocs, flux);
136  }
137  else
138  {
139  v_Solve(nDim, Fwd, Bwd, flux);
140  }
141  }
142 
143  /**
144  * @brief Rotate a vector field to trace normal.
145  *
146  * This function performs a rotation of a vector so that the first
147  * component aligns with the trace normal direction.
148  *
149  * The vectors components are stored in inarray. Their locations must
150  * be specified in the "vecLocs" array. vecLocs[0] contains the locations
151  * of the first vectors components, vecLocs[1] those of the second and
152  * so on.
153  *
154  * In 2D, this is accomplished through the transform:
155  *
156  * \f[ (u_x, u_y) = (n_x u_x + n_y u_y, -n_x v_x + n_y v_y) \f]
157  *
158  * In 3D, we generate a (non-unique) transformation using
159  * RiemannSolver::fromToRotation.
160  *
161  */
162  void RiemannSolver::rotateToNormal(
163  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
164  const Array<OneD, const Array<OneD, NekDouble> > &normals,
165  const Array<OneD, const Array<OneD, NekDouble> > &vecLocs,
166  Array<OneD, Array<OneD, NekDouble> > &outarray)
167  {
168  for (int i = 0; i < inarray.num_elements(); ++i)
169  {
170  Vmath::Vcopy(inarray[i].num_elements(), inarray[i], 1,
171  outarray[i], 1);
172  }
173 
174  for (int i = 0; i < vecLocs.num_elements(); i++)
175  {
176  ASSERTL1(vecLocs[i].num_elements() == normals.num_elements(),
177  "vecLocs[i] element count mismatch");
178 
179  switch (normals.num_elements())
180  {
181  case 1:
182  // do nothing
183  break;
184 
185  case 2:
186  {
187  const int nq = inarray[0].num_elements();
188  const int vx = (int)vecLocs[i][0];
189  const int vy = (int)vecLocs[i][1];
190 
191  Vmath::Vmul (nq, inarray [vx], 1, normals [0], 1,
192  outarray[vx], 1);
193  Vmath::Vvtvp(nq, inarray [vy], 1, normals [1], 1,
194  outarray[vx], 1, outarray[vx], 1);
195  Vmath::Vmul (nq, inarray [vx], 1, normals [1], 1,
196  outarray[vy], 1);
197  Vmath::Vvtvm(nq, inarray [vy], 1, normals [0], 1,
198  outarray[vy], 1, outarray[vy], 1);
199  break;
200  }
201 
202  case 3:
203  {
204  const int nq = inarray[0].num_elements();
205  const int vx = (int)vecLocs[i][0];
206  const int vy = (int)vecLocs[i][1];
207  const int vz = (int)vecLocs[i][2];
208 
209  // Generate matrices if they don't already exist.
210  if (m_rotMat.num_elements() == 0)
211  {
212  GenerateRotationMatrices(normals);
213  }
214 
215  // Apply rotation matrices.
216  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[0], 1,
217  inarray [vy], 1, m_rotMat[1], 1,
218  outarray[vx], 1);
219  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[2], 1,
220  outarray[vx], 1, outarray[vx], 1);
221  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[3], 1,
222  inarray [vy], 1, m_rotMat[4], 1,
223  outarray[vy], 1);
224  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[5], 1,
225  outarray[vy], 1, outarray[vy], 1);
226  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[6], 1,
227  inarray [vy], 1, m_rotMat[7], 1,
228  outarray[vz], 1);
229  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[8], 1,
230  outarray[vz], 1, outarray[vz], 1);
231  break;
232  }
233 
234  default:
235  ASSERTL1(false, "Invalid space dimension.");
236  break;
237  }
238  }
239  }
240 
241  /**
242  * @brief Rotate a vector field from trace normal.
243  *
244  * This function performs a rotation of the triad of vector components
245  * provided in inarray so that the first component aligns with the
246  * Cartesian components; it performs the inverse operation of
247  * RiemannSolver::rotateToNormal.
248  */
249  void RiemannSolver::rotateFromNormal(
250  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
251  const Array<OneD, const Array<OneD, NekDouble> > &normals,
252  const Array<OneD, const Array<OneD, NekDouble> > &vecLocs,
253  Array<OneD, Array<OneD, NekDouble> > &outarray)
254  {
255  for (int i = 0; i < inarray.num_elements(); ++i)
256  {
257  Vmath::Vcopy(inarray[i].num_elements(), inarray[i], 1,
258  outarray[i], 1);
259  }
260 
261  for (int i = 0; i < vecLocs.num_elements(); i++)
262  {
263  ASSERTL1(vecLocs[i].num_elements() == normals.num_elements(),
264  "vecLocs[i] element count mismatch");
265 
266  switch (normals.num_elements())
267  {
268  case 1:
269  // do nothing
270  break;
271 
272  case 2:
273  {
274  const int nq = normals[0].num_elements();
275  const int vx = (int)vecLocs[i][0];
276  const int vy = (int)vecLocs[i][1];
277 
278  Vmath::Vmul (nq, inarray [vy], 1, normals [1], 1,
279  outarray[vx], 1);
280  Vmath::Vvtvm(nq, inarray [vx], 1, normals [0], 1,
281  outarray[vx], 1, outarray[vx], 1);
282  Vmath::Vmul (nq, inarray [vx], 1, normals [1], 1,
283  outarray[vy], 1);
284  Vmath::Vvtvp(nq, inarray [vy], 1, normals [0], 1,
285  outarray[vy], 1, outarray[vy], 1);
286  break;
287  }
288 
289  case 3:
290  {
291  const int nq = normals[0].num_elements();
292  const int vx = (int)vecLocs[i][0];
293  const int vy = (int)vecLocs[i][1];
294  const int vz = (int)vecLocs[i][2];
295 
296  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[0], 1,
297  inarray [vy], 1, m_rotMat[3], 1,
298  outarray[vx], 1);
299  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[6], 1,
300  outarray[vx], 1, outarray[vx], 1);
301  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[1], 1,
302  inarray [vy], 1, m_rotMat[4], 1,
303  outarray[vy], 1);
304  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[7], 1,
305  outarray[vy], 1, outarray[vy], 1);
306  Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[2], 1,
307  inarray [vy], 1, m_rotMat[5], 1,
308  outarray[vz], 1);
309  Vmath::Vvtvp (nq, inarray [vz], 1, m_rotMat[8], 1,
310  outarray[vz], 1, outarray[vz], 1);
311  break;
312  }
313 
314  default:
315  ASSERTL1(false, "Invalid space dimension.");
316  break;
317  }
318  }
319  }
320 
321  /**
322  * @brief Determine whether a scalar has been defined in #m_scalars.
323  *
324  * @param name Scalar name.
325  */
326  bool RiemannSolver::CheckScalars(std::string name)
327  {
328  std::map<std::string, RSScalarFuncType>::iterator it =
329  m_scalars.find(name);
330 
331  return it != m_scalars.end();
332  }
333 
334  /**
335  * @brief Determine whether a vector has been defined in #m_vectors.
336  *
337  * @param name Vector name.
338  */
339  bool RiemannSolver::CheckVectors(std::string name)
340  {
341  std::map<std::string, RSVecFuncType>::iterator it =
342  m_vectors.find(name);
343 
344  return it != m_vectors.end();
345  }
346 
347  /**
348  * @brief Determine whether a parameter has been defined in #m_params.
349  *
350  * @param name Parameter name.
351  */
352  bool RiemannSolver::CheckParams(std::string name)
353  {
354  std::map<std::string, RSParamFuncType>::iterator it =
355  m_params.find(name);
356 
357  return it != m_params.end();
358  }
359 
360  /**
361  * @brief Determine whether a scalar has been defined in #m_auxScal.
362  *
363  * @param name Scalar name.
364  */
365  bool RiemannSolver::CheckAuxScal(std::string name)
366  {
367  std::map<std::string, RSScalarFuncType>::iterator it =
368  m_auxScal.find(name);
369 
370  return it != m_auxScal.end();
371  }
372 
373  /**
374  * @brief Determine whether a vector has been defined in #m_auxVec.
375  *
376  * @param name Vector name.
377  */
378  bool RiemannSolver::CheckAuxVec(std::string name)
379  {
380  std::map<std::string, RSVecFuncType>::iterator it =
381  m_auxVec.find(name);
382 
383  return it != m_auxVec.end();
384  }
385 
386  /**
387  * @brief Generate rotation matrices for 3D expansions.
388  */
389  void RiemannSolver::GenerateRotationMatrices(
390  const Array<OneD, const Array<OneD, NekDouble> > &normals)
391  {
392  Array<OneD, NekDouble> xdir(3,0.0);
393  Array<OneD, NekDouble> tn (3);
394  NekDouble tmp[9];
395  const int nq = normals[0].num_elements();
396  int i, j;
397  xdir[0] = 1.0;
398 
399  // Allocate storage for rotation matrices.
400  m_rotMat = Array<OneD, Array<OneD, NekDouble> >(9);
401 
402  for (i = 0; i < 9; ++i)
403  {
404  m_rotMat[i] = Array<OneD, NekDouble>(nq);
405  }
406  for (i = 0; i < normals[0].num_elements(); ++i)
407  {
408  // Generate matrix which takes us from (1,0,0) vector to trace
409  // normal.
410  tn[0] = normals[0][i];
411  tn[1] = normals[1][i];
412  tn[2] = normals[2][i];
413  FromToRotation(tn, xdir, tmp);
414 
415  for (j = 0; j < 9; ++j)
416  {
417  m_rotMat[j][i] = tmp[j];
418  }
419  }
420  }
421 
422  /**
423  * @brief A function for creating a rotation matrix that rotates a
424  * vector @a from into another vector @a to.
425  *
426  * Authors: Tomas Möller, John Hughes
427  * "Efficiently Building a Matrix to Rotate One Vector to
428  * Another" Journal of Graphics Tools, 4(4):1-4, 1999
429  *
430  * @param from Normalised 3-vector to rotate from.
431  * @param to Normalised 3-vector to rotate to.
432  * @param out Resulting 3x3 rotation matrix (row-major order).
433  */
434  void RiemannSolver::FromToRotation(
435  Array<OneD, const NekDouble> &from,
436  Array<OneD, const NekDouble> &to,
437  NekDouble *mat)
438  {
439  NekDouble v[3];
440  NekDouble e, h, f;
441 
442  CROSS(v, from, to);
443  e = DOT(from, to);
444  f = (e < 0)? -e:e;
445  if (f > 1.0 - EPSILON)
446  {
447  NekDouble u[3], v[3];
448  NekDouble x[3];
449  NekDouble c1, c2, c3;
450  int i, j;
451 
452  x[0] = (from[0] > 0.0)? from[0] : -from[0];
453  x[1] = (from[1] > 0.0)? from[1] : -from[1];
454  x[2] = (from[2] > 0.0)? from[2] : -from[2];
455 
456  if (x[0] < x[1])
457  {
458  if (x[0] < x[2])
459  {
460  x[0] = 1.0; x[1] = x[2] = 0.0;
461  }
462  else
463  {
464  x[2] = 1.0; x[0] = x[1] = 0.0;
465  }
466  }
467  else
468  {
469  if (x[1] < x[2])
470  {
471  x[1] = 1.0; x[0] = x[2] = 0.0;
472  }
473  else
474  {
475  x[2] = 1.0; x[0] = x[1] = 0.0;
476  }
477  }
478 
479  u[0] = x[0] - from[0];
480  u[1] = x[1] - from[1];
481  u[2] = x[2] - from[2];
482  v[0] = x[0] - to [0];
483  v[1] = x[1] - to [1];
484  v[2] = x[2] - to [2];
485 
486  c1 = 2.0 / DOT(u, u);
487  c2 = 2.0 / DOT(v, v);
488  c3 = c1 * c2 * DOT(u, v);
489 
490  for (i = 0; i < 3; i++) {
491  for (j = 0; j < 3; j++) {
492  mat[3*i+j] = - c1 * u[i] * u[j]
493  - c2 * v[i] * v[j]
494  + c3 * v[i] * u[j];
495  }
496  mat[i+3*i] += 1.0;
497  }
498  }
499  else
500  {
501  NekDouble hvx, hvz, hvxy, hvxz, hvyz;
502  h = 1.0/(1.0 + e);
503  hvx = h * v[0];
504  hvz = h * v[2];
505  hvxy = hvx * v[1];
506  hvxz = hvx * v[2];
507  hvyz = hvz * v[1];
508  mat[0] = e + hvx * v[0];
509  mat[1] = hvxy - v[2];
510  mat[2] = hvxz + v[1];
511  mat[3] = hvxy + v[2];
512  mat[4] = e + h * v[1] * v[1];
513  mat[5] = hvyz - v[0];
514  mat[6] = hvxz - v[1];
515  mat[7] = hvyz + v[0];
516  mat[8] = e + hvz * v[2];
517  }
518  }
519  }
520 }
Generated on Sun Mar 15 2015 00:11:34 for Nektar++ by   doxygen 1.8.1.2