Nektar++
Functions
Nektar::RiemannTests Namespace Reference

Functions

 BOOST_AUTO_TEST_CASE (RoeAlongXconstSolution)
 
 BOOST_AUTO_TEST_CASE (RoeAlongYconstSolution)
 
 BOOST_AUTO_TEST_CASE (RoeAlongZconstSolution)
 
 BOOST_AUTO_TEST_CASE (RoeAlongXdensityJump)
 
 BOOST_AUTO_TEST_CASE (RRR)
 

Function Documentation

◆ BOOST_AUTO_TEST_CASE() [1/5]

Nektar::RiemannTests::BOOST_AUTO_TEST_CASE ( RoeAlongXconstSolution  )

Definition at line 46 of file TestRiemann.cpp.

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.
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}
double NekDouble

References CellMLToNektar.cellml_metadata::p.

◆ BOOST_AUTO_TEST_CASE() [2/5]

Nektar::RiemannTests::BOOST_AUTO_TEST_CASE ( RoeAlongXdensityJump  )

Definition at line 341 of file TestRiemann.cpp.

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.
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}

References CellMLToNektar.cellml_metadata::p.

◆ BOOST_AUTO_TEST_CASE() [3/5]

Nektar::RiemannTests::BOOST_AUTO_TEST_CASE ( RoeAlongYconstSolution  )

Definition at line 149 of file TestRiemann.cpp.

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.
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}

References CellMLToNektar.cellml_metadata::p.

◆ BOOST_AUTO_TEST_CASE() [4/5]

Nektar::RiemannTests::BOOST_AUTO_TEST_CASE ( RoeAlongZconstSolution  )

Definition at line 245 of file TestRiemann.cpp.

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.
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}

References CellMLToNektar.cellml_metadata::p.

◆ BOOST_AUTO_TEST_CASE() [5/5]

Nektar::RiemannTests::BOOST_AUTO_TEST_CASE ( RRR  )

Definition at line 43 of file UnitTestRiemann.cpp.

44{
45
46 BOOST_CHECK_CLOSE(1., 1., 1e-10);
47}