1 // Copyright (c) 2017-2021, Lawrence Livermore National Security, LLC and
2 // other Axom Project Developers. See the top-level LICENSE file for details.
3 //
4 // SPDX-License-Identifier: (BSD-3-Clause)
5
6 /**************************************************************************
7 * Program: ShockTube.C
8 * Purpose: 1D shock tube, split flux Euler equations
9 *
10 * | m | | mv |
11 * Q = | mv | F = | mv^2 + P |
12 * | E | | v(E+P) |
13 *
14 * P = (gamma - 1.0)[E - 0.5 mv^2 ]
15 *
16 * Cp
17 * gamma = -- m = mass/volume v = velocity
18 * Cv
19 *
20 * All quantiites are non-dimensionalized.
21 *
22 * @Q @F @Q @F @Q
23 * -- + -- = -- + -- -- = 0
24 * @t @x @t @Q @x
25 *
26 *************************************************************************/
27
28 //#include "Vista.h"
29 //#include "View.h"
30
31 #include "axom/config.hpp"
32
33 // Sidre component headers
34 #include "axom/sidre.hpp"
35
36 // Standard library headers
37 #include <cstdio>
38 #include <cmath>
39 #include <cstdlib>
40
41 using axom::sidre::Buffer;
42 using axom::sidre::DataStore;
43 using axom::sidre::Group;
44 using axom::sidre::IndexType;
45 using axom::sidre::TypeID;
46 using axom::sidre::View;
47
48 using namespace conduit;
49
50 #define UPWIND 0
51 #define DOWNWIND 1
52
53 const double gammaa = M_SQRT2;
54 const double gammaaInverse = M_SQRT1_2;
55
CreateScalarIntViewAndSetVal(Group * const grp,const std::string & name,int32 const value)56 void CreateScalarIntViewAndSetVal(Group* const grp,
57 const std::string& name,
58 int32 const value)
59 {
60 grp->createViewScalar(name, value);
61 }
62
CreateScalarFloatBufferViewAndSetVal(Group * const grp,const std::string & name,float64 const value)63 void CreateScalarFloatBufferViewAndSetVal(Group* const grp,
64 const std::string& name,
65 float64 const value)
66 {
67 grp->createViewScalar(name, value);
68 }
69
70 /**************************************************************************
71 * Subroutine: GetUserInput
72 * Purpose : Ask for control and output information
73 *************************************************************************/
74
GetUserInput(Group * const prob)75 void GetUserInput(Group* const prob)
76 {
77 /**********************************/
78 /* Get mesh info, and create mesh */
79 /**********************************/
80 {
81 int numElems, numFaces;
82
83 //printf("How many zones for the 1D shock tube? ");
84 //scanf("%d", &numElems);
85 numElems = 10;
86
87 numElems += 2; /* add an inflow and outflow zone */
88 numFaces = numElems - 1;
89
90 // create buffer and view, and set value
91 CreateScalarIntViewAndSetVal(prob, "numElems", numElems);
92
93 CreateScalarIntViewAndSetVal(prob, "numFaces", numFaces);
94 }
95
96 /********************/
97 /* Get physics info */
98 /********************/
99 {
100 double pratio = -1.0, dratio = -1.0;
101
102 #if 0
103 printf("What cfl number would you like to use? ");
104 scanf("%lf", &cfl);
105 #endif
106
107 while(pratio < 0.0 || pratio > 1.0)
108 {
109 //printf("What pressure ratio would you like (0 <= x <= 1)? ");
110 //scanf("%lf", &pratio);
111 pratio = 0.5;
112 }
113
114 while(dratio < 0.0 || dratio > 1.0)
115 {
116 //printf("What density ratio would you like (0 <= x <= 1)? ");
117 //scanf("%lf", &dratio);
118 dratio = 0.5;
119 }
120
121 CreateScalarFloatBufferViewAndSetVal(prob, "pressureRatio", pratio);
122
123 CreateScalarFloatBufferViewAndSetVal(prob, "densityRatio", dratio);
124 }
125
126 /********************/
127 /* Get output info */
128 /********************/
129 {
130 int numUltraDumps, numCyclesPerDump;
131
132 //printf("How many Ultra dumps would you like? ");
133 // scanf("%d", &numUltraDumps);
134 numUltraDumps = 10;
135
136 //printf("How many cycles per Ultra dump would you like? ");
137 // scanf("%d", &numCyclesPerDump);
138 numCyclesPerDump = 10;
139
140 CreateScalarIntViewAndSetVal(prob, "numUltraDumps", numUltraDumps);
141
142 CreateScalarIntViewAndSetVal(prob, "numCyclesPerDump", numCyclesPerDump);
143
144 CreateScalarIntViewAndSetVal(prob,
145 "numTotalCycles",
146 numUltraDumps * numCyclesPerDump);
147 }
148
149 return;
150 }
151
152 /**************************************************************************
153 * Subroutine: CreateShockTubeMesh
154 * Purpose : Build an empty mesh for the shock tube
155
156
157 Gaps between elements are faces
158 |
159 -------------------------------
160 | | | | | |
161
162 ### ### ### ### ### ### ### ###
163 ### ### ### ### ... ### ### ### ### <--- 1D Shock tube model
164 ### ### ### ### ### ### ### ###
165
166 | | | |
167 | ----------------------------- |
168 Inflow | Outflow
169 Element Tube Elements Element
170
171
172 *************************************************************************/
173
CreateShockTubeMesh(Group * const prob)174 void CreateShockTubeMesh(Group* const prob)
175 {
176 int i;
177 int32 const numElems = prob->getView("numElems")->getData();
178 int32 const numFaces = prob->getView("numFaces")->getData();
179
180 //int inflow[1];
181 //int outflow[1];
182
183 /* create element and face classes */
184
185 Group* const elem = prob->createGroup("elem");
186 Group* const face = prob->createGroup("face");
187
188 /* set up some important views */
189
190 //inflow[0] = 0; /* identify inflow elements */
191 // elem->viewCreate("inflow", new IndexSet(1, inflow));
192 elem->createGroup("inflow");
193
194 //outflow[0] = numElems - 1; /* identify outflow elements */
195 // elem->viewCreate("outflow", new IndexSet(1, outflow));
196 elem->createGroup("outflow");
197
198 int32 numTubeElems = numElems - 2;
199 Group* const tube = elem->createGroup(
200 "tube"); //->SetDataShape(DataStoreNS::DataShape(numTubeElems));
201
202 int32* const mapToElems =
203 tube->createViewAndAllocate("mapToElems", DataType::int32(numTubeElems))
204 ->getData();
205
206 for(int k = 0u; k < numTubeElems; ++k)
207 {
208 mapToElems[k] = k + 1;
209 }
210
211 /* Set up some important data relations */
212
213 /* Each face connects to two elements */
214
215 int32* const faceToElem =
216 face->createViewAndAllocate("faceToElem", DataType::int32(2 * numFaces))
217 ->getData();
218
219 for(i = 0; i < numFaces; ++i)
220 {
221 faceToElem[i * 2 + UPWIND] = i;
222 faceToElem[i * 2 + DOWNWIND] = i + 1;
223 }
224
225 /* Each element connects to two faces */ //
226 // Relation &elemToFace = *tube->relationCreate("elemToFace", 2);
227 int32* elemToFace =
228 tube->createViewAndAllocate("elemToFace", DataType::int32(2 * numElems))
229 ->getData();
230
231 for(i = 0; i < numElems; ++i)
232 {
233 elemToFace[i * 2 + UPWIND] = i; /* same map as above by coincidence */
234 elemToFace[i * 2 + DOWNWIND] = i + 1;
235 }
236
237 return;
238 }
239
240 /**************************************************************************
241 * Subroutine: InitializeShockTube
242 * Purpose : Populate the mesh with values
243 *************************************************************************/
244
InitializeShockTube(Group * const prob)245 void InitializeShockTube(Group* const prob)
246 {
247 int i;
248
249 /* These were created in GetUserInput() */
250 Group* const elem = (prob->getGroup("elem"));
251 Group* const face = (prob->getGroup("face"));
252
253 /* Create element centered quantities */
254
255 int32 const numElems = prob->getView("numElems")->getData();
256 int32 const numFaces = prob->getView("numFaces")->getData();
257
258 float64* const mass =
259 elem->createViewAndAllocate("mass", DataType::float64(numElems))->getData();
260
261 float64* const momentum =
262 elem->createViewAndAllocate("momentum", DataType::float64(numElems))->getData();
263
264 float64* const energy =
265 elem->createViewAndAllocate("energy", DataType::float64(numElems))->getData();
266
267 float64* const pressure =
268 elem->createViewAndAllocate("pressure", DataType::float64(numElems))->getData();
269
270 /* Create face centered quantities */
271
272 face->createViewAndAllocate("F0", DataType::float64(numFaces));
273 face->createViewAndAllocate("F1", DataType::float64(numFaces));
274 face->createViewAndAllocate("F2", DataType::float64(numFaces));
275
276 // face->fieldCreateReal("F", 3); /* mv, mv^2+P, and v(E+P) */
277
278 /* Fill left half with high pressure, right half with low pressure */
279 int startTube = 0;
280 int endTube = numElems;
281 int midTube = endTube / 2;
282
283 /* Non-dimensionalized reference values */
284 double massInitial = 1.0;
285 double momentumInitial = 0.0;
286 double pressureInitial = gammaaInverse;
287 double energyInitial = pressureInitial / (gammaa - 1.0);
288
289 /* Initialize zonal quantities*/
290 for(i = startTube; i < midTube; ++i)
291 {
292 mass[i] = massInitial;
293 momentum[i] = momentumInitial;
294 pressure[i] = pressureInitial;
295 energy[i] = energyInitial;
296 }
297
298 /* adjust parameters for low pressure portion of tube */
299 double dratio = prob->getView("densityRatio")->getData();
300 double pratio = prob->getView("pressureRatio")->getData();
301
302 massInitial *= dratio;
303 pressureInitial *= pratio;
304 energyInitial = pressureInitial / (gammaa - 1.0);
305
306 for(i = midTube; i < endTube; ++i)
307 {
308 mass[i] = massInitial;
309 momentum[i] = momentumInitial;
310 pressure[i] = pressureInitial;
311 energy[i] = energyInitial;
312 }
313
314 /* Create needed time info */
315
316 CreateScalarFloatBufferViewAndSetVal(prob, "time", 0.0);
317 CreateScalarIntViewAndSetVal(prob, "cycle", 0);
318
319 CreateScalarFloatBufferViewAndSetVal(prob, "dx", (1.0 / ((double)endTube)));
320 double dx = prob->getView("dx")->getData();
321 CreateScalarFloatBufferViewAndSetVal(prob, "dt", 0.4 * dx);
322
323 return;
324 }
325
326 /**************************************************************************
327 * Subroutine: ComputeFaceInfo
328 * Purpose : Compute F quantities at faces.
329 *
330 * @F @F0 @F1 @F2
331 * -- = --- + --- + ---
332 * @x @x @x @x
333 *
334 * Calculate F0, F1 and F2 at the face centers.
335 *
336 *************************************************************************/
337
ComputeFaceInfo(Group * const prob)338 void ComputeFaceInfo(Group* const prob)
339 {
340 int i;
341 Group* const face = prob->getGroup("face");
342 // Relation &faceToElem = *face->relation("faceToElem");
343 int32 const* const faceToElem = face->getView("faceToElem")->getData();
344
345 float64* const F0 = face->getView("F0")->getData();
346 float64* const F1 = face->getView("F1")->getData();
347 float64* const F2 = face->getView("F2")->getData();
348 int numFaces = face->getView("F0")->getNumElements();
349
350 Group* const elem = prob->getGroup("elem");
351 float64* const mass = elem->getView("mass")->getData();
352 float64* const momentum = elem->getView("momentum")->getData();
353 float64* const energy = elem->getView("energy")->getData();
354
355 for(i = 0; i < numFaces; ++i)
356 {
357 /* each face has an upwind and downwind element. */
358 int upWind = faceToElem[i * 2 + UPWIND]; /* upwind element */
359 int downWind = faceToElem[i * 2 + DOWNWIND]; /* downwind element */
360
361 /* calculate face centered quantities */
362 double massf = 0.5 * (mass[upWind] + mass[downWind]);
363 double momentumf = 0.5 * (momentum[upWind] + momentum[downWind]);
364 double energyf = 0.5 * (energy[upWind] + energy[downWind]);
365 double pressuref =
366 (gammaa - 1.0) * (energyf - 0.5 * momentumf * momentumf / massf);
367 double c = sqrt(gammaa * pressuref / massf);
368 double v = momentumf / massf;
369 double ev;
370 double cLocal;
371 int contributor;
372
373 /* Now that we have the wave speeds, we might want to */
374 /* look for the max wave speed here, and update dt */
375 /* appropriately right before leaving this function. */
376 /* ... */
377
378 /* OK, calculate face quantities */
379
380 F0[i] = F1[i] = F2[i] = 0.0;
381
382 contributor = ((v >= 0.0) ? upWind : downWind);
383 massf = mass[contributor];
384 momentumf = momentum[contributor];
385 energyf = energy[contributor];
386 pressuref = energyf - 0.5 * momentumf * momentumf / massf;
387 ev = v * (gammaa - 1.0);
388
389 F0[i] += ev * massf;
390 F1[i] += ev * momentumf;
391 F2[i] += ev * (energyf - pressuref);
392
393 contributor = ((v + c >= 0.0) ? upWind : downWind);
394 massf = mass[contributor];
395 momentumf = momentum[contributor];
396 energyf = energy[contributor];
397 pressuref = (gammaa - 1.0) * (energyf - 0.5 * momentumf * momentumf / massf);
398 ev = 0.5 * (v + c);
399 cLocal = sqrt(gammaa * pressuref / massf);
400
401 F0[i] += ev * massf;
402 F1[i] += ev * (momentumf + massf * cLocal);
403 F2[i] += ev * (energyf + pressuref + momentumf * cLocal);
404
405 contributor = ((v - c >= 0.0) ? upWind : downWind);
406 massf = mass[contributor];
407 momentumf = momentum[contributor];
408 energyf = energy[contributor];
409 pressuref = (gammaa - 1.0) * (energyf - 0.5 * momentumf * momentumf / massf);
410 ev = 0.5 * (v - c);
411 cLocal = sqrt(gammaa * pressuref / massf);
412
413 F0[i] += ev * massf;
414 F1[i] += ev * (momentumf - massf * cLocal);
415 F2[i] += ev * (energyf + pressuref - momentumf * cLocal);
416 }
417
418 return;
419 }
420
421 /**************************************************************************
422 * Subroutine: UpdateElemInfo
423 * Purpose : Q(elem) = Q(elem) + deltaQ(elem)
424 *
425 * deltaQ(elem) = - (F(downWindFace) - F(upWindFace)) * dt / dx ;
426 *
427 *************************************************************************/
428
UpdateElemInfo(Group * const prob)429 void UpdateElemInfo(Group* const prob)
430 {
431 int i;
432
433 /* get the element quantities we want to update */
434 Group* const elem = prob->getGroup("elem");
435 float64* const mass = elem->getView("mass")->getData();
436 float64* const momentum = elem->getView("momentum")->getData();
437 float64* const energy = elem->getView("energy")->getData();
438 float64* const pressure = elem->getView("pressure")->getData();
439
440 /* focus on just the elements within the shock tube */
441 Group* const tube = elem->getGroup("tube");
442 int32* const elemToFace = tube->getView("elemToFace")->getData();
443
444 // Relation &elemToFace = *tube->relation("elemToFace");
445 int numTubeElems = tube->getView("mapToElems")->getNumElements();
446
447 // int *is = tube->map();
448 int32* const is = tube->getView("mapToElems")->getData();
449
450 /* The element update is calculated as the flux between faces */
451 Group* const face = prob->getGroup("face");
452 float64* const F0 = face->getView("F0")->getData();
453 float64* const F1 = face->getView("F1")->getData();
454 float64* const F2 = face->getView("F2")->getData();
455
456 double const dx = prob->getView("dx")->getData();
457 double const dt = prob->getView("dt")->getData();
458 float64& time = *prob->getView("time")->getData<float64*>();
459
460 for(i = 0; i < numTubeElems; ++i)
461 {
462 /* recalculate elements in the shocktube, don't touch inflow/outflow */
463 int elemIdx = is[i];
464
465 /* each element inside the tube has an upwind and downwind face */
466 int upWind = elemToFace[i * 2 + UPWIND]; /* upwind face */
467 int downWind = elemToFace[i * 2 + DOWNWIND]; /* downwind face */
468
469 mass[elemIdx] -= gammaaInverse * (F0[downWind] - F0[upWind]) * dt / dx;
470 momentum[elemIdx] -= gammaaInverse * (F1[downWind] - F1[upWind]) * dt / dx;
471 energy[elemIdx] -= gammaaInverse * (F2[downWind] - F2[upWind]) * dt / dx;
472 pressure[elemIdx] = (gammaa - 1.0) *
473 (energy[elemIdx] -
474 0.5 * momentum[elemIdx] * momentum[elemIdx] / mass[elemIdx]);
475 }
476
477 /* update the time */
478 time += dt;
479
480 return;
481 }
482
483 #include <stdio.h>
484 #include <string.h>
485 #include <ctype.h>
486
DumpUltra(Group * const prob)487 void DumpUltra(Group* const prob)
488 {
489 #if 1
490 FILE* fp;
491 char fname[100];
492 char* tail;
493
494 // VHashTraverse_t content ;
495
496 strcpy(fname, "problem");
497
498 /* Skip past the junk */
499 for(tail = fname; isalpha(*tail); ++tail)
500 ;
501
502 sprintf(tail, "_%04d.ult", prob->getView("cycle")->getData<int>());
503
504 if((fp = fopen(fname, "w")) == NULL)
505 {
506 printf("Could not open file %s. Aborting.\n", fname);
507 exit(-1);
508 }
509
510 fprintf(fp, "# Ultra Plot File\n");
511 fprintf(fp, "# Problem: %s\n", "problem");
512
513 for(IndexType i = 0; i < prob->getNumViews(); i++)
514 {
515 View* const view = prob->getView(i);
516 const int length = view->getNumElements();
517 const std::string& name = view->getName();
518 if(length <= 1)
519 {
520 if(view->getTypeID() == axom::sidre::INT32_ID)
521 {
522 fprintf(fp, "# %s = %d\n", name.c_str(), view->getData<int>());
523 }
524 else if(view->getTypeID() == axom::sidre::FLOAT64_ID)
525 {
526 fprintf(fp, "# %s = %f\n", name.c_str(), view->getData<double>());
527 }
528 }
529 }
530
531 Group* const elem = prob->getGroup("elem");
532
533 for(IndexType i = 0; i < elem->getNumViews(); i++)
534 {
535 View* const view = elem->getView(i);
536 const int length = view->getNumElements();
537 const std::string& name = view->getName();
538 fprintf(fp, "# %s\n", name.c_str());
539
540 if(view->getTypeID() == axom::sidre::INT32_ID)
541 {
542 int32 const* const data = view->getData();
543 for(int j = 0; j < length; ++j)
544 {
545 fprintf(fp, "%f %f\n", (double)j, (double)data[j]);
546 }
547 fprintf(fp, "\n");
548 }
549 else if(view->getTypeID() == axom::sidre::FLOAT64_ID)
550 {
551 float64 const* const data = view->getData();
552 for(int j = 0; j < length; ++j)
553 {
554 fprintf(fp, "%f %f\n", (double)j, (double)data[j]);
555 }
556 fprintf(fp, "\n");
557 }
558 }
559
560 fclose(fp);
561 #endif
562 return;
563 }
564 /**************************************************************************
565 * Subroutine: main
566 * Purpose : Simulate a 1D Shock Tube using split flux Euler formulation
567 *************************************************************************/
568
main(void)569 int main(void)
570 {
571 // extern void DumpUltra(Group * const prob);
572
573 /* We should be able to parallelize pretty easily by */
574 /* adding an MPI_Init() here, modifying the setup slightly, */
575 /* adding a communication subroutine in the main loop, */
576 /* and calling MPI_Finalize() at the end of main() */
577
578 DataStore DATASTORE;
579 DataStore* const dataStore = &DATASTORE;
580 Group* const rootGroup = dataStore->getRoot();
581 Group* const prob = rootGroup->createGroup("problem");
582
583 GetUserInput(prob);
584 CreateShockTubeMesh(prob);
585 InitializeShockTube(prob);
586
587 /* use a reference when you want to update the param directly */
588 int* const currCycle = prob->getView("cycle")->getData();
589
590 int numTotalCycles = prob->getView("numTotalCycles")->getData();
591 int dumpInterval = prob->getView("numCyclesPerDump")->getData();
592
593 std::cout << "Starting problem run." << std::endl;
594
595 for(*currCycle = 0; *currCycle < numTotalCycles; ++(*currCycle))
596 {
597 std::cout << " cycle " << *currCycle << std::endl;
598 /* dump the ultra file, based on the user chosen attribute mask */
599 if((*currCycle) % dumpInterval == 0)
600 {
601 DumpUltra(prob);
602 }
603
604 ComputeFaceInfo(prob);
605 UpdateElemInfo(prob);
606 }
607
608 // DumpUltra(prob); /* One last dump */
609
610 std::cout << "Finished problem run." << std::endl;
611 return 0;
612 }
613