Nektar++
RoeSolver.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: RoeSolver.h
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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Roe Riemann solver.
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #ifndef NEKTAR_SOLVERS_COMPRESSIBLEFLOWSOLVER_RIEMANNSOLVER_ROESOLVER
36 #define NEKTAR_SOLVERS_COMPRESSIBLEFLOWSOLVER_RIEMANNSOLVER_ROESOLVER
37 
39 
40 namespace Nektar
41 {
43 {
44 public:
47  {
48  return RiemannSolverSharedPtr(new RoeSolver(pSession));
49  }
50 
51  static std::string solverName;
52 
53  /// programmatic ctor
54  RoeSolver();
55 
56 protected:
58 
59  using ND = NekDouble;
60 
61  void v_PointSolve(ND rhoL, ND rhouL, ND rhovL, ND rhowL, ND EL, ND rhoR,
62  ND rhouR, ND rhovR, ND rhowR, ND ER, ND &rhof, ND &rhouf,
63  ND &rhovf, ND &rhowf, ND &Ef) final;
64 
65  void v_ArraySolve(const Array<OneD, const Array<OneD, ND>> &Fwd,
66  const Array<OneD, const Array<OneD, ND>> &Bwd,
67  Array<OneD, Array<OneD, ND>> &flux) final;
68 };
69 
70 template <class T, typename = typename std::enable_if<
71  std::is_floating_point<T>::value ||
73 inline void RoeKernel(T &rhoL, T &rhouL, T &rhovL, T &rhowL, T &EL, T &rhoR,
74  T &rhouR, T &rhovR, T &rhowR, T &ER, T &rhof, T &rhouf,
75  T &rhovf, T &rhowf, T &Ef, NekDouble gamma)
76 {
77  // Left and right velocities
78  T uL = rhouL / rhoL;
79  T vL = rhovL / rhoL;
80  T wL = rhowL / rhoL;
81  T uR = rhouR / rhoR;
82  T vR = rhovR / rhoR;
83  T wR = rhowR / rhoR;
84 
85  // Left and right pressures
86  T pL = (gamma - 1.0) * (EL - 0.5 * (rhouL * uL + rhovL * vL + rhowL * wL));
87  T pR = (gamma - 1.0) * (ER - 0.5 * (rhouR * uR + rhovR * vR + rhowR * wR));
88 
89  // Left and right enthalpy
90  T hL = (EL + pL) / rhoL;
91  T hR = (ER + pR) / rhoR;
92 
93  // Square root of rhoL and rhoR.
94  T srL = sqrt(rhoL);
95  T srR = sqrt(rhoR);
96  T srLR = srL + srR;
97 
98  // Velocity, enthalpy and sound speed Roe averages (equation 11.60).
99  T uRoe = (srL * uL + srR * uR) / srLR;
100  T vRoe = (srL * vL + srR * vR) / srLR;
101  T wRoe = (srL * wL + srR * wR) / srLR;
102  T hRoe = (srL * hL + srR * hR) / srLR;
103  T URoe = (uRoe * uRoe + vRoe * vRoe + wRoe * wRoe);
104  T cRoe = sqrt((gamma - 1.0) * (hRoe - 0.5 * URoe));
105 
106  // Compute eigenvectors (equation 11.59).
107  T k[5][5] = {{1., uRoe - cRoe, vRoe, wRoe, hRoe - uRoe * cRoe},
108  {1., uRoe, vRoe, wRoe, 0.5 * URoe},
109  {0., 0., 1., 0., vRoe},
110  {0., 0., 0., 1., wRoe},
111  {1., uRoe + cRoe, vRoe, wRoe, hRoe + uRoe * cRoe}};
112 
113  // Calculate jumps \Delta u_i (defined preceding equation 11.67).
114  T jump[5] = {rhoR - rhoL, rhouR - rhouL, rhovR - rhovL, rhowR - rhowL,
115  ER - EL};
116 
117  // Define \Delta u_5 (equation 11.70).
118  T jumpbar = jump[4] - (jump[2] - vRoe * jump[0]) * vRoe -
119  (jump[3] - wRoe * jump[0]) * wRoe;
120 
121  // Compute wave amplitudes (equations 11.68, 11.69).
122  T alpha[5];
123  alpha[1] = (gamma - 1.0) *
124  (jump[0] * (hRoe - uRoe * uRoe) + uRoe * jump[1] - jumpbar) /
125  (cRoe * cRoe);
126  alpha[0] =
127  (jump[0] * (uRoe + cRoe) - jump[1] - cRoe * alpha[1]) / (2.0 * cRoe);
128  alpha[4] = jump[0] - (alpha[0] + alpha[1]);
129  alpha[2] = jump[2] - vRoe * jump[0];
130  alpha[3] = jump[3] - wRoe * jump[0];
131 
132  // Compute average of left and right fluxes needed for equation 11.29.
133  rhof = 0.5 * (rhoL * uL + rhoR * uR);
134  rhouf = 0.5 * (pL + rhoL * uL * uL + pR + rhoR * uR * uR);
135  rhovf = 0.5 * (rhoL * uL * vL + rhoR * uR * vR);
136  rhowf = 0.5 * (rhoL * uL * wL + rhoR * uR * wR);
137  Ef = 0.5 * (uL * (EL + pL) + uR * (ER + pR));
138 
139  // Needed to get right overload resolution for std::abs
140  using std::abs;
141 
142  // Compute eigenvalues \lambda_i (equation 11.58).
143  T uRoeAbs = abs(uRoe);
144  T lambda[5] = {abs(uRoe - cRoe), uRoeAbs, uRoeAbs, uRoeAbs,
145  abs(uRoe + cRoe)};
146 
147  // Finally perform summation (11.29).
148  for (size_t i = 0; i < 5; ++i)
149  {
150  uRoeAbs = 0.5 * alpha[i] * lambda[i];
151 
152  rhof -= uRoeAbs * k[i][0];
153  rhouf -= uRoeAbs * k[i][1];
154  rhovf -= uRoeAbs * k[i][2];
155  rhowf -= uRoeAbs * k[i][3];
156  Ef -= uRoeAbs * k[i][4];
157  }
158 }
159 
160 } // namespace Nektar
161 
162 #endif
void v_PointSolve(ND rhoL, ND rhouL, ND rhovL, ND rhowL, ND EL, ND rhoR, ND rhouR, ND rhovR, ND rhowR, ND ER, ND &rhof, ND &rhouf, ND &rhovf, ND &rhowf, ND &Ef) final
Roe Riemann solver.
Definition: RoeSolver.cpp:83
void v_ArraySolve(const Array< OneD, const Array< OneD, ND >> &Fwd, const Array< OneD, const Array< OneD, ND >> &Bwd, Array< OneD, Array< OneD, ND >> &flux) final
Definition: RoeSolver.cpp:95
RoeSolver()
programmatic ctor
Definition: RoeSolver.cpp:52
static RiemannSolverSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Definition: RoeSolver.h:45
static std::string solverName
Definition: RoeSolver.h:51
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< RiemannSolver > RiemannSolverSharedPtr
A shared pointer to an EquationSystem object.
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
void RoeKernel(T &rhoL, T &rhouL, T &rhovL, T &rhowL, T &EL, T &rhoR, T &rhouR, T &rhovR, T &rhowR, T &ER, T &rhof, T &rhouf, T &rhovf, T &rhowf, T &Ef, NekDouble gamma)
Definition: RoeSolver.h:73
double NekDouble
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:295
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291