Nektar++
RoeSolver.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: RoeSolver.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: Roe Riemann solver.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 
38 namespace Nektar
39 {
40  std::string RoeSolver::solverName =
42  "Roe",
44  "Roe Riemann solver");
45 
47  {
48 
49  }
50 
51  /**
52  * @brief Roe Riemann solver.
53  *
54  * Stated equations numbers are from:
55  *
56  * "Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical
57  * Introduction", E. F. Toro (3rd edition, 2009).
58  *
59  * We follow the algorithm prescribed following equation 11.70.
60  *
61  * @param rhoL Density left state.
62  * @param rhoR Density right state.
63  * @param rhouL x-momentum component left state.
64  * @param rhouR x-momentum component right state.
65  * @param rhovL y-momentum component left state.
66  * @param rhovR y-momentum component right state.
67  * @param rhowL z-momentum component left state.
68  * @param rhowR z-momentum component right state.
69  * @param EL Energy left state.
70  * @param ER Energy right state.
71  * @param rhof Computed Riemann flux for density.
72  * @param rhouf Computed Riemann flux for x-momentum component
73  * @param rhovf Computed Riemann flux for y-momentum component
74  * @param rhowf Computed Riemann flux for z-momentum component
75  * @param Ef Computed Riemann flux for energy.
76  */
78  double rhoL, double rhouL, double rhovL, double rhowL, double EL,
79  double rhoR, double rhouR, double rhovR, double rhowR, double ER,
80  double &rhof, double &rhouf, double &rhovf, double &rhowf, double &Ef)
81  {
82  static NekDouble gamma = m_params["gamma"]();
83 
84  // Left and right velocities
85  NekDouble uL = rhouL / rhoL;
86  NekDouble vL = rhovL / rhoL;
87  NekDouble wL = rhowL / rhoL;
88  NekDouble uR = rhouR / rhoR;
89  NekDouble vR = rhovR / rhoR;
90  NekDouble wR = rhowR / rhoR;
91 
92  // Left and right pressures
93  NekDouble pL = (gamma - 1.0) *
94  (EL - 0.5 * (rhouL * uL + rhovL * vL + rhowL * wL));
95  NekDouble pR = (gamma - 1.0) *
96  (ER - 0.5 * (rhouR * uR + rhovR * vR + rhowR * wR));
97 
98  // Left and right enthalpy
99  NekDouble hL = (EL + pL) / rhoL;
100  NekDouble hR = (ER + pR) / rhoR;
101 
102  // Square root of rhoL and rhoR.
103  NekDouble srL = sqrt(rhoL);
104  NekDouble srR = sqrt(rhoR);
105  NekDouble srLR = srL + srR;
106 
107  // Velocity, enthalpy and sound speed Roe averages (equation 11.60).
108  NekDouble uRoe = (srL * uL + srR * uR) / srLR;
109  NekDouble vRoe = (srL * vL + srR * vR) / srLR;
110  NekDouble wRoe = (srL * wL + srR * wR) / srLR;
111  NekDouble hRoe = (srL * hL + srR * hR) / srLR;
112  NekDouble URoe = (uRoe * uRoe + vRoe * vRoe + wRoe * wRoe);
113  NekDouble cRoe = sqrt((gamma - 1.0)*(hRoe - 0.5 * URoe));
114 
115  // Compute eigenvectors (equation 11.59).
116  NekDouble k[5][5] = {
117  {1, uRoe - cRoe, vRoe, wRoe, hRoe - uRoe * cRoe},
118  {1, uRoe, vRoe, wRoe, 0.5 * URoe},
119  {0, 0, 1, 0, vRoe},
120  {0, 0, 0, 1, wRoe},
121  {1, uRoe+cRoe, vRoe, wRoe, hRoe + uRoe*cRoe}
122  };
123 
124  // Calculate jumps \Delta u_i (defined preceding equation 11.67).
125  NekDouble jump[5] = {
126  rhoR - rhoL,
127  rhouR - rhouL,
128  rhovR - rhovL,
129  rhowR - rhowL,
130  ER - EL
131  };
132 
133  // Define \Delta u_5 (equation 11.70).
134  NekDouble jumpbar = jump[4] - (jump[2]-vRoe*jump[0])*vRoe -
135  (jump[3]-wRoe*jump[0])*wRoe;
136 
137  // Compute wave amplitudes (equations 11.68, 11.69).
138  NekDouble alpha[5];
139  alpha[1] = (gamma-1.0)*(jump[0]*(hRoe - uRoe*uRoe) + uRoe*jump[1] -
140  jumpbar)/(cRoe*cRoe);
141  alpha[0] = (jump[0]*(uRoe + cRoe) - jump[1] - cRoe*alpha[1])/(2.0*cRoe);
142  alpha[4] = jump[0] - (alpha[0] + alpha[1]);
143  alpha[2] = jump[2] - vRoe * jump[0];
144  alpha[3] = jump[3] - wRoe * jump[0];
145 
146  // Compute average of left and right fluxes needed for equation 11.29.
147  rhof = 0.5*(rhoL*uL + rhoR*uR);
148  rhouf = 0.5*(pL + rhoL*uL*uL + pR + rhoR*uR*uR);
149  rhovf = 0.5*(rhoL*uL*vL + rhoR*uR*vR);
150  rhowf = 0.5*(rhoL*uL*wL + rhoR*uR*wR);
151  Ef = 0.5*(uL*(EL + pL) + uR*(ER + pR));
152 
153  // Compute eigenvalues \lambda_i (equation 11.58).
154  NekDouble uRoeAbs = fabs(uRoe);
155  NekDouble lambda[5] = {
156  fabs(uRoe - cRoe),
157  uRoeAbs,
158  uRoeAbs,
159  uRoeAbs,
160  fabs(uRoe + cRoe)
161  };
162 
163  // Finally perform summation (11.29).
164  for (int i = 0; i < 5; ++i)
165  {
166  NekDouble ahat = 0.5*alpha[i]*lambda[i];
167  rhof -= ahat*k[i][0];
168  rhouf -= ahat*k[i][1];
169  rhovf -= ahat*k[i][2];
170  rhowf -= ahat*k[i][3];
171  Ef -= ahat*k[i][4];
172  }
173  }
174 }
static std::string solverName
Definition: RoeSolver.h:52
static RiemannSolverSharedPtr create()
Definition: RoeSolver.h:46
virtual void v_PointSolve(double rhoL, double rhouL, double rhovL, double rhowL, double EL, double rhoR, double rhouR, double rhovR, double rhowR, double ER, double &rhof, double &rhouf, double &rhovf, double &rhowf, double &Ef)
Roe Riemann solver.
Definition: RoeSolver.cpp:77
RiemannSolverFactory & GetRiemannSolverFactory()
double NekDouble
std::map< std::string, RSParamFuncType > m_params
Map of parameter function types.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215