Peano 4
Loading...
Searching...
No Matches
ADERDGTest.cpp
Go to the documentation of this file.
1// This file is part of the Peano project. For conditions of distribution and
2// use, please see the copyright notice at www.peano-framework.org
3#include "ADERDGTest.h"
4
5#include <sstream>
6#include <iostream>
7#include <iomanip>
8
9#include <limits>
10#include <cmath>
11
12#include "tarch/la/Vector.h"
13
19
21 TestCase ("exahype2::aderdg::ADERDGTest"),
22 QuadraturePoints{0.069431844202973713731097404888714663684368133544921875,0.33000947820757187134432797392946667969226837158203125,0.66999052179242812865567202607053332030773162841796875,0.9305681557970262307577513638534583151340484619140625},
23 QuadratureWeights{0.1739274225687269248563637802362791262567043304443359375,0.32607257743127304738806060413480736315250396728515625,0.32607257743127304738806060413480736315250396728515625,0.1739274225687269248563637802362791262567043304443359375},
24 BarycentricWeights{-7.4205400680389477230392003548331558704376220703125,18.79544940755506132745722425170242786407470703125,-18.79544940755506132745722425170242786407470703125,7.42054006803894861121762005495838820934295654296875},
25 BasisFunctionValuesLeft{1.5267881254572663873858573424513451755046844482421875,-0.8136324494869271450880887641687877476215362548828125,0.40076152031165046540905905203544534742832183837890625,-0.11391719628198997138479597879268112592399120330810546875},
26 BasisFunctionValuesRight{-0.113917196281990040773735017864964902400970458984375,0.400761520311650742964815208324580453336238861083984375,-0.813632449486927811221903539262712001800537109375,1.5267881254572674976088819676078855991363525390625},
27 DerivativeOperator{-6.6640004727045631938153746887110173702239990234375,9.7203088313703940315235740854404866695404052734375,-4.21756469699035818621268845163285732269287109375,1.161256338324529568950538305216468870639801025390625,-1.5151152295984677831341969067580066621303558349609375,-0.7688287844464174458636307463166303932666778564453125,2.941340462561433444221847821609117090702056884765625,-0.6573964485165488813578349436284042894840240478515625,0.65739644851654832624632263105013407766819000244140625,-2.941340462561433444221847821609117090702056884765625,0.768828784446416335640606121160089969635009765625,1.515115229598468449268011681851930916309356689453125,-1.161256338324528680772118605091236531734466552734375,4.21756469699035729803426875150762498378753662109375,-9.720308831370392255166734685190021991729736328125,6.66400047270456497017221408896148204803466796875},
28 MassMatrix{0.1739274225687269248563637802362791262567043304443359375,0.0,0.0,0.0,0.0,0.32607257743127304738806060413480736315250396728515625,0.0,0.0,0.0,0.0,0.32607257743127304738806060413480736315250396728515625,0.0,0.0,0.0,0.0,0.1739274225687269248563637802362791262567043304443359375},
29 StiffnessOperator{-1.1590524262142825051569161587394773960113525390625,-0.494037528020547400675610560938366688787937164306640625,0.214358954361956122180998818294028751552104949951171875,-0.201974321866382811041518152705975808203220367431640625,1.69062826161228674237690938753075897693634033203125,-0.250693983347795967819848783619818277657032012939453125,-0.95909046573029943516530693159438669681549072265625,0.73355015726438665968345276269246824085712432861328125,-0.7335501572643867707057552252081222832202911376953125,0.95909046573029943516530693159438669681549072265625,0.250693983347795634752941396072856150567531585693359375,-1.690628261612286298287699537468142807483673095703125,0.2019743218663829775749718464794568717479705810546875,-0.21435895436195628871445251206750981509685516357421875,0.494037528020547622720215485969674773514270782470703125,1.159052426214282949246126008802093565464019775390625},
30 InvertedPredictorLhsOperator{0.546435362419644660267499525208023669620160944759845733642578125,-0.14432618329325673533526652736469486626447178423404693603515625,0.1014623598288631781773111092959105405952868750318884849548339844,-0.06687578310368458280271795196592066190532932523638010025024414063,1.0188533167713026317442726043083212061901576817035675048828125,0.5847599728573232642954426996340089317527599632740020751953125,-0.1702637242678432449784643037959952494020399171859025955200195313,0.1014623598288624077161422867843221240491402568295598030090332031,1.0240105066930918847022125017787175238481722772121429443359375,1.0007437785531960714770216558378024274134077131748199462890625,0.58475997285732400610656911421614267965196631848812103271484375,-0.1443261832932565697233846802038925716260564513504505157470703125,0.97400505826439649705202061813480440832790918648242950439453125,1.02401050669309141394162920857269227781216613948345184326171875,1.0188533167713030021077347253566358631360344588756561279296875,0.546435362419644295325048266587231182711548171937465667724609375},
31 EquidistantGridProjector{1.5267881254572663873858573424513451755046844482421875,-0.0049591884824765342099084364235750399529933929443359375,-0.0021913277260105974188209021491502426215447485446929931640625,-0.113917196281990040773735017864964902400970458984375,-0.8136324494869271450880887641687877476215362548828125,0.997304018973186767738070557243190705776214599609375,0.0098464972353006600946923043693459476344287395477294921875,0.400761520311650742964815208324580453336238861083984375,0.40076152031165046540905905203544534742832183837890625,0.00984649723530049529596208657267197850160300731658935546875,0.99730401897318643467116316969622857868671417236328125,-0.813632449486927811221903539262712001800537109375,-0.11391719628198997138479597879268112592399120330810546875,-0.00219132772601056229067051361880658078007400035858154296875,-0.004959188482476617476635283310315571725368499755859375,1.5267881254572674976088819676078855991363525390625},
32 FineGridProjector{{1.3365809435384148340375531915924511849880218505859375,0.75017378205554174908087361473008058965206146240234375,0.25006831351327185597455127208377234637737274169921875,0.03282922721362747930928804862560355104506015777587890625,-0.5106599509271045889136075857095420360565185546875,0.350399125663716670686653742450289428234100341796875,0.9137546431442424843538674394949339330196380615234375,1.01007139472141194147525311564095318317413330078125,0.2422582771987365768406874622087343595921993255615234375,-0.1376638598070825392216676164025557227432727813720703125,-0.2182390043731604889476471953457803465425968170166015625,-0.055641038587549929150810612554778344929218292236328125,-0.068179269810046794209057452462729997932910919189453125,0.037090952087824181904185394387241103686392307281494140625,0.054416047715646100046971156416475423611700534820556640625,0.01274041665251061418440148287345436983741819858551025390625},{-0.035350042596419377349814539002181845717132091522216796875,-0.0928683884767426970352488524440559558570384979248046875,-0.071267786528997401074292383782449178397655487060546875,-0.01767502129820965051099079801133484579622745513916015625,0.97104619867609065497759956997469998896121978759765625,0.7760907833371601949323803637525998055934906005859375,0.388045391668580041955038950618472881615161895751953125,0.08197886521853813002191913028582348488271236419677734375,0.081978865218538310433160631873761303722858428955078125,0.38804539166857987542158525684499181807041168212890625,0.776090783337159972887775438721291720867156982421875,0.97104619867609087702220449500600807368755340576171875,-0.0176750212982097025526950773155476781539618968963623046875,-0.07126778652899741495208019159690593369305133819580078125,-0.0928683884767427525464000837018829770386219024658203125,-0.035350042596419335716451115558811579830944538116455078125},{0.01274041665251053785656853989394221571274101734161376953125,0.054416047715646016780244309529734891839325428009033203125,0.037090952087824084759670739686043816618621349334716796875,-0.0681792698100467664534818368338164873421192169189453125,-0.055641038587549630778372744543958106078207492828369140625,-0.2182390043731603779253447328301263041794300079345703125,-0.1376638598070822616659114601134206168353557586669921875,0.242258277198736549085111846579820849001407623291015625,1.01007139472141194147525311564095318317413330078125,0.9137546431442427063984723645262420177459716796875,0.350399125663715838019385273582884110510349273681640625,-0.5106599509271045889136075857095420360565185546875,0.032829227213627291959152643130437354557216167449951171875,0.250068313513271689441097578310291282832622528076171875,0.7501737820555423041923859273083508014678955078125,1.3365809435384148340375531915924511849880218505859375}},
33 testEuler_UIn{
34 1.00000000000000, 0.100000000000000, 0.200000000000000, 0.300000000000000, 2.57000000000000},
35 testEuler_QOut{
36 1.00000000000000, 0.100000000000000, 0.200000000000000, 0.300000000000000, 2.57000000000000},
37 testAdvection_UIn{
38 1.00000000000000, 1.000000000000000, 1.123456789000000, 0.987654321000000}
39{
40}
41
43 std::function< void(
44 const double * const __restrict__ Q,
46 double t,
47 int normal,
48 double * __restrict__ F
49 ) > flux,
50 std::function< double(
51 const double * const __restrict__ Q,
53 double t,
54 int normal
55 ) > maxEigenvalue,
56 std::function< void(
57 const double* Q,
58 const int i_Q,
59 const int unknown
60 )> validatePicardLoopHullAndCorrectorResult,
62 const double dx,
63 const double t,
64 const double dt,
65 const int unknowns,
66 const int auxiliaryVariables,
67 const int order,
68 const double* test_UIn,
69 const bool verbose
70) {
71 std::ostringstream buffer;
72
73 const int strideQ = unknowns+auxiliaryVariables;
74 const int nodesPerAxis = (order+1);
75 int nodesPerFace = nodesPerAxis;
76 if ( Dimensions == 3 ) {
77 nodesPerFace *= nodesPerAxis;
78 }
79 const int nodesPerCell = nodesPerFace*nodesPerAxis;
80 const int spaceTimeNodesPerCell = nodesPerCell*nodesPerAxis;
81
82 // In-/Outputs:
83 buffer << "\ninput (U):\n\n";
84 double U[nodesPerCell*strideQ];
85 for (int i = 0; i < nodesPerCell; i++) {
86 for (int m=0; m < unknowns; m++) {
87 const int i_U = i*strideQ + m;
88 U[i_U] = test_UIn[m];
89 buffer << std::fixed << std::showpoint << std::setprecision(6) << U[i_U] << " ";
90 }
91 for (int m=unknowns; m < strideQ; m++) {
92 const int i_U = i*strideQ + m;
93 U[i_U] = 0.0;
94 buffer << std::fixed << std::showpoint << std::setprecision(6) << U[i_U] << " ";
95 }
96 buffer << "\n";
97 }
98 if ( verbose ) {
99 std::cout << buffer.str() << std::flush;
100 }
101
102 // Outputs:
103 double Q[spaceTimeNodesPerCell*strideQ];
104 double _QHull[Dimensions*2][2*nodesPerCell*strideQ];
105 double* QHull[Dimensions*2];
106 double* QHull_faceSwap[Dimensions*2]; // left becomes right
107 double riemannResult[Dimensions*2*nodesPerAxis];
108 bool atBoundary[Dimensions*2];
109 for (int face = 0; face < Dimensions*2; face++) {
110 const int orientation = face/Dimensions;
111 const int direction = face - Dimensions*orientation;
112 QHull [ face ] = &_QHull[face][0];
113 QHull_faceSwap[ face ] = &_QHull[Dimensions*(1-orientation)+direction][0];
114 atBoundary [ face ] = false;
115 }
116 double maxEigenvaluePerFace[Dimensions*2];
117
119 flux,
120 [&](
121 const double * const __restrict__ Q,
123 double t,
124 double * __restrict__ S
125 )->void {},
126 [&](
127 const double * const __restrict__ Q,
128 const double * const __restrict__ dQ_or_dQdn,
130 double t,
131 int direction,
132 double * __restrict__ BgradQ
133 )->void {},
134 Q, // Q
135 U, // U
136 QuadraturePoints,
137 QuadratureWeights,
138 StiffnessOperator, // Kxi,
139 InvertedPredictorLhsOperator, // iK1,
140 BasisFunctionValuesLeft, // FLCoeff,
141 DerivativeOperator, // dudx,
142 x,
143 dx, // we assume cubic/square cells
144 t,
145 dt,
146 order,
147 unknowns,
148 auxiliaryVariables,
149 1e-8, // atol,
150 true,
151 false,
152 false
153 );
154
155 // print result
156 buffer << "\nresult (Q):\n\n";
157 for (int i = 0; i < spaceTimeNodesPerCell; i++) {
158 for (int m=0; m < unknowns; m++) {
159 const int i_Q = i*unknowns + m;
160 validatePicardLoopHullAndCorrectorResult(Q,i_Q,m);
161 buffer << std::fixed << std::showpoint << std::setprecision(6) << Q[i_Q] << " ";
162 }
163 buffer << "\n";
164 }
165 if ( verbose ) {
166 std::cout << buffer.str() << std::flush;
167 }
168
170 [&] (
171 const double * const __restrict__ Q,
173 double t,
174 int direction,
175 double * __restrict__ F
176 ) -> void {
177 flux(Q,x,t,direction,F);
178 },
179 [&] (
180 const double * const __restrict__ Q,
182 double t,
183 double * __restrict__ S
184 ) -> void {
185 },
186 [&] (
187 const double * const __restrict__ Q,
188 double * __restrict__ dQ_or_deltaQ,
190 double t,
191 int direction,
192 double * __restrict__ BgradQ
193 ) -> void {
194 },
195 U, // UOut,
196 Q, // QIn,
197 QuadraturePoints, // nodes,
198 QuadratureWeights, // weights,
199 StiffnessOperator, // Kxi,
200 DerivativeOperator, // dudx,
201 x, // cellCentre,
202 dx, // dx,
203 t, // t,
204 dt, // dt,
205 order, // order,
206 unknowns, // unknowns,
207 auxiliaryVariables, // auxiliaryVariables,
208 true /*callFlux*/, // callFlux,
209 false /*callSource*/, // callSource,
210 false /*callNonconservativeProduct*/
211 ); // callNonconservativeProduct)
212
213
214 // TODO Using simple extrapolation leads to wrong result
215 //::exahype2::aderdg::spaceTimePredictor_extrapolateInTime_loop_AoS(
216 // U,
217 // Q,
218 // BasisFunctionValuesRight, // FRCoeff,
219 // order,
220 // unknowns,
221 // auxiliaryVariables);
222
223 // print result
224 buffer << "\nresult (U, post volume integral):\n\n";
225 for (int i = 0; i < nodesPerCell; i++) {
226 for (int m=0; m < unknowns; m++) {
227 const int i_U = i*unknowns + m;
228 buffer << std::fixed << std::showpoint << std::setprecision(6) << U[i_U] << " ";
229 }
230 buffer << "\n";
231 }
232 if ( verbose ) {
233 std::cout << buffer.str() << std::flush;
234 }
235
236 // prepare QHull
237 for (int face = 0; face < 2*Dimensions; face++) {
238 for (int i = 0; i < 2*nodesPerCell; i++) {
239 for (int m=0; m < unknowns; m++) {
240 const int i_QHull = i*unknowns + m;
241 QHull[face][i_QHull] = std::numeric_limits<double>::quiet_NaN();
242 }
243 }
244 }
245
247 QHull,
248 Q,
249 BasisFunctionValuesLeft, // FLCoeff,
250 BasisFunctionValuesRight, // FRCoeff,
251 order,
252 unknowns,
253 auxiliaryVariables);
254
255 // print result
256 buffer << "\nresult (QHull):\n\n";
257 for (int face = 0; face < 2*Dimensions; face++) {
258 for (int i = 0; i < 2*nodesPerCell; i++) {
259 for (int m=0; m < unknowns; m++) {
260 const int i_QHull = i*unknowns + m;
261 const double value = QHull[face][i_QHull];
262 buffer << std::fixed << std::showpoint << std::setprecision(6) << value << " ";
263 }
264 buffer << "\n";
265 }
266 buffer << "\n";
267 }
268 if ( verbose ) {
269 std::cout << buffer.str() << std::flush;
270 }
271
272 // extrapolate again with swapped face pointers (left<->right) to fill the neigbour blocks
274 QHull_faceSwap,
275 Q,
276 BasisFunctionValuesLeft, // FLCoeff,
277 BasisFunctionValuesRight, // FRCoeff,
278 order,
279 unknowns,
280 auxiliaryVariables);
281
282 // print result
283 buffer << "\nresult (QHull, post filling neighbour arrays):\n\n";
284 for (int face = 0; face < 2*Dimensions; face++) {
285 for (int i = 0; i < 2*nodesPerCell; i++) {
286 for (int m=0; m < unknowns; m++) {
287 const int i_QHull = i*unknowns + m;
288 const double value = QHull[face][i_QHull];
289 validatePicardLoopHullAndCorrectorResult(QHull[face],i_QHull,m);
290 buffer << std::fixed << std::showpoint << std::setprecision(6) << value << " ";
291 }
292 buffer << "\n";
293 }
294 buffer << "\n";
295 }
296 if ( verbose ) {
297 std::cout << buffer.str() << std::flush;
298 }
299
300 // max eigenvalue
302 maxEigenvalue,
303 maxEigenvaluePerFace,
304 QHull,
305 QuadraturePoints,
306 x,
307 dx,
308 t,
309 dt,
310 order,
311 unknowns,
312 auxiliaryVariables);
313
314 buffer << "\nresult (maxEigenvaluePerFace):\n\n";
315 for (int face = 0; face < 2*Dimensions; face++) {
316 buffer << std::fixed << std::showpoint << std::setprecision(6) << maxEigenvaluePerFace[face] << " ";
317 }
318 if ( verbose ) {
319 std::cout << buffer.str() << std::flush;
320 }
321 // Rusanov test
323 flux,
324 flux,
325 [&](
326 double * __restrict__ Q,
327 double * __restrict__ dQ_or_deltaQ,
329 double t,
330 int direction,
331 double * __restrict__ BgradQ
332 ) -> void {},
333 [&](
334 double * __restrict__ Q,
335 double * __restrict__ dQ_or_deltaQ,
337 double t,
338 int direction,
339 double * __restrict__ BgradQ
340 ) -> void {},
341 riemannResult,
342 QHull,
343 maxEigenvaluePerFace,
344 QuadraturePoints,
345 QuadratureWeights,
346 x,
347 dx,
348 t,
349 dt,
350 order,
351 unknowns,
352 auxiliaryVariables,
353 atBoundary,
354 true /*callFlux*/,
355 false/*callNonconservativeProduct*/);
356
357 // print result
358 buffer << "\nresult (riemannResult):\n\n";
359 for (int face = 0; face < Dimensions*2; face++) {
360 for (int i = 0; i < nodesPerFace; i++) {
361 for (int m=0; m < unknowns; m++) {
362 const int i_QFace = (face*nodesPerFace+i)*unknowns + m;
363 buffer << std::fixed << std::showpoint << std::setprecision(6) << riemannResult[i_QFace] << " ";
364 }
365 buffer << "\n";
366 }
367 buffer << "\n";
368 }
369 if ( verbose ) {
370 std::cout << buffer.str() << std::flush;
371 }
372
374 [&] (
375 double * __restrict__ Q,
377 double t
378 ) -> void {
379 },
380 maxEigenvalue,
381 U,
382 riemannResult,
383 QuadraturePoints,
384 QuadratureWeights,
385 BasisFunctionValuesLeft, // only left
386 x,
387 dx,
388 t,
389 dt,
390 order,
391 unknowns,
392 auxiliaryVariables,
393 false); // callMaxEigenvalue
394
395 buffer << "\nresult (U, post Riemann):\n\n";
396 for (int i = 0; i < nodesPerCell; i++) {
397 for (int m=0; m < unknowns; m++) {
398 const int i_U = i*unknowns + m;
399 validatePicardLoopHullAndCorrectorResult(U,i_U,m);
400 buffer << std::fixed << std::showpoint << std::setprecision(6) << U[i_U] << " ";
401 }
402 buffer << "\n";
403 }
404 if ( verbose ) {
405 std::cout << buffer.str() << std::flush;
406 }
407}
408
410 testMethod (testAdvection)
411 testMethod (testEuler)
412 testMethod (testInterpolate)
413}
414
416 std::ostringstream buffer;
417 const bool verbose = false;
418
419 constexpr int unknowns = 5;
420 constexpr int auxiliaryVariables = 0;
421 constexpr int order = 3; // order must match nodes, weights etc.
422
423 constexpr int strideQ = unknowns+auxiliaryVariables;
424 constexpr int nodesPerAxis = (order+1);
425 int nodesPerFace = nodesPerAxis;
426 if ( Dimensions == 3 ) {
427 nodesPerFace *= nodesPerAxis;
428 }
429 const int nodesPerCell = nodesPerFace*nodesPerAxis;
430
431 if ( verbose ) {
432 std::cout << "# INTERPOLATION:\n" << std::endl;
433 }
434
435 // In-/Outputs:
436 buffer << "\ninput (U):\n\n";
437 double U[nodesPerCell*strideQ];
438 for (int i = 0; i < nodesPerCell; i++) {
439 for (int m=0; m < unknowns; m++) {
440 const int i_U = i*strideQ + m;
441 U[i_U] = testEuler_UIn[m];
442 buffer << std::fixed << std::showpoint << std::setprecision(6) << U[i_U] << " ";
443 }
444 for (int m=unknowns; m < strideQ; m++) {
445 const int i_U = i*strideQ + m;
446 U[i_U] = 0.0;
447 buffer << std::fixed << std::showpoint << std::setprecision(6) << U[i_U] << " ";
448 }
449 buffer << "\n";
450 }
451 if ( verbose ) {
452 std::cout << buffer.str() << std::flush;
453 }
454
455 const tarch::la::Vector<Dimensions,double> refX(0.5); // must be in [0,1]^d
456
457 double pointwiseQOut[strideQ]; // unknowns + auxiliaryVariables
459 U,
460 QuadraturePoints,
461 BarycentricWeights,
462 refX,
463 nodesPerAxis,
464 strideQ,
465 pointwiseQOut);
466
467 buffer << "\noutput (pointwiseQOut):\n\n";
468 for (int m=0; m < unknowns; m++) {
469 const double value = pointwiseQOut[m];
470 const double eps = 1.0e-6;
472 value, testEuler_UIn[m], eps, m); // constant solution
473 buffer << std::fixed << std::showpoint << std::setprecision(6) << value << " ";
474 }
475 buffer << "\n";
476 if ( verbose ) {
477 std::cout << buffer.str() << std::flush;
478 }
479}
480
482 // geometry
483 const tarch::la::Vector<Dimensions, double> x({0.0, 0.0});
484 const double dx = 1;
485 const double t = 0.0;
486 const double dt = 0.001;
487
488 const int unknowns = 4;
489 const int auxiliaryVariables = 0;
490 const int order = 3; // order must match nodes, weights etc.
491
492 const bool verbose = false;
493
494 if ( verbose ) {
495 std::cout << "# LINEAR ADVECTION:\n" << std::endl;
496 }
497 runADERDGStep(
498 [&] (
499 const double * const __restrict__ Q,
501 double t,
502 int direction,
503 double * __restrict__ F
504 ) -> void {
505 const double irho = 1./Q[0];
506 const double velocity = irho*Q[direction+1];
507 F[0] = velocity*Q[0];
508 F[1] = 0.0;
509 F[2] = 0.0;
510 F[3] = 0.0;
511 },
512 [&] (
513 const double * const __restrict__ Q,
515 double t,
516 const int direction
517 ) -> double {
518 const double irho = 1./Q[0];
519 const double velocity = irho*Q[direction+1];
520 return std::abs(velocity);
521 },
522 [&] (
523 const double* Q,
524 const int i_Q,
525 const int m
526 ) -> void {
527 const double eps = 1.0e-6;
529 Q[i_Q], testAdvection_UIn[m], eps, i_Q); // constant solution
530 },
531 x,dx,t,dt,unknowns,auxiliaryVariables,order,
532 testAdvection_UIn,
533 verbose
534 );
535}
536
538 // geometry
539 const tarch::la::Vector<Dimensions, double> x({0.0, 0.0});
540 const double dx = 5e-02;
541 const double t = 0.0;
542 const double dt = 1.686854344081342E-003;
543
544 const int unknowns = 5;
545 const int auxiliaryVariables = 0;
546 const int order = 3; // order must match nodes, weights etc.
547
548 const bool verbose = false;
549
550 if ( verbose ) {
551 std::cout << "\n# COMPRESSIBLE EULER:\n" << std::endl;
552 }
553 runADERDGStep(
554 [&] (
555 const double * const __restrict__ Q,
557 double t,
558 int direction,
559 double * __restrict__ F
560 ) -> void {
561 constexpr double gamma = 1.4;
562 const double irho = 1./Q[0];
563 const double p = (gamma-1) * (Q[4] - 0.5*irho*(Q[1]*Q[1]+Q[2]*Q[2]));
564
565 const double velocity = irho*Q[direction+1];
566 F[0] = velocity*Q[0];
567 F[1] = velocity*Q[1];
568 F[2] = velocity*Q[2];
569 F[3] = velocity*Q[3];
570 F[direction+1] += p;
571 F[4] = velocity*(Q[4]+p);
572 },
573 [&] (
574 const double * const __restrict__ Q,
576 double t,
577 const int direction
578 ) -> double {
579 constexpr double gamma = 1.4;
580 const double irho = 1./Q[0];
581 #if Dimensions==3
582 const double p = (gamma-1) * (Q[4] - 0.5*irho*(Q[1]*Q[1]+Q[2]*Q[2]+Q[3]*Q[3]));
583 #else
584 const double p = (gamma-1) * (Q[4] - 0.5*irho*(Q[1]*Q[1]+Q[2]*Q[2]));
585 #endif
586
587 const double u_n = Q[direction + 1] * irho;
588 const double c = std::sqrt(gamma * p * irho);
589
590 return std::max( std::abs(u_n - c), std::abs(u_n + c) );
591 },
592 [&] (
593 const double* Q,
594 const int i_Q,
595 const int m
596 ) -> void {
597 const double eps = 1.0e-5;
599 Q[i_Q], testEuler_QOut[m], eps, i_Q);
600 },
601 x,dx,t,dt,unknowns,auxiliaryVariables,order,
602 testEuler_UIn,
603 verbose
604 );
605}
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ so we are integrating with f$ phi_k phi_k dx
#define testMethod(name)
Run a test method and check for errors.
Definition TestMacros.h:24
#define validateNumericalEqualsWithEpsWithParams1(actualValue, validValue, eps, param0)
Definition TestMacros.h:449
void testEuler()
Tests ADER-DG with Euler fluxes and eigenvalues.
void testInterpolate()
Test interpolation operator.
void testAdvection()
Tests ADER-DG with linear advection flux and eigenvalues.
void run() override
This routine is triggered by the TestCaseCollection.
void runADERDGStep(std::function< void(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ F) > flux, std::function< double(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, int normal) > maxEigenvalue, std::function< void(const double *Q, const int i_Q, const int unknown)> validatePicardLoopHullAndCorrectorResult, const tarch::la::Vector< Dimensions, double > &x, const double dx, const double t, const double dt, const int unknowns, const int auxiliaryVariables, const int order, const double *test_UIn, const bool verbose)
void rusanovNonlinear_loop_AoS(std::function< void(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ F) > flux, std::function< void(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ F) > boundaryFlux, std::function< void(double *__restrict__ Q, double *__restrict__ dQ_or_deltaQ, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ BgradQ) > nonconservativeProduct, std::function< void(double *__restrict__ Q, double *__restrict__ dQ_or_deltaQ, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ BgradQ) > boundaryNonconservativeProduct, double *__restrict__ riemannResultOut, const double *const __restrict__ QHullIn[Dimensions *2], const double maxEigenvaluePerFace[Dimensions *2], const double *const __restrict__ nodes, const double *const __restrict__ weights, const tarch::la::Vector< Dimensions, double > &cellCentre, const double dx, const double t, const double dt, const int order, const int unknowns, const int auxiliaryVariables, const bool atBoundary[Dimensions *2], const bool callFlux, const bool callNonconservativeProduct)
Apply a Rusanov Riemann solver to the pair of extrapolated predictor values from two neighbouring cel...
void spaceTimePredictor_PicardLoop_loop_AoS(std::function< void(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ F) > flux, std::function< void(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, double *__restrict__ S) > algebraicSource, std::function< void(const double *const __restrict__ Q, const double *const __restrict__ dQ_or_deltaQ, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ BgradQ) > nonconservativeProduct, double *__restrict__ QOut, const double *const __restrict__ UIn, const double *const __restrict__ nodes, const double *const __restrict__ weights, const double *const __restrict__ Kxi, const double *const __restrict__ iK1, const double *const __restrict__ FLCoeff, const double *const __restrict__ dudx, const tarch::la::Vector< Dimensions, double > &cellCentre, const double dx, const double t, const double dt, const int order, const int unknowns, const int auxiliaryVariables, const double atol, const bool callFlux, const bool callSource, const bool callNonconservativeProduct)
Compute the space-time predictor (Qout) from the current solution (UIn).
void corrector_addCellContributions_loop_AoS(std::function< void(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ F) > flux, std::function< void(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, double *__restrict__ S) > algebraicSource, std::function< void(const double *const __restrict__ Q, double *__restrict__ dQ_or_deltaQ, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ BgradQ) > nonconservativeProduct, double *__restrict__ UOut, const double *const __restrict__ QIn, const double *const __restrict__ nodes, const double *const __restrict__ weights, const double *const __restrict__ Kxi, const double *const __restrict__ dudx, const tarch::la::Vector< Dimensions, double > &cellCentre, const double dx, const double t, const double dt, const int order, const int unknowns, const int auxiliaryVariables, const bool callFlux, const bool callSource, const bool callNonconservativeProduct)
Add cell-local contributions to new solution or update vector.
void spaceTimePredictor_extrapolate_loop_AoS(double *__restrict__ QHullOut[Dimensions *2], const double *const __restrict__ QIn, const double *const __restrict__ FLCoeff, const double *const __restrict__ FRCoeff, const int order, const int unknowns, const int auxiliaryVariables)
Extrapolate all of cell's predictor DOF to the hull of the element.
double corrector_addRiemannContributions_loop_AoS(std::function< void(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t) > adjustSolution, std::function< double(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, const int normal) > maxAbsoluteEigenvalue, double *__restrict__ UOut, const double *const __restrict__ riemannResultIn, const double *const __restrict__ nodes, const double *const __restrict__ weights, const double *const __restrict__ FLCoeff, const tarch::la::Vector< Dimensions, double > &cellCentre, const double dx, const double t, const double dt, const int order, const int unknowns, const int auxiliaryVariables, const bool callMaxEigenvalue)
Add cell-local contributions to new solution or update vector.
GPUCallableMethod void interpolate_AoS(const double *__restrict__ const UIn, const double *__restrict__ const nodes, const double *__restrict__ const barycentricWeights, const tarch::la::Vector< Dimensions, double > &referenceCoodinates, const int nodesPerAxis, const int strideQ, double *__restrict__ pointwiseQOut)
Evalutes DG polynomial at arbitrary reference space coordinate in reference domain [0,...
void riemann_maxAbsoluteEigenvalue_loop_AoS(std::function< double(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, const int direction) > maxAbsoluteEigenvalue, double maxEigenvaluePerFaceOut[Dimensions *2], double *const QHullIn[Dimensions *2], const double *const __restrict__ nodes, const tarch::la::Vector< Dimensions, double > &cellCentre, const double dx, const double t, const double dt, const int order, const int unknowns, const int auxiliaryVariables)
Per face, determine the eigenvalue with the largest absolute value among the boundary-extrapolated pr...
Simple vector class.
Definition Vector.h:134