Nektar++
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ExactSolverToro.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ExactSolverToro.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 // 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: Exact Riemann solver (Toro 2009).
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
36 
37 #define TOL 1e-6
38 #define NRITER 20
39 
40 namespace Nektar
41 {
42 std::string ExactSolverToro::solverName =
44  "ExactToro", ExactSolverToro::create, "Exact Riemann solver");
45 
48  : CompressibleSolver(pSession)
49 {
50 }
51 
52 /**
53  * @brief Use either PVRS, two-rarefaction or two-shock Riemann solvers to
54  * calculate an initial pressure for the Newton-Raphson scheme.
55  *
56  * @param g Array of calculated gamma values.
57  * @param rhoL Density left state.
58  * @param rhoR Density right state.
59  * @param uL x-velocity component left state.
60  * @param uR x-velocity component right state.
61  * @param pL Pressure component left state.
62  * @param pR Pressure component right state.
63  * @param cL Sound speed component left state.
64  * @param cR Sound speed component right state.
65  * @return Computed initial guess for the Newton-Raphson scheme.
66  */
68  NekDouble pL, NekDouble cL, NekDouble rhoR,
69  NekDouble uR, NekDouble pR, NekDouble cR)
70 {
71  const NekDouble quser = 2.0;
72  NekDouble cup, ppv, pmin, pmax, qmax;
73 
74  cup = 0.25 * (rhoL + rhoR) * (cL + cR);
75  ppv = 0.5 * (pL + pR) + 0.5 * (uL - uR) * cup;
76  ppv = std::max(0.0, ppv);
77  pmin = std::min(pL, pR);
78  pmax = std::max(pL, pR);
79  qmax = pmax / pmin;
80 
81  if (qmax <= quser && pmin <= ppv && ppv <= pmax)
82  {
83  // Select PVRS Riemann solver.
84  return ppv;
85  }
86  else if (ppv < pmin)
87  {
88  // Select two-rarefaction Riemann solver.
89  NekDouble pq = pow(pL / pR, g[1]);
90  NekDouble um =
91  (pq * uL / cL + uR / cR + g[4] * (pq - 1.0)) / (pq / cL + 1.0 / cR);
92  NekDouble ptL = 1.0 + g[7] * (uL - um) / cL;
93  NekDouble ptR = 1.0 + g[7] * (um - uR) / cR;
94  return 0.5 * (pL * pow(ptL, g[3]) + pR * pow(ptR, g[3]));
95  }
96  else
97  {
98  // Select two-shock Riemann solver with PVRS as estimate.
99  NekDouble geL = sqrt((g[5] / rhoL) / (g[6] * pL + ppv));
100  NekDouble geR = sqrt((g[5] / rhoR) / (g[6] * pR + ppv));
101  return (geL * pL + geR * pR - (uR - uL)) / (geL + geR);
102  }
103 }
104 
105 /**
106  * @brief Evaluate pressure functions fL and fR in Newton iteration of
107  * Riemann solver (see equation 4.85 and 4.86 of Toro 2009).
108  *
109  * @param g Array of gamma parameters.
110  * @param p Pressure at current iteration.
111  * @param dk Density (left or right state)
112  * @param pk Pressure (left or right state)
113  * @param ck Sound speed (left or right state)
114  * @param f Computed pressure (shock).
115  * @param fd Computed pressure (rarefaction).
116  */
117 inline void prefun(NekDouble *g, NekDouble p, NekDouble dk, NekDouble pk,
118  NekDouble ck, NekDouble &f, NekDouble &fd)
119 {
120  if (p <= pk)
121  {
122  // Rarefaction wave
123  NekDouble prat = p / pk;
124  f = g[4] * ck * (pow(prat, g[1]) - 1.0);
125  fd = pow(prat, -g[2]) / (dk * ck);
126  }
127  else
128  {
129  // Shock wave
130  NekDouble ak = g[5] / dk;
131  NekDouble bk = g[6] * pk;
132  NekDouble qrt = sqrt(ak / (bk + p));
133  f = (p - pk) * qrt;
134  fd = (1.0 - 0.5 * (p - pk) / (bk + p)) * qrt;
135  }
136 }
137 
138 /**
139  * @brief Exact Riemann solver for the Euler equations.
140  *
141  * This algorithm is transcribed from:
142  *
143  * "Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical
144  * Introduction", E. F. Toro (3rd edition, 2009).
145  *
146  * The full Fortran 77 routine can be found at the end of chapter 4 (section
147  * 4.9). This transcription is essentially the functions STARPU and SAMPLE
148  * glued together, and variable names are kept mostly the same. See the
149  * preceding chapter which explains the derivation of the solver. The
150  * routines PREFUN and GUESSP are kept separate and are reproduced above.
151  *
152  * @param rhoL Density left state.
153  * @param rhoR Density right state.
154  * @param rhouL x-momentum component left state.
155  * @param rhouR x-momentum component right state.
156  * @param rhovL y-momentum component left state.
157  * @param rhovR y-momentum component right state.
158  * @param rhowL z-momentum component left state.
159  * @param rhowR z-momentum component right state.
160  * @param EL Energy left state.
161  * @param ER Energy right state.
162  * @param rhof Computed Riemann flux for density.
163  * @param rhouf Computed Riemann flux for x-momentum component
164  * @param rhovf Computed Riemann flux for y-momentum component
165  * @param rhowf Computed Riemann flux for z-momentum component
166  * @param Ef Computed Riemann flux for energy.
167  */
169  NekDouble rhoL, NekDouble rhouL, NekDouble rhovL, NekDouble rhowL,
170  NekDouble EL, NekDouble rhoR, NekDouble rhouR, NekDouble rhovR,
171  NekDouble rhowR, NekDouble ER, NekDouble &rhof, NekDouble &rhouf,
172  NekDouble &rhovf, NekDouble &rhowf, NekDouble &Ef)
173 {
174  static NekDouble gamma = m_params["gamma"]();
175 
176  // Left and right variables.
177  NekDouble uL = rhouL / rhoL;
178  NekDouble vL = rhovL / rhoL;
179  NekDouble wL = rhowL / rhoL;
180  NekDouble uR = rhouR / rhoR;
181  NekDouble vR = rhovR / rhoR;
182  NekDouble wR = rhowR / rhoR;
183 
184  // Left and right pressure.
185  NekDouble pL =
186  (gamma - 1.0) * (EL - 0.5 * (rhouL * uL + rhovL * vL + rhowL * wL));
187  NekDouble pR =
188  (gamma - 1.0) * (ER - 0.5 * (rhouR * uR + rhovR * vR + rhowR * wR));
189 
190  // Compute gammas.
191  NekDouble g[] = {gamma,
192  (gamma - 1.0) / (2.0 * gamma),
193  (gamma + 1.0) / (2.0 * gamma),
194  2.0 * gamma / (gamma - 1.0),
195  2.0 / (gamma - 1.0),
196  2.0 / (gamma + 1.0),
197  (gamma - 1.0) / (gamma + 1.0),
198  0.5 * (gamma - 1.0),
199  gamma - 1.0};
200 
201  // Compute sound speeds.
202  NekDouble cL = sqrt(gamma * pL / rhoL);
203  NekDouble cR = sqrt(gamma * pR / rhoR);
204 
205  ASSERTL0(g[4] * (cL + cR) > (uR - uL),
206  "Vacuum is generated by given data.");
207 
208  // Guess initial pressure.
209  NekDouble pOld = guessp(g, rhoL, uL, pL, cL, rhoR, uR, pR, cR);
210  NekDouble uDiff = uR - uL;
211  NekDouble p, fL, fR, fLd, fRd, change;
212  int k;
213 
214  // Newton-Raphson iteration for pressure in star region.
215  for (k = 0; k < NRITER; ++k)
216  {
217  prefun(g, pOld, rhoL, pL, cL, fL, fLd);
218  prefun(g, pOld, rhoR, pR, cR, fR, fRd);
219  p = pOld - (fL + fR + uDiff) / (fLd + fRd);
220  change = 2 * fabs((p - pOld) / (p + pOld));
221 
222  if (change <= TOL)
223  {
224  break;
225  }
226 
227  if (p < 0.0)
228  {
229  p = TOL;
230  }
231 
232  pOld = p;
233  }
234 
235  ASSERTL0(k < NRITER, "Divergence in Newton-Raphson scheme");
236 
237  // Compute velocity in star region.
238  NekDouble u = 0.5 * (uL + uR + fR - fL);
239 
240  // -- SAMPLE ROUTINE --
241  // The variable S aligns with the Fortran parameter S of the SAMPLE
242  // routine, but is hard-coded as 0 (and should be optimised out by the
243  // compiler). Since we are using a Godunov scheme we pick this as 0 (see
244  // chapter 6 of Toro 2009).
245  const NekDouble S = 0.0;
246 
247  // Computed primitive variables.
248  NekDouble outRho, outU, outV, outW, outP;
249 
250  if (S <= u)
251  {
252  if (p <= pL)
253  {
254  // Left rarefaction
255  NekDouble shL = uL - cL;
256  if (S <= shL)
257  {
258  // Sampled point is left data state.
259  outRho = rhoL;
260  outU = uL;
261  outV = vL;
262  outW = wL;
263  outP = pL;
264  }
265  else
266  {
267  NekDouble cmL = cL * pow(p / pL, g[1]);
268  NekDouble stL = u - cmL;
269 
270  if (S > stL)
271  {
272  // Sampled point is star left state
273  outRho = rhoL * pow(p / pL, 1.0 / gamma);
274  outU = u;
275  outV = vL;
276  outW = wL;
277  outP = p;
278  }
279  else
280  {
281  // Sampled point is inside left fan
282  NekDouble c = g[5] * (cL + g[7] * (uL - S));
283  outRho = rhoL * pow(c / cL, g[4]);
284  outU = g[5] * (cL + g[7] * uL + S);
285  outV = vL;
286  outW = wL;
287  outP = pL * pow(c / cL, g[3]);
288  }
289  }
290  }
291  else
292  {
293  // Left shock
294  NekDouble pmL = p / pL;
295  NekDouble SL = uL - cL * sqrt(g[2] * pmL + g[1]);
296  if (S <= SL)
297  {
298  // Sampled point is left data state
299  outRho = rhoL;
300  outU = uL;
301  outV = vL;
302  outW = wL;
303  outP = pL;
304  }
305  else
306  {
307  // Sampled point is star left state
308  outRho = rhoL * (pmL + g[6]) / (pmL * g[6] + 1.0);
309  outU = u;
310  outV = vL;
311  outW = wL;
312  outP = p;
313  }
314  }
315  }
316  else
317  {
318  if (p > pR)
319  {
320  // Right shock
321  NekDouble pmR = p / pR;
322  NekDouble SR = uR + cR * sqrt(g[2] * pmR + g[1]);
323  if (S >= SR)
324  {
325  // Sampled point is right data state
326  outRho = rhoR;
327  outU = uR;
328  outV = vR;
329  outW = wR;
330  outP = pR;
331  }
332  else
333  {
334  // Sampled point is star right state
335  outRho = rhoR * (pmR + g[6]) / (pmR * g[6] + 1.0);
336  outU = u;
337  outV = vR;
338  outW = wR;
339  outP = p;
340  }
341  }
342  else
343  {
344  // Right rarefaction
345  NekDouble shR = uR + cR;
346 
347  if (S >= shR)
348  {
349  // Sampled point is right data state
350  outRho = rhoR;
351  outU = uR;
352  outV = vR;
353  outW = wR;
354  outP = pR;
355  }
356  else
357  {
358  NekDouble cmR = cR * pow(p / pR, g[1]);
359  NekDouble stR = u + cmR;
360 
361  if (S <= stR)
362  {
363  // Sampled point is star right state
364  outRho = rhoR * pow(p / pR, 1.0 / gamma);
365  outU = u;
366  outV = vR;
367  outW = wR;
368  outP = p;
369  }
370  else
371  {
372  // Sampled point is inside left fan
373  NekDouble c = g[5] * (cR - g[7] * (uR - S));
374  outRho = rhoR * pow(c / cR, g[4]);
375  outU = g[5] * (-cR + g[7] * uR + S);
376  outV = vR;
377  outW = wR;
378  outP = pR * pow(c / cR, g[3]);
379  }
380  }
381  }
382  }
383 
384  // Transform computed primitive variables to fluxes.
385  rhof = outRho * outU;
386  rhouf = outP + outRho * outU * outU;
387  rhovf = outRho * outU * outV;
388  rhowf = outRho * outU * outW;
389  Ef = outU *
390  (outP / (gamma - 1.0) +
391  0.5 * outRho * (outU * outU + outV * outV + outW * outW) + outP);
392 }
393 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define TOL
#define NRITER
static std::string solverName
ExactSolverToro(const LibUtilities::SessionReaderSharedPtr &pSession)
static RiemannSolverSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
virtual void v_PointSolve(NekDouble rhoL, NekDouble rhouL, NekDouble rhovL, NekDouble rhowL, NekDouble EL, NekDouble rhoR, NekDouble rhouR, NekDouble rhovR, NekDouble rhowR, NekDouble ER, NekDouble &rhof, NekDouble &rhouf, NekDouble &rhovf, NekDouble &rhowf, NekDouble &Ef) override
Exact Riemann solver for the Euler equations.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
std::map< std::string, RSParamFuncType > m_params
Map of parameter function types.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
RiemannSolverFactory & GetRiemannSolverFactory()
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
NekDouble guessp(NekDouble g[], NekDouble rhoL, NekDouble uL, NekDouble pL, NekDouble cL, NekDouble rhoR, NekDouble uR, NekDouble pR, NekDouble cR)
Use either PVRS, two-rarefaction or two-shock Riemann solvers to calculate an initial pressure for th...
void prefun(NekDouble *g, NekDouble p, NekDouble dk, NekDouble pk, NekDouble ck, NekDouble &f, NekDouble &fd)
Evaluate pressure functions fL and fR in Newton iteration of Riemann solver (see equation 4....
double NekDouble
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294