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