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// 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
36
38
39namespace Nektar
40{
41
42std::string RoeSolver::solverName =
44 "Roe", RoeSolver::create, "Roe Riemann solver");
45
47 : CompressibleSolver(pSession)
48{
49 // m_pointSolve = false;
50}
51
52/// programmatic ctor
54{
55 // m_pointSolve = false;
56}
57
58/**
59 * @brief Roe Riemann solver.
60 *
61 * Stated equations numbers are from:
62 *
63 * "Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical
64 * Introduction", E. F. Toro (3rd edition, 2009).
65 *
66 * We follow the algorithm prescribed following equation 11.70.
67 *
68 * @param rhoL Density left state.
69 * @param rhoR Density right state.
70 * @param rhouL x-momentum component left state.
71 * @param rhouR x-momentum component right state.
72 * @param rhovL y-momentum component left state.
73 * @param rhovR y-momentum component right state.
74 * @param rhowL z-momentum component left state.
75 * @param rhowR z-momentum component right state.
76 * @param EL Energy left state.
77 * @param ER Energy right state.
78 * @param rhof Computed Riemann flux for density.
79 * @param rhouf Computed Riemann flux for x-momentum component
80 * @param rhovf Computed Riemann flux for y-momentum component
81 * @param rhowf Computed Riemann flux for z-momentum component
82 * @param Ef Computed Riemann flux for energy.
83 */
84void RoeSolver::v_PointSolve(double rhoL, double rhouL, double rhovL,
85 double rhowL, double EL, double rhoR, double rhouR,
86 double rhovR, double rhowR, double ER,
87 double &rhof, double &rhouf, double &rhovf,
88 double &rhowf, double &Ef)
89{
90 static NekDouble gamma = m_params["gamma"]();
91
92 RoeKernel(rhoL, rhouL, rhovL, rhowL, EL, rhoR, rhouR, rhovR, rhowR, ER,
93 rhof, rhouf, rhovf, rhowf, Ef, gamma);
94}
95
96/**
97 *
98 */
100 const Array<OneD, const Array<OneD, NekDouble>> &fwd,
101 const Array<OneD, const Array<OneD, NekDouble>> &bwd,
103{
104 static auto gamma = m_params["gamma"]();
105 static size_t nVars = fwd.size();
106 static size_t spaceDim = nVars - 2;
107
108 using namespace tinysimd;
109 using vec_t = simd<NekDouble>;
110
111 // get limit of vectorizable chunk
112 size_t sizeScalar = fwd[0].size();
113 size_t sizeVec = (sizeScalar / vec_t::width) * vec_t::width;
114
115 // SIMD loop
116 size_t i = 0;
117 for (; i < sizeVec; i += vec_t::width)
118 {
119 vec_t rhoL{}, rhouL{}, rhovL{}, rhowL{}, EL{};
120 vec_t rhoR{}, rhouR{}, rhovR{}, rhowR{}, ER{};
121
122 // load
123 rhoL.load(&(fwd[0][i]), is_not_aligned);
124 rhouL.load(&(fwd[1][i]), is_not_aligned);
125 EL.load(&(fwd[spaceDim + 1][i]), is_not_aligned);
126 rhoR.load(&(bwd[0][i]), is_not_aligned);
127 rhouR.load(&(bwd[1][i]), is_not_aligned);
128 ER.load(&(bwd[spaceDim + 1][i]), is_not_aligned);
129
130 if (spaceDim == 2)
131 {
132 rhovL.load(&(fwd[2][i]), is_not_aligned);
133 rhovR.load(&(bwd[2][i]), is_not_aligned);
134 }
135 else if (spaceDim == 3)
136 {
137 rhovL.load(&(fwd[2][i]), is_not_aligned);
138 rhowL.load(&(fwd[3][i]), is_not_aligned);
139 rhovR.load(&(bwd[2][i]), is_not_aligned);
140 rhowR.load(&(bwd[3][i]), is_not_aligned);
141 }
142
143 vec_t rhof{}, rhouf{}, rhovf{}, rhowf{}, Ef{};
144
145 RoeKernel(rhoL, rhouL, rhovL, rhowL, EL, rhoR, rhouR, rhovR, rhowR, ER,
146 rhof, rhouf, rhovf, rhowf, Ef, gamma);
147
148 // store
149 rhof.store(&(flux[0][i]), is_not_aligned);
150 rhouf.store(&(flux[1][i]), is_not_aligned);
151 Ef.store(&(flux[nVars - 1][i]), is_not_aligned);
152 if (spaceDim == 2)
153 {
154 rhovf.store(&(flux[2][i]), is_not_aligned);
155 }
156 else if (spaceDim == 3)
157 {
158 rhovf.store(&(flux[2][i]), is_not_aligned);
159 rhowf.store(&(flux[3][i]), is_not_aligned);
160 }
161
162 } // avx loop
163
164 // spillover loop
165 for (; i < sizeScalar; ++i)
166 {
167 NekDouble rhoL{}, rhouL{}, rhovL{}, rhowL{}, EL{};
168 NekDouble rhoR{}, rhouR{}, rhovR{}, rhowR{}, ER{};
169
170 // load
171 rhoL = fwd[0][i];
172 rhouL = fwd[1][i];
173 EL = fwd[spaceDim + 1][i];
174 rhoR = bwd[0][i];
175 rhouR = bwd[1][i];
176 ER = bwd[spaceDim + 1][i];
177
178 if (spaceDim == 2)
179 {
180 rhovL = fwd[2][i];
181 rhovR = bwd[2][i];
182 }
183 else if (spaceDim == 3)
184 {
185 rhovL = fwd[2][i];
186 rhowL = fwd[3][i];
187 rhovR = bwd[2][i];
188 rhowR = bwd[3][i];
189 }
190
191 NekDouble rhof{}, rhouf{}, rhovf{}, rhowf{}, Ef{};
192
193 RoeKernel(rhoL, rhouL, rhovL, rhowL, EL, rhoR, rhouR, rhovR, rhowR, ER,
194 rhof, rhouf, rhovf, rhowf, Ef, gamma);
195
196 // store
197 flux[0][i] = rhof;
198 flux[1][i] = rhouf;
199 flux[nVars - 1][i] = Ef;
200 if (spaceDim == 2)
201 {
202 flux[2][i] = rhovf;
203 }
204 else if (spaceDim == 3)
205 {
206 flux[2][i] = rhovf;
207 flux[3][i] = rhowf;
208 }
209
210 } // loop
211}
212
213} // namespace Nektar
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
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::map< std::string, RSParamFuncType > m_params
Map of parameter function types.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
RiemannSolverFactory & GetRiemannSolverFactory()
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
tinysimd::simd< NekDouble > vec_t
double NekDouble
static constexpr struct tinysimd::is_not_aligned_t is_not_aligned
typename abi< ScalarType, width >::type simd
Definition: tinysimd.hpp:80