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