Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ExactSolverToro.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ExactSolver.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: Exact Riemann solver (Toro 2009).
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <CompressibleFlowSolver/RiemannSolvers/ExactSolverToro.h>
37 
38 #define TOL 1e-6
39 #define NRITER 20
40 
41 namespace Nektar
42 {
43  std::string ExactSolverToro::solverName =
44  SolverUtils::GetRiemannSolverFactory().RegisterCreatorFunction(
45  "ExactToro",
46  ExactSolverToro::create,
47  "Exact Riemann solver");
48 
49  ExactSolverToro::ExactSolverToro() : CompressibleSolver()
50  {
51 
52  }
53 
54  /**
55  * @brief Use either PVRS, two-rarefaction or two-shock Riemann solvers to
56  * calculate an initial pressure for the Newton-Raphson scheme.
57  *
58  * @param g Array of calculated gamma values.
59  * @param rhoL Density left state.
60  * @param rhoR Density right state.
61  * @param uL x-velocity component left state.
62  * @param uR x-velocity component right state.
63  * @param pL Pressure component left state.
64  * @param pR Pressure component right state.
65  * @param cL Sound speed component left state.
66  * @param cR Sound speed component right state.
67  * @return Computed initial guess for the Newton-Raphson scheme.
68  */
69  inline NekDouble guessp(
70  NekDouble g[],
71  NekDouble rhoL, NekDouble uL, NekDouble pL, NekDouble cL,
72  NekDouble rhoR, NekDouble uR, NekDouble pR, NekDouble cR)
73  {
74  const NekDouble quser = 2.0;
75  NekDouble cup, ppv, pmin, pmax, qmax;
76 
77  cup = 0.25*(rhoL + rhoR)*(cL + cR);
78  ppv = 0.5 *(pL + pR) + 0.5*(uL - uR)*cup;
79  ppv = std::max(0.0, ppv);
80  pmin = std::min(pL, pR);
81  pmax = std::max(pL, pR);
82  qmax = pmax/pmin;
83 
84  if (qmax <= quser && pmin <= ppv && ppv <= pmax)
85  {
86  // Select PVRS Riemann solver.
87  return ppv;
88  }
89  else if (ppv < pmin)
90  {
91  // Select two-rarefaction Riemann solver.
92  NekDouble pq = pow(pL/pR, g[1]);
93  NekDouble um = (pq*uL/cL + uR/cR + g[4]*(pq - 1.0))/(pq/cL + 1.0/cR);
94  NekDouble ptL = 1.0 + g[7]*(uL - um)/cL;
95  NekDouble ptR = 1.0 + g[7]*(um - uR)/cR;
96  return 0.5*(pL*pow(ptL, g[3]) + pR*pow(ptR, g[3]));
97  }
98  else
99  {
100  // Select two-shock Riemann solver with PVRS as estimate.
101  NekDouble geL = sqrt((g[5]/rhoL)/(g[6]*pL + ppv));
102  NekDouble geR = sqrt((g[5]/rhoR)/(g[6]*pR + ppv));
103  return (geL*pL + geR*pR - (uR - uL))/(geL + geR);
104  }
105  }
106 
107  /**
108  * @brief Evaluate pressure functions fL and fR in Newton iteration of
109  * Riemann solver (see equation 4.85 and 4.86 of Toro 2009).
110  *
111  * @param g Array of gamma parameters.
112  * @param p Pressure at current iteration.
113  * @param dk Density (left or right state)
114  * @param pk Pressure (left or right state)
115  * @param ck Sound speed (left or right state)
116  * @param f Computed pressure (shock).
117  * @param fd Computed pressure (rarefaction).
118  */
119  inline void prefun(NekDouble *g, NekDouble p, NekDouble dk, NekDouble pk,
120  NekDouble ck, NekDouble &f, NekDouble &fd)
121  {
122  if (p <= pk)
123  {
124  // Rarefaction wave
125  NekDouble prat = p/pk;
126  f = g[4]*ck*(pow(prat, g[1])-1.0);
127  fd = pow(prat, -g[2])/(dk*ck);
128  }
129  else
130  {
131  // Shock wave
132  NekDouble ak = g[5]/dk;
133  NekDouble bk = g[6]*pk;
134  NekDouble qrt = sqrt(ak/(bk+p));
135  f = (p-pk)*qrt;
136  fd = (1.0-0.5*(p-pk)/(bk+p))*qrt;
137  }
138  }
139 
140  /**
141  * @brief Exact Riemann solver for the Euler equations.
142  *
143  * This algorithm is transcribed from:
144  *
145  * "Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical
146  * Introduction", E. F. Toro (3rd edition, 2009).
147  *
148  * The full Fortran 77 routine can be found at the end of chapter 4 (section
149  * 4.9). This transcription is essentially the functions STARPU and SAMPLE
150  * glued together, and variable names are kept mostly the same. See the
151  * preceding chapter which explains the derivation of the solver. The
152  * routines PREFUN and GUESSP are kept separate and are reproduced above.
153  *
154  * @param rhoL Density left state.
155  * @param rhoR Density right state.
156  * @param rhouL x-momentum component left state.
157  * @param rhouR x-momentum component right state.
158  * @param rhovL y-momentum component left state.
159  * @param rhovR y-momentum component right state.
160  * @param rhowL z-momentum component left state.
161  * @param rhowR z-momentum component right state.
162  * @param EL Energy left state.
163  * @param ER Energy right state.
164  * @param rhof Computed Riemann flux for density.
165  * @param rhouf Computed Riemann flux for x-momentum component
166  * @param rhovf Computed Riemann flux for y-momentum component
167  * @param rhowf Computed Riemann flux for z-momentum component
168  * @param Ef Computed Riemann flux for energy.
169  */
170  void ExactSolverToro::v_PointSolve(
171  NekDouble rhoL, NekDouble rhouL, NekDouble rhovL, NekDouble rhowL, NekDouble EL,
172  NekDouble rhoR, NekDouble rhouR, NekDouble rhovR, NekDouble rhowR, NekDouble ER,
173  NekDouble &rhof, NekDouble &rhouf, NekDouble &rhovf, NekDouble &rhowf, NekDouble &Ef)
174  {
175  static NekDouble gamma = m_params["gamma"]();
176 
177  // Left and right variables.
178  NekDouble uL = rhouL / rhoL;
179  NekDouble vL = rhovL / rhoL;
180  NekDouble wL = rhowL / rhoL;
181  NekDouble uR = rhouR / rhoR;
182  NekDouble vR = rhovR / rhoR;
183  NekDouble wR = rhowR / rhoR;
184 
185  // Left and right pressure.
186  NekDouble pL = (gamma - 1.0) *
187  (EL - 0.5 * (rhouL*uL + rhovL*vL + rhowL*wL));
188  NekDouble pR = (gamma - 1.0) *
189  (ER - 0.5 * (rhouR*uR + rhovR*vR + rhowR*wR));
190 
191  // Compute gammas.
192  NekDouble g[] = { gamma,
193  (gamma-1.0)/(2.0*gamma),
194  (gamma+1.0)/(2.0*gamma),
195  2.0*gamma/(gamma-1.0),
196  2.0/(gamma-1.0),
197  2.0/(gamma+1.0),
198  (gamma-1.0)/(gamma+1.0),
199  0.5*(gamma-1.0),
200  gamma-1.0 };
201 
202  // Compute sound speeds.
203  NekDouble cL = sqrt(gamma * pL / rhoL);
204  NekDouble cR = sqrt(gamma * pR / rhoR);
205 
206  ASSERTL0(g[4]*(cL+cR) > (uR-uL), "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*(outP/(gamma-1.0) + 0.5*outRho*
390  (outU*outU + outV*outV + outW*outW) + outP);
391  }
392 }
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2