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