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/tools/floating_point_comparison.hpp>
33 #include <boost/test/unit_test.hpp>
34 
35 // #include <SolverUtils/RiemannSolvers/RiemannSolver.h>
36 #include "../RoeSolver.h"
37 #include "../RoeSolverSIMD.h"
38 
39 namespace Nektar
40 {
41 namespace RiemannTests
42 {
43 BOOST_AUTO_TEST_CASE(RoeAlongXconstSolution)
44 {
45  // LibUtilities::SessionReaderSharedPtr dummySession;
46  // SolverUtils::RiemannSolverSharedPtr riemannSolver;
47  // std::string riemannName = "Roe";
48  // riemannSolver = SolverUtils::GetRiemannSolverFactory()
49  // .CreateInstance(riemName, m_session);
50  // auto riemannSolver = RoeSolver();
51  auto riemannSolver = RoeSolverSIMD();
52  // Setting up parameters for Riemann solver
53  NekDouble gamma = 1.4;
54  riemannSolver.SetParam("gamma",
55  [&gamma]() -> NekDouble & { return gamma; });
56 
57  size_t spaceDim = 3;
58  size_t npts = 5; // so that avx spillover loop is engaged
59 
60  // Set up locations of velocity vector.
61  Array<OneD, Array<OneD, NekDouble>> vecLocs(1);
62  vecLocs[0] = Array<OneD, NekDouble>(spaceDim);
63  for (size_t i = 0; i < spaceDim; ++i)
64  {
65  vecLocs[0][i] = 1 + i;
66  }
67  riemannSolver.SetAuxVec(
68  "vecLocs",
69  [&vecLocs]() -> const Array<OneD, const Array<OneD, NekDouble>> & {
70  return vecLocs;
71  });
72 
73  // setup normals
74  Array<OneD, Array<OneD, NekDouble>> normals(spaceDim);
75  for (size_t i = 0; i < spaceDim; ++i)
76  {
77  normals[i] = Array<OneD, NekDouble>(npts);
78  }
79  riemannSolver.SetVector(
80  "N", [&normals]() -> const Array<OneD, const Array<OneD, NekDouble>> & {
81  return normals;
82  });
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 =
120  rhoe + 0.5 * (rhou * rhou + rhov * rhov + rhow * rhow) / rho;
121  fwd[nFields - 1][i] = E;
122  bwd[nFields - 1][i] = E;
123  // set face normal along x
124  normals[0][i] = 1.0;
125  // Ref solution
126  flxRef[0][i] = rhou;
127  flxRef[1][i] = rhou * rhou / rho + p;
128  flxRef[2][i] = rhou * rhov / rho;
129  flxRef[3][i] = rhou * rhow / rho;
130  flxRef[nFields - 1][i] = (E + p) * rhou / rho;
131  }
132 
133  riemannSolver.Solve(spaceDim, fwd, bwd, flx);
134 
135  // check fluxes
136  for (size_t i = 0; i < npts; ++i)
137  {
138  BOOST_CHECK_CLOSE(flxRef[0][i], flx[0][i], 1e-10);
139  BOOST_CHECK_CLOSE(flxRef[1][i], flx[1][i], 1e-10);
140  BOOST_CHECK_CLOSE(flxRef[2][i], flx[2][i], 1e-10);
141  BOOST_CHECK_CLOSE(flxRef[3][i], flx[3][i], 1e-10);
142  BOOST_CHECK_CLOSE(flxRef[4][i], flx[4][i], 1e-10);
143  }
144 }
145 
146 BOOST_AUTO_TEST_CASE(RoeAlongYconstSolution)
147 {
148  // LibUtilities::SessionReaderSharedPtr dummySession;
149  // SolverUtils::RiemannSolverSharedPtr riemannSolver;
150  // std::string riemannName = "Roe";
151  // riemannSolver = SolverUtils::GetRiemannSolverFactory()
152  // .CreateInstance(riemName, m_session);
153  // auto riemannSolver = RoeSolver();
154  auto riemannSolver = RoeSolverSIMD();
155  // Setting up parameters for Riemann solver
156  NekDouble gamma = 1.4;
157  riemannSolver.SetParam("gamma",
158  [&gamma]() -> NekDouble & { return gamma; });
159 
160  size_t spaceDim = 3;
161  size_t npts = 5;
162 
163  // Set up locations of velocity vector.
164  Array<OneD, Array<OneD, NekDouble>> vecLocs(1);
165  vecLocs[0] = Array<OneD, NekDouble>(spaceDim);
166  for (size_t i = 0; i < spaceDim; ++i)
167  {
168  vecLocs[0][i] = 1 + i;
169  }
170  riemannSolver.SetAuxVec(
171  "vecLocs",
172  [&vecLocs]() -> const Array<OneD, const Array<OneD, NekDouble>> & {
173  return vecLocs;
174  });
175 
176  // setup normals
177  Array<OneD, Array<OneD, NekDouble>> normals(spaceDim);
178  for (size_t i = 0; i < spaceDim; ++i)
179  {
180  normals[i] = Array<OneD, NekDouble>(npts);
181  }
182  riemannSolver.SetVector(
183  "N", [&normals]() -> const Array<OneD, const Array<OneD, NekDouble>> & {
184  return normals;
185  });
186 
187  size_t nFields = spaceDim + 2;
188  Array<OneD, Array<OneD, NekDouble>> fwd(nFields), bwd(nFields),
189  flx(nFields), flxRef(nFields);
190  for (size_t i = 0; i < nFields; ++i)
191  {
192  fwd[i] = Array<OneD, NekDouble>(npts);
193  bwd[i] = Array<OneD, NekDouble>(npts);
194  flx[i] = Array<OneD, NekDouble>(npts);
195  flxRef[i] = Array<OneD, NekDouble>(npts);
196  }
197 
198  // up to now it is "boiler plate" code to set up the test
199  // below are the conditions for the test
200 
201  // density
202  NekDouble rho = 0.9;
203  fwd[0][0] = rho;
204  bwd[0][0] = rho;
205  // x-momentum
206  NekDouble rhou = rho * 1.0;
207  fwd[1][0] = rhou;
208  bwd[1][0] = rhou;
209  // y-momentum
210  NekDouble rhov = rho * 2.0;
211  fwd[2][0] = rhov;
212  bwd[2][0] = rhov;
213  // z-momentum
214  NekDouble rhow = rho * 3.0;
215  fwd[3][0] = rhow;
216  bwd[3][0] = rhow;
217  // energy
218  NekDouble p = 1.0;
219  NekDouble rhoe = p / (gamma - 1.0);
220  NekDouble E = rhoe + 0.5 * (rhou * rhou + rhov * rhov + rhow * rhow) / rho;
221  fwd[nFields - 1][0] = E;
222  bwd[nFields - 1][0] = E;
223  // set face normal along y
224  normals[1][0] = 1.0;
225  // Ref solution
226  flxRef[0][0] = rhov;
227  flxRef[1][0] = rhov * rhou / rho;
228  flxRef[2][0] = rhov * rhov / rho + p;
229  flxRef[3][0] = rhov * rhow / rho;
230  flxRef[nFields - 1][0] = (E + p) * rhov / rho;
231 
232  riemannSolver.Solve(spaceDim, fwd, bwd, flx);
233 
234  // check fluxes
235  BOOST_CHECK_CLOSE(flxRef[0][0], flx[0][0], 1e-10);
236  BOOST_CHECK_CLOSE(flxRef[1][0], flx[1][0], 1e-10);
237  BOOST_CHECK_CLOSE(flxRef[2][0], flx[2][0], 1e-10);
238  BOOST_CHECK_CLOSE(flxRef[3][0], flx[3][0], 1e-10);
239  BOOST_CHECK_CLOSE(flxRef[4][0], flx[4][0], 1e-10);
240 }
241 
242 BOOST_AUTO_TEST_CASE(RoeAlongZconstSolution)
243 {
244  // LibUtilities::SessionReaderSharedPtr dummySession;
245  // SolverUtils::RiemannSolverSharedPtr riemannSolver;
246  // std::string riemannName = "Roe";
247  // riemannSolver = SolverUtils::GetRiemannSolverFactory()
248  // .CreateInstance(riemName, m_session);
249  // auto riemannSolver = RoeSolver();
250  auto riemannSolver = RoeSolverSIMD();
251  // Setting up parameters for Riemann solver
252  NekDouble gamma = 1.4;
253  riemannSolver.SetParam("gamma",
254  [&gamma]() -> NekDouble & { return gamma; });
255 
256  size_t spaceDim = 3;
257  size_t npts = 1;
258 
259  // Set up locations of velocity vector.
260  Array<OneD, Array<OneD, NekDouble>> vecLocs(1);
261  vecLocs[0] = Array<OneD, NekDouble>(spaceDim);
262  for (size_t i = 0; i < spaceDim; ++i)
263  {
264  vecLocs[0][i] = 1 + i;
265  }
266  riemannSolver.SetAuxVec(
267  "vecLocs",
268  [&vecLocs]() -> const Array<OneD, const Array<OneD, NekDouble>> & {
269  return vecLocs;
270  });
271 
272  // setup normals
273  Array<OneD, Array<OneD, NekDouble>> normals(spaceDim);
274  for (size_t i = 0; i < spaceDim; ++i)
275  {
276  normals[i] = Array<OneD, NekDouble>(npts);
277  }
278  riemannSolver.SetVector(
279  "N", [&normals]() -> const Array<OneD, const Array<OneD, NekDouble>> & {
280  return normals;
281  });
282 
283  size_t nFields = spaceDim + 2;
284  Array<OneD, Array<OneD, NekDouble>> fwd(nFields), bwd(nFields),
285  flx(nFields), flxRef(nFields);
286  for (size_t i = 0; i < nFields; ++i)
287  {
288  fwd[i] = Array<OneD, NekDouble>(npts);
289  bwd[i] = Array<OneD, NekDouble>(npts);
290  flx[i] = Array<OneD, NekDouble>(npts);
291  flxRef[i] = Array<OneD, NekDouble>(npts);
292  }
293 
294  // up to now it is "boiler plate" code to set up the test
295  // below are the conditions for the test
296 
297  // density
298  NekDouble rho = 0.9;
299  fwd[0][0] = rho;
300  bwd[0][0] = rho;
301  // x-momentum
302  NekDouble rhou = rho * 1.0;
303  fwd[1][0] = rhou;
304  bwd[1][0] = rhou;
305  // y-momentum
306  NekDouble rhov = rho * 2.0;
307  fwd[2][0] = rhov;
308  bwd[2][0] = rhov;
309  // z-momentum
310  NekDouble rhow = rho * 3.0;
311  fwd[3][0] = rhow;
312  bwd[3][0] = rhow;
313  // energy
314  NekDouble p = 1.0;
315  NekDouble rhoe = p / (gamma - 1.0);
316  NekDouble E = rhoe + 0.5 * (rhou * rhou + rhov * rhov + rhow * rhow) / rho;
317  fwd[nFields - 1][0] = E;
318  bwd[nFields - 1][0] = E;
319  // set face normal along y
320  normals[2][0] = 1.0;
321  // Ref solution
322  flxRef[0][0] = rhow;
323  flxRef[1][0] = rhow * rhou / rho;
324  flxRef[2][0] = rhow * rhov / rho;
325  flxRef[3][0] = rhow * rhow / rho + p;
326  flxRef[nFields - 1][0] = (E + p) * rhow / rho;
327 
328  riemannSolver.Solve(spaceDim, fwd, bwd, flx);
329 
330  // check fluxes
331  BOOST_CHECK_CLOSE(flxRef[0][0], flx[0][0], 1e-10);
332  BOOST_CHECK_CLOSE(flxRef[1][0], flx[1][0], 1e-10);
333  BOOST_CHECK_CLOSE(flxRef[2][0], flx[2][0], 1e-10);
334  BOOST_CHECK_CLOSE(flxRef[3][0], flx[3][0], 1e-10);
335  BOOST_CHECK_CLOSE(flxRef[4][0], flx[4][0], 1e-10);
336 }
337 
338 BOOST_AUTO_TEST_CASE(RoeAlongXdensityJump)
339 {
340  // LibUtilities::SessionReaderSharedPtr dummySession;
341  // SolverUtils::RiemannSolverSharedPtr riemannSolver;
342  // std::string riemannName = "Roe";
343  // riemannSolver = SolverUtils::GetRiemannSolverFactory()
344  // .CreateInstance(riemName, m_session);
345  // auto riemannSolver = RoeSolver();
346  auto riemannSolver = RoeSolverSIMD();
347  // Setting up parameters for Riemann solver
348  NekDouble gamma = 1.4;
349  riemannSolver.SetParam("gamma",
350  [&gamma]() -> NekDouble & { return gamma; });
351 
352  size_t spaceDim = 3;
353  size_t npts = 11; // so that avx spillover loop is engaged
354 
355  // Set up locations of velocity vector.
356  Array<OneD, Array<OneD, NekDouble>> vecLocs(1);
357  vecLocs[0] = Array<OneD, NekDouble>(spaceDim);
358  for (size_t i = 0; i < spaceDim; ++i)
359  {
360  vecLocs[0][i] = 1 + i;
361  }
362  riemannSolver.SetAuxVec(
363  "vecLocs",
364  [&vecLocs]() -> const Array<OneD, const Array<OneD, NekDouble>> & {
365  return vecLocs;
366  });
367 
368  // setup normals
369  Array<OneD, Array<OneD, NekDouble>> normals(spaceDim);
370  for (size_t i = 0; i < spaceDim; ++i)
371  {
372  normals[i] = Array<OneD, NekDouble>(npts);
373  }
374  riemannSolver.SetVector(
375  "N", [&normals]() -> const Array<OneD, const Array<OneD, NekDouble>> & {
376  return normals;
377  });
378 
379  size_t nFields = spaceDim + 2;
380  Array<OneD, Array<OneD, NekDouble>> fwd(nFields), bwd(nFields),
381  flx(nFields), flxRef(nFields);
382  for (size_t i = 0; i < nFields; ++i)
383  {
384  fwd[i] = Array<OneD, NekDouble>(npts);
385  bwd[i] = Array<OneD, NekDouble>(npts);
386  flx[i] = Array<OneD, NekDouble>(npts);
387  flxRef[i] = Array<OneD, NekDouble>(npts);
388  }
389 
390  // up to now it is "boiler plate" code to set up the test
391  // below are the conditions for the test
392 
393  for (size_t i = 0; i < npts; ++i)
394  {
395  // density
396  NekDouble rhoL = 1.0;
397  NekDouble rhoR = 2 * rhoL;
398  fwd[0][i] = rhoL;
399  bwd[0][i] = rhoR;
400  // x-momentum
401  NekDouble rhou = rhoL * 1.0;
402  fwd[1][i] = rhou;
403  bwd[1][i] = rhou;
404  // y-momentum
405  NekDouble rhov = rhoL * 2.0;
406  fwd[2][i] = rhov;
407  bwd[2][i] = rhov;
408  // z-momentum
409  NekDouble rhow = rhoL * 3.0;
410  fwd[3][i] = rhow;
411  bwd[3][i] = rhow;
412  // energy
413  NekDouble p = 1.0;
414  NekDouble rhoe = p / (gamma - 1.0);
415  NekDouble EL =
416  rhoe + 0.5 * (rhou * rhou + rhov * rhov + rhow * rhow) / rhoL;
417  NekDouble ER =
418  rhoe + 0.5 * (rhou * rhou + rhov * rhov + rhow * rhow) / rhoR;
419  fwd[nFields - 1][i] = EL;
420  bwd[nFields - 1][i] = ER;
421  // set face normal along x
422  normals[0][i] = 1.0;
423  // Ref solution (from point solve)
424  flxRef[0][i] = 0.87858599768171342;
425  flxRef[1][i] = 2.0449028304431223;
426  flxRef[2][i] = 1.8282946712594808;
427  flxRef[3][i] = 2.7424420068892208;
428  flxRef[nFields - 1][i] = 9.8154698039903128;
429  }
430 
431  riemannSolver.Solve(spaceDim, fwd, bwd, flx);
432 
433  // check fluxes
434  for (size_t i = 0; i < npts; ++i)
435  {
436  BOOST_CHECK_CLOSE(flxRef[0][i], flx[0][i], 1e-10);
437  BOOST_CHECK_CLOSE(flxRef[1][i], flx[1][i], 1e-10);
438  BOOST_CHECK_CLOSE(flxRef[2][i], flx[2][i], 1e-10);
439  BOOST_CHECK_CLOSE(flxRef[3][i], flx[3][i], 1e-10);
440  BOOST_CHECK_CLOSE(flxRef[4][i], flx[4][i], 1e-10);
441  }
442 }
443 
444 } // namespace RiemannTests
445 } // namespace Nektar
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::RiemannTests::BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(RoeAlongXconstSolution)
Definition: TestRiemann.cpp:43
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 Tue Aug 23 2022 07:27:15 for Nektar++ by doxygen 1.9.1