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
40namespace Nektar
41{
42
44{
45public:
48 {
49 return RiemannSolverSharedPtr(new RoeSolver(pSession));
50 }
51
52 static std::string solverName;
53
54 /// programmatic ctor
55 RoeSolver();
56
57protected:
59
60 using ND = NekDouble;
61
62 void v_PointSolve(ND rhoL, ND rhouL, ND rhovL, ND rhowL, ND EL, ND rhoR,
63 ND rhouR, ND rhovR, ND rhowR, ND ER, ND &rhof, ND &rhouf,
64 ND &rhovf, ND &rhowf, ND &Ef) final;
65
66 void v_ArraySolve(const Array<OneD, const Array<OneD, ND>> &Fwd,
67 const Array<OneD, const Array<OneD, ND>> &Bwd,
68 Array<OneD, Array<OneD, ND>> &flux) final;
69};
70
71template <class T, typename = typename std::enable_if<
72 std::is_floating_point_v<T> ||
73 tinysimd::is_vector_floating_point_v<T>>::type>
74inline void RoeKernel(T &rhoL, T &rhouL, T &rhovL, T &rhowL, T &EL, T &rhoR,
75 T &rhouR, T &rhovR, T &rhowR, T &ER, T &rhof, T &rhouf,
76 T &rhovf, T &rhowf, T &Ef, NekDouble gamma)
77{
78 // Left and right velocities
79 T uL = rhouL / rhoL;
80 T vL = rhovL / rhoL;
81 T wL = rhowL / rhoL;
82 T uR = rhouR / rhoR;
83 T vR = rhovR / rhoR;
84 T wR = rhowR / rhoR;
85
86 // Left and right pressures
87 T pL = (gamma - 1.0) * (EL - 0.5 * (rhouL * uL + rhovL * vL + rhowL * wL));
88 T pR = (gamma - 1.0) * (ER - 0.5 * (rhouR * uR + rhovR * vR + rhowR * wR));
89
90 // Left and right enthalpy
91 T hL = (EL + pL) / rhoL;
92 T hR = (ER + pR) / rhoR;
93
94 // Square root of rhoL and rhoR.
95 T srL = sqrt(rhoL);
96 T srR = sqrt(rhoR);
97 T srLR = srL + srR;
98
99 // Velocity, enthalpy and sound speed Roe averages (equation 11.60).
100 T uRoe = (srL * uL + srR * uR) / srLR;
101 T vRoe = (srL * vL + srR * vR) / srLR;
102 T wRoe = (srL * wL + srR * wR) / srLR;
103 T hRoe = (srL * hL + srR * hR) / srLR;
104 T URoe = (uRoe * uRoe + vRoe * vRoe + wRoe * wRoe);
105 T cRoe = sqrt((gamma - 1.0) * (hRoe - 0.5 * URoe));
106
107 // Compute eigenvectors (equation 11.59).
108 T k[5][5] = {{1., uRoe - cRoe, vRoe, wRoe, hRoe - uRoe * cRoe},
109 {1., uRoe, vRoe, wRoe, 0.5 * URoe},
110 {0., 0., 1., 0., vRoe},
111 {0., 0., 0., 1., wRoe},
112 {1., uRoe + cRoe, vRoe, wRoe, hRoe + uRoe * cRoe}};
113
114 // Calculate jumps \Delta u_i (defined preceding equation 11.67).
115 T jump[5] = {rhoR - rhoL, rhouR - rhouL, rhovR - rhovL, rhowR - rhowL,
116 ER - EL};
117
118 // Define \Delta u_5 (equation 11.70).
119 T jumpbar = jump[4] - (jump[2] - vRoe * jump[0]) * vRoe -
120 (jump[3] - wRoe * jump[0]) * wRoe;
121
122 // Compute wave amplitudes (equations 11.68, 11.69).
123 T alpha[5];
124 alpha[1] = (gamma - 1.0) *
125 (jump[0] * (hRoe - uRoe * uRoe) + uRoe * jump[1] - jumpbar) /
126 (cRoe * cRoe);
127 alpha[0] =
128 (jump[0] * (uRoe + cRoe) - jump[1] - cRoe * alpha[1]) / (2.0 * cRoe);
129 alpha[4] = jump[0] - (alpha[0] + alpha[1]);
130 alpha[2] = jump[2] - vRoe * jump[0];
131 alpha[3] = jump[3] - wRoe * jump[0];
132
133 // Compute average of left and right fluxes needed for equation 11.29.
134 rhof = 0.5 * (rhoL * uL + rhoR * uR);
135 rhouf = 0.5 * (pL + rhoL * uL * uL + pR + rhoR * uR * uR);
136 rhovf = 0.5 * (rhoL * uL * vL + rhoR * uR * vR);
137 rhowf = 0.5 * (rhoL * uL * wL + rhoR * uR * wR);
138 Ef = 0.5 * (uL * (EL + pL) + uR * (ER + pR));
139
140 // Needed to get right overload resolution for std::abs
141 using std::abs;
142
143 // Compute eigenvalues \lambda_i (equation 11.58).
144 T uRoeAbs = abs(uRoe);
145 T lambda[5] = {abs(uRoe - cRoe), uRoeAbs, uRoeAbs, uRoeAbs,
146 abs(uRoe + cRoe)};
147
148 // Finally perform summation (11.29).
149 for (size_t i = 0; i < 5; ++i)
150 {
151 uRoeAbs = 0.5 * alpha[i] * lambda[i];
152
153 rhof -= uRoeAbs * k[i][0];
154 rhouf -= uRoeAbs * k[i][1];
155 rhovf -= uRoeAbs * k[i][2];
156 rhowf -= uRoeAbs * k[i][3];
157 Ef -= uRoeAbs * k[i][4];
158 }
159}
160
161} // namespace Nektar
162
163#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:84
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:99
RoeSolver()
programmatic ctor
Definition: RoeSolver.cpp:53
static RiemannSolverSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Definition: RoeSolver.h:46
static std::string solverName
Definition: RoeSolver.h:52
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< RiemannSolver > RiemannSolverSharedPtr
A shared pointer to an EquationSystem object.
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:74
double NekDouble
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:289
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285