Peano
Loading...
Searching...
No Matches
CellIntegralTest.cpp
Go to the documentation of this file.
1#include "CellIntegralTest.h"
2
3#include <cmath>
4#include <iostream>
5#include <fstream>
6
7#include "../DGUtils.h"
8#include "../CellIntegral.h"
9#include "../Riemann.h"
11#include "../rusanov/Rusanov.h"
12#include "TestUtils.h"
13#include "exahype2/CellData.h"
14
15/*
16 * We use precomputed values for lagrange polynomials as basis function
17 * on Gauss-Legendre nodes
18 */
20 TestCase ("exahype2::dg::tests::CellIntegralTest"),
21 _Order4_QuadraturePoints1d { 0.04691007703066802,0.23076534494715845,0.5,0.7692346550528415,0.9530899229693319},
22 _Order4_MassMatrixDiagonal1d { 0.11846344252809464,0.23931433524968335,0.2844444444444444,0.2393143352496833,0.11846344252809476},
23 _Order4_StiffnessOperator1d {
24 -1.2005181447816875017764459698810242116451263427734375,-0.459606063929972219472830374797922559082508087158203125,0.1713312197562774918946360003246809355914592742919921875,-0.11698479961911900648630080468137748539447784423828125,0.130728401276001104935886587554705329239368438720703125,
25 1.824799516391880427335081549244932830333709716796875,-0.362969592101977378550969888237887062132358551025390625,-0.81657642230598626031934372804244048893451690673828125,0.4444344937740536405357261173776350915431976318359375,-0.464471255981286024194787387386895716190338134765625,
26 -0.95802422631547734521717529787565581500530242919921875,1.15002535018688423207322557573206722736358642578125,-0.0000000000000003789561257387200741506709440429594952952620993784360603484628882,-1.150025350186884676162435425794683396816253662109375,0.95802422631547756726178022290696389973163604736328125,
27 0.464471255981286079705938618644722737371921539306640625,-0.44443449377405397360263350492459721863269805908203125,0.81657642230598648236394865307374857366085052490234375,0.362969592101977489573272350753541104495525360107421875,-1.824799516391880427335081549244932830333709716796875,
28 -0.13072840127600116044703781881253235042095184326171875,0.11698479961911904811966422812474775128066539764404296875,-0.1713312197562775196502116159535944461822509765625,0.459606063929972219472830374797922559082508087158203125,1.2005181447816877238210508949123322963714599609375
29 },
30 _Order4_DerivativeOperator1d {
31 -10.1340811913090842466544927447102963924407958984375,15.4039041703444912201348415692336857318878173828125,-8.0870875087753280041624748264439404010772705078125,3.920798231666558830710300753707997500896453857421875,-1.1035337019266335811806811761925928294658660888671875,
32 -1.9205120472639156670169313656515441834926605224609375,-1.516706434335750142139431773102842271327972412109375,4.80550130432859834428427348029799759387969970703125,-1.8571160532876695992143822877551428973674774169921875,0.488833230558736009374598552312818355858325958251953125,
33 0.602336319455663016242397134192287921905517578125,-2.870776484669482986333832741365768015384674072265625,-0.000000000000001332267629550187848508358001708984375,2.870776484669483874512252441491000354290008544921875,-0.60233631945566312726469959670794196426868438720703125,
34 -0.48883323055873584284114485853933729231357574462890625,1.8571160532876682669467527375672943890094757080078125,-4.80550130432860012064111288054846227169036865234375,1.5167064343357505862286416231654584407806396484375,1.9205120472639156670169313656515441834926605224609375,
35 1.1035337019266331370914713261299766600131988525390625,-3.920798231666557942531881053582765161991119384765625,8.087087508775329780519314226694405078887939453125,-15.4039041703444912201348415692336857318878173828125,10.13408119130908602301133214496076107025146484375
36 },
37 _Order4_BasisFunctionValuesLeft1d { 1.5514080490943134,-0.8931583920000723,0.5333333333333335,-0.2679416522233877,0.07635866179581295},
38 _Order2_QuadraturePoints1d{0.1127016653792583,0.5,0.8872983346207417},
39 _Order2_MassMatrixDiagonal1d{0.2777777777777777,0.44444444444444436,0.27777777777777773},
40 _Order2_StiffnessOperator1d{-1.0758287072798378,-0.5737753105492469,0.35860956909327923,1.4344382763731172,6.459277139373494e-17,-1.4344382763731172,-0.35860956909327923,0.5737753105492469,1.0758287072798378},
41 _Order2_DerivativeOperator1d{-3.872983346207417,5.163977794943222,-1.2909944487358056,-1.2909944487358056,0.0,1.2909944487358056,1.2909944487358056,-5.163977794943222,3.872983346207417} {}
42
44 testMethod(runEulerOrder2OnStationarySetup);
45 testMethod(runEulerOrder4OnStationarySetup);
46}
47
49 //using euler equations as example
50 constexpr int unknowns = 5;
51 constexpr int auxiliaryVariables = 0;
52 constexpr int order = 2;
53
54 #if Dimensions==2
55 const tarch::la::Vector<2, double> x({0.0, 0.0});
56 const tarch::la::Vector<2, double> cellSize({1.0, 1.0});
57 constexpr int numberOfDoFs = (order+1)*(order+1);
58 #else
59 const tarch::la::Vector<3, double> x({0.0, 0.0, 0.0});
60 const tarch::la::Vector<3, double> cellSize({1.0, 1.0, 1.0});
61 constexpr int numberOfDoFs = (order+1)*(order+1)*(order+1);
62 #endif
63 //const double dx = 0.01;
64 const double t = 0.0;
65 const double dt = 0.01;
66
67 double QIn[numberOfDoFs * (unknowns+auxiliaryVariables)];
68 double QOut[numberOfDoFs * (unknowns+auxiliaryVariables)];
69
70 for (int dof=0; dof<numberOfDoFs; dof++) {
71 QIn[dof*(unknowns+auxiliaryVariables)+0] = 0.1;
72 QIn[dof*(unknowns+auxiliaryVariables)+1] = 0.0;
73 QIn[dof*(unknowns+auxiliaryVariables)+2] = 0.0;
74 QIn[dof*(unknowns+auxiliaryVariables)+3] = 0.0;
75 QIn[dof*(unknowns+auxiliaryVariables)+4] = 0.0;
76 }
77
78 ::exahype2::CellData cellData(
79 QIn,
80 x,
81 cellSize,
82 t,
83 dt,
84 QOut
85 );
86
88 cellData,
89 order,
90 unknowns,
91 auxiliaryVariables,
92 [&](
93 const double * __restrict__ Q,
95 double t,
96 double dt,
97 int normal,
98 double * __restrict__ F
99 ) -> void {
100 eulerFlux(Q,x,t,dt,normal,F);
101 },
102 nullptr, //ncp
103 nullptr, //source
104 nullptr, //pointSources
105 _Order2_QuadraturePoints1d,
106 _Order2_MassMatrixDiagonal1d, //massMatrixDiagonal
107 _Order2_StiffnessOperator1d,
108 _Order2_DerivativeOperator1d,
109 true, //flux
110 false, //ncp
111 false, //source
112 false); //pointSources
113
114
115 for (int dof=0; dof<numberOfDoFs; dof++) {
116 validateNumericalEqualsWithParams2( QOut[dof*(unknowns+auxiliaryVariables)+0], 0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
117 validateNumericalEqualsWithParams2( QOut[dof*(unknowns+auxiliaryVariables)+1], 0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
118 validateNumericalEqualsWithParams2( QOut[dof*(unknowns+auxiliaryVariables)+2], 0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
119 validateNumericalEqualsWithParams2( QOut[dof*(unknowns+auxiliaryVariables)+3], 0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
120 validateNumericalEqualsWithParams2( QOut[dof*(unknowns+auxiliaryVariables)+4], 0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
121 }
122
123 for (int dof=0; dof<numberOfDoFs; dof++) {
124 QIn[dof*(unknowns+auxiliaryVariables)+0] = 0.1;
125 QIn[dof*(unknowns+auxiliaryVariables)+1] = 0.0;
126 QIn[dof*(unknowns+auxiliaryVariables)+2] = 0.0;
127 QIn[dof*(unknowns+auxiliaryVariables)+3] = 0.0;
128 QIn[dof*(unknowns+auxiliaryVariables)+4] = 4.0;
129 }
130
132 cellData,
133 order,
134 unknowns,
135 auxiliaryVariables,
136 [&](
137 const double * __restrict__ Q,
139 double t,
140 double dt,
141 int normal,
142 double * __restrict__ F
143 ) -> void {
144 eulerFlux(Q,x,t,dt,normal,F);
145 },
146 nullptr, //ncp
147 nullptr, //source
148 nullptr, //pointSources
149 _Order2_QuadraturePoints1d,
150 _Order2_MassMatrixDiagonal1d, //massMatrixDiagonal
151 _Order2_StiffnessOperator1d,
152 _Order2_DerivativeOperator1d,
153 true, //flux
154 false, //ncp
155 false, //source
156 false); //pointSources
157
158
159 dfor(dof,order+1) {
160 int dofLinearised = peano4::utils::dLinearised(dof,order+1);
161 validateNumericalEqualsWithParams2( QOut[dofLinearised*(unknowns+auxiliaryVariables)+0], 0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
162 if (dof(0)<(order+1)/2) {
163 validateWithParams2( QOut[dofLinearised*(unknowns+auxiliaryVariables)+1]<0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
164 }
165 else if (dof(0)>(order+1)/2) {
166 validateWithParams2( QOut[dofLinearised*(unknowns+auxiliaryVariables)+1]>0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
167 }
168 if (dof(1)<(order+1)/2) {
169 validateWithParams2( QOut[dofLinearised*(unknowns+auxiliaryVariables)+2]<0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
170 }
171 else if (dof(1)>(order+1)/2) {
172 validateWithParams2( QOut[dofLinearised*(unknowns+auxiliaryVariables)+2]>0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
173 }
174 validateNumericalEqualsWithParams2( QOut[dofLinearised*(unknowns+auxiliaryVariables)+4], 0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
175 }
176}
177
178
180 //using euler equations as example
181 constexpr int unknowns = 5;
182 constexpr int auxiliaryVariables = 0;
183 constexpr int order = 4;
184
185 #if Dimensions==2
186 const tarch::la::Vector<2, double> x({0.0, 0.0});
187 const tarch::la::Vector<2, double> cellSize({1.0, 1.0});
188 constexpr int numberOfDoFs = (order+1)*(order+1);
189 #else
190 const tarch::la::Vector<3, double> x({0.0, 0.0, 0.0});
191 const tarch::la::Vector<3, double> cellSize({1.0, 1.0, 1.0});
192 constexpr int numberOfDoFs = (order+1)*(order+1)*(order+1);
193 #endif
194 //const double dx = 0.01;
195 const double t = 0.0;
196 const double dt = 0.01;
197
198 double QIn[numberOfDoFs * (unknowns+auxiliaryVariables)];
199 double QOut[numberOfDoFs * (unknowns+auxiliaryVariables)];
200
201 for (int dof=0; dof<numberOfDoFs; dof++) {
202 QIn[dof*(unknowns+auxiliaryVariables)+0] = 0.1;
203 QIn[dof*(unknowns+auxiliaryVariables)+1] = 0.0;
204 QIn[dof*(unknowns+auxiliaryVariables)+2] = 0.0;
205 QIn[dof*(unknowns+auxiliaryVariables)+3] = 0.0;
206 QIn[dof*(unknowns+auxiliaryVariables)+4] = 0.0;
207 }
208
209 ::exahype2::CellData cellData(
210 QIn,
211 x,
212 cellSize,
213 t,
214 dt,
215 QOut
216 );
217
219 cellData,
220 order,
221 unknowns,
222 auxiliaryVariables,
223 [&](
224 const double * __restrict__ Q,
226 double t,
227 double dt,
228 int normal,
229 double * __restrict__ F
230 ) -> void {
231 eulerFlux(Q,x,t,dt,normal,F);
232 },
233 nullptr, //ncp
234 nullptr, //source
235 nullptr, //pointSources
236 _Order4_QuadraturePoints1d,
237 _Order4_MassMatrixDiagonal1d, //massMatrixDiagonal
238 _Order4_StiffnessOperator1d,
239 _Order4_DerivativeOperator1d,
240 true, //flux
241 false, //ncp
242 false, //source
243 false); //pointSources
244
245
246 for (int dof=0; dof<numberOfDoFs; dof++) {
247 validateNumericalEqualsWithParams2( QOut[dof*(unknowns+auxiliaryVariables)+0], 0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
248 validateNumericalEqualsWithParams2( QOut[dof*(unknowns+auxiliaryVariables)+1], 0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
249 validateNumericalEqualsWithParams2( QOut[dof*(unknowns+auxiliaryVariables)+2], 0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
250 validateNumericalEqualsWithParams2( QOut[dof*(unknowns+auxiliaryVariables)+3], 0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
251 validateNumericalEqualsWithParams2( QOut[dof*(unknowns+auxiliaryVariables)+4], 0.0, dof, plotCell( QOut, order, unknowns, auxiliaryVariables ) );
252 }
253}
254
255
256/*
257void exahype2::dg::tests::CellIntegralTest::runDGEuler(){
258 validate(false);
259
260 file_stream << "Starting Euler test in " << Dimensions << " dimensions\n\n";
261
262 // geometry
263#if Dimensions==2
264 const tarch::la::Vector<2, double> x({0.0, 0.0});
265 const tarch::la::Vector<2, double> cellSize({1.0, 1.0});
266#else
267 const tarch::la::Vector<3, double> x({0.0, 0.0, 0.0});
268 const tarch::la::Vector<3, double> cellSize({1.0, 1.0, 1.0});
269#endif
270 const double dx = 0.01;
271 const double t = 0.0;
272 const double dt = 0.01;
273
274 //using euler equations as example
275 const int unknowns = 5;
276 const int auxiliaryVariables = 0;
277 const int order = 4;
278
279 //Parameter initialisation and memory allocation
280 const int nodesPerAxis = order+1;
281 const int nodesPerCell = getNodesPerCell(nodesPerAxis);
282 const int nodesPerFace = nodesPerCell/nodesPerAxis;
283 const int strideQ = unknowns+auxiliaryVariables;
284
285 //values in the cell at the beginning and end
286 double* QIn = new double[nodesPerCell*strideQ];
287 double* QOut = new double[nodesPerCell*strideQ];
288 //values on the faces
289 double* QLeft = new double[2*nodesPerFace*strideQ];
290 double* QRight = new double[2*nodesPerFace*strideQ];
291 double* QDown = new double[2*nodesPerFace*strideQ];
292 double* QUp = new double[2*nodesPerFace*strideQ];
293#if Dimensions==3
294 double* QFront = new double[2*nodesPerFace*strideQ];
295 double* QBack = new double[2*nodesPerFace*strideQ];
296#endif
297
298 //values of the rusanov fluxes on the faces
299 double* QHullForRiemann = new double[2*Dimensions*2*nodesPerFace*unknowns];
300
301
302 ::exahype2::CellData cellData(
303 QIn,
304 x,
305 cellSize,
306 t,
307 dt,
308 QOut
309 );
310
311
312 //setting initial conditions
313
314 file_stream << "Initial conditions: \n";
315 for(int i=0; i<nodesPerCell; i++){
316 eulerInitial(&QIn[i*strideQ]);
317 file_stream << "node " << i << " [";
318 for(int var=0; var<strideQ; var++){
319 file_stream << QIn[i*strideQ+var] << ", ";
320 }
321 file_stream << "]\n";
322 }
323
324 //for later comparisons
325 double QInit[strideQ];
326 eulerInitial(QInit);
327
328 cellIntegral_patchwise_in_situ_GaussLegendre_functors(
329 cellData,
330 order,
331 unknowns,
332 auxiliaryVariables,
333 [&](
334 const double * __restrict__ Q,
335 const tarch::la::Vector<Dimensions,double>& x,
336 double t,
337 double dt,
338 int normal,
339 double * __restrict__ F
340 ) -> void {
341 eulerFlux(Q,x,t,dt,normal,F);
342 },
343 nullptr, //ncp
344 nullptr, //source
345 nullptr, //pointSources
346 QuadraturePoints1d,
347 MassMatrixDiagonal1d, //massMatrixDiagonal
348 StiffnessOperator1d,
349 DerivativeOperator1d,
350 true, //flux
351 false, //ncp
352 false, //source
353 false); //pointSources
354
355
356 file_stream << "\nComputed volume integral values are: \n";
357 for(int i=0; i<nodesPerCell; i++){
358 file_stream << "node " << i << ": [";
359 for(int var=0; var<unknowns;var++){
360 file_stream << QOut[i*unknowns+var] << ", ";
361 }
362 file_stream << "]\n";
363 }
364 file_stream << "\n";
365
366 // Step 2), the values computed in step 1) are extrapolated to the surface boundary
367 projectVolumetricDataOntoFaces(
368 QIn,
369 order,
370 unknowns,
371 auxiliaryVariables,
372 BasisFunctionValuesLeft1d,
373 QLeft,
374 QRight,
375 QDown,
376 QUp
377#if Dimensions==3
378 ,QFront,
379 QBack
380#endif
381 );
382
383 double* F[2*Dimensions];
384
385 F[0] = QLeft;
386 F[1] = QRight;
387 F[2] = QDown;
388 F[3] = QUp;
389#if Dimensions==3
390 F[4] = QFront;
391 F[5] = QBack;
392#endif
393
394
395 // * project data onto opposite face as well to trade with "identical" cell,
396 // * this essentially simulates the exchange of surface data with the
397 // * neighbouring cells which Peano would take care of for us normally
398 projectVolumetricDataOntoFaces(
399 QIn,
400 order,
401 unknowns,
402 auxiliaryVariables,
403 BasisFunctionValuesLeft1d,
404 QRight,
405 QLeft,
406 QUp,
407 QDown
408#if Dimensions==3
409 ,QBack,
410 QFront
411#endif
412 );
413
414
415 for(int f=0; f<2*Dimensions; f++){
416
417 double* currentFace = F[f];
418 file_stream << "\nFace " << f << ": \n";
419
420 for(int node=0; node<2*nodesPerFace; node++){
421 file_stream << "node " << node << ": [";
422 for(int var=0; var<strideQ; var++){
423 file_stream << (std::abs(currentFace[node*strideQ+var])<0.00001 ? 0.0 : currentFace[node*strideQ+var]) << ", ";
424 }
425 file_stream << "]\n";
426 }
427
428 }
429
430
431 // The faces are assembled into one pointer construct for ease of access in the
432 // next steps. In addition information about whether the faces are at the boundary
433 // is provided since at the boundary the flux must be computed with boundary values
434 // instead of using the extrapolated solution from the neighbour.
435 double* QHull[Dimensions*2];
436#if Dimensions==2
437 QHull[0] = QLeft;
438 QHull[2] = QRight;
439 QHull[1] = QDown;
440 QHull[3] = QUp;
441#elif Dimensions==3
442 QHull[0] = QLeft;
443 QHull[3] = QRight;
444 QHull[1] = QDown;
445 QHull[4] = QUp;
446 QHull[2] = QBack;
447 QHull[5] = QFront;
448#endif
449
450
451 for(int i=0; i<2*Dimensions*2*nodesPerFace*strideQ; i++){
452 QHullForRiemann[i] = 0.0;
453 }
454
455
456 for(int face=0; face<2*Dimensions; face++){
457
458 rusanov::solveRiemannProblem_pointwise_in_situ(
459 [&](
460 const double * __restrict__ Q,
461 const tarch::la::Vector<Dimensions,double>& x,
462 double t,
463 double dt,
464 int normal,
465 double * __restrict__ F
466 ) -> void {
467 eulerFlux(Q,x,t,dt,normal,F);
468 },
469 nullptr,
470 [&](
471 const double * __restrict__ Q,
472 const tarch::la::Vector<Dimensions,double>& x,
473 double t,
474 double dt,
475 int normal
476 ) -> double {
477 return eulerEigenvalue(Q,x,t,dt,normal);
478 },
479 x,
480 cellSize,
481 t,
482 dt,
483 order,
484 unknowns,
485 auxiliaryVariables,
486 face, //faceNumber
487 QuadraturePoints1d,
488 true, //useFlux
489 false, //useNCP
490 QHull[face], //TODO
491 &QHullForRiemann[face*2*nodesPerFace*unknowns]
492 );
493
494 }
495
496
497 file_stream << "\nRusanov Fluxes per Face:\n";
498 for(int f=0; f<2*Dimensions; f++){
499 file_stream << "At face " << f << "\n";
500 for(int i=0; i<2*nodesPerFace; i++){
501 file_stream << "at node " << i << ": [";
502 for(int var=0; var<unknowns; var++){
503 file_stream << (std::abs(QHullForRiemann[(f*2*nodesPerFace+i)*unknowns+var])<0.000001 ? 0.0 : QHullForRiemann[(f*2*nodesPerFace+i)*unknowns+var]) << ", ";
504 }
505 file_stream << "]\n";
506 }
507 }
508
509
510 // The contributions from the computed rusanov flux are then added to the nodes
511 // within the cell.
512 integrateOverRiemannSolutionsAndAddToVolume_GaussLegendre(
513#if Dimensions==2
514 &QHullForRiemann[(0*2+0)*nodesPerFace*unknowns], //left
515 &QHullForRiemann[(2*2+0)*nodesPerFace*unknowns], //right
516 &QHullForRiemann[(1*2+0)*nodesPerFace*unknowns], //down
517 &QHullForRiemann[(3*2+0)*nodesPerFace*unknowns], //up
518#else
519 &QHullForRiemann[(0*2+0)*nodesPerFace*unknowns], //left
520 &QHullForRiemann[(3*2+0)*nodesPerFace*unknowns], //right
521 &QHullForRiemann[(1*2+0)*nodesPerFace*unknowns], //down
522 &QHullForRiemann[(4*2+0)*nodesPerFace*unknowns], //up
523 &QHullForRiemann[(2*2+0)*nodesPerFace*unknowns], //front
524 &QHullForRiemann[(5*2+0)*nodesPerFace*unknowns], //back
525#endif
526 order,
527 unknowns,
528 auxiliaryVariables,
529 cellSize,
530 BasisFunctionValuesLeft1d,
531 MassMatrixDiagonal1d,
532 QOut
533 );
534
535
536
537 file_stream << "\nValues after Riemann Integral are: \n";
538 for(int i=0; i<nodesPerCell; i++){
539 file_stream << "node " << i << ": [";
540 for(int var=0; var<unknowns;var++){
541 file_stream << (std::abs(QOut[i*unknowns+var]) <0.0000001 ? 0.0 : QOut[i*unknowns+var]) << ", ";
542 }
543 file_stream << "]\n";
544 }
545
546
547 multiplyWithInvertedMassMatrix_GaussLegendre(
548 cellData,
549 order,
550 unknowns,
551 auxiliaryVariables,
552 MassMatrixDiagonal1d
553 );
554
555 //Finally we check that after all these steps have been taken QOut, which should be the sum of the volume and face contributions, is zero
556
557 file_stream << "\nFinal values are: \n";
558 for(int i=0; i<nodesPerCell; i++){
559 file_stream << "node " << i << ": [";
560 for(int var=0; var<unknowns;var++){
561 file_stream << (std::abs(QOut[i*unknowns+var]) <0.0000001 ? 0.0 : QOut[i*unknowns+var]) << ", ";
562// validateNumericalEquals(QOut[i*strideQ+var], 0.0);
563 }
564 file_stream << "]\n";
565 }
566
567
568 delete[] QIn;
569 delete[] QOut;
570 delete[] QLeft;
571 delete[] QRight;
572 delete[] QDown;
573 delete[] QUp;
574#if Dimensions==3
575 delete[] QFront;
576 delete[] QBack;
577#endif
578 delete[] QHullForRiemann;
579
580}//runDGEuler
581
582
583
584void exahype2::dg::tests::CellIntegralTest::runDGElastic(){
585
586
587 file_stream << "\n\nStarting Elastic test in " << Dimensions << " dimensions\n\n";
588
589 // geometry
590#if Dimensions==2
591 const tarch::la::Vector<2, double> x({0.0, 0.0});
592 const tarch::la::Vector<2, double> cellSize({1.0, 1.0});
593#else
594 const tarch::la::Vector<3, double> x({0.0, 0.0, 0.0});
595 const tarch::la::Vector<3, double> cellSize({1.0, 1.0, 1.0});
596#endif
597 const double dx = 0.01;
598 const double t = 0.0;
599 const double dt = 0.01;
600
601 //using euler equations as example
602 const int unknowns = 9;
603 const int auxiliaryVariables = 3;
604 const int order = 4;
605
606 //Parameter initialisation and memory allocation
607 const int nodesPerAxis = order+1;
608 const int nodesPerCell = getNodesPerCell(nodesPerAxis);
609 const int nodesPerFace = nodesPerCell/nodesPerAxis;
610 const int strideQ = unknowns+auxiliaryVariables;
611
612 //values in the cell at the beginning and end
613 double* QIn = new double[nodesPerCell*strideQ];
614 double* QOut = new double[nodesPerCell*unknowns];
615 //values on the faces
616 double* QLeft = new double[2*nodesPerFace*strideQ];
617 double* QRight = new double[2*nodesPerFace*strideQ];
618 double* QDown = new double[2*nodesPerFace*strideQ];
619 double* QUp = new double[2*nodesPerFace*strideQ];
620#if Dimensions==3
621 double* QFront = new double[2*nodesPerFace*strideQ];
622 double* QBack = new double[2*nodesPerFace*strideQ];
623#endif
624
625 //values of the rusanov fluxes on the faces
626 double* QHullForRiemann = new double[2*Dimensions*2*nodesPerFace*strideQ];
627
628
629 ::exahype2::CellData cellData(
630 QIn,
631 x,
632 cellSize,
633 t,
634 dt,
635 QOut
636 );
637
638
639 //setting initial conditions
640
641 file_stream << "Initial conditions: \n";
642 for(int i=0; i<nodesPerCell; i++){
643 elasticInitial(&QIn[i*strideQ]);
644 file_stream << "node " << i << " [";
645 for(int var=0; var<strideQ; var++){
646 file_stream << QIn[i*strideQ+var] << ", ";
647 }
648 file_stream << "]\n";
649 }
650
651 //for later comparisons
652 double QInit[strideQ];
653 elasticInitial(QInit);
654
655 cellIntegral_patchwise_in_situ_GaussLegendre_functors(
656 cellData,
657 order,
658 unknowns,
659 auxiliaryVariables,
660 [&](
661 const double * __restrict__ Q,
662 const tarch::la::Vector<Dimensions,double>& x,
663 double t,
664 double dt,
665 int normal,
666 double * __restrict__ F
667 ) -> void {
668 elasticFlux(Q,x,t,dt,normal,F);
669 },
670 [&](
671 const double * __restrict__ Q,
672 const double * __restrict__ dQdx,
673 const tarch::la::Vector<Dimensions,double>& x,
674 double t,
675 double dt,
676 int normal,
677 double * __restrict__ BgradQ
678 ) -> void{
679 elasticNonConservativeProduct(Q, dQdx, x, t, dt, normal, BgradQ);
680 },
681 nullptr, //source
682 nullptr, //pointSources
683 QuadraturePoints1d,
684 MassMatrixDiagonal1d, //massMatrixDiagonal
685 StiffnessOperator1d,
686 DerivativeOperator1d,
687 true, //flux
688 true, //ncp
689 false, //source
690 false); //pointSources
691
692 file_stream << "\nComputed volume integral values are: \n";
693 for(int i=0; i<nodesPerCell; i++){
694 file_stream << "node " << i << ": [";
695 for(int var=0; var<unknowns;var++){
696 file_stream << (std::abs(QOut[i*unknowns+var]) <0.0000001 ? 0.0 : QOut[i*unknowns+var]) << ", ";
697 }
698 file_stream << "]\n";
699 }
700 file_stream << "\n";
701
702 // Step 2), the values computed in step 1) are extrapolated to the surface boundary
703 projectVolumetricDataOntoFaces(
704 QIn,
705 order,
706 unknowns,
707 auxiliaryVariables,
708 BasisFunctionValuesLeft1d,
709 QLeft,
710 QRight,
711 QDown,
712 QUp
713#if Dimensions==3
714 ,QFront,
715 QBack
716#endif
717 );
718
719 double* F[2*Dimensions];
720
721 F[0] = QLeft;
722 F[1] = QRight;
723 F[2] = QDown;
724 F[3] = QUp;
725#if Dimensions==3
726 F[4] = QFront;
727 F[5] = QBack;
728#endif
729
730
731 // project data onto opposite face as well to trade with "identical" cell,
732 // this essentially simulates the exchange of surface data with the
733 // neighbouring cells which Peano would take care of for us normally
734 projectVolumetricDataOntoFaces(
735 QIn,
736 order,
737 unknowns,
738 auxiliaryVariables,
739 BasisFunctionValuesLeft1d,
740 QRight,
741 QLeft,
742 QUp,
743 QDown
744#if Dimensions==3
745 ,QBack,
746 QFront
747#endif
748 );
749
750
751 for(int f=0; f<2*Dimensions; f++){
752
753 double* currentFace = F[f];
754 file_stream << "\nFace " << f << ": \n";
755
756 for(int node=0; node<2*nodesPerFace; node++){
757 file_stream << "node " << node << ": [";
758 for(int var=0; var<strideQ; var++){
759 file_stream << (std::abs(currentFace[node*strideQ+var])<0.00001 ? 0.0 : currentFace[node*strideQ+var]) << ", ";
760 }
761 file_stream << "]\n";
762 }
763
764 }
765
766
767 // The faces are assembled into one pointer construct for ease of access in the
768 // next steps. In addition information about whether the faces are at the boundary
769 // is provided since at the boundary the flux must be computed with boundary values
770 // instead of using the extrapolated solution from the neighbour.
771 double* QHull[Dimensions*2];
772#if Dimensions==2
773 QHull[0] = QLeft;
774 QHull[2] = QRight;
775 QHull[1] = QDown;
776 QHull[3] = QUp;
777#elif Dimensions==3
778 QHull[0] = QLeft;
779 QHull[3] = QRight;
780 QHull[1] = QDown;
781 QHull[4] = QUp;
782 QHull[2] = QBack;
783 QHull[5] = QFront;
784#endif
785
786
787 for(int i=0; i<2*Dimensions*2*nodesPerFace*unknowns; i++){
788 QHullForRiemann[i] = 0.0;
789 }
790
791
792 for(int face=0; face<2*Dimensions; face++){
793
794 rusanov::solveRiemannProblem_pointwise_in_situ(
795 [&](
796 const double * __restrict__ Q,
797 const tarch::la::Vector<Dimensions,double>& x,
798 double t,
799 double dt,
800 int normal,
801 double * __restrict__ F
802 ) -> void {
803 elasticFlux(Q,x,t,dt,normal,F);
804 },
805 [&](
806 const double * __restrict__ Q,
807 const double * __restrict__ dQdx,
808 const tarch::la::Vector<Dimensions,double>& x,
809 double t,
810 double dt,
811 int normal,
812 double * __restrict__ BgradQ
813 ) -> void{
814 elasticNonConservativeProduct(Q, dQdx, x, t, dt, normal, BgradQ);
815 },
816 [&](
817 const double * __restrict__ Q,
818 const tarch::la::Vector<Dimensions,double>& x,
819 double t,
820 double dt,
821 int normal
822 ) -> double {
823 return elasticEigenvalue(Q,x,t,dt,normal);
824 },
825 x,
826 cellSize,
827 t,
828 dt,
829 order,
830 unknowns,
831 auxiliaryVariables,
832 face, //faceNumber
833 QuadraturePoints1d,
834 true, //useFlux
835 true, //useNCP
836 QHull[face], //TODO
837 &QHullForRiemann[face*2*nodesPerFace*unknowns]
838 );
839
840 }
841
842
843 file_stream << "\nRusanov Fluxes per Face:\n";
844 for(int f=0; f<2*Dimensions; f++){
845 file_stream << "At face " << f << "\n";
846 for(int i=0; i<2*nodesPerFace; i++){
847 file_stream << "at node " << i << ": [";
848 for(int var=0; var<unknowns; var++){
849 file_stream << (std::abs(QHullForRiemann[(f*2*nodesPerFace+i)*unknowns+var])<0.000001 ? 0.0 : QHullForRiemann[(f*2*nodesPerFace+i)*unknowns+var]) << ", ";
850 }
851 file_stream << "]\n";
852 }
853 }
854
855
856 // The contributions from the computed rusanov flux are then added to the nodes
857 // within the cell.
858 integrateOverRiemannSolutionsAndAddToVolume_GaussLegendre(
859#if Dimensions==2
860 &QHullForRiemann[(0*2+0)*nodesPerFace*unknowns], //left
861 &QHullForRiemann[(2*2+0)*nodesPerFace*unknowns], //right
862 &QHullForRiemann[(1*2+0)*nodesPerFace*unknowns], //down
863 &QHullForRiemann[(3*2+0)*nodesPerFace*unknowns], //up
864#else
865 &QHullForRiemann[(0*2+0)*nodesPerFace*unknowns], //left
866 &QHullForRiemann[(3*2+0)*nodesPerFace*unknowns], //right
867 &QHullForRiemann[(1*2+0)*nodesPerFace*unknowns], //down
868 &QHullForRiemann[(4*2+0)*nodesPerFace*unknowns], //up
869 &QHullForRiemann[(2*2+0)*nodesPerFace*unknowns], //front
870 &QHullForRiemann[(5*2+0)*nodesPerFace*unknowns], //back
871#endif
872 order,
873 unknowns,
874 auxiliaryVariables,
875 cellSize,
876 BasisFunctionValuesLeft1d,
877 MassMatrixDiagonal1d,
878 QOut
879 );
880
881
882
883 file_stream << "\nValues after Riemann Integral are: \n";
884 for(int i=0; i<nodesPerCell; i++){
885 file_stream << "node " << i << ": [";
886 for(int var=0; var<unknowns;var++){
887 file_stream << (std::abs(QOut[i*unknowns+var]) <0.0000001 ? 0.0 : QOut[i*unknowns+var]) << ", ";
888 }
889 file_stream << "]\n";
890 }
891
892
893 multiplyWithInvertedMassMatrix_GaussLegendre(
894 cellData,
895 order,
896 unknowns,
897 auxiliaryVariables,
898 MassMatrixDiagonal1d
899 );
900
901 //Finally we check that after all these steps have been taken QOut, which should be the sum of the volume and face contributions, is zero
902
903 file_stream << "\nFinal values are: \n";
904 for(int i=0; i<nodesPerCell; i++){
905 file_stream << "node " << i << ": [";
906 for(int var=0; var<unknowns;var++){
907 file_stream << (std::abs(QOut[i*unknowns+var]) <0.0000001 ? 0.0 : QOut[i*unknowns+var]) << ", ";
908// validateNumericalEquals(QOut[i*strideQ+var], 0.0);
909 }
910 file_stream << "]\n";
911 }
912
913
914 delete[] QIn;
915 delete[] QOut;
916 delete[] QLeft;
917 delete[] QRight;
918 delete[] QDown;
919 delete[] QUp;
920#if Dimensions==3
921 delete[] QFront;
922 delete[] QBack;
923#endif
924
925 delete[] QHullForRiemann;
926
927
928} //runDGElastic
929*/
930
#define dfor(counter, max)
d-dimensional Loop
Definition Loop.h:313
#define testMethod(name)
Run a test method and check for errors.
Definition TestMacros.h:24
#define validateWithParams2(booleanExpr, param0, param1)
Definition TestMacros.h:69
#define validateNumericalEqualsWithParams2(actualValue, validValue, param0, param1)
Definition TestMacros.h:497
void runEulerOrder2OnStationarySetup()
I initialise all velocity data with zero, and set density to 1.
virtual void run() override
This routine is triggered by the TestCaseCollection.
void eulerFlux(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, double dt, int normal, double *__restrict__ F)
Definition TestUtils.cpp:3
void cellIntegral_patchwise_in_situ_GaussLegendre_functors(::exahype2::CellData< double, double > &cellData, const int order, const int unknowns, const int auxiliaryVariables, Flux flux, NonConservativeProduct nonconservativeProduct, Source source, PointSources pointSources, const double *__restrict__ QuadratureNodes1d, const double *__restrict__ MassMatrixDiagonal1d, const double *__restrict__ StiffnessMatrix1d, const double *__restrict__ DerivativeOperator1d, bool evaluateFlux, bool evaluateNonconservativeProduct, bool evaluateSource, bool evaluatePointSources)
std::string plotCell(const double *__restrict__ Q, const int order, const int unknowns, const int auxiliaryVariables)
Definition DGUtils.cpp:160
CPUGPUMethod int dLinearised(const tarch::la::Vector< Dimensions, int > &counter, int max)
Map d-dimensional vector onto integer index.
Definition Loop.cpp:106
Representation of a number of cells which contains all information that's required to process the sto...
Definition CellData.h:77
Simple vector class.
Definition Vector.h:134