Nektar++
TestRiemann.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 // For more information, please see: http://www.nektar.info
3 //
4 // The MIT License
5 //
6 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
7 // Department of Aeronautics, Imperial College London (UK), and Scientific
8 // Computing and Imaging Institute, University of Utah (USA).
9 //
10 // Permission is hereby granted, free of charge, to any person obtaining a
11 // copy of this software and associated documentation files (the "Software"),
12 // to deal in the Software without restriction, including without limitation
13 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
14 // and/or sell copies of the Software, and to permit persons to whom the
15 // Software is furnished to do so, subject to the following conditions:
16 //
17 /// The above copyright notice and this permission notice shall be included
18 // in all copies or substantial portions of the Software.
19 //
20 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26 // DEALINGS IN THE SOFTWARE.
27 //
28 // Description:
29 //
30 ///////////////////////////////////////////////////////////////////////////////
31 
32 #include <boost/test/auto_unit_test.hpp>
33 #include <boost/test/test_case_template.hpp>
34 #include <boost/test/floating_point_comparison.hpp>
35 #include <boost/test/unit_test.hpp>
36 
37 // #include <SolverUtils/RiemannSolvers/RiemannSolver.h>
38 #include "../RoeSolver.h"
39 #include "../RoeSolverSIMD.h"
40 
41 
42 namespace Nektar
43 {
44 namespace RiemannTests
45 {
46  BOOST_AUTO_TEST_CASE(RoeAlongXconstSolution)
47  {
48  // LibUtilities::SessionReaderSharedPtr dummySession;
49  // SolverUtils::RiemannSolverSharedPtr riemannSolver;
50  // std::string riemannName = "Roe";
51  // riemannSolver = SolverUtils::GetRiemannSolverFactory()
52  // .CreateInstance(riemName, m_session);
53  // auto riemannSolver = RoeSolver();
54  auto riemannSolver = RoeSolverSIMD();
55  // Setting up parameters for Riemann solver
56  NekDouble gamma = 1.4;
57  riemannSolver.SetParam("gamma", [&gamma]()
58  -> NekDouble& {return gamma;});
59 
60  size_t spaceDim = 3;
61  size_t npts = 5; // so that avx spillover loop is engaged
62 
63  // Set up locations of velocity vector.
64  Array<OneD, Array<OneD, NekDouble>> vecLocs(1);
65  vecLocs[0] = Array<OneD, NekDouble>(spaceDim);
66  for (size_t i = 0; i < spaceDim; ++i)
67  {
68  vecLocs[0][i] = 1 + i;
69  }
70  riemannSolver.SetAuxVec("vecLocs", [&vecLocs]()
71  -> const Array<OneD, const Array<OneD, NekDouble>>&
72  {return vecLocs;});
73 
74  // setup normals
75  Array<OneD, Array<OneD, NekDouble>> normals(spaceDim);
76  for (size_t i = 0; i < spaceDim; ++i)
77  {
78  normals[i] = Array<OneD, NekDouble>(npts);
79  }
80  riemannSolver.SetVector("N", [&normals]()
81  -> const Array<OneD, const Array<OneD, NekDouble>>&
82  {return normals;});
83 
84  size_t nFields = spaceDim + 2;
85  Array<OneD, Array<OneD, NekDouble>> fwd(nFields), bwd(nFields),
86  flx(nFields), flxRef(nFields);
87  for (size_t i = 0; i < nFields; ++i)
88  {
89  fwd[i] = Array<OneD, NekDouble>(npts);
90  bwd[i] = Array<OneD, NekDouble>(npts);
91  flx[i] = Array<OneD, NekDouble>(npts);
92  flxRef[i] = Array<OneD, NekDouble>(npts);
93  }
94 
95  // up to now it is "boiler plate" code to set up the test
96  // below are the conditions for the test
97 
98  for (size_t i = 0; i < npts; ++i)
99  {
100  // density
101  NekDouble rho = 0.9;
102  fwd[0][i] = rho;
103  bwd[0][i] = rho;
104  // x-momentum
105  NekDouble rhou = rho*1.0;
106  fwd[1][i] = rhou;
107  bwd[1][i] = rhou;
108  // y-momentum
109  NekDouble rhov = rho*2.0;
110  fwd[2][i] = rhov;
111  bwd[2][i] = rhov;
112  // z-momentum
113  NekDouble rhow = rho*3.0;
114  fwd[3][i] = rhow;
115  bwd[3][i] = rhow;
116  // energy
117  NekDouble p = 1.0;
118  NekDouble rhoe = p / (gamma - 1.0);
119  NekDouble E = rhoe + 0.5*(rhou*rhou + rhov*rhov + rhow*rhow) / rho;
120  fwd[nFields-1][i] = E;
121  bwd[nFields-1][i] = E;
122  // set face normal along x
123  normals[0][i] = 1.0;
124  // Ref solution
125  flxRef[0][i] = rhou;
126  flxRef[1][i] = rhou*rhou/rho + p;
127  flxRef[2][i] = rhou*rhov/rho;
128  flxRef[3][i] = rhou*rhow/rho;
129  flxRef[nFields-1][i] = (E+p)*rhou/rho;
130  }
131 
132  riemannSolver.Solve(spaceDim, fwd, bwd, flx);
133 
134  // check fluxes
135  for (size_t i = 0; i < npts; ++i)
136  {
137  BOOST_CHECK_CLOSE(flxRef[0][i], flx[0][i], 1e-10);
138  BOOST_CHECK_CLOSE(flxRef[1][i], flx[1][i], 1e-10);
139  BOOST_CHECK_CLOSE(flxRef[2][i], flx[2][i], 1e-10);
140  BOOST_CHECK_CLOSE(flxRef[3][i], flx[3][i], 1e-10);
141  BOOST_CHECK_CLOSE(flxRef[4][i], flx[4][i], 1e-10);
142  }
143 
144  }
145 
146 
147  BOOST_AUTO_TEST_CASE(RoeAlongYconstSolution)
148  {
149  // LibUtilities::SessionReaderSharedPtr dummySession;
150  // SolverUtils::RiemannSolverSharedPtr riemannSolver;
151  // std::string riemannName = "Roe";
152  // riemannSolver = SolverUtils::GetRiemannSolverFactory()
153  // .CreateInstance(riemName, m_session);
154  // auto riemannSolver = RoeSolver();
155  auto riemannSolver = RoeSolverSIMD();
156  // Setting up parameters for Riemann solver
157  NekDouble gamma = 1.4;
158  riemannSolver.SetParam("gamma", [&gamma]()
159  -> NekDouble& {return gamma;});
160 
161  size_t spaceDim = 3;
162  size_t npts = 5;
163 
164  // Set up locations of velocity vector.
165  Array<OneD, Array<OneD, NekDouble>> vecLocs(1);
166  vecLocs[0] = Array<OneD, NekDouble>(spaceDim);
167  for (size_t i = 0; i < spaceDim; ++i)
168  {
169  vecLocs[0][i] = 1 + i;
170  }
171  riemannSolver.SetAuxVec("vecLocs", [&vecLocs]()
172  -> const Array<OneD, const Array<OneD, NekDouble>>&
173  {return vecLocs;});
174 
175  // setup normals
176  Array<OneD, Array<OneD, NekDouble>> normals(spaceDim);
177  for (size_t i = 0; i < spaceDim; ++i)
178  {
179  normals[i] = Array<OneD, NekDouble>(npts);
180  }
181  riemannSolver.SetVector("N", [&normals]()
182  -> const Array<OneD, const Array<OneD, NekDouble>>&
183  {return normals;});
184 
185  size_t nFields = spaceDim + 2;
186  Array<OneD, Array<OneD, NekDouble>> fwd(nFields), bwd(nFields),
187  flx(nFields), flxRef(nFields);
188  for (size_t i = 0; i < nFields; ++i)
189  {
190  fwd[i] = Array<OneD, NekDouble>(npts);
191  bwd[i] = Array<OneD, NekDouble>(npts);
192  flx[i] = Array<OneD, NekDouble>(npts);
193  flxRef[i] = Array<OneD, NekDouble>(npts);
194  }
195 
196  // up to now it is "boiler plate" code to set up the test
197  // below are the conditions for the test
198 
199  // density
200  NekDouble rho = 0.9;
201  fwd[0][0] = rho;
202  bwd[0][0] = rho;
203  // x-momentum
204  NekDouble rhou = rho*1.0;
205  fwd[1][0] = rhou;
206  bwd[1][0] = rhou;
207  // y-momentum
208  NekDouble rhov = rho*2.0;
209  fwd[2][0] = rhov;
210  bwd[2][0] = rhov;
211  // z-momentum
212  NekDouble rhow = rho*3.0;
213  fwd[3][0] = rhow;
214  bwd[3][0] = rhow;
215  // energy
216  NekDouble p = 1.0;
217  NekDouble rhoe = p / (gamma - 1.0);
218  NekDouble E = rhoe + 0.5*(rhou*rhou + rhov*rhov + rhow*rhow) / rho;
219  fwd[nFields-1][0] = E;
220  bwd[nFields-1][0] = E;
221  // set face normal along y
222  normals[1][0] = 1.0;
223  // Ref solution
224  flxRef[0][0] = rhov;
225  flxRef[1][0] = rhov*rhou/rho;
226  flxRef[2][0] = rhov*rhov/rho + p;
227  flxRef[3][0] = rhov*rhow/rho;
228  flxRef[nFields-1][0] = (E+p)*rhov/rho;
229 
230  riemannSolver.Solve(spaceDim, fwd, bwd, flx);
231 
232  // check fluxes
233  BOOST_CHECK_CLOSE(flxRef[0][0], flx[0][0], 1e-10);
234  BOOST_CHECK_CLOSE(flxRef[1][0], flx[1][0], 1e-10);
235  BOOST_CHECK_CLOSE(flxRef[2][0], flx[2][0], 1e-10);
236  BOOST_CHECK_CLOSE(flxRef[3][0], flx[3][0], 1e-10);
237  BOOST_CHECK_CLOSE(flxRef[4][0], flx[4][0], 1e-10);
238 
239  }
240 
241  BOOST_AUTO_TEST_CASE(RoeAlongZconstSolution)
242  {
243  // LibUtilities::SessionReaderSharedPtr dummySession;
244  // SolverUtils::RiemannSolverSharedPtr riemannSolver;
245  // std::string riemannName = "Roe";
246  // riemannSolver = SolverUtils::GetRiemannSolverFactory()
247  // .CreateInstance(riemName, m_session);
248  // auto riemannSolver = RoeSolver();
249  auto riemannSolver = RoeSolverSIMD();
250  // Setting up parameters for Riemann solver
251  NekDouble gamma = 1.4;
252  riemannSolver.SetParam("gamma", [&gamma]()
253  -> NekDouble& {return gamma;});
254 
255  size_t spaceDim = 3;
256  size_t npts = 1;
257 
258  // Set up locations of velocity vector.
259  Array<OneD, Array<OneD, NekDouble>> vecLocs(1);
260  vecLocs[0] = Array<OneD, NekDouble>(spaceDim);
261  for (size_t i = 0; i < spaceDim; ++i)
262  {
263  vecLocs[0][i] = 1 + i;
264  }
265  riemannSolver.SetAuxVec("vecLocs", [&vecLocs]()
266  -> const Array<OneD, const Array<OneD, NekDouble>>&
267  {return vecLocs;});
268 
269  // setup normals
270  Array<OneD, Array<OneD, NekDouble>> normals(spaceDim);
271  for (size_t i = 0; i < spaceDim; ++i)
272  {
273  normals[i] = Array<OneD, NekDouble>(npts);
274  }
275  riemannSolver.SetVector("N", [&normals]()
276  -> const Array<OneD, const Array<OneD, NekDouble>>&
277  {return normals;});
278 
279  size_t nFields = spaceDim + 2;
280  Array<OneD, Array<OneD, NekDouble>> fwd(nFields), bwd(nFields),
281  flx(nFields), flxRef(nFields);
282  for (size_t i = 0; i < nFields; ++i)
283  {
284  fwd[i] = Array<OneD, NekDouble>(npts);
285  bwd[i] = Array<OneD, NekDouble>(npts);
286  flx[i] = Array<OneD, NekDouble>(npts);
287  flxRef[i] = Array<OneD, NekDouble>(npts);
288  }
289 
290  // up to now it is "boiler plate" code to set up the test
291  // below are the conditions for the test
292 
293  // density
294  NekDouble rho = 0.9;
295  fwd[0][0] = rho;
296  bwd[0][0] = rho;
297  // x-momentum
298  NekDouble rhou = rho*1.0;
299  fwd[1][0] = rhou;
300  bwd[1][0] = rhou;
301  // y-momentum
302  NekDouble rhov = rho*2.0;
303  fwd[2][0] = rhov;
304  bwd[2][0] = rhov;
305  // z-momentum
306  NekDouble rhow = rho*3.0;
307  fwd[3][0] = rhow;
308  bwd[3][0] = rhow;
309  // energy
310  NekDouble p = 1.0;
311  NekDouble rhoe = p / (gamma - 1.0);
312  NekDouble E = rhoe + 0.5*(rhou*rhou + rhov*rhov + rhow*rhow) / rho;
313  fwd[nFields-1][0] = E;
314  bwd[nFields-1][0] = E;
315  // set face normal along y
316  normals[2][0] = 1.0;
317  // Ref solution
318  flxRef[0][0] = rhow;
319  flxRef[1][0] = rhow*rhou/rho;
320  flxRef[2][0] = rhow*rhov/rho;
321  flxRef[3][0] = rhow*rhow/rho + p;
322  flxRef[nFields-1][0] = (E+p)*rhow/rho;
323 
324  riemannSolver.Solve(spaceDim, fwd, bwd, flx);
325 
326  // check fluxes
327  BOOST_CHECK_CLOSE(flxRef[0][0], flx[0][0], 1e-10);
328  BOOST_CHECK_CLOSE(flxRef[1][0], flx[1][0], 1e-10);
329  BOOST_CHECK_CLOSE(flxRef[2][0], flx[2][0], 1e-10);
330  BOOST_CHECK_CLOSE(flxRef[3][0], flx[3][0], 1e-10);
331  BOOST_CHECK_CLOSE(flxRef[4][0], flx[4][0], 1e-10);
332 
333  }
334 
335 
336  BOOST_AUTO_TEST_CASE(RoeAlongXdensityJump)
337  {
338  // LibUtilities::SessionReaderSharedPtr dummySession;
339  // SolverUtils::RiemannSolverSharedPtr riemannSolver;
340  // std::string riemannName = "Roe";
341  // riemannSolver = SolverUtils::GetRiemannSolverFactory()
342  // .CreateInstance(riemName, m_session);
343  // auto riemannSolver = RoeSolver();
344  auto riemannSolver = RoeSolverSIMD();
345  // Setting up parameters for Riemann solver
346  NekDouble gamma = 1.4;
347  riemannSolver.SetParam("gamma", [&gamma]()
348  -> NekDouble& {return gamma;});
349 
350  size_t spaceDim = 3;
351  size_t npts = 11; // so that avx spillover loop is engaged
352 
353  // Set up locations of velocity vector.
354  Array<OneD, Array<OneD, NekDouble>> vecLocs(1);
355  vecLocs[0] = Array<OneD, NekDouble>(spaceDim);
356  for (size_t i = 0; i < spaceDim; ++i)
357  {
358  vecLocs[0][i] = 1 + i;
359  }
360  riemannSolver.SetAuxVec("vecLocs", [&vecLocs]()
361  -> const Array<OneD, const Array<OneD, NekDouble>>&
362  {return vecLocs;});
363 
364  // setup normals
365  Array<OneD, Array<OneD, NekDouble>> normals(spaceDim);
366  for (size_t i = 0; i < spaceDim; ++i)
367  {
368  normals[i] = Array<OneD, NekDouble>(npts);
369  }
370  riemannSolver.SetVector("N", [&normals]()
371  -> const Array<OneD, const Array<OneD, NekDouble>>&
372  {return normals;});
373 
374  size_t nFields = spaceDim + 2;
375  Array<OneD, Array<OneD, NekDouble>> fwd(nFields), bwd(nFields),
376  flx(nFields), flxRef(nFields);
377  for (size_t i = 0; i < nFields; ++i)
378  {
379  fwd[i] = Array<OneD, NekDouble>(npts);
380  bwd[i] = Array<OneD, NekDouble>(npts);
381  flx[i] = Array<OneD, NekDouble>(npts);
382  flxRef[i] = Array<OneD, NekDouble>(npts);
383  }
384 
385  // up to now it is "boiler plate" code to set up the test
386  // below are the conditions for the test
387 
388  for (size_t i = 0; i < npts; ++i)
389  {
390  // density
391  NekDouble rhoL = 1.0;
392  NekDouble rhoR = 2*rhoL;
393  fwd[0][i] = rhoL;
394  bwd[0][i] = rhoR;
395  // x-momentum
396  NekDouble rhou = rhoL*1.0;
397  fwd[1][i] = rhou;
398  bwd[1][i] = rhou;
399  // y-momentum
400  NekDouble rhov = rhoL*2.0;
401  fwd[2][i] = rhov;
402  bwd[2][i] = rhov;
403  // z-momentum
404  NekDouble rhow = rhoL*3.0;
405  fwd[3][i] = rhow;
406  bwd[3][i] = rhow;
407  // energy
408  NekDouble p = 1.0;
409  NekDouble rhoe = p / (gamma - 1.0);
410  NekDouble EL = rhoe + 0.5*(rhou*rhou + rhov*rhov + rhow*rhow) / rhoL;
411  NekDouble ER = rhoe + 0.5*(rhou*rhou + rhov*rhov + rhow*rhow) / rhoR;
412  fwd[nFields-1][i] = EL;
413  bwd[nFields-1][i] = ER;
414  // set face normal along x
415  normals[0][i] = 1.0;
416  // Ref solution (from point solve)
417  flxRef[0][i] = 0.87858599768171342;
418  flxRef[1][i] = 2.0449028304431223;
419  flxRef[2][i] = 1.8282946712594808;
420  flxRef[3][i] = 2.7424420068892208;
421  flxRef[nFields-1][i] = 9.8154698039903128;
422  }
423 
424  riemannSolver.Solve(spaceDim, fwd, bwd, flx);
425 
426  // check fluxes
427  for (size_t i = 0; i < npts; ++i)
428  {
429  BOOST_CHECK_CLOSE(flxRef[0][i], flx[0][i], 1e-10);
430  BOOST_CHECK_CLOSE(flxRef[1][i], flx[1][i], 1e-10);
431  BOOST_CHECK_CLOSE(flxRef[2][i], flx[2][i], 1e-10);
432  BOOST_CHECK_CLOSE(flxRef[3][i], flx[3][i], 1e-10);
433  BOOST_CHECK_CLOSE(flxRef[4][i], flx[4][i], 1e-10);
434  }
435 
436  }
437 
438 
439 
440 }
441 }
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::RiemannTests::BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(RoeAlongXconstSolution)
Definition: TestRiemann.cpp:46
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Generated on Wed Nov 24 2021 10:33:29 for Nektar++ by doxygen 1.9.1