Nektar++
RoeSolverSIMD.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: RoeSolverSIMD.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 using simd types.
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
37 
39 
40 namespace Nektar
41 {
42 std::string RoeSolverSIMD::solverName =
44  "RoeOpt", RoeSolverSIMD::create, "Roe Riemann solver opt");
45 
48  : CompressibleSolver(pSession)
49 {
50  m_requiresRotation = false;
51 }
52 
53 /// programmatic ctor
55 {
56  m_requiresRotation = false;
57 }
58 
60  const int nDim, const Array<OneD, const Array<OneD, NekDouble>> &fwd,
61  const Array<OneD, const Array<OneD, NekDouble>> &bwd,
63 {
64  static auto gamma = m_params["gamma"]();
65  static auto nVars = fwd.size();
66  static auto spaceDim = nVars - 2;
67 
68  // 3D case only so far
69  ASSERTL0(spaceDim == 3, "SIMD Roe implemented only for 3D case...");
70 
71  using namespace tinysimd;
72  using vec_t = simd<NekDouble>;
73  // using vec_t = typename tinysimd::abi::scalar<NekDouble>::type;
74 
75  // get limit of vectorizable chunk
76  size_t sizeScalar = fwd[0].size();
77  size_t sizeVec = (sizeScalar / vec_t::width) * vec_t::width;
78 
79  // get normal, vellocs
80  ASSERTL1(CheckVectors("N"), "N not defined.");
81  // ASSERTL1(CheckAuxVec("vecLocs"), "vecLocs not defined.");
83  // const Array<OneD, const Array<OneD, NekDouble> > vecLocs =
84  // m_auxVec["vecLocs"]();
85 
86  // const unsigned int vx = (int)vecLocs[0][0];
87  // const unsigned int vy = (int)vecLocs[0][1];
88  // const unsigned int vz = (int)vecLocs[0][2];
89 
90  // Generate matrices if they don't already exist.
91  if (m_rotMat.size() == 0)
92  {
93  GenerateRotationMatrices(normals);
94  }
95 
96  // SIMD loop
97  size_t i = 0;
98  for (; i < sizeVec; i += vec_t::width)
99  {
100  // load scalars
101  vec_t rhoL, rhoR, ER, EL;
102  rhoL.load(&(fwd[0][i]), is_not_aligned);
103  rhoR.load(&(bwd[0][i]), is_not_aligned);
104  ER.load(&(bwd[spaceDim + 1][i]), is_not_aligned);
105  EL.load(&(fwd[spaceDim + 1][i]), is_not_aligned);
106 
107  // load vectors left
108  vec_t tmpIn[3], tmpOut[3];
109  tmpIn[0].load(&(fwd[1][i]), is_not_aligned);
110  tmpIn[1].load(&(fwd[2][i]), is_not_aligned);
111  tmpIn[2].load(&(fwd[3][i]), is_not_aligned);
112 
113  // load rotation matrix
114  vec_t rotMat[9];
115  for (size_t j = 0; j < 9; ++j)
116  {
117  rotMat[j].load(&(m_rotMat[j][i]), is_not_aligned);
118  }
119 
120  // rotateTo kernel Fwd
121  rotateToNormalKernel(tmpIn, rotMat, tmpOut);
122 
123  vec_t rhouL = tmpOut[0];
124  vec_t rhovL = tmpOut[1];
125  vec_t rhowL = tmpOut[2];
126 
127  // load vectors right
128  tmpIn[0].load(&(bwd[1][i]), is_not_aligned);
129  tmpIn[1].load(&(bwd[2][i]), is_not_aligned);
130  tmpIn[2].load(&(bwd[3][i]), is_not_aligned);
131 
132  // rotateTo kernel Bwd
133  rotateToNormalKernel(tmpIn, rotMat, tmpOut);
134 
135  vec_t rhouR = tmpOut[0];
136  vec_t rhovR = tmpOut[1];
137  vec_t rhowR = tmpOut[2];
138 
139  // Roe kernel
140  vec_t rhof{}, Ef{};
141  RoeKernel(rhoL, rhouL, rhovL, rhowL, EL, rhoR, rhouR, rhovR, rhowR, ER,
142  rhof, tmpIn[0], tmpIn[1], tmpIn[2], Ef, gamma);
143 
144  // rotateFrom kernel
145  rotateFromNormalKernel(tmpIn, rotMat, tmpOut);
146 
147  // store scalar
148  rhof.store(&(flux[0][i]), is_not_aligned);
149  Ef.store(&(flux[nVars - 1][i]), is_not_aligned);
150 
151  // store vector 3D only
152  tmpOut[0].store(&(flux[1][i]), is_not_aligned);
153  tmpOut[1].store(&(flux[2][i]), is_not_aligned);
154  tmpOut[2].store(&(flux[3][i]), is_not_aligned);
155  }
156 
157  // spillover loop
158  for (; i < sizeScalar; ++i)
159  {
160  // load scalars
161  NekDouble rhoL = fwd[0][i];
162  NekDouble rhoR = bwd[0][i];
163  NekDouble EL = fwd[spaceDim + 1][i];
164  NekDouble ER = bwd[spaceDim + 1][i];
165 
166  // 3D case only
167  // load vectors left
168  NekDouble tmpIn[3], tmpOut[3];
169  tmpIn[0] = fwd[1][i];
170  tmpIn[1] = fwd[2][i];
171  tmpIn[2] = fwd[3][i];
172 
173  // load rotation matrix
174  NekDouble rotMat[9];
175  for (size_t j = 0; j < 9; ++j)
176  {
177  rotMat[j] = m_rotMat[j][i];
178  }
179 
180  // rotateTo kernel Fwd
181  rotateToNormalKernel(tmpIn, rotMat, tmpOut);
182 
183  NekDouble rhouL = tmpOut[0];
184  NekDouble rhovL = tmpOut[1];
185  NekDouble rhowL = tmpOut[2];
186 
187  // load vectors right
188  tmpIn[0] = bwd[1][i];
189  tmpIn[1] = bwd[2][i];
190  tmpIn[2] = bwd[3][i];
191 
192  // rotateTo kernel Bwd
193  rotateToNormalKernel(tmpIn, rotMat, tmpOut);
194 
195  NekDouble rhouR = tmpOut[0];
196  NekDouble rhovR = tmpOut[1];
197  NekDouble rhowR = tmpOut[2];
198 
199  // Roe kernel
200  NekDouble rhof{}, Ef{};
201  RoeKernel(rhoL, rhouL, rhovL, rhowL, EL, rhoR, rhouR, rhovR, rhowR, ER,
202  rhof, tmpIn[0], tmpIn[1], tmpIn[2], Ef, gamma);
203 
204  // rotateFrom kernel
205  rotateFromNormalKernel(tmpIn, rotMat, tmpOut);
206 
207  // store scalar
208  flux[0][i] = rhof;
209  flux[nVars - 1][i] = Ef;
210 
211  // store vector 3D only
212  flux[1][i] = tmpOut[0];
213  flux[2][i] = tmpOut[1];
214  flux[3][i] = tmpOut[2];
215  }
216 }
217 
218 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
void v_Solve(const int nDim, const Array< OneD, const Array< OneD, ND >> &Fwd, const Array< OneD, const Array< OneD, ND >> &Bwd, Array< OneD, Array< OneD, ND >> &flux) final
static RiemannSolverSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Definition: RoeSolverSIMD.h:46
static std::string solverName
Definition: RoeSolverSIMD.h:52
RoeSolverSIMD()
programmatic ctor
bool m_requiresRotation
Indicates whether the Riemann solver requires a rotation to be applied to the velocity fields.
SOLVER_UTILS_EXPORT void GenerateRotationMatrices(const Array< OneD, const Array< OneD, NekDouble >> &normals)
Generate rotation matrices for 3D expansions.
Array< OneD, Array< OneD, NekDouble > > m_rotMat
Rotation matrices for each trace quadrature point.
std::map< std::string, RSVecFuncType > m_vectors
Map of vector function types.
SOLVER_UTILS_EXPORT bool CheckVectors(std::string name)
Determine whether a vector has been defined in m_vectors.
std::map< std::string, RSParamFuncType > m_params
Map of parameter function types.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
void rotateFromNormalKernel(T *in, T *rotMat, T *out)
void rotateToNormalKernel(T *in, T *rotMat, T *out)
RiemannSolverFactory & GetRiemannSolverFactory()
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
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:73
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