PlasCom2  1.0
XPACC Multi-physics simluation application
TestIntegrated.C
Go to the documentation of this file.
1 #include "Testing.H"
2 #include "PlasCom2.H"
3 #include "Simulation.H"
4 #include "RKTestFixtures.H"
5 #include "EulerRHS.H"
6 #include "Stencil.H"
7 #include "PCPPCommUtil.H"
8 #include "PCPPReport.H"
9 #include "PCPPIntervalUtils.H"
10 #include "TestFixtures.H"
11 #include "EulerTestFixtures.H"
12 #include "ViscidTestFixtures.H"
13 #include "PC2Util.H"
14 
15 #include <cmath>
16 
18 
22 
24  pcpp::CommunicatorType &testComm)
25 {
26 
27  bool PoiseuilleTestResult = true;
28  int myRank = testComm.Rank();
29  int numProc = testComm.NProc();
30 
31  const std::string testFunctionName("TestIntegrated_Poiseuille2DX");
32  const std::string testSuiteName("ViscidTest");
33  const std::string testCaseName("Poiseuille2DX");
34  const std::string testDirectory(testSuiteName+"/"+testCaseName);
35  const std::string originalDirectory(ix::sys::CWD());
36  std::ostringstream compareStream;
37  double errorTolerancePC1_2 = 1.e-4;
38  double errorToleranceExact = 2.e-5;
39  if(myRank == 0)
40  std::cout << "Running " << testFunctionName << std::endl;
41 
42  int argc = 4;
43  const char *argv[] = {"plascom2x","-c","poiseuille2d.config",NULL};
44 
45  pcpp::CommunicatorType nonDimenComm2(testComm);
46  plascom2::application simulationApplicationNonDimen2(argc,const_cast<char **>(argv),nonDimenComm2.Comm());
47 
48  int returnValue = 0;
49 
50  if(myRank == 0) {
51  if(!ix::sys::FILEEXISTS(testDirectory)){
52  PoiseuilleTestResult = false;
53  std::cout << "TestIntegrated_Poiseuille2D: testing directory ("
54  << testDirectory << ") did not exist. Aborting test."
55  << std::endl;
56  }
57  }
58  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
59  if(!PoiseuilleTestResult)
60  goto endTest;
61  ix::sys::ChDir(testDirectory);
62 
63  // nonDimensional run using PlasComCM (rho, C, T) non-dimensionalization
64  {
65  ix::sys::ChDir("nonDimenPC1");
66 
67  if(myRank == 0)
68  std::cout << "Running nonDimensional PlasCom2... (PC1ND)" << std::endl;
69  returnValue = application::ApplicationDriver(simulationApplicationNonDimen2);
70  if(myRank == 0)
71  std::cout << "Done running nonDimensional PlasCom2 (PC1ND)." << std::endl;
72 
73  if(returnValue > 0){
74  std::cout << "TestIntegrated_Poiseuille2DX: PlasCom2 returned error code ("
75  << returnValue << "). Test failed." << std::endl;
76  PoiseuilleTestResult = false;
77  goto endTest;
78  }
79 
80  if(myRank == 0){
81  if(!ix::sys::FILEEXISTS("PlasCom2_000006000.h5")){
82  std::cout << "TestIntegrated_Poiseuille2DX: PlasCom2 failed to produce expected" << std::endl
83  << "output file (PlasCom2_000006000.h5). Test failed." << std::endl;
84  PoiseuilleTestResult = false;
85  }
86  }
87  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
88  if(!PoiseuilleTestResult)
89  goto endTest;
90  }
91 
92  // get the domain from the application and copy the state
93  // make a new scope for variable declaration, needed for goto usage :-(
94  {
95  DomainVector &simulationDomain(simulationApplicationNonDimen2.GetDomains());
96 
97  // generate a state containing the exact solution and write it to a file
98  int direction=0;
99  // single grid
100  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
101  state_t exactState(simulationDomain[0].State(0));
102  if(testfixtures::viscid::GeneratePoiseuilleExact(exactGrid,exactState,direction)){
103  PoiseuilleTestResult = false;
104  std::cout << "Fail to calculated exact poiseuille solution" << std::endl;
105  }
106  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
107  if(!PoiseuilleTestResult)
108  goto endTest;
109 
110  size_t numPointsBuffer = exactGrid.BufferSize();
111 
112  int rhoHandle = exactState.GetDataIndex("rho");
113  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
114  double *rhoPtr = rhoData.Data<double>();
115 
116  int rhoVHandle = exactState.GetDataIndex("rhoV");
117  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
118  double *rhoVPtr = rhoVData.Data<double>();
119 
120  int rhoEHandle = exactState.GetDataIndex("rhoE");
121  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
122  double *rhoEPtr = rhoEData.Data<double>();
123 
124  int velocityHandle = exactState.GetDataIndex("velocity");
125  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
126  double *velocityPtr = velocityData.Data<double>();
127 
128  int pressureHandle = exactState.GetDataIndex("pressure");
129  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
130  double *pressurePtr = pressureData.Data<double>();
131 
132  int temperatureHandle = exactState.GetDataIndex("temperature");
133  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
134  double *temperaturePtr = temperatureData.Data<double>();
135 
136  {
137  if(myRank == 0) {
138  std::cout << "Writing exact solution for nonDimensional PlasCom2"
139  << " with legacy PlasComCM nonDimensionalization" << std::endl;
140  }
141  //
142  // output exact nonDimensional solution in the PlasComCM dimensionalization
143  //
144  // we are using a non-unity length reference, which is embedded in the mesh
145  // modify the exact solution for velocity to account for this
146  //
147  double gamma = 1.4;
148  double refRho = 1.225;
149  double refPress = 101325;
150  double refSndSpd = sqrt(gamma*refPress/refRho);
151  double refTemp = (gamma-1)*288.15;
152  double refLength = 0.5;
153  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
154  // subtract out kinetic energy
155  rhoEPtr[iPoint] -= 0.5*rhoVPtr[iPoint]*rhoVPtr[iPoint]/rhoPtr[iPoint];
156 
157  // modify velocity and momentum by reference length
158  rhoPtr[iPoint] /= refRho;
159  rhoVPtr[iPoint] /= refRho*refSndSpd/refLength;
160 
161  // scale internal energy and add kinetic energy back
162  rhoEPtr[iPoint] /= refSndSpd*refSndSpd*refRho;
163  rhoEPtr[iPoint] += 0.5*rhoVPtr[iPoint]*rhoVPtr[iPoint]/rhoPtr[iPoint];
164 
165  velocityPtr[iPoint] /= refSndSpd;
166  pressurePtr[iPoint] /= refRho*refSndSpd*refSndSpd;
167  temperaturePtr[iPoint] /= refTemp;
168  }
169 
170  std::ostringstream Ostr;
171  Ostr << "exact_Poiseuille_PC1ND_6000";
172  const std::string domainName(Ostr.str());
173  const std::string hdfFileName(domainName+".h5");
174  const std::string xdmfFileName(domainName+".xdmf");
175  const std::string geometryName("poiseuille2d");
176  const std::string gridName("grid1");
177  if(ix::sys::FILEEXISTS(hdfFileName)){
178  ix::sys::Remove(hdfFileName);
179  ix::sys::Remove(xdmfFileName);
180  }
181  fixtures::ConfigurationType testConfig;
182  state_t paramState;
183  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
184  gridName,exactGrid,exactState,paramState,
185  testConfig,compareStream)){
186  std::cout << compareStream.str();
187  std::cout << "Output to HDF5 failed:\n";
188  std::cout << "Continuing....\n";
189 
190  }
191  pcpp::io::RenewStream(compareStream);
192  }
193  ix::sys::ChDir("../");
194  }
195 
196  if(myRank == 0) {
197 
198  if(!ix::sys::FILEEXISTS("nonDimenPC1/exact_Poiseuille_PC1ND_6000.h5")){
199  std::cout << "TestIntegrated_Poiseuille2DX: failed to produce expected" << std::endl
200  << "output file (exact_Poiseuille_PC1ND_6000.h5). Test failed." << std::endl;
201  PoiseuilleTestResult = false;
202  }
203  }
204  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
205  if(!PoiseuilleTestResult)
206  goto endTest;
207 
208  {
209  if(myRank == 0) {
210  std::cout << "Comparing exact solution against nonDimensional PlasCom2"
211  << " with legacy PlasComCM nonDimensionalization" << std::endl;
212  ix::sys::ChDir("nonDimenPC1");
213  returnValue = plascom2::util::PC2Compare("PlasCom2_000006000.h5",
214  "exact_Poiseuille_PC1ND_6000.h5",
215  errorToleranceExact,compareStream);
216  if(returnValue > 0){
217  std::cout << "TestIntegrated_Poiseuille2DX: Resulting state comparison yielded negative result:"
218  << std::endl
219  << compareStream.str()
220  << std::endl;
221  PoiseuilleTestResult = false;
222  } else {
223  std::cout << "Passed!" << std::endl;
224  }
225  pcpp::io::RenewStream(compareStream);
226  ix::sys::ChDir("../");
227  }
228  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
229  testComm.Barrier();
230  }
231 
232  endTest:
233  {
234  ix::sys::ChDir(originalDirectory);
235  parallelUnitResults.UpdateResult("Integrated:Viscid:Poiseuille2DX",PoiseuilleTestResult);
236  }
237 
238  return;
239 }
240 
242  pcpp::CommunicatorType &testComm)
243 {
244 
245  bool PoiseuilleTestResult = true;
246  int myRank = testComm.Rank();
247  int numProc = testComm.NProc();
248 
249  const std::string testFunctionName("TestIntegrated_PFRectilinear2DX");
250  const std::string testSuiteName("ViscidTest");
251  const std::string testCaseName("Poiseuille2DX");
252  const std::string testDirectory(testSuiteName+"/"+testCaseName);
253  const std::string originalDirectory(ix::sys::CWD());
254  std::ostringstream compareStream;
255  double errorTolerancePC1_2 = 1.e-4;
256  double errorToleranceExact = 2.e-5;
257  if(myRank == 0)
258  std::cout << "Running " << testFunctionName << std::endl;
259 
260  int argc = 4;
261  const char *argv[] = {"plascom2x","-c","rectilinear.config",NULL};
262 
263  pcpp::CommunicatorType nonDimenComm2(testComm);
264  plascom2::application simulationApplicationNonDimen2(argc,const_cast<char **>(argv),nonDimenComm2.Comm());
265 
266  int returnValue = 0;
267 
268  if(myRank == 0) {
269  if(!ix::sys::FILEEXISTS(testDirectory)){
270  PoiseuilleTestResult = false;
271  std::cout << testFunctionName << ": testing directory ("
272  << testDirectory << ") did not exist. Aborting test."
273  << std::endl;
274  }
275  }
276  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
277  if(!PoiseuilleTestResult)
278  goto endTest;
279  ix::sys::ChDir(testDirectory);
280 
281  // nonDimensional run using PlasComCM (rho, C, T) non-dimensionalization
282  {
283  ix::sys::ChDir("nonDimenPC1");
284 
285  if(myRank == 0){
286  if(ix::sys::FILEEXISTS("PlasCom2_000006000.h5"))
287  ix::sys::Remove("PlasCom2_000006000.h5");
288  std::cout << "Running nonDimensional PlasCom2... (PC1ND)" << std::endl;
289  }
290  testComm.Barrier();
291 
292  returnValue = application::ApplicationDriver(simulationApplicationNonDimen2);
293 
294  testComm.Barrier();
295 
296  if(myRank == 0)
297  std::cout << "Done running nonDimensional PlasCom2 (PC1ND)." << std::endl;
298 
299  if(returnValue > 0){
300  std::cout << testFunctionName << ": PlasCom2 returned error code ("
301  << returnValue << "). Test failed." << std::endl;
302  PoiseuilleTestResult = false;
303  goto endTest;
304  }
305 
306  if(myRank == 0){
307  if(!ix::sys::FILEEXISTS("PlasCom2_000006000.h5")){
308  std::cout << testFunctionName << ": PlasCom2 failed to produce expected" << std::endl
309  << "output file (PlasCom2_000006000.h5). Test failed." << std::endl;
310  PoiseuilleTestResult = false;
311  }
312  }
313  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
314  if(!PoiseuilleTestResult)
315  goto endTest;
316  }
317 
318  // get the domain from the application and copy the state
319  // make a new scope for variable declaration, needed for goto usage :-(
320  {
321  DomainVector &simulationDomain(simulationApplicationNonDimen2.GetDomains());
322 
323  // generate a state containing the exact solution and write it to a file
324  int direction=0;
325  // single grid
326  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
327  state_t exactState(simulationDomain[0].State(0));
328  if(testfixtures::viscid::GeneratePoiseuilleExact(exactGrid,exactState,direction)){
329  PoiseuilleTestResult = false;
330  std::cout << "Fail to calculated exact poiseuille solution" << std::endl;
331  }
332  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
333  if(!PoiseuilleTestResult)
334  goto endTest;
335 
336  size_t numPointsBuffer = exactGrid.BufferSize();
337 
338  int rhoHandle = exactState.GetDataIndex("rho");
339  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
340  double *rhoPtr = rhoData.Data<double>();
341 
342  int rhoVHandle = exactState.GetDataIndex("rhoV");
343  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
344  double *rhoVPtr = rhoVData.Data<double>();
345 
346  int rhoEHandle = exactState.GetDataIndex("rhoE");
347  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
348  double *rhoEPtr = rhoEData.Data<double>();
349 
350  int velocityHandle = exactState.GetDataIndex("velocity");
351  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
352  double *velocityPtr = velocityData.Data<double>();
353 
354  int pressureHandle = exactState.GetDataIndex("pressure");
355  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
356  double *pressurePtr = pressureData.Data<double>();
357 
358  int temperatureHandle = exactState.GetDataIndex("temperature");
359  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
360  double *temperaturePtr = temperatureData.Data<double>();
361 
362  {
363  if(myRank == 0) {
364  std::cout << "Writing exact solution for nonDimensional PlasCom2"
365  << " with legacy PlasComCM nonDimensionalization" << std::endl;
366  }
367  //
368  // output exact nonDimensional solution in the PlasComCM dimensionalization
369  //
370  // we are using a non-unity length reference, which is embedded in the mesh
371  // modify the exact solution for velocity to account for this
372  //
373  double gamma = 1.4;
374  double refRho = 1.225;
375  double refPress = 101325;
376  double refSndSpd = sqrt(gamma*refPress/refRho);
377  double refTemp = (gamma-1)*288.15;
378  double refLength = 0.5;
379  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
380  // subtract out kinetic energy
381  rhoEPtr[iPoint] -= 0.5*rhoVPtr[iPoint]*rhoVPtr[iPoint]/rhoPtr[iPoint];
382 
383  // modify velocity and momentum by reference length
384  rhoPtr[iPoint] /= refRho;
385  rhoVPtr[iPoint] /= refRho*refSndSpd/refLength;
386 
387  // scale internal energy and add kinetic energy back
388  rhoEPtr[iPoint] /= refSndSpd*refSndSpd*refRho;
389  rhoEPtr[iPoint] += 0.5*rhoVPtr[iPoint]*rhoVPtr[iPoint]/rhoPtr[iPoint];
390 
391  velocityPtr[iPoint] /= refSndSpd;
392  pressurePtr[iPoint] /= refRho*refSndSpd*refSndSpd;
393  temperaturePtr[iPoint] /= refTemp;
394  }
395 
396  std::ostringstream Ostr;
397  Ostr << "exact_Poiseuille_PC1ND_6000";
398  const std::string domainName(Ostr.str());
399  const std::string hdfFileName(domainName+".h5");
400  const std::string xdmfFileName(domainName+".xdmf");
401  const std::string geometryName("poiseuille2d");
402  const std::string gridName("grid1");
403  if(ix::sys::FILEEXISTS(hdfFileName)){
404  ix::sys::Remove(hdfFileName);
405  ix::sys::Remove(xdmfFileName);
406  }
407  fixtures::ConfigurationType testConfig;
408  state_t paramState;
409  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
410  gridName,exactGrid,exactState,paramState,
411  testConfig,compareStream)){
412  std::cout << compareStream.str();
413  std::cout << "Output to HDF5 failed:\n";
414  std::cout << "Continuing....\n";
415 
416  }
417  pcpp::io::RenewStream(compareStream);
418  }
419  ix::sys::ChDir("../");
420  }
421 
422  if(myRank == 0) {
423 
424  if(!ix::sys::FILEEXISTS("nonDimenPC1/exact_Poiseuille_PC1ND_6000.h5")){
425  std::cout << testFunctionName << ": failed to produce expected" << std::endl
426  << "output file (exact_Poiseuille_PC1ND_6000.h5). Test failed." << std::endl;
427  PoiseuilleTestResult = false;
428  }
429  }
430  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
431  if(!PoiseuilleTestResult)
432  goto endTest;
433 
434  {
435  if(myRank == 0) {
436  std::cout << "Comparing exact solution against nonDimensional PlasCom2"
437  << " with legacy PlasComCM nonDimensionalization" << std::endl;
438  ix::sys::ChDir("nonDimenPC1");
439  returnValue = plascom2::util::PC2Compare("PlasCom2_000006000.h5",
440  "exact_Poiseuille_PC1ND_6000.h5",
441  errorToleranceExact,compareStream);
442  if(returnValue > 0){
443  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
444  << std::endl
445  << compareStream.str()
446  << std::endl;
447  PoiseuilleTestResult = false;
448  } else {
449  std::cout << "Passed!" << std::endl;
450  }
451  pcpp::io::RenewStream(compareStream);
452  ix::sys::ChDir("../");
453  }
454  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
455  testComm.Barrier();
456  }
457 
458  endTest:
459  {
460  ix::sys::ChDir(originalDirectory);
461  parallelUnitResults.UpdateResult("Integrated:Viscid:PFRectilinear2DX",PoiseuilleTestResult);
462  }
463 
464  return;
465 }
466 
468  pcpp::CommunicatorType &testComm)
469 {
470 
471  bool PoiseuilleTestResult = true;
472  int myRank = testComm.Rank();
473  int numProc = testComm.NProc();
474 
475  const std::string testFunctionName("TestIntegrated_PFCurvilinear2DX");
476  const std::string testSuiteName("ViscidTest");
477  const std::string testCaseName("Poiseuille2DX");
478  const std::string testDirectory(testSuiteName+"/"+testCaseName);
479  const std::string originalDirectory(ix::sys::CWD());
480  std::ostringstream compareStream;
481  double errorToleranceExact = 9.e-3;
482  // MJA
483  //double errorTolerancePC1_2 = 1.e-4;
484  //double errorToleranceExact = 2.e-5;
485  if(myRank == 0)
486  std::cout << "Running " << testFunctionName << std::endl;
487 
488  int argc = 4;
489  const char *argv[] = {"plascom2x","-c","curvilinear.config",NULL};
490 
491  pcpp::CommunicatorType nonDimenComm2(testComm);
492  plascom2::application simulationApplicationNonDimen2(argc,const_cast<char **>(argv),nonDimenComm2.Comm());
493 
494  int returnValue = 0;
495 
496  if(myRank == 0) {
497  if(!ix::sys::FILEEXISTS(testDirectory)){
498  PoiseuilleTestResult = false;
499  std::cout << testFunctionName << ": testing directory ("
500  << testDirectory << ") did not exist. Aborting test."
501  << std::endl;
502  }
503  }
504  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
505  if(!PoiseuilleTestResult)
506  goto endTest;
507  ix::sys::ChDir(testDirectory);
508 
509  // nonDimensional run using PlasComCM (rho, C, T) non-dimensionalization
510  {
511  ix::sys::ChDir("nonDimenPC1");
512 
513  if(myRank == 0){
514  if(ix::sys::FILEEXISTS("PlasCom2_000006000.h5"))
515  ix::sys::Remove("PlasCom2_000006000.h5");
516  std::cout << "Running nonDimensional PlasCom2... (PC1ND)" << std::endl;
517  }
518  testComm.Barrier();
519 
520  returnValue = application::ApplicationDriver(simulationApplicationNonDimen2);
521 
522  testComm.Barrier();
523 
524  if(myRank == 0)
525  std::cout << "Done running nonDimensional PlasCom2 (PC1ND)." << std::endl;
526 
527  if(returnValue > 0){
528  std::cout << testFunctionName << ": PlasCom2 returned error code ("
529  << returnValue << "). Test failed." << std::endl;
530  PoiseuilleTestResult = false;
531  goto endTest;
532  }
533 
534  if(myRank == 0){
535  if(!ix::sys::FILEEXISTS("PlasCom2_000006000.h5")){
536  std::cout << testFunctionName << ": PlasCom2 failed to produce expected" << std::endl
537  << "output file (PlasCom2_000006000.h5). Test failed." << std::endl;
538  PoiseuilleTestResult = false;
539  }
540  }
541  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
542  if(!PoiseuilleTestResult)
543  goto endTest;
544  }
545 
546  // get the domain from the application and copy the state
547  // make a new scope for variable declaration, needed for goto usage :-(
548  {
549  DomainVector &simulationDomain(simulationApplicationNonDimen2.GetDomains());
550 
551  // generate a state containing the exact solution and write it to a file
552  int direction=0;
553  // single grid
554  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
555  state_t exactState(simulationDomain[0].State(0));
556  if(testfixtures::viscid::GeneratePoiseuilleExact(exactGrid,exactState,direction)){
557  PoiseuilleTestResult = false;
558  std::cout << "Fail to calculated exact poiseuille solution" << std::endl;
559  }
560  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
561  if(!PoiseuilleTestResult)
562  goto endTest;
563 
564  size_t numPointsBuffer = exactGrid.BufferSize();
565 
566  int rhoHandle = exactState.GetDataIndex("rho");
567  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
568  double *rhoPtr = rhoData.Data<double>();
569 
570  int rhoVHandle = exactState.GetDataIndex("rhoV");
571  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
572  double *rhoVPtr = rhoVData.Data<double>();
573 
574  int rhoEHandle = exactState.GetDataIndex("rhoE");
575  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
576  double *rhoEPtr = rhoEData.Data<double>();
577 
578  int velocityHandle = exactState.GetDataIndex("velocity");
579  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
580  double *velocityPtr = velocityData.Data<double>();
581 
582  int pressureHandle = exactState.GetDataIndex("pressure");
583  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
584  double *pressurePtr = pressureData.Data<double>();
585 
586  int temperatureHandle = exactState.GetDataIndex("temperature");
587  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
588  double *temperaturePtr = temperatureData.Data<double>();
589 
590  {
591  if(myRank == 0) {
592  std::cout << "Writing exact solution for nonDimensional PlasCom2"
593  << " with legacy PlasComCM nonDimensionalization" << std::endl;
594  }
595  //
596  // output exact nonDimensional solution in the PlasComCM dimensionalization
597  //
598  // we are using a non-unity length reference, which is embedded in the mesh
599  // modify the exact solution for velocity to account for this
600  //
601  double gamma = 1.4;
602  double refRho = 1.225;
603  double refPress = 101325;
604  double refSndSpd = sqrt(gamma*refPress/refRho);
605  double refTemp = (gamma-1)*288.15;
606  double refLength = 0.5;
607  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
608  // subtract out kinetic energy
609  rhoEPtr[iPoint] -= 0.5*rhoVPtr[iPoint]*rhoVPtr[iPoint]/rhoPtr[iPoint];
610 
611  // modify velocity and momentum by reference length
612  rhoPtr[iPoint] /= refRho;
613  rhoVPtr[iPoint] /= refRho*refSndSpd/refLength;
614 
615  // scale internal energy and add kinetic energy back
616  rhoEPtr[iPoint] /= refSndSpd*refSndSpd*refRho;
617  rhoEPtr[iPoint] += 0.5*rhoVPtr[iPoint]*rhoVPtr[iPoint]/rhoPtr[iPoint];
618 
619  velocityPtr[iPoint] /= refSndSpd;
620  pressurePtr[iPoint] /= refRho*refSndSpd*refSndSpd;
621  temperaturePtr[iPoint] /= refTemp;
622  }
623 
624  std::ostringstream Ostr;
625  Ostr << "exact_Poiseuille_PC1ND_6000";
626  const std::string domainName(Ostr.str());
627  const std::string hdfFileName(domainName+".h5");
628  const std::string xdmfFileName(domainName+".xdmf");
629  const std::string geometryName("poiseuille2d");
630  const std::string gridName("grid1");
631  if(ix::sys::FILEEXISTS(hdfFileName)){
632  ix::sys::Remove(hdfFileName);
633  ix::sys::Remove(xdmfFileName);
634  }
635  fixtures::ConfigurationType testConfig;
636  state_t paramState;
637  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
638  gridName,exactGrid,exactState,paramState,
639  testConfig,compareStream)){
640  std::cout << compareStream.str();
641  std::cout << "Output to HDF5 failed:\n";
642  std::cout << "Continuing....\n";
643 
644  }
645  pcpp::io::RenewStream(compareStream);
646  }
647  ix::sys::ChDir("../");
648  }
649 
650  if(myRank == 0) {
651 
652  if(!ix::sys::FILEEXISTS("nonDimenPC1/exact_Poiseuille_PC1ND_6000.h5")){
653  std::cout << testFunctionName << ": failed to produce expected" << std::endl
654  << "output file (exact_Poiseuille_PC1ND_6000.h5). Test failed." << std::endl;
655  PoiseuilleTestResult = false;
656  }
657  }
658  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
659  if(!PoiseuilleTestResult)
660  goto endTest;
661 
662  {
663  if(myRank == 0) {
664  std::cout << "Comparing exact solution against nonDimensional PlasCom2"
665  << " with legacy PlasComCM nonDimensionalization" << std::endl;
666  ix::sys::ChDir("nonDimenPC1");
667  returnValue = plascom2::util::PC2Compare("PlasCom2_000006000.h5",
668  "exact_Poiseuille_PC1ND_6000.h5",
669  errorToleranceExact,compareStream);
670  if(returnValue > 0){
671  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
672  << std::endl
673  << compareStream.str()
674  << std::endl;
675  PoiseuilleTestResult = false;
676  } else {
677  std::cout << "Passed!" << std::endl;
678  }
679  pcpp::io::RenewStream(compareStream);
680  ix::sys::ChDir("../");
681  }
682  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
683  testComm.Barrier();
684  }
685 
686  endTest:
687  {
688  ix::sys::ChDir(originalDirectory);
689  parallelUnitResults.UpdateResult("Integrated:Viscid:PFCurvilinear2DX",PoiseuilleTestResult);
690  }
691 
692  return;
693 }
694 
695 
697  pcpp::CommunicatorType &testComm)
698 {
699 
700  bool PoiseuilleTestResult = true;
701  int myRank = testComm.Rank();
702  int numProc = testComm.NProc();
703 
704  const std::string testFunctionName("TestIntegrated_Poiseuille2DY");
705  const std::string testSuiteName("ViscidTest");
706  const std::string testCaseName("Poiseuille2DY");
707  const std::string testDirectory(testSuiteName+"/"+testCaseName);
708  const std::string originalDirectory(ix::sys::CWD());
709  std::ostringstream compareStream;
710  double errorToleranceExact = 2.e-5;
711  if(myRank == 0)
712  std::cout << "Running " << testFunctionName << std::endl;
713 
714  int argc = 4;
715  const char *argv[] = {"plascom2x","-c","poiseuille2d.config",NULL};
716 
717  pcpp::CommunicatorType dimenComm(testComm);
718  plascom2::application simulationApplicationDimen(argc,const_cast<char **>(argv),dimenComm.Comm());
719 
720  int returnValue = 0;
721 
722  if(myRank == 0) {
723  if(!ix::sys::FILEEXISTS(testDirectory)){
724  PoiseuilleTestResult = false;
725  std::cout << "TestIntegrated_Poiseuille2D: testing directory ("
726  << testDirectory << ") did not exist. Aborting test."
727  << std::endl;
728  }
729  }
730  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
731  if(!PoiseuilleTestResult)
732  goto endTest;
733  ix::sys::ChDir(testDirectory);
734 
735  // simulation in dimensional coordinates
736  {
737  ix::sys::ChDir("dimen");
738 
739  if(myRank == 0)
740  std::cout << "Running dimensional PlasCom2..." << std::endl;
741  returnValue = application::ApplicationDriver(simulationApplicationDimen);
742  if(myRank == 0)
743  std::cout << "Done running dimensional PlasCom2." << std::endl;
744 
745  if(returnValue > 0){
746  std::cout << "TestIntegrated_Poiseuille2DY: PlasCom2 returned error code ("
747  << returnValue << "). Test failed." << std::endl;
748  PoiseuilleTestResult = false;
749  goto endTest;
750  }
751 
752  if(myRank == 0){
753  if(!ix::sys::FILEEXISTS("PlasCom2_000003000.h5")){
754  std::cout << "TestIntegrated_Poiseuille2DY: PlasCom2 failed to produce expected"
755  << "output file (PlasCom2_000003000.h5). Test failed." << std::endl;
756  PoiseuilleTestResult = false;
757  }
758  }
759  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
760  if(!PoiseuilleTestResult)
761  goto endTest;
762  }
763 
764  // get the domain from the application and copy the state
765  // make a new scope for variable declaration, needed for goto usage :-(
766  {
767  DomainVector &simulationDomain(simulationApplicationDimen.GetDomains());
768 
769  // generate a state containing the exact solution and write it to a file
770  int direction=1;
771  // single grid
772  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
773  state_t exactState(simulationDomain[0].State(0));
774  if(testfixtures::viscid::GeneratePoiseuilleExact(exactGrid,exactState,direction)){
775  PoiseuilleTestResult = false;
776  std::cout << "Fail to calculate exact poiseuille solution" << std::endl;
777  }
778  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
779  if(!PoiseuilleTestResult)
780  goto endTest;
781 
782  size_t numPointsBuffer = exactGrid.BufferSize();
783 
784  int rhoHandle = exactState.GetDataIndex("rho");
785  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
786  double *rhoPtr = rhoData.Data<double>();
787 
788  int rhoVHandle = exactState.GetDataIndex("rhoV");
789  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
790  double *rhoVPtr = rhoVData.Data<double>();
791 
792  int rhoEHandle = exactState.GetDataIndex("rhoE");
793  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
794  double *rhoEPtr = rhoEData.Data<double>();
795 
796  int velocityHandle = exactState.GetDataIndex("velocity");
797  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
798  double *velocityPtr = velocityData.Data<double>();
799 
800  int pressureHandle = exactState.GetDataIndex("pressure");
801  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
802  double *pressurePtr = pressureData.Data<double>();
803 
804  int temperatureHandle = exactState.GetDataIndex("temperature");
805  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
806  double *temperaturePtr = temperatureData.Data<double>();
807 
808  {
809  // output exact dimensional solution
810  if(myRank == 0)
811  std::cout << "Writing exact solution for dimensional PlasCom2" << std::endl;
812  std::ostringstream Ostr;
813  Ostr << "exact_Poiseuille_dimen_3000";
814  const std::string domainName(Ostr.str());
815  const std::string hdfFileName(domainName+".h5");
816  const std::string xdmfFileName(domainName+".xdmf");
817  const std::string geometryName("poiseuille");
818  const std::string gridName("grid1");
819  if(ix::sys::FILEEXISTS(hdfFileName)){
820  ix::sys::Remove(hdfFileName);
821  ix::sys::Remove(xdmfFileName);
822  }
823  fixtures::ConfigurationType testConfig;
824  state_t paramState;
825  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
826  gridName,exactGrid,exactState,paramState,
827  testConfig,compareStream)){
828  std::cout << compareStream.str();
829  std::cout << "Output to HDF5 failed:\n";
830  std::cout << "Continuing....\n";
831 
832  }
833  pcpp::io::RenewStream(compareStream);
834  }
835  ix::sys::ChDir("../");
836  }
837 
838  if(myRank == 0) {
839  if(!ix::sys::FILEEXISTS("dimen/exact_Poiseuille_dimen_3000.h5")){
840  std::cout << "TestIntegrated_Poiseuille2DY: failed to produce expected" << std::endl
841  << "output file (exact_Poiseuille_dimen_3000.h5). Test failed." << std::endl;
842  PoiseuilleTestResult = false;
843  }
844 
845  }
846  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
847  if(!PoiseuilleTestResult)
848  goto endTest;
849 
850  {
851  if(myRank == 0) {
852  std::cout << "Comparing exact solution against dimensional PlasCom2" << std::endl;
853  ix::sys::ChDir("dimen");
854  returnValue = plascom2::util::PC2Compare("PlasCom2_000003000.h5",
855  "exact_Poiseuille_dimen_3000.h5",
856  errorToleranceExact,compareStream);
857 
858  if(returnValue > 0){
859  std::cout << "TestIntegrated_Poiseuille2DY: Resulting state comparison yielded negative result:"
860  << std::endl
861  << compareStream.str()
862  << std::endl;
863  PoiseuilleTestResult = false;
864  } else {
865  std::cout << "Passed!" << std::endl;
866  }
867 
868  pcpp::io::RenewStream(compareStream);
869  ix::sys::ChDir("../");
870  }
871  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
872  }
873 
874  endTest:
875  {
876  ix::sys::ChDir(originalDirectory);
877  parallelUnitResults.UpdateResult("Integrated:Viscid:Poiseuille2DY",PoiseuilleTestResult);
878  }
879 
880  return;
881 }
882 
884  pcpp::CommunicatorType &testComm)
885 {
886 
887  bool PoiseuilleTestResult = true;
888  int myRank = testComm.Rank();
889  int numProc = testComm.NProc();
890 
891  const std::string testFunctionName("TestIntegrated_PFRectilinear2DY");
892  const std::string testSuiteName("ViscidTest");
893  const std::string testCaseName("Poiseuille2DY");
894  const std::string testDirectory(testSuiteName+"/"+testCaseName);
895  const std::string originalDirectory(ix::sys::CWD());
896  std::ostringstream compareStream;
897  double errorToleranceExact = 2.e-5;
898  if(myRank == 0)
899  std::cout << "Running " << testFunctionName << std::endl;
900 
901  int argc = 4;
902  const char *argv[] = {"plascom2x","-c","rectilinear.config",NULL};
903 
904  pcpp::CommunicatorType dimenComm(testComm);
905  plascom2::application simulationApplicationDimen(argc,const_cast<char **>(argv),dimenComm.Comm());
906 
907  int returnValue = 0;
908 
909  if(myRank == 0) {
910  if(!ix::sys::FILEEXISTS(testDirectory)){
911  PoiseuilleTestResult = false;
912  std::cout << testFunctionName << ": testing directory ("
913  << testDirectory << ") did not exist. Aborting test."
914  << std::endl;
915  }
916  }
917  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
918  if(!PoiseuilleTestResult)
919  goto endTest;
920  ix::sys::ChDir(testDirectory);
921 
922  // simulation in dimensional coordinates
923  {
924  ix::sys::ChDir("dimen");
925 
926  if(myRank == 0)
927  std::cout << "Running dimensional PlasCom2..." << std::endl;
928  returnValue = application::ApplicationDriver(simulationApplicationDimen);
929  if(myRank == 0)
930  std::cout << "Done running dimensional PlasCom2." << std::endl;
931 
932  if(returnValue > 0){
933  std::cout << testFunctionName << ": PlasCom2 returned error code ("
934  << returnValue << "). Test failed." << std::endl;
935  PoiseuilleTestResult = false;
936  goto endTest;
937  }
938 
939  if(myRank == 0){
940  if(!ix::sys::FILEEXISTS("PlasCom2_000003000.h5")){
941  std::cout << testFunctionName << ": PlasCom2 failed to produce expected"
942  << "output file (PlasCom2_000003000.h5). Test failed." << std::endl;
943  PoiseuilleTestResult = false;
944  }
945  }
946  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
947  if(!PoiseuilleTestResult)
948  goto endTest;
949  }
950 
951  // get the domain from the application and copy the state
952  // make a new scope for variable declaration, needed for goto usage :-(
953  {
954  DomainVector &simulationDomain(simulationApplicationDimen.GetDomains());
955 
956  // generate a state containing the exact solution and write it to a file
957  int direction=1;
958  // single grid
959  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
960  state_t exactState(simulationDomain[0].State(0));
961  if(testfixtures::viscid::GeneratePoiseuilleExact(exactGrid,exactState,direction)){
962  PoiseuilleTestResult = false;
963  std::cout << "Fail to calculate exact poiseuille solution" << std::endl;
964  }
965  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
966  if(!PoiseuilleTestResult)
967  goto endTest;
968 
969  size_t numPointsBuffer = exactGrid.BufferSize();
970 
971  int rhoHandle = exactState.GetDataIndex("rho");
972  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
973  double *rhoPtr = rhoData.Data<double>();
974 
975  int rhoVHandle = exactState.GetDataIndex("rhoV");
976  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
977  double *rhoVPtr = rhoVData.Data<double>();
978 
979  int rhoEHandle = exactState.GetDataIndex("rhoE");
980  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
981  double *rhoEPtr = rhoEData.Data<double>();
982 
983  int velocityHandle = exactState.GetDataIndex("velocity");
984  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
985  double *velocityPtr = velocityData.Data<double>();
986 
987  int pressureHandle = exactState.GetDataIndex("pressure");
988  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
989  double *pressurePtr = pressureData.Data<double>();
990 
991  int temperatureHandle = exactState.GetDataIndex("temperature");
992  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
993  double *temperaturePtr = temperatureData.Data<double>();
994 
995  {
996  // output exact dimensional solution
997  if(myRank == 0)
998  std::cout << "Writing exact solution for dimensional PlasCom2" << std::endl;
999  std::ostringstream Ostr;
1000  Ostr << "exact_Poiseuille_dimen_3000";
1001  const std::string domainName(Ostr.str());
1002  const std::string hdfFileName(domainName+".h5");
1003  const std::string xdmfFileName(domainName+".xdmf");
1004  const std::string geometryName("poiseuille");
1005  const std::string gridName("grid1");
1006  if(ix::sys::FILEEXISTS(hdfFileName)){
1007  ix::sys::Remove(hdfFileName);
1008  ix::sys::Remove(xdmfFileName);
1009  }
1010  fixtures::ConfigurationType testConfig;
1011  state_t paramState;
1012  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
1013  gridName,exactGrid,exactState,paramState,
1014  testConfig,compareStream)){
1015  std::cout << compareStream.str();
1016  std::cout << "Output to HDF5 failed:\n";
1017  std::cout << "Continuing....\n";
1018 
1019  }
1020  pcpp::io::RenewStream(compareStream);
1021  }
1022  ix::sys::ChDir("../");
1023  }
1024 
1025  if(myRank == 0) {
1026  if(!ix::sys::FILEEXISTS("dimen/exact_Poiseuille_dimen_3000.h5")){
1027  std::cout << testFunctionName << ": failed to produce expected" << std::endl
1028  << "output file (exact_Poiseuille_dimen_3000.h5). Test failed." << std::endl;
1029  PoiseuilleTestResult = false;
1030  }
1031 
1032  }
1033  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1034  if(!PoiseuilleTestResult)
1035  goto endTest;
1036 
1037  {
1038  if(myRank == 0) {
1039  std::cout << "Comparing exact solution against dimensional PlasCom2" << std::endl;
1040  ix::sys::ChDir("dimen");
1041  returnValue = plascom2::util::PC2Compare("PlasCom2_000003000.h5",
1042  "exact_Poiseuille_dimen_3000.h5",
1043  errorToleranceExact,compareStream);
1044 
1045  if(returnValue > 0){
1046  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
1047  << std::endl
1048  << compareStream.str()
1049  << std::endl;
1050  PoiseuilleTestResult = false;
1051  } else {
1052  std::cout << "Passed!" << std::endl;
1053  }
1054 
1055  pcpp::io::RenewStream(compareStream);
1056  ix::sys::ChDir("../");
1057  }
1058  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1059  }
1060 
1061  endTest:
1062  {
1063  ix::sys::ChDir(originalDirectory);
1064  parallelUnitResults.UpdateResult("Integrated:Viscid:PFRectilinear2DY",PoiseuilleTestResult);
1065  }
1066 
1067  return;
1068 }
1069 
1071  pcpp::CommunicatorType &testComm)
1072 {
1073 
1074  bool PoiseuilleTestResult = true;
1075  int myRank = testComm.Rank();
1076  int numProc = testComm.NProc();
1077 
1078  const std::string testFunctionName("TestIntegrated_PFCurvilinear2DY");
1079  const std::string testSuiteName("ViscidTest");
1080  const std::string testCaseName("Poiseuille2DY");
1081  const std::string testDirectory(testSuiteName+"/"+testCaseName);
1082  const std::string originalDirectory(ix::sys::CWD());
1083  std::ostringstream compareStream;
1084  double errorToleranceExact = 2.e-5;
1085  if(myRank == 0)
1086  std::cout << "Running " << testFunctionName << std::endl;
1087 
1088  int argc = 4;
1089  const char *argv[] = {"plascom2x","-c","curvilinear.config",NULL};
1090 
1091  pcpp::CommunicatorType dimenComm(testComm);
1092  plascom2::application simulationApplicationDimen(argc,const_cast<char **>(argv),dimenComm.Comm());
1093 
1094  int returnValue = 0;
1095 
1096  if(myRank == 0) {
1097  if(!ix::sys::FILEEXISTS(testDirectory)){
1098  PoiseuilleTestResult = false;
1099  std::cout << testFunctionName << ": testing directory ("
1100  << testDirectory << ") did not exist. Aborting test."
1101  << std::endl;
1102  }
1103  }
1104  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1105  if(!PoiseuilleTestResult)
1106  goto endTest;
1107  ix::sys::ChDir(testDirectory);
1108 
1109  // simulation in dimensional coordinates
1110  {
1111  ix::sys::ChDir("dimen");
1112 
1113  if(myRank == 0)
1114  std::cout << "Running dimensional PlasCom2..." << std::endl;
1115  returnValue = application::ApplicationDriver(simulationApplicationDimen);
1116  if(myRank == 0)
1117  std::cout << "Done running dimensional PlasCom2." << std::endl;
1118 
1119  if(returnValue > 0){
1120  std::cout << testFunctionName << ": PlasCom2 returned error code ("
1121  << returnValue << "). Test failed." << std::endl;
1122  PoiseuilleTestResult = false;
1123  goto endTest;
1124  }
1125 
1126  if(myRank == 0){
1127  if(!ix::sys::FILEEXISTS("PlasCom2_000003000.h5")){
1128  std::cout << testFunctionName << ": PlasCom2 failed to produce expected"
1129  << "output file (PlasCom2_000003000.h5). Test failed." << std::endl;
1130  PoiseuilleTestResult = false;
1131  }
1132  }
1133  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1134  if(!PoiseuilleTestResult)
1135  goto endTest;
1136  }
1137 
1138  // get the domain from the application and copy the state
1139  // make a new scope for variable declaration, needed for goto usage :-(
1140  {
1141  DomainVector &simulationDomain(simulationApplicationDimen.GetDomains());
1142 
1143  // generate a state containing the exact solution and write it to a file
1144  int direction=1;
1145  // single grid
1146  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
1147  state_t exactState(simulationDomain[0].State(0));
1148  if(testfixtures::viscid::GeneratePoiseuilleExact(exactGrid,exactState,direction)){
1149  PoiseuilleTestResult = false;
1150  std::cout << "Fail to calculate exact poiseuille solution" << std::endl;
1151  }
1152  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1153  if(!PoiseuilleTestResult)
1154  goto endTest;
1155 
1156  size_t numPointsBuffer = exactGrid.BufferSize();
1157 
1158  int rhoHandle = exactState.GetDataIndex("rho");
1159  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
1160  double *rhoPtr = rhoData.Data<double>();
1161 
1162  int rhoVHandle = exactState.GetDataIndex("rhoV");
1163  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
1164  double *rhoVPtr = rhoVData.Data<double>();
1165 
1166  int rhoEHandle = exactState.GetDataIndex("rhoE");
1167  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
1168  double *rhoEPtr = rhoEData.Data<double>();
1169 
1170  int velocityHandle = exactState.GetDataIndex("velocity");
1171  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
1172  double *velocityPtr = velocityData.Data<double>();
1173 
1174  int pressureHandle = exactState.GetDataIndex("pressure");
1175  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
1176  double *pressurePtr = pressureData.Data<double>();
1177 
1178  int temperatureHandle = exactState.GetDataIndex("temperature");
1179  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
1180  double *temperaturePtr = temperatureData.Data<double>();
1181 
1182  {
1183  // output exact dimensional solution
1184  if(myRank == 0)
1185  std::cout << "Writing exact solution for dimensional PlasCom2" << std::endl;
1186  std::ostringstream Ostr;
1187  Ostr << "exact_Poiseuille_dimen_3000";
1188  const std::string domainName(Ostr.str());
1189  const std::string hdfFileName(domainName+".h5");
1190  const std::string xdmfFileName(domainName+".xdmf");
1191  const std::string geometryName("poiseuille");
1192  const std::string gridName("grid1");
1193  if(ix::sys::FILEEXISTS(hdfFileName)){
1194  ix::sys::Remove(hdfFileName);
1195  ix::sys::Remove(xdmfFileName);
1196  }
1197  fixtures::ConfigurationType testConfig;
1198  state_t paramState;
1199  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
1200  gridName,exactGrid,exactState,paramState,
1201  testConfig,compareStream)){
1202  std::cout << compareStream.str();
1203  std::cout << "Output to HDF5 failed:\n";
1204  std::cout << "Continuing....\n";
1205 
1206  }
1207  pcpp::io::RenewStream(compareStream);
1208  }
1209  ix::sys::ChDir("../");
1210  }
1211 
1212  if(myRank == 0) {
1213  if(!ix::sys::FILEEXISTS("dimen/exact_Poiseuille_dimen_3000.h5")){
1214  std::cout << testFunctionName << ": failed to produce expected" << std::endl
1215  << "output file (exact_Poiseuille_dimen_3000.h5). Test failed." << std::endl;
1216  PoiseuilleTestResult = false;
1217  }
1218 
1219  }
1220  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1221  if(!PoiseuilleTestResult)
1222  goto endTest;
1223 
1224  {
1225  if(myRank == 0) {
1226  std::cout << "Comparing exact solution against dimensional PlasCom2" << std::endl;
1227  ix::sys::ChDir("dimen");
1228  returnValue = plascom2::util::PC2Compare("PlasCom2_000003000.h5",
1229  "exact_Poiseuille_dimen_3000.h5",
1230  errorToleranceExact,compareStream);
1231 
1232  if(returnValue > 0){
1233  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
1234  << std::endl
1235  << compareStream.str()
1236  << std::endl;
1237  PoiseuilleTestResult = false;
1238  } else {
1239  std::cout << "Passed!" << std::endl;
1240  }
1241 
1242  pcpp::io::RenewStream(compareStream);
1243  ix::sys::ChDir("../");
1244  }
1245  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1246  }
1247 
1248  endTest:
1249  {
1250  ix::sys::ChDir(originalDirectory);
1251  parallelUnitResults.UpdateResult("Integrated:Viscid:PFCurvilinear2DY",PoiseuilleTestResult);
1252  }
1253 
1254  return;
1255 }
1256 
1258  pcpp::CommunicatorType &testComm)
1259 {
1260 
1261  bool PoiseuilleTestResult = true;
1262  int myRank = testComm.Rank();
1263  int numProc = testComm.NProc();
1264 
1265  const std::string testFunctionName("TestIntegrated_Poiseuille3DZ");
1266  const std::string testSuiteName("ViscidTest");
1267  const std::string testCaseName("Poiseuille3DZ");
1268  const std::string testDirectory(testSuiteName+"/"+testCaseName);
1269  const std::string originalDirectory(ix::sys::CWD());
1270  std::ostringstream compareStream;
1271  double errorToleranceExact = 2.e-4;
1272  if(myRank == 0)
1273  std::cout << "Running " << testFunctionName << std::endl;
1274 
1275  int argc = 4;
1276  const char *argv[] = {"plascom2x","-c","poiseuille3d.config",NULL};
1277 
1278  pcpp::CommunicatorType nonDimenComm(testComm);
1279  plascom2::application simulationApplicationNonDimen(argc,const_cast<char **>(argv),nonDimenComm.Comm());
1280 
1281  int returnValue = 0;
1282 
1283  if(myRank == 0) {
1284  if(!ix::sys::FILEEXISTS(testDirectory)){
1285  PoiseuilleTestResult = false;
1286  std::cout << "TestIntegrated_poiseuille3d: testing directory ("
1287  << testDirectory << ") did not exist. Aborting test."
1288  << std::endl;
1289  }
1290  }
1291  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1292  if(!PoiseuilleTestResult)
1293  goto endTest;
1294  ix::sys::ChDir(testDirectory);
1295 
1296  // nonDimensional run using PlasComCM (rho, C, T) non-dimensionalization
1297  {
1298  ix::sys::ChDir("nonDimen");
1299 
1300  if(myRank == 0)
1301  std::cout << "Running nonDimensional PlasCom2..." << std::endl;
1302  returnValue = application::ApplicationDriver(simulationApplicationNonDimen);
1303  if(myRank == 0)
1304  std::cout << "Done running nonDimensional PlasCom2." << std::endl;
1305 
1306  if(returnValue > 0){
1307  std::cout << "TestIntegrated_Poiseuille3DZ: PlasCom2 returned error code ("
1308  << returnValue << "). Test failed." << std::endl;
1309  PoiseuilleTestResult = false;
1310  goto endTest;
1311  }
1312 
1313  if(myRank == 0){
1314  if(!ix::sys::FILEEXISTS("PlasCom2_000003000.h5")){
1315  std::cout << "TestIntegrated_Poiseuille3DZ: PlasCom2 failed to produce expected" << std::endl
1316  << "output file (PlasCom2_000003000.h5). Test failed." << std::endl;
1317  PoiseuilleTestResult = false;
1318  }
1319  }
1320  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1321  if(!PoiseuilleTestResult)
1322  goto endTest;
1323  }
1324 
1325  // get the domain from the application and copy the state
1326  // make a new scope for variable declaration, needed for goto usage :-(
1327  {
1328  DomainVector &simulationDomain(simulationApplicationNonDimen.GetDomains());
1329 
1330  // generate a state containing the exact solution and write it to a file
1331  int direction=2;
1332  // single grid
1333  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
1334  state_t exactState(simulationDomain[0].State(0));
1335  if(testfixtures::viscid::GeneratePoiseuilleExact(exactGrid,exactState,direction)){
1336  PoiseuilleTestResult = false;
1337  std::cout << "Fail to calculated exact poiseuille solution" << std::endl;
1338  }
1339  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1340  if(!PoiseuilleTestResult)
1341  goto endTest;
1342 
1343  size_t numPointsBuffer = exactGrid.BufferSize();
1344 
1345  int rhoHandle = exactState.GetDataIndex("rho");
1346  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
1347  double *rhoPtr = rhoData.Data<double>();
1348 
1349  int rhoVHandle = exactState.GetDataIndex("rhoV");
1350  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
1351  double *rhoVPtr = rhoVData.Data<double>();
1352 
1353  int rhoEHandle = exactState.GetDataIndex("rhoE");
1354  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
1355  double *rhoEPtr = rhoEData.Data<double>();
1356 
1357  int velocityHandle = exactState.GetDataIndex("velocity");
1358  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
1359  double *velocityPtr = velocityData.Data<double>();
1360 
1361  int pressureHandle = exactState.GetDataIndex("pressure");
1362  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
1363  double *pressurePtr = pressureData.Data<double>();
1364 
1365  int temperatureHandle = exactState.GetDataIndex("temperature");
1366  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
1367  double *temperaturePtr = temperatureData.Data<double>();
1368 
1369  {
1370  if(myRank == 0) {
1371  std::cout << "Writing exact solution for nonDimensional PlasCom2" << std::endl;
1372  }
1373  //
1374  // output exact nonDimensional solution in the PlasComCM dimensionalization
1375  //
1376  // we are using a non-unity length reference, which is embedded in the mesh
1377  // modify the exact solution for velocity to account for this
1378  //
1379  double gamma = 1.4;
1380  double refRho = 1.225;
1381  double refPress = 101325;
1382  double refSndSpd = sqrt(gamma*refPress/refRho);
1383  double refVel = sqrt(refPress/refRho);
1384  double refTemp = 288.15;
1385  double refLength = 1.0;
1386  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
1387  // subtract out kinetic energy
1388  rhoEPtr[iPoint] -= 0.5*rhoVPtr[iPoint]*rhoVPtr[iPoint]/rhoPtr[iPoint];
1389 
1390  // modify velocity and momentum by reference length
1391  rhoPtr[iPoint] /= refRho;
1392  rhoVPtr[iPoint+2*numPointsBuffer] /= refRho*refVel/refLength;
1393 
1394  // scale internal energy and add kinetic energy back
1395  rhoEPtr[iPoint] /= refVel*refVel*refRho;
1396  rhoEPtr[iPoint] += 0.5*rhoVPtr[iPoint]*rhoVPtr[iPoint]/rhoPtr[iPoint];
1397 
1398  velocityPtr[iPoint+2*numPointsBuffer] /= refVel;
1399  pressurePtr[iPoint] /= refPress;
1400  temperaturePtr[iPoint] /= refTemp;
1401  }
1402 
1403  std::ostringstream Ostr;
1404  Ostr << "exact_Poiseuille_ND_3000";
1405  const std::string domainName(Ostr.str());
1406  const std::string hdfFileName(domainName+".h5");
1407  const std::string xdmfFileName(domainName+".xdmf");
1408  const std::string geometryName("poiseuille3d");
1409  const std::string gridName("grid1");
1410  if(ix::sys::FILEEXISTS(hdfFileName)){
1411  ix::sys::Remove(hdfFileName);
1412  ix::sys::Remove(xdmfFileName);
1413  }
1414  fixtures::ConfigurationType testConfig;
1415  state_t paramState;
1416  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
1417  gridName,exactGrid,exactState,paramState,
1418  testConfig,compareStream)){
1419  std::cout << compareStream.str();
1420  std::cout << "Output to HDF5 failed:\n";
1421  std::cout << "Continuing....\n";
1422 
1423  }
1424  pcpp::io::RenewStream(compareStream);
1425  }
1426  ix::sys::ChDir("../");
1427  }
1428 
1429  if(myRank == 0) {
1430 
1431  if(!ix::sys::FILEEXISTS("nonDimen/exact_Poiseuille_ND_3000.h5")){
1432  std::cout << "TestIntegrated_Poiseuille3DZ: failed to produce expected" << std::endl
1433  << "output file (exact_Poiseuille_ND_3000.h5). Test failed." << std::endl;
1434  PoiseuilleTestResult = false;
1435  }
1436  }
1437  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1438  if(!PoiseuilleTestResult)
1439  goto endTest;
1440 
1441  {
1442  if(myRank == 0) {
1443  std::cout << "Comparing exact solution against nonDimensional PlasCom2" << std::endl;
1444  ix::sys::ChDir("nonDimen");
1445  returnValue = plascom2::util::PC2Compare("PlasCom2_000003000.h5",
1446  "exact_Poiseuille_ND_3000.h5",
1447  errorToleranceExact,compareStream);
1448  if(returnValue > 0){
1449  std::cout << "TestIntegrated_Poiseuille3DZ: Resulting state comparison yielded negative result:"
1450  << std::endl
1451  << compareStream.str()
1452  << std::endl;
1453  PoiseuilleTestResult = false;
1454  } else {
1455  std::cout << "Passed!" << std::endl;
1456  }
1457  pcpp::io::RenewStream(compareStream);
1458  ix::sys::ChDir("../");
1459  }
1460  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1461  }
1462 
1463  endTest:
1464  {
1465  ix::sys::ChDir(originalDirectory);
1466  parallelUnitResults.UpdateResult("Integrated:Viscid:Poiseuille3DZ",PoiseuilleTestResult);
1467  }
1468 
1469  return;
1470 }
1471 
1473  pcpp::CommunicatorType &testComm)
1474 {
1475 
1476  bool PoiseuilleTestResult = true;
1477  int myRank = testComm.Rank();
1478  int numProc = testComm.NProc();
1479 
1480  const std::string testFunctionName("TestIntegrated_PFRectilinear3DZ");
1481  const std::string testSuiteName("ViscidTest");
1482  const std::string testCaseName("Poiseuille3DZ");
1483  const std::string testDirectory(testSuiteName+"/"+testCaseName);
1484  const std::string originalDirectory(ix::sys::CWD());
1485  std::ostringstream compareStream;
1486  double errorToleranceExact = 2.e-4;
1487  if(myRank == 0)
1488  std::cout << "Running " << testFunctionName << std::endl;
1489 
1490  int argc = 4;
1491  const char *argv[] = {"plascom2x","-c","rectilinear.config",NULL};
1492 
1493  pcpp::CommunicatorType nonDimenComm(testComm);
1494  plascom2::application simulationApplicationNonDimen(argc,const_cast<char **>(argv),nonDimenComm.Comm());
1495 
1496  int returnValue = 0;
1497 
1498  if(myRank == 0) {
1499  if(!ix::sys::FILEEXISTS(testDirectory)){
1500  PoiseuilleTestResult = false;
1501  std::cout << "TestIntegrated_poiseuille3d: testing directory ("
1502  << testDirectory << ") did not exist. Aborting test."
1503  << std::endl;
1504  }
1505  }
1506  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1507  if(!PoiseuilleTestResult)
1508  goto endTest;
1509  ix::sys::ChDir(testDirectory);
1510 
1511  // nonDimensional run using PlasComCM (rho, C, T) non-dimensionalization
1512  {
1513  ix::sys::ChDir("nonDimen");
1514 
1515  if(myRank == 0)
1516  std::cout << "Running nonDimensional PlasCom2..." << std::endl;
1517  returnValue = application::ApplicationDriver(simulationApplicationNonDimen);
1518  if(myRank == 0)
1519  std::cout << "Done running nonDimensional PlasCom2." << std::endl;
1520 
1521  if(returnValue > 0){
1522  std::cout << testFunctionName << ": PlasCom2 returned error code ("
1523  << returnValue << "). Test failed." << std::endl;
1524  PoiseuilleTestResult = false;
1525  goto endTest;
1526  }
1527 
1528  if(myRank == 0){
1529  if(!ix::sys::FILEEXISTS("PlasCom2_000003000.h5")){
1530  std::cout << testFunctionName << ": PlasCom2 failed to produce expected" << std::endl
1531  << "output file (PlasCom2_000003000.h5). Test failed." << std::endl;
1532  PoiseuilleTestResult = false;
1533  }
1534  }
1535  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1536  if(!PoiseuilleTestResult)
1537  goto endTest;
1538  }
1539 
1540  // get the domain from the application and copy the state
1541  // make a new scope for variable declaration, needed for goto usage :-(
1542  {
1543  DomainVector &simulationDomain(simulationApplicationNonDimen.GetDomains());
1544 
1545  // generate a state containing the exact solution and write it to a file
1546  int direction=2;
1547  // single grid
1548  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
1549  state_t exactState(simulationDomain[0].State(0));
1550  if(testfixtures::viscid::GeneratePoiseuilleExact(exactGrid,exactState,direction)){
1551  PoiseuilleTestResult = false;
1552  std::cout << "Fail to calculated exact poiseuille solution" << std::endl;
1553  }
1554  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1555  if(!PoiseuilleTestResult)
1556  goto endTest;
1557 
1558  size_t numPointsBuffer = exactGrid.BufferSize();
1559 
1560  int rhoHandle = exactState.GetDataIndex("rho");
1561  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
1562  double *rhoPtr = rhoData.Data<double>();
1563 
1564  int rhoVHandle = exactState.GetDataIndex("rhoV");
1565  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
1566  double *rhoVPtr = rhoVData.Data<double>();
1567 
1568  int rhoEHandle = exactState.GetDataIndex("rhoE");
1569  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
1570  double *rhoEPtr = rhoEData.Data<double>();
1571 
1572  int velocityHandle = exactState.GetDataIndex("velocity");
1573  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
1574  double *velocityPtr = velocityData.Data<double>();
1575 
1576  int pressureHandle = exactState.GetDataIndex("pressure");
1577  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
1578  double *pressurePtr = pressureData.Data<double>();
1579 
1580  int temperatureHandle = exactState.GetDataIndex("temperature");
1581  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
1582  double *temperaturePtr = temperatureData.Data<double>();
1583 
1584  {
1585  if(myRank == 0) {
1586  std::cout << "Writing exact solution for nonDimensional PlasCom2" << std::endl;
1587  }
1588  //
1589  // output exact nonDimensional solution in the PlasComCM dimensionalization
1590  //
1591  // we are using a non-unity length reference, which is embedded in the mesh
1592  // modify the exact solution for velocity to account for this
1593  //
1594  double gamma = 1.4;
1595  double refRho = 1.225;
1596  double refPress = 101325;
1597  double refSndSpd = sqrt(gamma*refPress/refRho);
1598  double refVel = sqrt(refPress/refRho);
1599  double refTemp = 288.15;
1600  double refLength = 1.0;
1601  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
1602  // subtract out kinetic energy
1603  rhoEPtr[iPoint] -= 0.5*rhoVPtr[iPoint]*rhoVPtr[iPoint]/rhoPtr[iPoint];
1604 
1605  // modify velocity and momentum by reference length
1606  rhoPtr[iPoint] /= refRho;
1607  rhoVPtr[iPoint+2*numPointsBuffer] /= refRho*refVel/refLength;
1608 
1609  // scale internal energy and add kinetic energy back
1610  rhoEPtr[iPoint] /= refVel*refVel*refRho;
1611  rhoEPtr[iPoint] += 0.5*rhoVPtr[iPoint]*rhoVPtr[iPoint]/rhoPtr[iPoint];
1612 
1613  velocityPtr[iPoint+2*numPointsBuffer] /= refVel;
1614  pressurePtr[iPoint] /= refPress;
1615  temperaturePtr[iPoint] /= refTemp;
1616  }
1617 
1618  std::ostringstream Ostr;
1619  Ostr << "exact_Poiseuille_ND_3000";
1620  const std::string domainName(Ostr.str());
1621  const std::string hdfFileName(domainName+".h5");
1622  const std::string xdmfFileName(domainName+".xdmf");
1623  const std::string geometryName("poiseuille3d");
1624  const std::string gridName("grid1");
1625  if(ix::sys::FILEEXISTS(hdfFileName)){
1626  ix::sys::Remove(hdfFileName);
1627  ix::sys::Remove(xdmfFileName);
1628  }
1629  fixtures::ConfigurationType testConfig;
1630  state_t paramState;
1631  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
1632  gridName,exactGrid,exactState,paramState,
1633  testConfig,compareStream)){
1634  std::cout << compareStream.str();
1635  std::cout << "Output to HDF5 failed:\n";
1636  std::cout << "Continuing....\n";
1637 
1638  }
1639  pcpp::io::RenewStream(compareStream);
1640  }
1641  ix::sys::ChDir("../");
1642  }
1643 
1644  if(myRank == 0) {
1645 
1646  if(!ix::sys::FILEEXISTS("nonDimen/exact_Poiseuille_ND_3000.h5")){
1647  std::cout << "TestIntegrated_Poiseuille3DZ: failed to produce expected" << std::endl
1648  << "output file (exact_Poiseuille_ND_3000.h5). Test failed." << std::endl;
1649  PoiseuilleTestResult = false;
1650  }
1651  }
1652  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1653  if(!PoiseuilleTestResult)
1654  goto endTest;
1655 
1656  {
1657  if(myRank == 0) {
1658  std::cout << "Comparing exact solution against nonDimensional PlasCom2" << std::endl;
1659  ix::sys::ChDir("nonDimen");
1660  returnValue = plascom2::util::PC2Compare("PlasCom2_000003000.h5",
1661  "exact_Poiseuille_ND_3000.h5",
1662  errorToleranceExact,compareStream);
1663  if(returnValue > 0){
1664  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
1665  << std::endl
1666  << compareStream.str()
1667  << std::endl;
1668  PoiseuilleTestResult = false;
1669  } else {
1670  std::cout << "Passed!" << std::endl;
1671  }
1672  pcpp::io::RenewStream(compareStream);
1673  ix::sys::ChDir("../");
1674  }
1675  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1676  }
1677 
1678  endTest:
1679  {
1680  ix::sys::ChDir(originalDirectory);
1681  parallelUnitResults.UpdateResult("Integrated:Viscid:PFRectilinear3DZ",PoiseuilleTestResult);
1682  }
1683 
1684  return;
1685 }
1686 
1688  pcpp::CommunicatorType &testComm)
1689 {
1690 
1691  bool PoiseuilleTestResult = true;
1692  int myRank = testComm.Rank();
1693  int numProc = testComm.NProc();
1694 
1695  const std::string testFunctionName("TestIntegrated_PFCurvilinear3DZ");
1696  const std::string testSuiteName("ViscidTest");
1697  const std::string testCaseName("Poiseuille3DZ");
1698  const std::string testDirectory(testSuiteName+"/"+testCaseName);
1699  const std::string originalDirectory(ix::sys::CWD());
1700  std::ostringstream compareStream;
1701  double errorToleranceExact = 2.e-4;
1702  if(myRank == 0)
1703  std::cout << "Running " << testFunctionName << std::endl;
1704 
1705  int argc = 4;
1706  const char *argv[] = {"plascom2x","-c","curvilinear.config",NULL};
1707 
1708  pcpp::CommunicatorType nonDimenComm(testComm);
1709  plascom2::application simulationApplicationNonDimen(argc,const_cast<char **>(argv),nonDimenComm.Comm());
1710 
1711  int returnValue = 0;
1712 
1713  if(myRank == 0) {
1714  if(!ix::sys::FILEEXISTS(testDirectory)){
1715  PoiseuilleTestResult = false;
1716  std::cout << "TestIntegrated_poiseuille3d: testing directory ("
1717  << testDirectory << ") did not exist. Aborting test."
1718  << std::endl;
1719  }
1720  }
1721  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1722  if(!PoiseuilleTestResult)
1723  goto endTest;
1724  ix::sys::ChDir(testDirectory);
1725 
1726  // nonDimensional run using PlasComCM (rho, C, T) non-dimensionalization
1727  {
1728  ix::sys::ChDir("nonDimen");
1729 
1730  if(myRank == 0)
1731  std::cout << "Running nonDimensional PlasCom2..." << std::endl;
1732  returnValue = application::ApplicationDriver(simulationApplicationNonDimen);
1733  if(myRank == 0)
1734  std::cout << "Done running nonDimensional PlasCom2." << std::endl;
1735 
1736  if(returnValue > 0){
1737  std::cout << testFunctionName << ": PlasCom2 returned error code ("
1738  << returnValue << "). Test failed." << std::endl;
1739  PoiseuilleTestResult = false;
1740  goto endTest;
1741  }
1742 
1743  if(myRank == 0){
1744  if(!ix::sys::FILEEXISTS("PlasCom2_000003000.h5")){
1745  std::cout << testFunctionName << ": PlasCom2 failed to produce expected" << std::endl
1746  << "output file (PlasCom2_000003000.h5). Test failed." << std::endl;
1747  PoiseuilleTestResult = false;
1748  }
1749  }
1750  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1751  if(!PoiseuilleTestResult)
1752  goto endTest;
1753  }
1754 
1755  // get the domain from the application and copy the state
1756  // make a new scope for variable declaration, needed for goto usage :-(
1757  {
1758  DomainVector &simulationDomain(simulationApplicationNonDimen.GetDomains());
1759 
1760  // generate a state containing the exact solution and write it to a file
1761  int direction=2;
1762  // single grid
1763  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
1764  state_t exactState(simulationDomain[0].State(0));
1765  if(testfixtures::viscid::GeneratePoiseuilleExact(exactGrid,exactState,direction)){
1766  PoiseuilleTestResult = false;
1767  std::cout << "Fail to calculated exact poiseuille solution" << std::endl;
1768  }
1769  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1770  if(!PoiseuilleTestResult)
1771  goto endTest;
1772 
1773  size_t numPointsBuffer = exactGrid.BufferSize();
1774 
1775  int rhoHandle = exactState.GetDataIndex("rho");
1776  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
1777  double *rhoPtr = rhoData.Data<double>();
1778 
1779  int rhoVHandle = exactState.GetDataIndex("rhoV");
1780  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
1781  double *rhoVPtr = rhoVData.Data<double>();
1782 
1783  int rhoEHandle = exactState.GetDataIndex("rhoE");
1784  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
1785  double *rhoEPtr = rhoEData.Data<double>();
1786 
1787  int velocityHandle = exactState.GetDataIndex("velocity");
1788  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
1789  double *velocityPtr = velocityData.Data<double>();
1790 
1791  int pressureHandle = exactState.GetDataIndex("pressure");
1792  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
1793  double *pressurePtr = pressureData.Data<double>();
1794 
1795  int temperatureHandle = exactState.GetDataIndex("temperature");
1796  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
1797  double *temperaturePtr = temperatureData.Data<double>();
1798 
1799  {
1800  if(myRank == 0) {
1801  std::cout << "Writing exact solution for nonDimensional PlasCom2" << std::endl;
1802  }
1803  //
1804  // output exact nonDimensional solution in the PlasComCM dimensionalization
1805  //
1806  // we are using a non-unity length reference, which is embedded in the mesh
1807  // modify the exact solution for velocity to account for this
1808  //
1809  double gamma = 1.4;
1810  double refRho = 1.225;
1811  double refPress = 101325;
1812  double refSndSpd = sqrt(gamma*refPress/refRho);
1813  double refVel = sqrt(refPress/refRho);
1814  double refTemp = 288.15;
1815  double refLength = 1.0;
1816  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
1817  // subtract out kinetic energy
1818  rhoEPtr[iPoint] -= 0.5*rhoVPtr[iPoint]*rhoVPtr[iPoint]/rhoPtr[iPoint];
1819 
1820  // modify velocity and momentum by reference length
1821  rhoPtr[iPoint] /= refRho;
1822  rhoVPtr[iPoint+2*numPointsBuffer] /= refRho*refVel/refLength;
1823 
1824  // scale internal energy and add kinetic energy back
1825  rhoEPtr[iPoint] /= refVel*refVel*refRho;
1826  rhoEPtr[iPoint] += 0.5*rhoVPtr[iPoint]*rhoVPtr[iPoint]/rhoPtr[iPoint];
1827 
1828  velocityPtr[iPoint+2*numPointsBuffer] /= refVel;
1829  pressurePtr[iPoint] /= refPress;
1830  temperaturePtr[iPoint] /= refTemp;
1831  }
1832 
1833  std::ostringstream Ostr;
1834  Ostr << "exact_Poiseuille_ND_3000";
1835  const std::string domainName(Ostr.str());
1836  const std::string hdfFileName(domainName+".h5");
1837  const std::string xdmfFileName(domainName+".xdmf");
1838  const std::string geometryName("poiseuille3d");
1839  const std::string gridName("grid1");
1840  if(ix::sys::FILEEXISTS(hdfFileName)){
1841  ix::sys::Remove(hdfFileName);
1842  ix::sys::Remove(xdmfFileName);
1843  }
1844  fixtures::ConfigurationType testConfig;
1845  state_t paramState;
1846  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
1847  gridName,exactGrid,exactState,paramState,
1848  testConfig,compareStream)){
1849  std::cout << compareStream.str();
1850  std::cout << "Output to HDF5 failed:\n";
1851  std::cout << "Continuing....\n";
1852 
1853  }
1854  pcpp::io::RenewStream(compareStream);
1855  }
1856  ix::sys::ChDir("../");
1857  }
1858 
1859  if(myRank == 0) {
1860 
1861  if(!ix::sys::FILEEXISTS("nonDimen/exact_Poiseuille_ND_3000.h5")){
1862  std::cout << "TestIntegrated_Poiseuille3DZ: failed to produce expected" << std::endl
1863  << "output file (exact_Poiseuille_ND_3000.h5). Test failed." << std::endl;
1864  PoiseuilleTestResult = false;
1865  }
1866  }
1867  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1868  if(!PoiseuilleTestResult)
1869  goto endTest;
1870 
1871  {
1872  if(myRank == 0) {
1873  std::cout << "Comparing exact solution against nonDimensional PlasCom2" << std::endl;
1874  ix::sys::ChDir("nonDimen");
1875  returnValue = plascom2::util::PC2Compare("PlasCom2_000003000.h5",
1876  "exact_Poiseuille_ND_3000.h5",
1877  errorToleranceExact,compareStream);
1878  if(returnValue > 0){
1879  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
1880  << std::endl
1881  << compareStream.str()
1882  << std::endl;
1883  PoiseuilleTestResult = false;
1884  } else {
1885  std::cout << "Passed!" << std::endl;
1886  }
1887  pcpp::io::RenewStream(compareStream);
1888  ix::sys::ChDir("../");
1889  }
1890  PoiseuilleTestResult = CheckResult(PoiseuilleTestResult,testComm);
1891  }
1892 
1893  endTest:
1894  {
1895  ix::sys::ChDir(originalDirectory);
1896  parallelUnitResults.UpdateResult("Integrated:Viscid:PFCurvilinear3DZ",PoiseuilleTestResult);
1897  }
1898 
1899  return;
1900 }
1901 
1903  pcpp::CommunicatorType &testComm)
1904 {
1905 
1906  bool viscousShockTestResult = true;
1907  int myRank = testComm.Rank();
1908  int numProc = testComm.NProc();
1909 
1910  const std::string testFunctionName("TestIntegrated_ViscousShock2DX");
1911  const std::string testSuiteName("ViscidTest");
1912  const std::string testCaseName("ViscousShock2DX");
1913  const std::string testDirectory(testSuiteName+"/"+testCaseName);
1914  const std::string originalDirectory(ix::sys::CWD());
1915  std::ostringstream compareStream;
1916  double errorTolerancePC1_2 = 1.e-4;
1917  double errorToleranceExact = 5.e-2;
1918  if(myRank == 0)
1919  std::cout << "Running " << testFunctionName << std::endl;
1920 
1921  int argc = 4;
1922  const char *argv[] = {"plascom2x","-c","shock.test.config",NULL};
1923 
1924  pcpp::CommunicatorType dimenComm(testComm);
1925  pcpp::CommunicatorType nonDimenComm(testComm);
1926  pcpp::CommunicatorType nonDimenComm2(testComm);
1927 
1928  plascom2::application simulationApplicationDimen(argc,const_cast<char **>(argv),dimenComm.Comm());
1929  plascom2::application simulationApplicationNonDimen(argc,const_cast<char **>(argv),nonDimenComm.Comm());
1930  plascom2::application simulationApplicationNonDimen2(argc,const_cast<char **>(argv),nonDimenComm2.Comm());
1931 
1932  int returnValue = 0;
1933 
1934  if(myRank == 0) {
1935  if(!ix::sys::FILEEXISTS(testDirectory)){
1936  viscousShockTestResult = false;
1937  std::cout << "TestIntegrated_viscousShock2D: testing directory ("
1938  << testDirectory << ") did not exist. Aborting test."
1939  << std::endl;
1940  }
1941  }
1942  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
1943  if(!viscousShockTestResult)
1944  goto endTest;
1945  ix::sys::ChDir(testDirectory);
1946 
1947 
1948  // simulation in dimensional coordinates
1949  {
1950  ix::sys::ChDir("dimenX");
1951 
1952  if(myRank == 0)
1953  std::cout << "Running dimensional PlasCom2..." << std::endl;
1954  returnValue = application::ApplicationDriver(simulationApplicationDimen);
1955  if(myRank == 0)
1956  std::cout << "Done running dimensional PlasCom2." << std::endl;
1957 
1958  if(returnValue > 0){
1959  std::cout << "TestIntegrated_ViscousShock2DX: PlasCom2 returned error code ("
1960  << returnValue << "). Test failed." << std::endl;
1961  viscousShockTestResult = false;
1962  goto endTest;
1963  }
1964 
1965  if(myRank == 0){
1966  if(!ix::sys::FILEEXISTS("PlasCom2_000040000.h5")){
1967  std::cout << "TestIntegrated_ViscousShock2DX: PlasCom2 failed to produce expected"
1968  << "output file (PlasCom2_000040000.h5). Test failed." << std::endl;
1969  viscousShockTestResult = false;
1970  }
1971  }
1972  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
1973  if(!viscousShockTestResult)
1974  goto endTest;
1975  }
1976 
1977  // nonDimensional run using standard (rho, P, T) non-dimensionalization
1978  {
1979  ix::sys::ChDir("../nonDimenX");
1980 
1981  if(myRank == 0)
1982  std::cout << "Running nonDimensional PlasCom2..." << std::endl;
1983  returnValue = application::ApplicationDriver(simulationApplicationNonDimen);
1984  if(myRank == 0)
1985  std::cout << "Done running nonDimensional PlasCom2." << std::endl;
1986 
1987  if(returnValue > 0){
1988  std::cout << "TestIntegrated_ViscousShock2DX: PlasCom2 returned error code ("
1989  << returnValue << "). Test failed." << std::endl;
1990  viscousShockTestResult = false;
1991  goto endTest;
1992  }
1993 
1994  if(myRank == 0){
1995  if(!ix::sys::FILEEXISTS("PlasCom2_000040000.h5")){
1996  std::cout << "TestIntegrated_ViscousShock2DX: PlasCom2 failed to produce expected" << std::endl
1997  << "output file (PlasCom2_000002000.h5). Test failed." << std::endl;
1998  viscousShockTestResult = false;
1999  }
2000  }
2001  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2002  if(!viscousShockTestResult)
2003  goto endTest;
2004  }
2005 
2006  // nonDimensional run using PlasComCM (rho, C, T) non-dimensionalization
2007  {
2008  ix::sys::ChDir("../nonDimenPC1X");
2009 
2010  if(myRank == 0)
2011  std::cout << "Running nonDimensional PlasCom2... (PC1ND)" << std::endl;
2012  returnValue = application::ApplicationDriver(simulationApplicationNonDimen2);
2013  if(myRank == 0)
2014  std::cout << "Done running nonDimensional PlasCom2 (PC1ND)." << std::endl;
2015 
2016  if(returnValue > 0){
2017  std::cout << "TestIntegrated_ViscousShock2DX: PlasCom2 returned error code ("
2018  << returnValue << "). Test failed." << std::endl;
2019  viscousShockTestResult = false;
2020  goto endTest;
2021  }
2022 
2023  if(myRank == 0){
2024  if(!ix::sys::FILEEXISTS("PlasCom2_000040000.h5")){
2025  std::cout << "TestIntegrated_ViscousShock2DX: PlasCom2 failed to produce expected" << std::endl
2026  << "output file (PlasCom2_000002000.h5). Test failed." << std::endl;
2027  viscousShockTestResult = false;
2028  }
2029  }
2030  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2031  if(!viscousShockTestResult)
2032  goto endTest;
2033  }
2034 
2035  // get the domain from the application and copy the state
2036  // make a new scope for variable declaration, needed for goto usage :-(
2037  {
2038  DomainVector &simulationDomain(simulationApplicationDimen.GetDomains());
2039 
2040  // generate a state containing the exact solution and write it to a file
2041  int direction=0;
2042  // single grid
2043  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
2044  state_t exactState(simulationDomain[0].State(0));
2045  if(testfixtures::viscid::GenerateViscidShockExact(exactGrid,exactState,direction)){
2046  viscousShockTestResult = false;
2047  std::cout << "Fail to calculated exact viscid shock solution" << std::endl;
2048  }
2049  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2050  if(!viscousShockTestResult)
2051  goto endTest;
2052 
2053  size_t numPointsBuffer = exactGrid.BufferSize();
2054 
2055  int rhoHandle = exactState.GetDataIndex("rho");
2056  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
2057  double *rhoPtr = rhoData.Data<double>();
2058 
2059  int rhoVHandle = exactState.GetDataIndex("rhoV");
2060  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
2061  double *rhoVPtr = rhoVData.Data<double>();
2062 
2063  int rhoEHandle = exactState.GetDataIndex("rhoE");
2064  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
2065  double *rhoEPtr = rhoEData.Data<double>();
2066 
2067  int velocityHandle = exactState.GetDataIndex("velocity");
2068  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
2069  double *velocityPtr = velocityData.Data<double>();
2070 
2071  int pressureHandle = exactState.GetDataIndex("pressure");
2072  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
2073  double *pressurePtr = pressureData.Data<double>();
2074 
2075  int temperatureHandle = exactState.GetDataIndex("temperature");
2076  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
2077  double *temperaturePtr = temperatureData.Data<double>();
2078 
2079  {
2080  // output exact dimensional solution
2081  if(myRank == 0)
2082  std::cout << "Writing exact solution for dimensional PlasCom2" << std::endl;
2083  ix::sys::ChDir("../dimenX");
2084  std::ostringstream Ostr;
2085  Ostr << "exact_ViscousShock_dimen_40000";
2086  const std::string domainName(Ostr.str());
2087  const std::string hdfFileName(domainName+".h5");
2088  const std::string xdmfFileName(domainName+".xdmf");
2089  const std::string geometryName("shockTube");
2090  const std::string gridName("grid1");
2091  if(ix::sys::FILEEXISTS(hdfFileName)){
2092  ix::sys::Remove(hdfFileName);
2093  ix::sys::Remove(xdmfFileName);
2094  }
2095  fixtures::ConfigurationType testConfig;
2096  state_t paramState;
2097  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
2098  gridName,exactGrid,exactState,paramState,
2099  testConfig,compareStream)){
2100  std::cout << compareStream.str();
2101  std::cout << "Output to HDF5 failed:\n";
2102  std::cout << "Continuing....\n";
2103 
2104  }
2105  pcpp::io::RenewStream(compareStream);
2106  }
2107 
2108  {
2109  if(myRank == 0)
2110  std::cout << "Writing exact solution for nonDimensional PlasCom2"
2111  << std::endl;
2112  ix::sys::ChDir("../nonDimenX");
2113  // output exact nonDimensional solution in the PlasComCM dimensionalization
2114  // we set the reference length to 1 so the grid doesn't need to change
2115  double refRho = 1.225;
2116  double refPress = 101325;
2117  double refVel = sqrt(refPress/refRho);
2118  double refTemp = 288.15;
2119  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
2120 
2121  rhoPtr[iPoint] /= refRho;
2122  rhoVPtr[iPoint] /= refRho*refVel;
2123  rhoEPtr[iPoint] /= refVel*refVel*refRho;
2124 
2125  velocityPtr[iPoint] /= refVel;
2126  pressurePtr[iPoint] /= refRho*refVel*refVel;
2127  temperaturePtr[iPoint] /= refTemp;
2128  }
2129 
2130  std::ostringstream Ostr;
2131  Ostr << "exact_ViscousShock_ND_40000";
2132  const std::string domainName(Ostr.str());
2133  const std::string hdfFileName(domainName+".h5");
2134  const std::string xdmfFileName(domainName+".xdmf");
2135  const std::string geometryName("shockTube");
2136  const std::string gridName("grid1");
2137  if(ix::sys::FILEEXISTS(hdfFileName)){
2138  ix::sys::Remove(hdfFileName);
2139  ix::sys::Remove(xdmfFileName);
2140  }
2141  fixtures::ConfigurationType testConfig;
2142  state_t paramState;
2143  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
2144  gridName,exactGrid,exactState,paramState,
2145  testConfig,compareStream)){
2146  std::cout << compareStream.str();
2147  std::cout << "Output to HDF5 failed:\n";
2148  std::cout << "Continuing....\n";
2149 
2150  }
2151  pcpp::io::RenewStream(compareStream);
2152  }
2153 
2154  {
2155  if(myRank == 0) {
2156  std::cout << "Writing exact solution for nonDimensional PlasCom2"
2157  << " with legacy PlasComCM nonDimensionalization" << std::endl;
2158  }
2159  ix::sys::ChDir("../nonDimenPC1X");
2160  // output exact nonDimensional solution in the PlasComCM dimensionalization
2161  // we set the reference length to 1 so the grid doesn't need to change
2162  //
2163  // since we already modified the solution to be in the default PlasCom2
2164  // nonDimensionalization, we only need to multiply by the converseion
2165  // factor between the two
2166  double gamma = 1.4;
2167  double refRho = 1.0;
2168  double refPress = 101325;
2169  double refSndSpd = sqrt(gamma);
2170  double refTemp = gamma-1;
2171  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
2172  rhoPtr[iPoint] /= refRho;
2173  rhoVPtr[iPoint] /= refRho*refSndSpd;
2174  rhoEPtr[iPoint] /= refSndSpd*refSndSpd*refRho;
2175 
2176  velocityPtr[iPoint] /= refSndSpd;
2177  pressurePtr[iPoint] /= refRho*refSndSpd*refSndSpd;
2178  temperaturePtr[iPoint] /= refTemp;
2179  }
2180 
2181  std::ostringstream Ostr;
2182  Ostr << "exact_ViscousShock_PC1ND_40000";
2183  const std::string domainName(Ostr.str());
2184  const std::string hdfFileName(domainName+".h5");
2185  const std::string xdmfFileName(domainName+".xdmf");
2186  const std::string geometryName("shockTube");
2187  const std::string gridName("grid1");
2188  if(ix::sys::FILEEXISTS(hdfFileName)){
2189  ix::sys::Remove(hdfFileName);
2190  ix::sys::Remove(xdmfFileName);
2191  }
2192  fixtures::ConfigurationType testConfig;
2193  state_t paramState;
2194  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
2195  gridName,exactGrid,exactState,paramState,
2196  testConfig,compareStream)){
2197  std::cout << compareStream.str();
2198  std::cout << "Output to HDF5 failed:\n";
2199  std::cout << "Continuing....\n";
2200 
2201  }
2202  pcpp::io::RenewStream(compareStream);
2203  }
2204  ix::sys::ChDir("../");
2205  }
2206 
2207  if(myRank == 0) {
2208  if(!ix::sys::FILEEXISTS("dimenX/exact_ViscousShock_dimen_40000.h5")){
2209  std::cout << "TestIntegrated_ViscousShock2DX: failed to produce expected" << std::endl
2210  << "output file (exact_ViscousShock_dimen_40000.h5). Test failed." << std::endl;
2211  viscousShockTestResult = false;
2212  }
2213 
2214  if(!ix::sys::FILEEXISTS("nonDimenX/exact_ViscousShock_ND_40000.h5")){
2215  std::cout << "TestIntegrated_ViscousShock2DX: failed to produce expected" << std::endl
2216  << "output file (exact_ViscousShock_ND_40000.h5). Test failed." << std::endl;
2217  viscousShockTestResult = false;
2218  }
2219 
2220  if(!ix::sys::FILEEXISTS("nonDimenPC1X/exact_ViscousShock_PC1ND_40000.h5")){
2221  std::cout << "TestIntegrated_ViscousShock2DX: failed to produce expected" << std::endl
2222  << "output file (exact_ViscousShock_PC1ND_40000.h5). Test failed." << std::endl;
2223  viscousShockTestResult = false;
2224  }
2225  }
2226  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2227  if(!viscousShockTestResult)
2228  goto endTest;
2229 
2230  {
2231  if(myRank == 0) {
2232  std::cout << "Comparing exact solution against dimensional PlasCom2" << std::endl;
2233  ix::sys::ChDir("dimenX");
2234  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
2235  "exact_ViscousShock_dimen_40000.h5",
2236  errorToleranceExact,compareStream);
2237 
2238  if(returnValue > 0){
2239  std::cout << "TestIntegrated_ViscousShock2DX: Resulting state comparison yielded negative result:"
2240  << std::endl
2241  << compareStream.str()
2242  << std::endl;
2243  viscousShockTestResult = false;
2244  } else {
2245  std::cout << "Passed!" << std::endl;
2246  }
2247 
2248  pcpp::io::RenewStream(compareStream);
2249  ix::sys::ChDir("../");
2250  }
2251  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2252  }
2253 
2254  {
2255  if(myRank == 0) {
2256  std::cout << "Comparing exact solution against nonDimensional PlasCom2" << std::endl;
2257  ix::sys::ChDir("nonDimenX");
2258  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
2259  "exact_ViscousShock_ND_40000.h5",
2260  errorToleranceExact,compareStream);
2261  if(returnValue > 0){
2262  std::cout << "TestIntegrated_ViscousShock2DX: Resulting state comparison yielded negative result:"
2263  << std::endl
2264  << compareStream.str()
2265  << std::endl;
2266  viscousShockTestResult = false;
2267  } else {
2268  std::cout << "Passed!" << std::endl;
2269  }
2270  pcpp::io::RenewStream(compareStream);
2271  ix::sys::ChDir("../");
2272  }
2273  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2274  }
2275 
2276  {
2277  if(myRank == 0) {
2278  std::cout << "Comparing exact solution against nonDimensional PlasCom2"
2279  << " with legacy PlasComCM nonDimensionalization" << std::endl;
2280  ix::sys::ChDir("nonDimenPC1X");
2281  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
2282  "exact_ViscousShock_PC1ND_40000.h5",
2283  errorToleranceExact,compareStream);
2284  if(returnValue > 0){
2285  std::cout << "TestIntegrated_ViscousShock2DX: Resulting state comparison yielded negative result:"
2286  << std::endl
2287  << compareStream.str()
2288  << std::endl;
2289  viscousShockTestResult = false;
2290  } else {
2291  std::cout << "Passed!" << std::endl;
2292  }
2293  pcpp::io::RenewStream(compareStream);
2294 
2295  std::cout << "Comparing PlasComCM results against nonDimensional PlasCom2"
2296  << " with legacy PlasComCM nonDimensionalization" << std::endl;
2297  ix::sys::ChDir("nonDimenPC1X");
2298  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
2299  "PlasComCM/RocFlo-CM.00040000.h5",
2300  errorTolerancePC1_2,compareStream);
2301  if(returnValue > 0){
2302  std::cout << "TestIntegrated_ViscousShock2DX: Resulting state comparison yielded negative result:"
2303  << std::endl
2304  << compareStream.str()
2305  << std::endl;
2306  viscousShockTestResult = false;
2307  } else {
2308  std::cout << "Passed!" << std::endl;
2309  }
2310  pcpp::io::RenewStream(compareStream);
2311  ix::sys::ChDir("../");
2312  }
2313  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2314  testComm.Barrier();
2315  }
2316 
2317  endTest:
2318  {
2319  ix::sys::ChDir(originalDirectory);
2320  parallelUnitResults.UpdateResult("Integrated:Viscid:ViscousShock2DX",viscousShockTestResult);
2321  }
2322 
2323  return;
2324 }
2325 
2327  pcpp::CommunicatorType &testComm)
2328 {
2329 
2330  bool viscousShockTestResult = true;
2331  int myRank = testComm.Rank();
2332  int numProc = testComm.NProc();
2333 
2334  const std::string testFunctionName("TestIntegrated_VSRectilinear2DX");
2335  const std::string testSuiteName("ViscidTest");
2336  const std::string testCaseName("ViscousShock2DX");
2337  const std::string testDirectory(testSuiteName+"/"+testCaseName);
2338  const std::string originalDirectory(ix::sys::CWD());
2339  std::ostringstream compareStream;
2340  double errorTolerancePC1_2 = 1.e-4;
2341  double errorToleranceExact = 5.e-2;
2342  if(myRank == 0)
2343  std::cout << "Running " << testFunctionName << std::endl;
2344 
2345  int argc = 4;
2346  const char *argv[] = {"plascom2x","-c","rectilinear.config",NULL};
2347 
2348  pcpp::CommunicatorType dimenComm(testComm);
2349  pcpp::CommunicatorType nonDimenComm(testComm);
2350  pcpp::CommunicatorType nonDimenComm2(testComm);
2351 
2352  plascom2::application simulationApplicationDimen(argc,const_cast<char **>(argv),dimenComm.Comm());
2353  plascom2::application simulationApplicationNonDimen(argc,const_cast<char **>(argv),nonDimenComm.Comm());
2354  plascom2::application simulationApplicationNonDimen2(argc,const_cast<char **>(argv),nonDimenComm2.Comm());
2355 
2356  int returnValue = 0;
2357 
2358  if(myRank == 0) {
2359  if(!ix::sys::FILEEXISTS(testDirectory)){
2360  viscousShockTestResult = false;
2361  std::cout << ": testing directory ("
2362  << testDirectory << ") did not exist. Aborting test."
2363  << std::endl;
2364  }
2365  }
2366  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2367  if(!viscousShockTestResult)
2368  goto endTest;
2369  ix::sys::ChDir(testDirectory);
2370 
2371 
2372  // simulation in dimensional coordinates
2373  {
2374  ix::sys::ChDir("dimenX");
2375 
2376  if(myRank == 0)
2377  std::cout << "Running dimensional PlasCom2..." << std::endl;
2378  returnValue = application::ApplicationDriver(simulationApplicationDimen);
2379  if(myRank == 0)
2380  std::cout << "Done running dimensional PlasCom2." << std::endl;
2381 
2382  if(returnValue > 0){
2383  std::cout << testFunctionName << ": PlasCom2 returned error code ("
2384  << returnValue << "). Test failed." << std::endl;
2385  viscousShockTestResult = false;
2386  goto endTest;
2387  }
2388 
2389  if(myRank == 0){
2390  if(!ix::sys::FILEEXISTS("PlasCom2_000040000.h5")){
2391  std::cout << testFunctionName << ": PlasCom2 failed to produce expected"
2392  << "output file (PlasCom2_000040000.h5). Test failed." << std::endl;
2393  viscousShockTestResult = false;
2394  }
2395  }
2396  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2397  if(!viscousShockTestResult)
2398  goto endTest;
2399  }
2400 
2401  // nonDimensional run using standard (rho, P, T) non-dimensionalization
2402  {
2403  ix::sys::ChDir("../nonDimenX");
2404 
2405  if(myRank == 0)
2406  std::cout << "Running nonDimensional PlasCom2..." << std::endl;
2407  returnValue = application::ApplicationDriver(simulationApplicationNonDimen);
2408  if(myRank == 0)
2409  std::cout << "Done running nonDimensional PlasCom2." << std::endl;
2410 
2411  if(returnValue > 0){
2412  std::cout << testFunctionName << ": PlasCom2 returned error code ("
2413  << returnValue << "). Test failed." << std::endl;
2414  viscousShockTestResult = false;
2415  goto endTest;
2416  }
2417 
2418  if(myRank == 0){
2419  if(!ix::sys::FILEEXISTS("PlasCom2_000040000.h5")){
2420  std::cout << testFunctionName << ": PlasCom2 failed to produce expected" << std::endl
2421  << "output file (PlasCom2_000002000.h5). Test failed." << std::endl;
2422  viscousShockTestResult = false;
2423  }
2424  }
2425  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2426  if(!viscousShockTestResult)
2427  goto endTest;
2428  }
2429 
2430  // nonDimensional run using PlasComCM (rho, C, T) non-dimensionalization
2431  {
2432  ix::sys::ChDir("../nonDimenPC1X");
2433 
2434  if(myRank == 0)
2435  std::cout << "Running nonDimensional PlasCom2... (PC1ND)" << std::endl;
2436  returnValue = application::ApplicationDriver(simulationApplicationNonDimen2);
2437  if(myRank == 0)
2438  std::cout << "Done running nonDimensional PlasCom2 (PC1ND)." << std::endl;
2439 
2440  if(returnValue > 0){
2441  std::cout << testFunctionName << ": PlasCom2 returned error code ("
2442  << returnValue << "). Test failed." << std::endl;
2443  viscousShockTestResult = false;
2444  goto endTest;
2445  }
2446 
2447  if(myRank == 0){
2448  if(!ix::sys::FILEEXISTS("PlasCom2_000040000.h5")){
2449  std::cout << testFunctionName << ": PlasCom2 failed to produce expected" << std::endl
2450  << "output file (PlasCom2_000002000.h5). Test failed." << std::endl;
2451  viscousShockTestResult = false;
2452  }
2453  }
2454  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2455  if(!viscousShockTestResult)
2456  goto endTest;
2457  }
2458 
2459  // get the domain from the application and copy the state
2460  // make a new scope for variable declaration, needed for goto usage :-(
2461  {
2462  DomainVector &simulationDomain(simulationApplicationDimen.GetDomains());
2463 
2464  // generate a state containing the exact solution and write it to a file
2465  int direction=0;
2466  // single grid
2467  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
2468  state_t exactState(simulationDomain[0].State(0));
2469  if(testfixtures::viscid::GenerateViscidShockExact(exactGrid,exactState,direction)){
2470  viscousShockTestResult = false;
2471  std::cout << "Fail to calculated exact viscid shock solution" << std::endl;
2472  }
2473  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2474  if(!viscousShockTestResult)
2475  goto endTest;
2476 
2477  size_t numPointsBuffer = exactGrid.BufferSize();
2478 
2479  int rhoHandle = exactState.GetDataIndex("rho");
2480  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
2481  double *rhoPtr = rhoData.Data<double>();
2482 
2483  int rhoVHandle = exactState.GetDataIndex("rhoV");
2484  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
2485  double *rhoVPtr = rhoVData.Data<double>();
2486 
2487  int rhoEHandle = exactState.GetDataIndex("rhoE");
2488  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
2489  double *rhoEPtr = rhoEData.Data<double>();
2490 
2491  int velocityHandle = exactState.GetDataIndex("velocity");
2492  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
2493  double *velocityPtr = velocityData.Data<double>();
2494 
2495  int pressureHandle = exactState.GetDataIndex("pressure");
2496  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
2497  double *pressurePtr = pressureData.Data<double>();
2498 
2499  int temperatureHandle = exactState.GetDataIndex("temperature");
2500  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
2501  double *temperaturePtr = temperatureData.Data<double>();
2502 
2503  {
2504  // output exact dimensional solution
2505  if(myRank == 0)
2506  std::cout << "Writing exact solution for dimensional PlasCom2" << std::endl;
2507  ix::sys::ChDir("../dimenX");
2508  std::ostringstream Ostr;
2509  Ostr << "exact_ViscousShock_dimen_40000";
2510  const std::string domainName(Ostr.str());
2511  const std::string hdfFileName(domainName+".h5");
2512  const std::string xdmfFileName(domainName+".xdmf");
2513  const std::string geometryName("shockTube");
2514  const std::string gridName("grid1");
2515  if(ix::sys::FILEEXISTS(hdfFileName)){
2516  ix::sys::Remove(hdfFileName);
2517  ix::sys::Remove(xdmfFileName);
2518  }
2519  fixtures::ConfigurationType testConfig;
2520  state_t paramState;
2521  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
2522  gridName,exactGrid,exactState,paramState,
2523  testConfig,compareStream)){
2524  std::cout << compareStream.str();
2525  std::cout << "Output to HDF5 failed:\n";
2526  std::cout << "Continuing....\n";
2527 
2528  }
2529  pcpp::io::RenewStream(compareStream);
2530  }
2531 
2532  {
2533  if(myRank == 0)
2534  std::cout << "Writing exact solution for nonDimensional PlasCom2"
2535  << std::endl;
2536  ix::sys::ChDir("../nonDimenX");
2537  // output exact nonDimensional solution in the PlasComCM dimensionalization
2538  // we set the reference length to 1 so the grid doesn't need to change
2539  double refRho = 1.225;
2540  double refPress = 101325;
2541  double refVel = sqrt(refPress/refRho);
2542  double refTemp = 288.15;
2543  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
2544 
2545  rhoPtr[iPoint] /= refRho;
2546  rhoVPtr[iPoint] /= refRho*refVel;
2547  rhoEPtr[iPoint] /= refVel*refVel*refRho;
2548 
2549  velocityPtr[iPoint] /= refVel;
2550  pressurePtr[iPoint] /= refRho*refVel*refVel;
2551  temperaturePtr[iPoint] /= refTemp;
2552  }
2553 
2554  std::ostringstream Ostr;
2555  Ostr << "exact_ViscousShock_ND_40000";
2556  const std::string domainName(Ostr.str());
2557  const std::string hdfFileName(domainName+".h5");
2558  const std::string xdmfFileName(domainName+".xdmf");
2559  const std::string geometryName("shockTube");
2560  const std::string gridName("grid1");
2561  if(ix::sys::FILEEXISTS(hdfFileName)){
2562  ix::sys::Remove(hdfFileName);
2563  ix::sys::Remove(xdmfFileName);
2564  }
2565  fixtures::ConfigurationType testConfig;
2566  state_t paramState;
2567  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
2568  gridName,exactGrid,exactState,paramState,
2569  testConfig,compareStream)){
2570  std::cout << compareStream.str();
2571  std::cout << "Output to HDF5 failed:\n";
2572  std::cout << "Continuing....\n";
2573 
2574  }
2575  pcpp::io::RenewStream(compareStream);
2576  }
2577 
2578  {
2579  if(myRank == 0) {
2580  std::cout << "Writing exact solution for nonDimensional PlasCom2"
2581  << " with legacy PlasComCM nonDimensionalization" << std::endl;
2582  }
2583  ix::sys::ChDir("../nonDimenPC1X");
2584  // output exact nonDimensional solution in the PlasComCM dimensionalization
2585  // we set the reference length to 1 so the grid doesn't need to change
2586  //
2587  // since we already modified the solution to be in the default PlasCom2
2588  // nonDimensionalization, we only need to multiply by the converseion
2589  // factor between the two
2590  double gamma = 1.4;
2591  double refRho = 1.0;
2592  double refPress = 101325;
2593  double refSndSpd = sqrt(gamma);
2594  double refTemp = gamma-1;
2595  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
2596  rhoPtr[iPoint] /= refRho;
2597  rhoVPtr[iPoint] /= refRho*refSndSpd;
2598  rhoEPtr[iPoint] /= refSndSpd*refSndSpd*refRho;
2599 
2600  velocityPtr[iPoint] /= refSndSpd;
2601  pressurePtr[iPoint] /= refRho*refSndSpd*refSndSpd;
2602  temperaturePtr[iPoint] /= refTemp;
2603  }
2604 
2605  std::ostringstream Ostr;
2606  Ostr << "exact_ViscousShock_PC1ND_40000";
2607  const std::string domainName(Ostr.str());
2608  const std::string hdfFileName(domainName+".h5");
2609  const std::string xdmfFileName(domainName+".xdmf");
2610  const std::string geometryName("shockTube");
2611  const std::string gridName("grid1");
2612  if(ix::sys::FILEEXISTS(hdfFileName)){
2613  ix::sys::Remove(hdfFileName);
2614  ix::sys::Remove(xdmfFileName);
2615  }
2616  fixtures::ConfigurationType testConfig;
2617  state_t paramState;
2618  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
2619  gridName,exactGrid,exactState,paramState,
2620  testConfig,compareStream)){
2621  std::cout << compareStream.str();
2622  std::cout << "Output to HDF5 failed:\n";
2623  std::cout << "Continuing....\n";
2624 
2625  }
2626  pcpp::io::RenewStream(compareStream);
2627  }
2628  ix::sys::ChDir("../");
2629  }
2630 
2631  if(myRank == 0) {
2632  if(!ix::sys::FILEEXISTS("dimenX/exact_ViscousShock_dimen_40000.h5")){
2633  std::cout << testFunctionName << ": failed to produce expected" << std::endl
2634  << "output file (exact_ViscousShock_dimen_40000.h5). Test failed." << std::endl;
2635  viscousShockTestResult = false;
2636  }
2637 
2638  if(!ix::sys::FILEEXISTS("nonDimenX/exact_ViscousShock_ND_40000.h5")){
2639  std::cout << testFunctionName << ": failed to produce expected" << std::endl
2640  << "output file (exact_ViscousShock_ND_40000.h5). Test failed." << std::endl;
2641  viscousShockTestResult = false;
2642  }
2643 
2644  if(!ix::sys::FILEEXISTS("nonDimenPC1X/exact_ViscousShock_PC1ND_40000.h5")){
2645  std::cout << testFunctionName << ": failed to produce expected" << std::endl
2646  << "output file (exact_ViscousShock_PC1ND_40000.h5). Test failed." << std::endl;
2647  viscousShockTestResult = false;
2648  }
2649  }
2650  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2651  if(!viscousShockTestResult)
2652  goto endTest;
2653 
2654  {
2655  if(myRank == 0) {
2656  std::cout << "Comparing exact solution against dimensional PlasCom2" << std::endl;
2657  ix::sys::ChDir("dimenX");
2658  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
2659  "exact_ViscousShock_dimen_40000.h5",
2660  errorToleranceExact,compareStream);
2661 
2662  if(returnValue > 0){
2663  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
2664  << std::endl
2665  << compareStream.str()
2666  << std::endl;
2667  viscousShockTestResult = false;
2668  } else {
2669  std::cout << "Passed!" << std::endl;
2670  }
2671 
2672  pcpp::io::RenewStream(compareStream);
2673  ix::sys::ChDir("../");
2674  }
2675  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2676  }
2677 
2678  {
2679  if(myRank == 0) {
2680  std::cout << "Comparing exact solution against nonDimensional PlasCom2" << std::endl;
2681  ix::sys::ChDir("nonDimenX");
2682  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
2683  "exact_ViscousShock_ND_40000.h5",
2684  errorToleranceExact,compareStream);
2685  if(returnValue > 0){
2686  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
2687  << std::endl
2688  << compareStream.str()
2689  << std::endl;
2690  viscousShockTestResult = false;
2691  } else {
2692  std::cout << "Passed!" << std::endl;
2693  }
2694  pcpp::io::RenewStream(compareStream);
2695  ix::sys::ChDir("../");
2696  }
2697  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2698  }
2699 
2700  {
2701  if(myRank == 0) {
2702  std::cout << "Comparing exact solution against nonDimensional PlasCom2"
2703  << " with legacy PlasComCM nonDimensionalization" << std::endl;
2704  ix::sys::ChDir("nonDimenPC1X");
2705  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
2706  "exact_ViscousShock_PC1ND_40000.h5",
2707  errorToleranceExact,compareStream);
2708  if(returnValue > 0){
2709  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
2710  << std::endl
2711  << compareStream.str()
2712  << std::endl;
2713  viscousShockTestResult = false;
2714  } else {
2715  std::cout << "Passed!" << std::endl;
2716  }
2717  pcpp::io::RenewStream(compareStream);
2718 
2719  std::cout << "Comparing PlasComCM results against nonDimensional PlasCom2"
2720  << " with legacy PlasComCM nonDimensionalization" << std::endl;
2721  ix::sys::ChDir("nonDimenPC1X");
2722  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
2723  "PlasComCM/RocFlo-CM.00040000.h5",
2724  errorTolerancePC1_2,compareStream);
2725  if(returnValue > 0){
2726  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
2727  << std::endl
2728  << compareStream.str()
2729  << std::endl;
2730  viscousShockTestResult = false;
2731  } else {
2732  std::cout << "Passed!" << std::endl;
2733  }
2734  pcpp::io::RenewStream(compareStream);
2735  ix::sys::ChDir("../");
2736  }
2737  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2738  testComm.Barrier();
2739  }
2740 
2741  endTest:
2742  {
2743  ix::sys::ChDir(originalDirectory);
2744  parallelUnitResults.UpdateResult("Integrated:Viscid:VSRectilinear2DX",viscousShockTestResult);
2745  }
2746 
2747  return;
2748 }
2749 
2751  pcpp::CommunicatorType &testComm)
2752 {
2753 
2754  bool viscousShockTestResult = true;
2755  int myRank = testComm.Rank();
2756  int numProc = testComm.NProc();
2757 
2758  const std::string testFunctionName("TestIntegrated_VSCurvilinear2DX");
2759  const std::string testSuiteName("ViscidTest");
2760  const std::string testCaseName("ViscousShock2DX");
2761  const std::string testDirectory(testSuiteName+"/"+testCaseName);
2762  const std::string originalDirectory(ix::sys::CWD());
2763  std::ostringstream compareStream;
2764  double errorTolerancePC1_2 = 1.e-4;
2765  double errorToleranceExact = 5.e-2;
2766  if(myRank == 0)
2767  std::cout << "Running " << testFunctionName << std::endl;
2768 
2769  int argc = 4;
2770  const char *argv[] = {"plascom2x","-c","curvilinear.config",NULL};
2771 
2772  pcpp::CommunicatorType dimenComm(testComm);
2773  pcpp::CommunicatorType nonDimenComm(testComm);
2774  pcpp::CommunicatorType nonDimenComm2(testComm);
2775 
2776  plascom2::application simulationApplicationDimen(argc,const_cast<char **>(argv),dimenComm.Comm());
2777  plascom2::application simulationApplicationNonDimen(argc,const_cast<char **>(argv),nonDimenComm.Comm());
2778  plascom2::application simulationApplicationNonDimen2(argc,const_cast<char **>(argv),nonDimenComm2.Comm());
2779 
2780  int returnValue = 0;
2781 
2782  if(myRank == 0) {
2783  if(!ix::sys::FILEEXISTS(testDirectory)){
2784  viscousShockTestResult = false;
2785  std::cout << ": testing directory ("
2786  << testDirectory << ") did not exist. Aborting test."
2787  << std::endl;
2788  }
2789  }
2790  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2791  if(!viscousShockTestResult)
2792  goto endTest;
2793  ix::sys::ChDir(testDirectory);
2794 
2795 
2796  // simulation in dimensional coordinates
2797  {
2798  ix::sys::ChDir("dimenX");
2799 
2800  if(myRank == 0)
2801  std::cout << "Running dimensional PlasCom2..." << std::endl;
2802  returnValue = application::ApplicationDriver(simulationApplicationDimen);
2803  if(myRank == 0)
2804  std::cout << "Done running dimensional PlasCom2." << std::endl;
2805 
2806  if(returnValue > 0){
2807  std::cout << testFunctionName << ": PlasCom2 returned error code ("
2808  << returnValue << "). Test failed." << std::endl;
2809  viscousShockTestResult = false;
2810  goto endTest;
2811  }
2812 
2813  if(myRank == 0){
2814  if(!ix::sys::FILEEXISTS("PlasCom2_000040000.h5")){
2815  std::cout << testFunctionName << ": PlasCom2 failed to produce expected"
2816  << "output file (PlasCom2_000040000.h5). Test failed." << std::endl;
2817  viscousShockTestResult = false;
2818  }
2819  }
2820  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2821  if(!viscousShockTestResult)
2822  goto endTest;
2823  }
2824 
2825  // nonDimensional run using standard (rho, P, T) non-dimensionalization
2826  {
2827  ix::sys::ChDir("../nonDimenX");
2828 
2829  if(myRank == 0)
2830  std::cout << "Running nonDimensional PlasCom2..." << std::endl;
2831  returnValue = application::ApplicationDriver(simulationApplicationNonDimen);
2832  if(myRank == 0)
2833  std::cout << "Done running nonDimensional PlasCom2." << std::endl;
2834 
2835  if(returnValue > 0){
2836  std::cout << testFunctionName << ": PlasCom2 returned error code ("
2837  << returnValue << "). Test failed." << std::endl;
2838  viscousShockTestResult = false;
2839  goto endTest;
2840  }
2841 
2842  if(myRank == 0){
2843  if(!ix::sys::FILEEXISTS("PlasCom2_000040000.h5")){
2844  std::cout << testFunctionName << ": PlasCom2 failed to produce expected" << std::endl
2845  << "output file (PlasCom2_000002000.h5). Test failed." << std::endl;
2846  viscousShockTestResult = false;
2847  }
2848  }
2849  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2850  if(!viscousShockTestResult)
2851  goto endTest;
2852  }
2853 
2854  // nonDimensional run using PlasComCM (rho, C, T) non-dimensionalization
2855  {
2856  ix::sys::ChDir("../nonDimenPC1X");
2857 
2858  if(myRank == 0)
2859  std::cout << "Running nonDimensional PlasCom2... (PC1ND)" << std::endl;
2860  returnValue = application::ApplicationDriver(simulationApplicationNonDimen2);
2861  if(myRank == 0)
2862  std::cout << "Done running nonDimensional PlasCom2 (PC1ND)." << std::endl;
2863 
2864  if(returnValue > 0){
2865  std::cout << testFunctionName << ": PlasCom2 returned error code ("
2866  << returnValue << "). Test failed." << std::endl;
2867  viscousShockTestResult = false;
2868  goto endTest;
2869  }
2870 
2871  if(myRank == 0){
2872  if(!ix::sys::FILEEXISTS("PlasCom2_000040000.h5")){
2873  std::cout << testFunctionName << ": PlasCom2 failed to produce expected" << std::endl
2874  << "output file (PlasCom2_000002000.h5). Test failed." << std::endl;
2875  viscousShockTestResult = false;
2876  }
2877  }
2878  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2879  if(!viscousShockTestResult)
2880  goto endTest;
2881  }
2882 
2883  // get the domain from the application and copy the state
2884  // make a new scope for variable declaration, needed for goto usage :-(
2885  {
2886  DomainVector &simulationDomain(simulationApplicationDimen.GetDomains());
2887 
2888  // generate a state containing the exact solution and write it to a file
2889  int direction=0;
2890  // single grid
2891  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
2892  state_t exactState(simulationDomain[0].State(0));
2893  if(testfixtures::viscid::GenerateViscidShockExact(exactGrid,exactState,direction)){
2894  viscousShockTestResult = false;
2895  std::cout << "Fail to calculated exact viscid shock solution" << std::endl;
2896  }
2897  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
2898  if(!viscousShockTestResult)
2899  goto endTest;
2900 
2901  size_t numPointsBuffer = exactGrid.BufferSize();
2902 
2903  int rhoHandle = exactState.GetDataIndex("rho");
2904  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
2905  double *rhoPtr = rhoData.Data<double>();
2906 
2907  int rhoVHandle = exactState.GetDataIndex("rhoV");
2908  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
2909  double *rhoVPtr = rhoVData.Data<double>();
2910 
2911  int rhoEHandle = exactState.GetDataIndex("rhoE");
2912  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
2913  double *rhoEPtr = rhoEData.Data<double>();
2914 
2915  int velocityHandle = exactState.GetDataIndex("velocity");
2916  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
2917  double *velocityPtr = velocityData.Data<double>();
2918 
2919  int pressureHandle = exactState.GetDataIndex("pressure");
2920  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
2921  double *pressurePtr = pressureData.Data<double>();
2922 
2923  int temperatureHandle = exactState.GetDataIndex("temperature");
2924  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
2925  double *temperaturePtr = temperatureData.Data<double>();
2926 
2927  {
2928  // output exact dimensional solution
2929  if(myRank == 0)
2930  std::cout << "Writing exact solution for dimensional PlasCom2" << std::endl;
2931  ix::sys::ChDir("../dimenX");
2932  std::ostringstream Ostr;
2933  Ostr << "exact_ViscousShock_dimen_40000";
2934  const std::string domainName(Ostr.str());
2935  const std::string hdfFileName(domainName+".h5");
2936  const std::string xdmfFileName(domainName+".xdmf");
2937  const std::string geometryName("shockTube");
2938  const std::string gridName("grid1");
2939  if(ix::sys::FILEEXISTS(hdfFileName)){
2940  ix::sys::Remove(hdfFileName);
2941  ix::sys::Remove(xdmfFileName);
2942  }
2943  fixtures::ConfigurationType testConfig;
2944  state_t paramState;
2945  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
2946  gridName,exactGrid,exactState,paramState,
2947  testConfig,compareStream)){
2948  std::cout << compareStream.str();
2949  std::cout << "Output to HDF5 failed:\n";
2950  std::cout << "Continuing....\n";
2951 
2952  }
2953  pcpp::io::RenewStream(compareStream);
2954  }
2955 
2956  {
2957  if(myRank == 0)
2958  std::cout << "Writing exact solution for nonDimensional PlasCom2"
2959  << std::endl;
2960  ix::sys::ChDir("../nonDimenX");
2961  // output exact nonDimensional solution in the PlasComCM dimensionalization
2962  // we set the reference length to 1 so the grid doesn't need to change
2963  double refRho = 1.225;
2964  double refPress = 101325;
2965  double refVel = sqrt(refPress/refRho);
2966  double refTemp = 288.15;
2967  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
2968 
2969  rhoPtr[iPoint] /= refRho;
2970  rhoVPtr[iPoint] /= refRho*refVel;
2971  rhoEPtr[iPoint] /= refVel*refVel*refRho;
2972 
2973  velocityPtr[iPoint] /= refVel;
2974  pressurePtr[iPoint] /= refRho*refVel*refVel;
2975  temperaturePtr[iPoint] /= refTemp;
2976  }
2977 
2978  std::ostringstream Ostr;
2979  Ostr << "exact_ViscousShock_ND_40000";
2980  const std::string domainName(Ostr.str());
2981  const std::string hdfFileName(domainName+".h5");
2982  const std::string xdmfFileName(domainName+".xdmf");
2983  const std::string geometryName("shockTube");
2984  const std::string gridName("grid1");
2985  if(ix::sys::FILEEXISTS(hdfFileName)){
2986  ix::sys::Remove(hdfFileName);
2987  ix::sys::Remove(xdmfFileName);
2988  }
2989  fixtures::ConfigurationType testConfig;
2990  state_t paramState;
2991  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
2992  gridName,exactGrid,exactState,paramState,
2993  testConfig,compareStream)){
2994  std::cout << compareStream.str();
2995  std::cout << "Output to HDF5 failed:\n";
2996  std::cout << "Continuing....\n";
2997 
2998  }
2999  pcpp::io::RenewStream(compareStream);
3000  }
3001 
3002  {
3003  if(myRank == 0) {
3004  std::cout << "Writing exact solution for nonDimensional PlasCom2"
3005  << " with legacy PlasComCM nonDimensionalization" << std::endl;
3006  }
3007  ix::sys::ChDir("../nonDimenPC1X");
3008  // output exact nonDimensional solution in the PlasComCM dimensionalization
3009  // we set the reference length to 1 so the grid doesn't need to change
3010  //
3011  // since we already modified the solution to be in the default PlasCom2
3012  // nonDimensionalization, we only need to multiply by the converseion
3013  // factor between the two
3014  double gamma = 1.4;
3015  double refRho = 1.0;
3016  double refPress = 101325;
3017  double refSndSpd = sqrt(gamma);
3018  double refTemp = gamma-1;
3019  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
3020  rhoPtr[iPoint] /= refRho;
3021  rhoVPtr[iPoint] /= refRho*refSndSpd;
3022  rhoEPtr[iPoint] /= refSndSpd*refSndSpd*refRho;
3023 
3024  velocityPtr[iPoint] /= refSndSpd;
3025  pressurePtr[iPoint] /= refRho*refSndSpd*refSndSpd;
3026  temperaturePtr[iPoint] /= refTemp;
3027  }
3028 
3029  std::ostringstream Ostr;
3030  Ostr << "exact_ViscousShock_PC1ND_40000";
3031  const std::string domainName(Ostr.str());
3032  const std::string hdfFileName(domainName+".h5");
3033  const std::string xdmfFileName(domainName+".xdmf");
3034  const std::string geometryName("shockTube");
3035  const std::string gridName("grid1");
3036  if(ix::sys::FILEEXISTS(hdfFileName)){
3037  ix::sys::Remove(hdfFileName);
3038  ix::sys::Remove(xdmfFileName);
3039  }
3040  fixtures::ConfigurationType testConfig;
3041  state_t paramState;
3042  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
3043  gridName,exactGrid,exactState,paramState,
3044  testConfig,compareStream)){
3045  std::cout << compareStream.str();
3046  std::cout << "Output to HDF5 failed:\n";
3047  std::cout << "Continuing....\n";
3048 
3049  }
3050  pcpp::io::RenewStream(compareStream);
3051  }
3052  ix::sys::ChDir("../");
3053  }
3054 
3055  if(myRank == 0) {
3056  if(!ix::sys::FILEEXISTS("dimenX/exact_ViscousShock_dimen_40000.h5")){
3057  std::cout << testFunctionName << ": failed to produce expected" << std::endl
3058  << "output file (exact_ViscousShock_dimen_40000.h5). Test failed." << std::endl;
3059  viscousShockTestResult = false;
3060  }
3061 
3062  if(!ix::sys::FILEEXISTS("nonDimenX/exact_ViscousShock_ND_40000.h5")){
3063  std::cout << testFunctionName << ": failed to produce expected" << std::endl
3064  << "output file (exact_ViscousShock_ND_40000.h5). Test failed." << std::endl;
3065  viscousShockTestResult = false;
3066  }
3067 
3068  if(!ix::sys::FILEEXISTS("nonDimenPC1X/exact_ViscousShock_PC1ND_40000.h5")){
3069  std::cout << testFunctionName << ": failed to produce expected" << std::endl
3070  << "output file (exact_ViscousShock_PC1ND_40000.h5). Test failed." << std::endl;
3071  viscousShockTestResult = false;
3072  }
3073  }
3074  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3075  if(!viscousShockTestResult)
3076  goto endTest;
3077 
3078  {
3079  if(myRank == 0) {
3080  std::cout << "Comparing exact solution against dimensional PlasCom2" << std::endl;
3081  ix::sys::ChDir("dimenX");
3082  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
3083  "exact_ViscousShock_dimen_40000.h5",
3084  errorToleranceExact,compareStream);
3085 
3086  if(returnValue > 0){
3087  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
3088  << std::endl
3089  << compareStream.str()
3090  << std::endl;
3091  viscousShockTestResult = false;
3092  } else {
3093  std::cout << "Passed!" << std::endl;
3094  }
3095 
3096  pcpp::io::RenewStream(compareStream);
3097  ix::sys::ChDir("../");
3098  }
3099  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3100  }
3101 
3102  {
3103  if(myRank == 0) {
3104  std::cout << "Comparing exact solution against nonDimensional PlasCom2" << std::endl;
3105  ix::sys::ChDir("nonDimenX");
3106  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
3107  "exact_ViscousShock_ND_40000.h5",
3108  errorToleranceExact,compareStream);
3109  if(returnValue > 0){
3110  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
3111  << std::endl
3112  << compareStream.str()
3113  << std::endl;
3114  viscousShockTestResult = false;
3115  } else {
3116  std::cout << "Passed!" << std::endl;
3117  }
3118  pcpp::io::RenewStream(compareStream);
3119  ix::sys::ChDir("../");
3120  }
3121  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3122  }
3123 
3124  {
3125  if(myRank == 0) {
3126  std::cout << "Comparing exact solution against nonDimensional PlasCom2"
3127  << " with legacy PlasComCM nonDimensionalization" << std::endl;
3128  ix::sys::ChDir("nonDimenPC1X");
3129  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
3130  "exact_ViscousShock_PC1ND_40000.h5",
3131  errorToleranceExact,compareStream);
3132  if(returnValue > 0){
3133  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
3134  << std::endl
3135  << compareStream.str()
3136  << std::endl;
3137  viscousShockTestResult = false;
3138  } else {
3139  std::cout << "Passed!" << std::endl;
3140  }
3141  pcpp::io::RenewStream(compareStream);
3142 
3143  std::cout << "Comparing PlasComCM results against nonDimensional PlasCom2"
3144  << " with legacy PlasComCM nonDimensionalization" << std::endl;
3145  ix::sys::ChDir("nonDimenPC1X");
3146  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
3147  "PlasComCM/RocFlo-CM.00040000.h5",
3148  errorTolerancePC1_2,compareStream);
3149  if(returnValue > 0){
3150  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
3151  << std::endl
3152  << compareStream.str()
3153  << std::endl;
3154  viscousShockTestResult = false;
3155  } else {
3156  std::cout << "Passed!" << std::endl;
3157  }
3158  pcpp::io::RenewStream(compareStream);
3159  ix::sys::ChDir("../");
3160  }
3161  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3162  testComm.Barrier();
3163  }
3164 
3165  endTest:
3166  {
3167  ix::sys::ChDir(originalDirectory);
3168  parallelUnitResults.UpdateResult("Integrated:Viscid:VSCurvilinear2DX",viscousShockTestResult);
3169  }
3170 
3171  return;
3172 }
3173 
3175  pcpp::CommunicatorType &testComm)
3176 {
3177 
3178  bool viscousShockTestResult = true;
3179  int myRank = testComm.Rank();
3180  int numProc = testComm.NProc();
3181 
3182 
3183  const std::string testFunctionName("TestIntegrated_ViscousShock2DY");
3184  const std::string testSuiteName("ViscidTest");
3185  const std::string testCaseName("ViscousShock2DY");
3186  const std::string testDirectory(testSuiteName+"/"+testCaseName);
3187  const std::string originalDirectory(ix::sys::CWD());
3188  std::ostringstream compareStream;
3189 // double errorTolerancePC1_2 = 1.e-4;
3190  double errorToleranceExact = 5.e-2;
3191  if(myRank == 0)
3192  std::cout << "Running " << testFunctionName << std::endl;
3193 
3194  int argc = 4;
3195  const char *argv[] = {"plascom2x","-c","shock.test.config",NULL};
3196 
3197  pcpp::CommunicatorType nonDimenComm(testComm);
3198  plascom2::application simulationApplicationNonDimen(argc,const_cast<char **>(argv),nonDimenComm.Comm());
3199 
3200  int returnValue = 0;
3201 
3202  if(myRank == 0) {
3203  if(!ix::sys::FILEEXISTS(testDirectory)){
3204  viscousShockTestResult = false;
3205  std::cout << "TestIntegrated_viscousShock2D: testing directory ("
3206  << testDirectory << ") did not exist. Aborting test."
3207  << std::endl;
3208  }
3209  }
3210  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3211  if(!viscousShockTestResult)
3212  goto endTest;
3213  ix::sys::ChDir(testDirectory);
3214 
3215  // nonDimensional run using standard (rho, P, T) non-dimensionalization
3216  {
3217  ix::sys::ChDir("nonDimen");
3218 
3219  if(myRank == 0)
3220  std::cout << "Running nonDimensional PlasCom2..." << std::endl;
3221  returnValue = application::ApplicationDriver(simulationApplicationNonDimen);
3222  if(myRank == 0)
3223  std::cout << "Done running nonDimensional PlasCom2." << std::endl;
3224 
3225  if(returnValue > 0){
3226  std::cout << "TestIntegrated_ViscousShock2DY: PlasCom2 returned error code ("
3227  << returnValue << "). Test failed." << std::endl;
3228  viscousShockTestResult = false;
3229  goto endTest;
3230  }
3231 
3232  if(myRank == 0){
3233  if(!ix::sys::FILEEXISTS("PlasCom2_000040000.h5")){
3234  std::cout << "TestIntegrated_ViscousShock2DY: PlasCom2 failed to produce expected" << std::endl
3235  << "output file (PlasCom2_000040000.h5). Test failed." << std::endl;
3236  viscousShockTestResult = false;
3237  }
3238  }
3239  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3240  if(!viscousShockTestResult)
3241  goto endTest;
3242  }
3243 
3244  // get the domain from the application and copy the state
3245  // make a new scope for variable declaration, needed for goto usage :-(
3246  {
3247  DomainVector &simulationDomain(simulationApplicationNonDimen.GetDomains());
3248 
3249  // generate a state containing the exact solution and write it to a file
3250  int direction=1;
3251  // single grid
3252  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
3253  state_t exactState(simulationDomain[0].State(0));
3254  if(testfixtures::viscid::GenerateViscidShockExact(exactGrid,exactState,direction)){
3255  viscousShockTestResult = false;
3256  std::cout << "Fail to calculated exact viscid shock solution" << std::endl;
3257  }
3258  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3259  if(!viscousShockTestResult)
3260  goto endTest;
3261 
3262  size_t numPointsBuffer = exactGrid.BufferSize();
3263 
3264  int rhoHandle = exactState.GetDataIndex("rho");
3265  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
3266  double *rhoPtr = rhoData.Data<double>();
3267 
3268  int rhoVHandle = exactState.GetDataIndex("rhoV");
3269  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
3270  double *rhoVPtr = rhoVData.Data<double>();
3271 
3272  int rhoEHandle = exactState.GetDataIndex("rhoE");
3273  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
3274  double *rhoEPtr = rhoEData.Data<double>();
3275 
3276  int velocityHandle = exactState.GetDataIndex("velocity");
3277  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
3278  double *velocityPtr = velocityData.Data<double>();
3279 
3280  int pressureHandle = exactState.GetDataIndex("pressure");
3281  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
3282  double *pressurePtr = pressureData.Data<double>();
3283 
3284  int temperatureHandle = exactState.GetDataIndex("temperature");
3285  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
3286  double *temperaturePtr = temperatureData.Data<double>();
3287 
3288  {
3289  if(myRank == 0)
3290  std::cout << "Writing exact solution for nonDimensional PlasCom2"
3291  << std::endl;
3292  //ix::sys::ChDir("../nonDimenX");
3293  // output exact nonDimensional solution in the PlasComCM dimensionalization
3294  // we set the reference length to 1 so the grid doesn't need to change
3295  double refRho = 1.225;
3296  double refPress = 101325;
3297  double refVel = sqrt(refPress/refRho);
3298  double refTemp = 288.15;
3299  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
3300 
3301  rhoPtr[iPoint] /= refRho;
3302  rhoVPtr[iPoint+numPointsBuffer] /= refRho*refVel;
3303  rhoEPtr[iPoint] /= refVel*refVel*refRho;
3304 
3305  velocityPtr[iPoint+numPointsBuffer] /= refVel;
3306  pressurePtr[iPoint] /= refRho*refVel*refVel;
3307  temperaturePtr[iPoint] /= refTemp;
3308  }
3309 
3310  std::ostringstream Ostr;
3311  Ostr << "exact_ViscousShock_ND_40000";
3312  const std::string domainName(Ostr.str());
3313  const std::string hdfFileName(domainName+".h5");
3314  const std::string xdmfFileName(domainName+".xdmf");
3315  const std::string geometryName("shockTube");
3316  const std::string gridName("grid1");
3317  if(ix::sys::FILEEXISTS(hdfFileName)){
3318  ix::sys::Remove(hdfFileName);
3319  ix::sys::Remove(xdmfFileName);
3320  }
3321  fixtures::ConfigurationType testConfig;
3322  state_t paramState;
3323  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
3324  gridName,exactGrid,exactState,paramState,
3325  testConfig,compareStream)){
3326  std::cout << compareStream.str();
3327  std::cout << "Output to HDF5 failed:\n";
3328  std::cout << "Continuing....\n";
3329 
3330  }
3331  pcpp::io::RenewStream(compareStream);
3332  }
3333  }
3334 
3335  if(myRank == 0) {
3336  if(!ix::sys::FILEEXISTS("exact_ViscousShock_ND_40000.h5")){
3337  std::cout << "TestIntegrated_ViscousShock2DY: failed to produce expected" << std::endl
3338  << "output file (exact_ViscousShock_ND_40000.h5). Test failed." << std::endl;
3339  viscousShockTestResult = false;
3340  }
3341  }
3342  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3343  if(!viscousShockTestResult)
3344  goto endTest;
3345 
3346  {
3347  if(myRank == 0) {
3348  std::cout << "Comparing exact solution against nonDimensional PlasCom2" << std::endl;
3349  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
3350  "exact_ViscousShock_ND_40000.h5",
3351  errorToleranceExact,compareStream);
3352  if(returnValue > 0){
3353  std::cout << "TestIntegrated_ViscousShock2DY: Resulting state comparison yielded negative result:"
3354  << std::endl
3355  << compareStream.str()
3356  << std::endl;
3357  viscousShockTestResult = false;
3358  } else {
3359  std::cout << "Passed!" << std::endl;
3360  }
3361  pcpp::io::RenewStream(compareStream);
3362  }
3363  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3364  }
3365 
3366  endTest:
3367  {
3368  ix::sys::ChDir(originalDirectory);
3369  parallelUnitResults.UpdateResult("Integrated:Viscid:ViscousShock2DY",viscousShockTestResult);
3370  }
3371 
3372  return;
3373 }
3374 
3376  pcpp::CommunicatorType &testComm)
3377 {
3378 
3379  bool viscousShockTestResult = true;
3380  int myRank = testComm.Rank();
3381  int numProc = testComm.NProc();
3382 
3383 
3384  const std::string testFunctionName("TestIntegrated_VSRectilinear2DY");
3385  const std::string testSuiteName("ViscidTest");
3386  const std::string testCaseName("ViscousShock2DY");
3387  const std::string testDirectory(testSuiteName+"/"+testCaseName);
3388  const std::string originalDirectory(ix::sys::CWD());
3389  std::ostringstream compareStream;
3390 // double errorTolerancePC1_2 = 1.e-4;
3391  double errorToleranceExact = 5.e-2;
3392  if(myRank == 0)
3393  std::cout << "Running " << testFunctionName << std::endl;
3394 
3395  int argc = 4;
3396  const char *argv[] = {"plascom2x","-c","rectilinear.config",NULL};
3397 
3398  pcpp::CommunicatorType nonDimenComm(testComm);
3399  plascom2::application simulationApplicationNonDimen(argc,const_cast<char **>(argv),nonDimenComm.Comm());
3400 
3401  int returnValue = 0;
3402 
3403  if(myRank == 0) {
3404  if(!ix::sys::FILEEXISTS(testDirectory)){
3405  viscousShockTestResult = false;
3406  std::cout << testFunctionName << ": testing directory ("
3407  << testDirectory << ") did not exist. Aborting test."
3408  << std::endl;
3409  }
3410  }
3411  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3412  if(!viscousShockTestResult)
3413  goto endTest;
3414  ix::sys::ChDir(testDirectory);
3415 
3416  // nonDimensional run using standard (rho, P, T) non-dimensionalization
3417  {
3418  ix::sys::ChDir("nonDimen");
3419 
3420  if(myRank == 0)
3421  std::cout << "Running nonDimensional PlasCom2..." << std::endl;
3422  returnValue = application::ApplicationDriver(simulationApplicationNonDimen);
3423  if(myRank == 0)
3424  std::cout << "Done running nonDimensional PlasCom2." << std::endl;
3425 
3426  if(returnValue > 0){
3427  std::cout << testFunctionName << ": PlasCom2 returned error code ("
3428  << returnValue << "). Test failed." << std::endl;
3429  viscousShockTestResult = false;
3430  goto endTest;
3431  }
3432 
3433  if(myRank == 0){
3434  if(!ix::sys::FILEEXISTS("PlasCom2_000040000.h5")){
3435  std::cout << testFunctionName << ": PlasCom2 failed to produce expected" << std::endl
3436  << "output file (PlasCom2_000040000.h5). Test failed." << std::endl;
3437  viscousShockTestResult = false;
3438  }
3439  }
3440  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3441  if(!viscousShockTestResult)
3442  goto endTest;
3443  }
3444 
3445  // get the domain from the application and copy the state
3446  // make a new scope for variable declaration, needed for goto usage :-(
3447  {
3448  DomainVector &simulationDomain(simulationApplicationNonDimen.GetDomains());
3449 
3450  // generate a state containing the exact solution and write it to a file
3451  int direction=1;
3452  // single grid
3453  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
3454  state_t exactState(simulationDomain[0].State(0));
3455  if(testfixtures::viscid::GenerateViscidShockExact(exactGrid,exactState,direction)){
3456  viscousShockTestResult = false;
3457  std::cout << "Fail to calculated exact viscid shock solution" << std::endl;
3458  }
3459  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3460  if(!viscousShockTestResult)
3461  goto endTest;
3462 
3463  size_t numPointsBuffer = exactGrid.BufferSize();
3464 
3465  int rhoHandle = exactState.GetDataIndex("rho");
3466  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
3467  double *rhoPtr = rhoData.Data<double>();
3468 
3469  int rhoVHandle = exactState.GetDataIndex("rhoV");
3470  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
3471  double *rhoVPtr = rhoVData.Data<double>();
3472 
3473  int rhoEHandle = exactState.GetDataIndex("rhoE");
3474  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
3475  double *rhoEPtr = rhoEData.Data<double>();
3476 
3477  int velocityHandle = exactState.GetDataIndex("velocity");
3478  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
3479  double *velocityPtr = velocityData.Data<double>();
3480 
3481  int pressureHandle = exactState.GetDataIndex("pressure");
3482  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
3483  double *pressurePtr = pressureData.Data<double>();
3484 
3485  int temperatureHandle = exactState.GetDataIndex("temperature");
3486  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
3487  double *temperaturePtr = temperatureData.Data<double>();
3488 
3489  {
3490  if(myRank == 0)
3491  std::cout << "Writing exact solution for nonDimensional PlasCom2"
3492  << std::endl;
3493  //ix::sys::ChDir("../nonDimenX");
3494  // output exact nonDimensional solution in the PlasComCM dimensionalization
3495  // we set the reference length to 1 so the grid doesn't need to change
3496  double refRho = 1.225;
3497  double refPress = 101325;
3498  double refVel = sqrt(refPress/refRho);
3499  double refTemp = 288.15;
3500  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
3501 
3502  rhoPtr[iPoint] /= refRho;
3503  rhoVPtr[iPoint+numPointsBuffer] /= refRho*refVel;
3504  rhoEPtr[iPoint] /= refVel*refVel*refRho;
3505 
3506  velocityPtr[iPoint+numPointsBuffer] /= refVel;
3507  pressurePtr[iPoint] /= refRho*refVel*refVel;
3508  temperaturePtr[iPoint] /= refTemp;
3509  }
3510 
3511  std::ostringstream Ostr;
3512  Ostr << "exact_ViscousShock_ND_40000";
3513  const std::string domainName(Ostr.str());
3514  const std::string hdfFileName(domainName+".h5");
3515  const std::string xdmfFileName(domainName+".xdmf");
3516  const std::string geometryName("shockTube");
3517  const std::string gridName("grid1");
3518  if(ix::sys::FILEEXISTS(hdfFileName)){
3519  ix::sys::Remove(hdfFileName);
3520  ix::sys::Remove(xdmfFileName);
3521  }
3522  fixtures::ConfigurationType testConfig;
3523  state_t paramState;
3524  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
3525  gridName,exactGrid,exactState,paramState,
3526  testConfig,compareStream)){
3527  std::cout << compareStream.str();
3528  std::cout << "Output to HDF5 failed:\n";
3529  std::cout << "Continuing....\n";
3530 
3531  }
3532  pcpp::io::RenewStream(compareStream);
3533  }
3534  }
3535 
3536  if(myRank == 0) {
3537  if(!ix::sys::FILEEXISTS("exact_ViscousShock_ND_40000.h5")){
3538  std::cout << testFunctionName << ": failed to produce expected" << std::endl
3539  << "output file (exact_ViscousShock_ND_40000.h5). Test failed." << std::endl;
3540  viscousShockTestResult = false;
3541  }
3542  }
3543  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3544  if(!viscousShockTestResult)
3545  goto endTest;
3546 
3547  {
3548  if(myRank == 0) {
3549  std::cout << "Comparing exact solution against nonDimensional PlasCom2" << std::endl;
3550  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
3551  "exact_ViscousShock_ND_40000.h5",
3552  errorToleranceExact,compareStream);
3553  if(returnValue > 0){
3554  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
3555  << std::endl
3556  << compareStream.str()
3557  << std::endl;
3558  viscousShockTestResult = false;
3559  } else {
3560  std::cout << "Passed!" << std::endl;
3561  }
3562  pcpp::io::RenewStream(compareStream);
3563  }
3564  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3565  }
3566 
3567  endTest:
3568  {
3569  ix::sys::ChDir(originalDirectory);
3570  parallelUnitResults.UpdateResult("Integrated:Viscid:VSRectilinear2DY",viscousShockTestResult);
3571  }
3572 
3573  return;
3574 }
3575 
3577  pcpp::CommunicatorType &testComm)
3578 {
3579 
3580  bool viscousShockTestResult = true;
3581  int myRank = testComm.Rank();
3582  int numProc = testComm.NProc();
3583 
3584 
3585  const std::string testFunctionName("TestIntegrated_VSCurvilinear2DY");
3586  const std::string testSuiteName("ViscidTest");
3587  const std::string testCaseName("ViscousShock2DY");
3588  const std::string testDirectory(testSuiteName+"/"+testCaseName);
3589  const std::string originalDirectory(ix::sys::CWD());
3590  std::ostringstream compareStream;
3591 // double errorTolerancePC1_2 = 1.e-4;
3592  double errorToleranceExact = 5.e-2;
3593  if(myRank == 0)
3594  std::cout << "Running " << testFunctionName << std::endl;
3595 
3596  int argc = 4;
3597  const char *argv[] = {"plascom2x","-c","curvilinear.config",NULL};
3598 
3599  pcpp::CommunicatorType nonDimenComm(testComm);
3600  plascom2::application simulationApplicationNonDimen(argc,const_cast<char **>(argv),nonDimenComm.Comm());
3601 
3602  int returnValue = 0;
3603 
3604  if(myRank == 0) {
3605  if(!ix::sys::FILEEXISTS(testDirectory)){
3606  viscousShockTestResult = false;
3607  std::cout << testFunctionName << ": testing directory ("
3608  << testDirectory << ") did not exist. Aborting test."
3609  << std::endl;
3610  }
3611  }
3612  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3613  if(!viscousShockTestResult)
3614  goto endTest;
3615  ix::sys::ChDir(testDirectory);
3616 
3617  // nonDimensional run using standard (rho, P, T) non-dimensionalization
3618  {
3619  ix::sys::ChDir("nonDimen");
3620 
3621  if(myRank == 0)
3622  std::cout << "Running nonDimensional PlasCom2..." << std::endl;
3623  returnValue = application::ApplicationDriver(simulationApplicationNonDimen);
3624  if(myRank == 0)
3625  std::cout << "Done running nonDimensional PlasCom2." << std::endl;
3626 
3627  if(returnValue > 0){
3628  std::cout << testFunctionName << ": PlasCom2 returned error code ("
3629  << returnValue << "). Test failed." << std::endl;
3630  viscousShockTestResult = false;
3631  goto endTest;
3632  }
3633 
3634  if(myRank == 0){
3635  if(!ix::sys::FILEEXISTS("PlasCom2_000040000.h5")){
3636  std::cout << testFunctionName << ": PlasCom2 failed to produce expected" << std::endl
3637  << "output file (PlasCom2_000040000.h5). Test failed." << std::endl;
3638  viscousShockTestResult = false;
3639  }
3640  }
3641  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3642  if(!viscousShockTestResult)
3643  goto endTest;
3644  }
3645 
3646  // get the domain from the application and copy the state
3647  // make a new scope for variable declaration, needed for goto usage :-(
3648  {
3649  DomainVector &simulationDomain(simulationApplicationNonDimen.GetDomains());
3650 
3651  // generate a state containing the exact solution and write it to a file
3652  int direction=1;
3653  // single grid
3654  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
3655  state_t exactState(simulationDomain[0].State(0));
3656  if(testfixtures::viscid::GenerateViscidShockExact(exactGrid,exactState,direction)){
3657  viscousShockTestResult = false;
3658  std::cout << "Fail to calculated exact viscid shock solution" << std::endl;
3659  }
3660  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3661  if(!viscousShockTestResult)
3662  goto endTest;
3663 
3664  size_t numPointsBuffer = exactGrid.BufferSize();
3665 
3666  int rhoHandle = exactState.GetDataIndex("rho");
3667  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
3668  double *rhoPtr = rhoData.Data<double>();
3669 
3670  int rhoVHandle = exactState.GetDataIndex("rhoV");
3671  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
3672  double *rhoVPtr = rhoVData.Data<double>();
3673 
3674  int rhoEHandle = exactState.GetDataIndex("rhoE");
3675  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
3676  double *rhoEPtr = rhoEData.Data<double>();
3677 
3678  int velocityHandle = exactState.GetDataIndex("velocity");
3679  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
3680  double *velocityPtr = velocityData.Data<double>();
3681 
3682  int pressureHandle = exactState.GetDataIndex("pressure");
3683  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
3684  double *pressurePtr = pressureData.Data<double>();
3685 
3686  int temperatureHandle = exactState.GetDataIndex("temperature");
3687  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
3688  double *temperaturePtr = temperatureData.Data<double>();
3689 
3690  {
3691  if(myRank == 0)
3692  std::cout << "Writing exact solution for nonDimensional PlasCom2"
3693  << std::endl;
3694  //ix::sys::ChDir("../nonDimenX");
3695  // output exact nonDimensional solution in the PlasComCM dimensionalization
3696  // we set the reference length to 1 so the grid doesn't need to change
3697  double refRho = 1.225;
3698  double refPress = 101325;
3699  double refVel = sqrt(refPress/refRho);
3700  double refTemp = 288.15;
3701  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
3702 
3703  rhoPtr[iPoint] /= refRho;
3704  rhoVPtr[iPoint+numPointsBuffer] /= refRho*refVel;
3705  rhoEPtr[iPoint] /= refVel*refVel*refRho;
3706 
3707  velocityPtr[iPoint+numPointsBuffer] /= refVel;
3708  pressurePtr[iPoint] /= refRho*refVel*refVel;
3709  temperaturePtr[iPoint] /= refTemp;
3710  }
3711 
3712  std::ostringstream Ostr;
3713  Ostr << "exact_ViscousShock_ND_40000";
3714  const std::string domainName(Ostr.str());
3715  const std::string hdfFileName(domainName+".h5");
3716  const std::string xdmfFileName(domainName+".xdmf");
3717  const std::string geometryName("shockTube");
3718  const std::string gridName("grid1");
3719  if(ix::sys::FILEEXISTS(hdfFileName)){
3720  ix::sys::Remove(hdfFileName);
3721  ix::sys::Remove(xdmfFileName);
3722  }
3723  fixtures::ConfigurationType testConfig;
3724  state_t paramState;
3725  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
3726  gridName,exactGrid,exactState,paramState,
3727  testConfig,compareStream)){
3728  std::cout << compareStream.str();
3729  std::cout << "Output to HDF5 failed:\n";
3730  std::cout << "Continuing....\n";
3731 
3732  }
3733  pcpp::io::RenewStream(compareStream);
3734  }
3735  }
3736 
3737  if(myRank == 0) {
3738  if(!ix::sys::FILEEXISTS("exact_ViscousShock_ND_40000.h5")){
3739  std::cout << testFunctionName << ": failed to produce expected" << std::endl
3740  << "output file (exact_ViscousShock_ND_40000.h5). Test failed." << std::endl;
3741  viscousShockTestResult = false;
3742  }
3743  }
3744  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3745  if(!viscousShockTestResult)
3746  goto endTest;
3747 
3748  {
3749  if(myRank == 0) {
3750  std::cout << "Comparing exact solution against nonDimensional PlasCom2" << std::endl;
3751  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
3752  "exact_ViscousShock_ND_40000.h5",
3753  errorToleranceExact,compareStream);
3754  if(returnValue > 0){
3755  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
3756  << std::endl
3757  << compareStream.str()
3758  << std::endl;
3759  viscousShockTestResult = false;
3760  } else {
3761  std::cout << "Passed!" << std::endl;
3762  }
3763  pcpp::io::RenewStream(compareStream);
3764  }
3765  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3766  }
3767 
3768  endTest:
3769  {
3770  ix::sys::ChDir(originalDirectory);
3771  parallelUnitResults.UpdateResult("Integrated:Viscid:VSCurvilinear2DY",viscousShockTestResult);
3772  }
3773 
3774  return;
3775 }
3776 
3778  pcpp::CommunicatorType &testComm)
3779 {
3780 
3781  bool viscousShockTestResult = true;
3782  int myRank = testComm.Rank();
3783  int numProc = testComm.NProc();
3784 
3785 
3786  const std::string testFunctionName("TestIntegrated_ViscousShock3DZ");
3787  const std::string testSuiteName("ViscidTest");
3788  const std::string testCaseName("ViscousShock3DZ");
3789  const std::string testDirectory(testSuiteName+"/"+testCaseName);
3790  const std::string originalDirectory(ix::sys::CWD());
3791  std::ostringstream compareStream;
3792 // double errorTolerancePC1_2 = 1.e-4;
3793  double errorToleranceExact = 5.e-2;
3794  if(myRank == 0)
3795  std::cout << "Running " << testFunctionName << std::endl;
3796 
3797  int argc = 4;
3798  const char *argv[] = {"plascom2x","-c","shock.test.config",NULL};
3799 
3800  pcpp::CommunicatorType nonDimenComm(testComm);
3801  plascom2::application simulationApplicationNonDimen(argc,const_cast<char **>(argv),nonDimenComm.Comm());
3802 
3803  int returnValue = 0;
3804 
3805  if(myRank == 0) {
3806  if(!ix::sys::FILEEXISTS(testDirectory)){
3807  viscousShockTestResult = false;
3808  std::cout << "TestIntegrated_viscousShock3DZ: testing directory ("
3809  << testDirectory << ") did not exist. Aborting test."
3810  << std::endl;
3811  }
3812  }
3813  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3814  if(!viscousShockTestResult)
3815  goto endTest;
3816  ix::sys::ChDir(testDirectory);
3817 
3818  // nonDimensional run using standard (rho, P, T) non-dimensionalization
3819  {
3820  ix::sys::ChDir("nonDimen");
3821 
3822  if(myRank == 0)
3823  std::cout << "Running nonDimensional PlasCom2..." << std::endl;
3824  returnValue = application::ApplicationDriver(simulationApplicationNonDimen);
3825  if(myRank == 0)
3826  std::cout << "Done running nonDimensional PlasCom2." << std::endl;
3827 
3828  if(returnValue > 0){
3829  std::cout << "TestIntegrated_ViscousShock3DZ: PlasCom2 returned error code ("
3830  << returnValue << "). Test failed." << std::endl;
3831  viscousShockTestResult = false;
3832  goto endTest;
3833  }
3834 
3835  if(myRank == 0){
3836  if(!ix::sys::FILEEXISTS("PlasCom2_000040000.h5")){
3837  std::cout << "TestIntegrated_ViscousShock3DZ: PlasCom2 failed to produce expected" << std::endl
3838  << "output file (PlasCom2_000040000.h5). Test failed." << std::endl;
3839  viscousShockTestResult = false;
3840  }
3841  }
3842  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3843  if(!viscousShockTestResult)
3844  goto endTest;
3845  }
3846 
3847  // get the domain from the application and copy the state
3848  // make a new scope for variable declaration, needed for goto usage :-(
3849  {
3850  DomainVector &simulationDomain(simulationApplicationNonDimen.GetDomains());
3851 
3852  // generate a state containing the exact solution and write it to a file
3853  int direction=2;
3854  // single grid
3855  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
3856  state_t exactState(simulationDomain[0].State(0));
3857  if(testfixtures::viscid::GenerateViscidShockExact(exactGrid,exactState,direction)){
3858  viscousShockTestResult = false;
3859  std::cout << "Fail to calculated exact viscid shock solution" << std::endl;
3860  }
3861  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3862  if(!viscousShockTestResult)
3863  goto endTest;
3864 
3865  size_t numPointsBuffer = exactGrid.BufferSize();
3866 
3867  int rhoHandle = exactState.GetDataIndex("rho");
3868  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
3869  double *rhoPtr = rhoData.Data<double>();
3870 
3871  int rhoVHandle = exactState.GetDataIndex("rhoV");
3872  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
3873  double *rhoVPtr = rhoVData.Data<double>();
3874 
3875  int rhoEHandle = exactState.GetDataIndex("rhoE");
3876  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
3877  double *rhoEPtr = rhoEData.Data<double>();
3878 
3879  int velocityHandle = exactState.GetDataIndex("velocity");
3880  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
3881  double *velocityPtr = velocityData.Data<double>();
3882 
3883  int pressureHandle = exactState.GetDataIndex("pressure");
3884  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
3885  double *pressurePtr = pressureData.Data<double>();
3886 
3887  int temperatureHandle = exactState.GetDataIndex("temperature");
3888  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
3889  double *temperaturePtr = temperatureData.Data<double>();
3890 
3891  {
3892  if(myRank == 0)
3893  std::cout << "Writing exact solution for nonDimensional PlasCom2"
3894  << std::endl;
3895  //ix::sys::ChDir("../nonDimenX");
3896  // output exact nonDimensional solution in the PlasComCM dimensionalization
3897  // we set the reference length to 1 so the grid doesn't need to change
3898  double refRho = 1.225;
3899  double refPress = 101325;
3900  double refVel = sqrt(refPress/refRho);
3901  double refTemp = 288.15;
3902  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
3903 
3904  rhoPtr[iPoint] /= refRho;
3905  rhoVPtr[iPoint+2*numPointsBuffer] /= refRho*refVel;
3906  rhoEPtr[iPoint] /= refVel*refVel*refRho;
3907 
3908  velocityPtr[iPoint+2*numPointsBuffer] /= refVel;
3909  pressurePtr[iPoint] /= refRho*refVel*refVel;
3910  temperaturePtr[iPoint] /= refTemp;
3911  }
3912 
3913  std::ostringstream Ostr;
3914  Ostr << "exact_ViscousShock_ND_40000";
3915  const std::string domainName(Ostr.str());
3916  const std::string hdfFileName(domainName+".h5");
3917  const std::string xdmfFileName(domainName+".xdmf");
3918  const std::string geometryName("shockTube");
3919  const std::string gridName("grid1");
3920  if(ix::sys::FILEEXISTS(hdfFileName)){
3921  ix::sys::Remove(hdfFileName);
3922  ix::sys::Remove(xdmfFileName);
3923  }
3924  fixtures::ConfigurationType testConfig;
3925  state_t paramState;
3926  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
3927  gridName,exactGrid,exactState,paramState,
3928  testConfig,compareStream)){
3929  std::cout << compareStream.str();
3930  std::cout << "Output to HDF5 failed:\n";
3931  std::cout << "Continuing....\n";
3932 
3933  }
3934  pcpp::io::RenewStream(compareStream);
3935  }
3936  }
3937 
3938  if(myRank == 0) {
3939  if(!ix::sys::FILEEXISTS("exact_ViscousShock_ND_40000.h5")){
3940  std::cout << "TestIntegrated_ViscousShock3DZ: failed to produce expected" << std::endl
3941  << "output file (exact_ViscousShock_ND_40000.h5). Test failed." << std::endl;
3942  viscousShockTestResult = false;
3943  }
3944  }
3945  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3946  if(!viscousShockTestResult)
3947  goto endTest;
3948 
3949  {
3950  if(myRank == 0) {
3951  std::cout << "Comparing exact solution against nonDimensional PlasCom2" << std::endl;
3952  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
3953  "exact_ViscousShock_ND_40000.h5",
3954  errorToleranceExact,compareStream);
3955  if(returnValue > 0){
3956  std::cout << "TestIntegrated_ViscousShock3DZ: Resulting state comparison yielded negative result:"
3957  << std::endl
3958  << compareStream.str()
3959  << std::endl;
3960  viscousShockTestResult = false;
3961  } else {
3962  std::cout << "Passed!" << std::endl;
3963  }
3964  pcpp::io::RenewStream(compareStream);
3965  }
3966  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
3967  }
3968 
3969  endTest:
3970  {
3971  ix::sys::ChDir(originalDirectory);
3972  parallelUnitResults.UpdateResult("Integrated:Viscid:ViscousShock3DZ",viscousShockTestResult);
3973  }
3974 
3975  return;
3976 }
3977 
3978 
3979 
3981  pcpp::CommunicatorType &testComm)
3982 {
3983 
3984  bool viscousShockTestResult = true;
3985  int myRank = testComm.Rank();
3986  int numProc = testComm.NProc();
3987 
3988 
3989  const std::string testFunctionName("TestIntegrated_VSRectilinear3DZ");
3990  const std::string testSuiteName("ViscidTest");
3991  const std::string testCaseName("ViscousShock3DZ");
3992  const std::string testDirectory(testSuiteName+"/"+testCaseName);
3993  const std::string originalDirectory(ix::sys::CWD());
3994  std::ostringstream compareStream;
3995 // double errorTolerancePC1_2 = 1.e-4;
3996  double errorToleranceExact = 5.e-2;
3997  if(myRank == 0)
3998  std::cout << "Running " << testFunctionName << std::endl;
3999 
4000  int argc = 4;
4001  const char *argv[] = {"plascom2x","-c","rectilinear.config",NULL};
4002 
4003  pcpp::CommunicatorType nonDimenComm(testComm);
4004  plascom2::application simulationApplicationNonDimen(argc,const_cast<char **>(argv),nonDimenComm.Comm());
4005 
4006  int returnValue = 0;
4007 
4008  if(myRank == 0) {
4009  if(!ix::sys::FILEEXISTS(testDirectory)){
4010  viscousShockTestResult = false;
4011  std::cout << testFunctionName << ": testing directory ("
4012  << testDirectory << ") did not exist. Aborting test."
4013  << std::endl;
4014  }
4015  }
4016  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
4017  if(!viscousShockTestResult)
4018  goto endTest;
4019  ix::sys::ChDir(testDirectory);
4020 
4021  // nonDimensional run using standard (rho, P, T) non-dimensionalization
4022  {
4023  ix::sys::ChDir("nonDimen");
4024 
4025  if(myRank == 0)
4026  std::cout << "Running nonDimensional PlasCom2..." << std::endl;
4027  returnValue = application::ApplicationDriver(simulationApplicationNonDimen);
4028  if(myRank == 0)
4029  std::cout << "Done running nonDimensional PlasCom2." << std::endl;
4030 
4031  if(returnValue > 0){
4032  std::cout << testFunctionName << ": PlasCom2 returned error code ("
4033  << returnValue << "). Test failed." << std::endl;
4034  viscousShockTestResult = false;
4035  goto endTest;
4036  }
4037 
4038  if(myRank == 0){
4039  if(!ix::sys::FILEEXISTS("PlasCom2_000040000.h5")){
4040  std::cout << testFunctionName << ": PlasCom2 failed to produce expected" << std::endl
4041  << "output file (PlasCom2_000040000.h5). Test failed." << std::endl;
4042  viscousShockTestResult = false;
4043  }
4044  }
4045  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
4046  if(!viscousShockTestResult)
4047  goto endTest;
4048  }
4049 
4050  // get the domain from the application and copy the state
4051  // make a new scope for variable declaration, needed for goto usage :-(
4052  {
4053  DomainVector &simulationDomain(simulationApplicationNonDimen.GetDomains());
4054 
4055  // generate a state containing the exact solution and write it to a file
4056  int direction=2;
4057  // single grid
4058  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
4059  state_t exactState(simulationDomain[0].State(0));
4060  if(testfixtures::viscid::GenerateViscidShockExact(exactGrid,exactState,direction)){
4061  viscousShockTestResult = false;
4062  std::cout << "Fail to calculated exact viscid shock solution" << std::endl;
4063  }
4064  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
4065  if(!viscousShockTestResult)
4066  goto endTest;
4067 
4068  size_t numPointsBuffer = exactGrid.BufferSize();
4069 
4070  int rhoHandle = exactState.GetDataIndex("rho");
4071  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
4072  double *rhoPtr = rhoData.Data<double>();
4073 
4074  int rhoVHandle = exactState.GetDataIndex("rhoV");
4075  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
4076  double *rhoVPtr = rhoVData.Data<double>();
4077 
4078  int rhoEHandle = exactState.GetDataIndex("rhoE");
4079  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
4080  double *rhoEPtr = rhoEData.Data<double>();
4081 
4082  int velocityHandle = exactState.GetDataIndex("velocity");
4083  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
4084  double *velocityPtr = velocityData.Data<double>();
4085 
4086  int pressureHandle = exactState.GetDataIndex("pressure");
4087  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
4088  double *pressurePtr = pressureData.Data<double>();
4089 
4090  int temperatureHandle = exactState.GetDataIndex("temperature");
4091  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
4092  double *temperaturePtr = temperatureData.Data<double>();
4093 
4094  {
4095  if(myRank == 0)
4096  std::cout << "Writing exact solution for nonDimensional PlasCom2"
4097  << std::endl;
4098  //ix::sys::ChDir("../nonDimenX");
4099  // output exact nonDimensional solution in the PlasComCM dimensionalization
4100  // we set the reference length to 1 so the grid doesn't need to change
4101  double refRho = 1.225;
4102  double refPress = 101325;
4103  double refVel = sqrt(refPress/refRho);
4104  double refTemp = 288.15;
4105  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
4106 
4107  rhoPtr[iPoint] /= refRho;
4108  rhoVPtr[iPoint+2*numPointsBuffer] /= refRho*refVel;
4109  rhoEPtr[iPoint] /= refVel*refVel*refRho;
4110 
4111  velocityPtr[iPoint+2*numPointsBuffer] /= refVel;
4112  pressurePtr[iPoint] /= refRho*refVel*refVel;
4113  temperaturePtr[iPoint] /= refTemp;
4114  }
4115 
4116  std::ostringstream Ostr;
4117  Ostr << "exact_ViscousShock_ND_40000";
4118  const std::string domainName(Ostr.str());
4119  const std::string hdfFileName(domainName+".h5");
4120  const std::string xdmfFileName(domainName+".xdmf");
4121  const std::string geometryName("shockTube");
4122  const std::string gridName("grid1");
4123  if(ix::sys::FILEEXISTS(hdfFileName)){
4124  ix::sys::Remove(hdfFileName);
4125  ix::sys::Remove(xdmfFileName);
4126  }
4127  fixtures::ConfigurationType testConfig;
4128  state_t paramState;
4129  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
4130  gridName,exactGrid,exactState,paramState,
4131  testConfig,compareStream)){
4132  std::cout << compareStream.str();
4133  std::cout << "Output to HDF5 failed:\n";
4134  std::cout << "Continuing....\n";
4135 
4136  }
4137  pcpp::io::RenewStream(compareStream);
4138  }
4139  }
4140 
4141  if(myRank == 0) {
4142  if(!ix::sys::FILEEXISTS("exact_ViscousShock_ND_40000.h5")){
4143  std::cout << testFunctionName << ": failed to produce expected" << std::endl
4144  << "output file (exact_ViscousShock_ND_40000.h5). Test failed." << std::endl;
4145  viscousShockTestResult = false;
4146  }
4147  }
4148  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
4149  if(!viscousShockTestResult)
4150  goto endTest;
4151 
4152  {
4153  if(myRank == 0) {
4154  std::cout << "Comparing exact solution against nonDimensional PlasCom2" << std::endl;
4155  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
4156  "exact_ViscousShock_ND_40000.h5",
4157  errorToleranceExact,compareStream);
4158  if(returnValue > 0){
4159  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
4160  << std::endl
4161  << compareStream.str()
4162  << std::endl;
4163  viscousShockTestResult = false;
4164  } else {
4165  std::cout << "Passed!" << std::endl;
4166  }
4167  pcpp::io::RenewStream(compareStream);
4168  }
4169  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
4170  }
4171 
4172  endTest:
4173  {
4174  ix::sys::ChDir(originalDirectory);
4175  parallelUnitResults.UpdateResult("Integrated:Viscid:VSRectilinear3DZ",viscousShockTestResult);
4176  }
4177 
4178  return;
4179 }
4180 
4182  pcpp::CommunicatorType &testComm)
4183 {
4184 
4185  bool viscousShockTestResult = true;
4186  int myRank = testComm.Rank();
4187  int numProc = testComm.NProc();
4188 
4189 
4190  const std::string testFunctionName("TestIntegrated_VSCurvilinear3DZ");
4191  const std::string testSuiteName("ViscidTest");
4192  const std::string testCaseName("ViscousShock3DZ");
4193  const std::string testDirectory(testSuiteName+"/"+testCaseName);
4194  const std::string originalDirectory(ix::sys::CWD());
4195  std::ostringstream compareStream;
4196 // double errorTolerancePC1_2 = 1.e-4;
4197  double errorToleranceExact = 5.e-2;
4198  if(myRank == 0)
4199  std::cout << "Running " << testFunctionName << std::endl;
4200 
4201  int argc = 4;
4202  const char *argv[] = {"plascom2x","-c","curvilinear.config",NULL};
4203 
4204  pcpp::CommunicatorType nonDimenComm(testComm);
4205  plascom2::application simulationApplicationNonDimen(argc,const_cast<char **>(argv),nonDimenComm.Comm());
4206 
4207  int returnValue = 0;
4208 
4209  if(myRank == 0) {
4210  if(!ix::sys::FILEEXISTS(testDirectory)){
4211  viscousShockTestResult = false;
4212  std::cout << testFunctionName << ": testing directory ("
4213  << testDirectory << ") did not exist. Aborting test."
4214  << std::endl;
4215  }
4216  }
4217  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
4218  if(!viscousShockTestResult)
4219  goto endTest;
4220  ix::sys::ChDir(testDirectory);
4221 
4222  // nonDimensional run using standard (rho, P, T) non-dimensionalization
4223  {
4224  ix::sys::ChDir("nonDimen");
4225 
4226  if(myRank == 0)
4227  std::cout << "Running nonDimensional PlasCom2..." << std::endl;
4228  returnValue = application::ApplicationDriver(simulationApplicationNonDimen);
4229  if(myRank == 0)
4230  std::cout << "Done running nonDimensional PlasCom2." << std::endl;
4231 
4232  if(returnValue > 0){
4233  std::cout << testFunctionName << ": PlasCom2 returned error code ("
4234  << returnValue << "). Test failed." << std::endl;
4235  viscousShockTestResult = false;
4236  goto endTest;
4237  }
4238 
4239  if(myRank == 0){
4240  if(!ix::sys::FILEEXISTS("PlasCom2_000040000.h5")){
4241  std::cout << testFunctionName << ": PlasCom2 failed to produce expected" << std::endl
4242  << "output file (PlasCom2_000040000.h5). Test failed." << std::endl;
4243  viscousShockTestResult = false;
4244  }
4245  }
4246  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
4247  if(!viscousShockTestResult)
4248  goto endTest;
4249  }
4250 
4251  // get the domain from the application and copy the state
4252  // make a new scope for variable declaration, needed for goto usage :-(
4253  {
4254  DomainVector &simulationDomain(simulationApplicationNonDimen.GetDomains());
4255 
4256  // generate a state containing the exact solution and write it to a file
4257  int direction=2;
4258  // single grid
4259  pbsgrid_t exactGrid(simulationDomain[0].Grid(0));
4260  state_t exactState(simulationDomain[0].State(0));
4261  if(testfixtures::viscid::GenerateViscidShockExact(exactGrid,exactState,direction)){
4262  viscousShockTestResult = false;
4263  std::cout << "Fail to calculated exact viscid shock solution" << std::endl;
4264  }
4265  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
4266  if(!viscousShockTestResult)
4267  goto endTest;
4268 
4269  size_t numPointsBuffer = exactGrid.BufferSize();
4270 
4271  int rhoHandle = exactState.GetDataIndex("rho");
4272  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
4273  double *rhoPtr = rhoData.Data<double>();
4274 
4275  int rhoVHandle = exactState.GetDataIndex("rhoV");
4276  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
4277  double *rhoVPtr = rhoVData.Data<double>();
4278 
4279  int rhoEHandle = exactState.GetDataIndex("rhoE");
4280  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
4281  double *rhoEPtr = rhoEData.Data<double>();
4282 
4283  int velocityHandle = exactState.GetDataIndex("velocity");
4284  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
4285  double *velocityPtr = velocityData.Data<double>();
4286 
4287  int pressureHandle = exactState.GetDataIndex("pressure");
4288  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
4289  double *pressurePtr = pressureData.Data<double>();
4290 
4291  int temperatureHandle = exactState.GetDataIndex("temperature");
4292  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
4293  double *temperaturePtr = temperatureData.Data<double>();
4294 
4295  {
4296  if(myRank == 0)
4297  std::cout << "Writing exact solution for nonDimensional PlasCom2"
4298  << std::endl;
4299  //ix::sys::ChDir("../nonDimenX");
4300  // output exact nonDimensional solution in the PlasComCM dimensionalization
4301  // we set the reference length to 1 so the grid doesn't need to change
4302  double refRho = 1.225;
4303  double refPress = 101325;
4304  double refVel = sqrt(refPress/refRho);
4305  double refTemp = 288.15;
4306  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
4307 
4308  rhoPtr[iPoint] /= refRho;
4309  rhoVPtr[iPoint+2*numPointsBuffer] /= refRho*refVel;
4310  rhoEPtr[iPoint] /= refVel*refVel*refRho;
4311 
4312  velocityPtr[iPoint+2*numPointsBuffer] /= refVel;
4313  pressurePtr[iPoint] /= refRho*refVel*refVel;
4314  temperaturePtr[iPoint] /= refTemp;
4315  }
4316 
4317  std::ostringstream Ostr;
4318  Ostr << "exact_ViscousShock_ND_40000";
4319  const std::string domainName(Ostr.str());
4320  const std::string hdfFileName(domainName+".h5");
4321  const std::string xdmfFileName(domainName+".xdmf");
4322  const std::string geometryName("shockTube");
4323  const std::string gridName("grid1");
4324  if(ix::sys::FILEEXISTS(hdfFileName)){
4325  ix::sys::Remove(hdfFileName);
4326  ix::sys::Remove(xdmfFileName);
4327  }
4328  fixtures::ConfigurationType testConfig;
4329  state_t paramState;
4330  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
4331  gridName,exactGrid,exactState,paramState,
4332  testConfig,compareStream)){
4333  std::cout << compareStream.str();
4334  std::cout << "Output to HDF5 failed:\n";
4335  std::cout << "Continuing....\n";
4336 
4337  }
4338  pcpp::io::RenewStream(compareStream);
4339  }
4340  }
4341 
4342  if(myRank == 0) {
4343  if(!ix::sys::FILEEXISTS("exact_ViscousShock_ND_40000.h5")){
4344  std::cout << testFunctionName << ": failed to produce expected" << std::endl
4345  << "output file (exact_ViscousShock_ND_40000.h5). Test failed." << std::endl;
4346  viscousShockTestResult = false;
4347  }
4348  }
4349  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
4350  if(!viscousShockTestResult)
4351  goto endTest;
4352 
4353  {
4354  if(myRank == 0) {
4355  std::cout << "Comparing exact solution against nonDimensional PlasCom2" << std::endl;
4356  returnValue = plascom2::util::PC2Compare("PlasCom2_000040000.h5",
4357  "exact_ViscousShock_ND_40000.h5",
4358  errorToleranceExact,compareStream);
4359  if(returnValue > 0){
4360  std::cout << testFunctionName << ": Resulting state comparison yielded negative result:"
4361  << std::endl
4362  << compareStream.str()
4363  << std::endl;
4364  viscousShockTestResult = false;
4365  } else {
4366  std::cout << "Passed!" << std::endl;
4367  }
4368  pcpp::io::RenewStream(compareStream);
4369  }
4370  viscousShockTestResult = CheckResult(viscousShockTestResult,testComm);
4371  }
4372 
4373  endTest:
4374  {
4375  ix::sys::ChDir(originalDirectory);
4376  parallelUnitResults.UpdateResult("Integrated:Viscid:VSCurvilinear3DZ",viscousShockTestResult);
4377  }
4378 
4379  return;
4380 }
4381 
4382 //
4383 // Generate the exact solution for the viscous shock problem
4384 // takes a grid and state and populates them based on the dimension and direction
4385 // direction = 0 (x), 1 (y), 2(z) -- must be 3D for z
4386 //
4387 int testfixtures::viscid::GenerateViscidShockExact(pbsgrid_t &exactGrid, state_t &exactState, const int direction){
4388 
4390  &partitionBufferInterval(exactGrid.PartitionBufferInterval());
4391  const pcpp::IndexIntervalType &partitionInterval(exactGrid.PartitionInterval());
4392  const std::vector<double> &dX(exactGrid.DX());
4393 
4394  std::vector<size_t> gridSizes = exactGrid.GridSizes();
4395  int numDim = gridSizes.size();
4396  const std::vector<double> &coordinateData(exactGrid.CoordinateData());
4397  size_t numPointsBuffer = exactGrid.BufferSize();
4398  const std::vector<size_t> &bufferSizes(exactGrid.BufferSizes());
4399 
4400  int rhoHandle = exactState.GetDataIndex("rho");
4401  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
4402  double *rhoPtr = rhoData.Data<double>();
4403 
4404  int rhoVHandle = exactState.GetDataIndex("rhoV");
4405  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
4406  double *rhoVPtr = rhoVData.Data<double>();
4407 
4408  int rhoEHandle = exactState.GetDataIndex("rhoE");
4409  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
4410  double *rhoEPtr = rhoEData.Data<double>();
4411 
4412  int velocityHandle = exactState.GetDataIndex("velocity");
4413  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
4414  double *velocityPtr = velocityData.Data<double>();
4415 
4416  int pressureHandle = exactState.GetDataIndex("pressure");
4417  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
4418  double *pressurePtr = pressureData.Data<double>();
4419 
4420  int temperatureHandle = exactState.GetDataIndex("temperature");
4421  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
4422  double *temperaturePtr = temperatureData.Data<double>();
4423 
4424  // clear solution data
4425  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
4426  rhoPtr[iPoint] = 0.0;
4427  rhoEPtr[iPoint] = 0.0;
4428  pressurePtr[iPoint] = 0.0;
4429  temperaturePtr[iPoint] = 0.0;
4430  }
4431 
4432  for(int iDim = 0;iDim < numDim;iDim++){
4433  size_t dimOffset = iDim*numPointsBuffer;
4434  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
4435  rhoVPtr[iPoint+dimOffset] = 0.0;
4436  velocityPtr[iPoint+dimOffset]=0.0;
4437  }
4438  }
4439 
4440  // solve for exact solution iteratively
4441  //std::cout << "Solving for exact solution" << std::endl;
4442 
4443  // some constants, we could read these from the input file
4444  // but we know what is in there, so why bother?
4445  double mach = 2.0;
4446  double p1 = 101325;
4447  double rho1 = 1.225;
4448  double T1 = 288.15;
4449  double R = p1/rho1/T1;
4450  double gamma = 1.4;
4451  double a1 = sqrt(gamma*p1/rho1);
4452  double u1 = mach*a1;
4453  double E1 = p1/(gamma-1)/rho1 + 0.5*u1*u1;
4454 
4455  // perfect gas normal shock conditions
4456  // shock stationary reference frame
4457  double pressureRatio = (2.*gamma*mach*mach-(gamma-1.))/(gamma+1.); // p2/p1
4458  double densityRatio = (gamma+1.)*mach*mach/((gamma-1.)*mach*mach+2.); // rho2/rho1
4459  double rho2 = rho1*densityRatio;
4460  double u2 = u1/densityRatio;
4461  double p2 = p1*pressureRatio;
4462  double T2 = p2/rho2/R;
4463 
4464  double mu = 20.0;
4465  double aCoeff = 3*rho1*u1*(gamma+1)/(8*mu*gamma);
4466  double bCoeff = u1/(u1-u2);
4467  double cCoeff = u2/(u1-u2);
4468 
4469  // find the transition region bounds
4470  //double x0 = 0.75;
4471  //double x0 = 0.73-0.0025;
4472  double x0 = 0.7272167;
4473  double nu1 = u1*.999;
4474  double nu2 = u2*1.001;
4475  double xMin = (cCoeff*log(nu2-u2)-bCoeff*log(u1-nu2))/aCoeff+x0;
4476  double xMax = (cCoeff*log(nu1-u2)-bCoeff*log(u1-nu1))/aCoeff+x0;
4477 
4478  size_t iStart = partitionBufferInterval[0].first;
4479  size_t iEnd = partitionBufferInterval[0].second;
4480  size_t jStart = partitionBufferInterval[1].first;
4481  size_t jEnd = partitionBufferInterval[1].second;
4482  size_t kStart = 0;
4483  size_t kEnd = 0;
4484  size_t nPlane = 0;
4485 
4486  if (numDim > 2) {
4487  kStart = partitionBufferInterval[2].first;
4488  kEnd = partitionBufferInterval[2].second;
4489  nPlane = bufferSizes[0]*bufferSizes[1];
4490  }
4491 
4492  for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
4493  size_t kBufferIndex = kIndex*nPlane;
4494  size_t kPartIndex = (kIndex - kStart) + partitionInterval[2].first;
4495  double z = kPartIndex*dX[2];
4496  for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
4497  size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes[0];
4498  size_t jPartIndex = (jIndex - jStart) + partitionInterval[1].first;
4499  double y = jPartIndex*dX[1];
4500  for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
4501  size_t bufferIndex = jkBufferIndex + iIndex;
4502  size_t iPartIndex = (iIndex - iStart) + partitionInterval[0].first;
4503  double x = iPartIndex*dX[0];
4504 
4505  double nuGuess = 0.;
4506  double position = 0.;
4507  size_t velOffSet = 0;
4508  if(direction == 0) {
4509  position = x;
4510  velOffSet = 0;
4511  } else if(direction == 1){
4512  position = y;
4513  velOffSet = numPointsBuffer;
4514  } else if(direction == 2){
4515  position = z;
4516  velOffSet = 2*numPointsBuffer;
4517  }
4518 
4519  if (position > xMax) {
4520  rhoPtr[bufferIndex] = rho1;
4521  rhoVPtr[bufferIndex+velOffSet] = -rho1*u1;
4522  rhoEPtr[bufferIndex] = p1/(gamma-1)+0.5*rho1*u1*u1;
4523  velocityPtr[bufferIndex+velOffSet] = -u1;
4524  pressurePtr[bufferIndex] = p1;
4525  temperaturePtr[bufferIndex] = T1;
4526  } else if (position < xMin) {
4527  rhoPtr[bufferIndex] = rho2;
4528  rhoVPtr[bufferIndex+velOffSet] = -rho2*u2;
4529  rhoEPtr[bufferIndex] = p2/(gamma-1)+0.5*rho2*u2*u2;
4530  velocityPtr[bufferIndex+velOffSet] = -u2;
4531  pressurePtr[bufferIndex] = p2;
4532  temperaturePtr[bufferIndex] = T2;
4533  } else {
4534  // only need to do this for one j value
4535  {
4536  // solver iteratively for u at this x in the transition region
4537  nuGuess = nu2+(nu1-nu2)/(xMax-xMin)*(position-xMin);
4538  double xGuess = (cCoeff*log(nuGuess-u2)-bCoeff*log(u1-nuGuess))/aCoeff+x0;
4539  double error = std::abs((position-xGuess)/position);
4540  double toler = 1.e-9;
4541  int count = 0;
4542  while(error > toler) {
4543  double dNudX = (u1-nuGuess)*(nuGuess-u2)/nuGuess*aCoeff;
4544  double nuGuessNew = nuGuess+dNudX*(position-xGuess);
4545 
4546  // we can't go past xMin/xMax as the function is undefined
4547  // use a linear approximation
4548  if(nuGuessNew < nu2){
4549  nuGuess = nu2+(nuGuess-nu2)/(xGuess-xMin)*(position-xMin);
4550  } else if (nuGuessNew > nu1) {
4551  nuGuess = nu1-(nu1-nuGuess)/(xMax-xGuess)*(xMax-position);
4552  } else {
4553  nuGuess = nuGuessNew;
4554  }
4555 
4556  xGuess = (cCoeff*log(nuGuess-u2)-bCoeff*log(u1-nuGuess))/aCoeff+x0;
4557  error = std::abs((position-xGuess)/position);
4558  count++;
4559  if(count > 30){
4560  std::cout << "TestIntegrated_ViscousShockPulse2D: Newton iteration failed ("
4561  << "error=" << error << "). Test failed." << std::endl;
4562  return 1;
4563  }
4564  }
4565  }
4566 
4567  rhoPtr[bufferIndex] = rho1*u1/nuGuess;
4568  rhoVPtr[bufferIndex+velOffSet] = -rhoPtr[bufferIndex]*nuGuess;
4569  velocityPtr[bufferIndex+velOffSet] = -nuGuess;
4570 
4571  // total enthalpy constant through the shock (e+p/rho+0.5*u^2)
4572  double E = (E1+p1/rho1)/gamma+0.5*nuGuess*nuGuess*(gamma-1)/gamma;
4573  rhoEPtr[bufferIndex] = E*rhoPtr[bufferIndex];
4574  pressurePtr[bufferIndex] = (gamma-1)*rhoPtr[bufferIndex]*
4575  (E - 0.5*nuGuess*nuGuess);
4576  temperaturePtr[bufferIndex] = pressurePtr[bufferIndex]/rhoPtr[bufferIndex]/R;
4577  }
4578  }
4579  }
4580  }
4581  return 0;
4582 }
4583 
4584 //
4585 // Generate the exact solution for the viscous shock problem
4586 // takes a grid and state and populates them based on the dimension and direction
4587 // direction = 0 (x), 1 (y), 2(z) -- must be 3D for z
4588 //
4589 int testfixtures::viscid::GeneratePoiseuilleExact(pbsgrid_t &exactGrid, state_t &exactState, const int direction){
4590 
4592  &partitionBufferInterval(exactGrid.PartitionBufferInterval());
4593  const pcpp::IndexIntervalType &partitionInterval(exactGrid.PartitionInterval());
4594  const std::vector<double> &dX(exactGrid.DX());
4595 
4596  std::vector<size_t> gridSizes = exactGrid.GridSizes();
4597  int numDim = gridSizes.size();
4598  const std::vector<double> &coordinateData(exactGrid.CoordinateData());
4599  const std::vector<double> &readGridExtent(exactGrid.PhysicalExtent());
4600  size_t numPointsBuffer = exactGrid.BufferSize();
4601  const std::vector<size_t> &bufferSizes(exactGrid.BufferSizes());
4602 
4603  int rhoHandle = exactState.GetDataIndex("rho");
4604  pcpp::field::dataset::DataBufferType &rhoData(exactState.Field(rhoHandle));
4605  double *rhoPtr = rhoData.Data<double>();
4606 
4607  int rhoVHandle = exactState.GetDataIndex("rhoV");
4608  pcpp::field::dataset::DataBufferType &rhoVData(exactState.Field(rhoVHandle));
4609  double *rhoVPtr = rhoVData.Data<double>();
4610 
4611  int rhoEHandle = exactState.GetDataIndex("rhoE");
4612  pcpp::field::dataset::DataBufferType &rhoEData(exactState.Field(rhoEHandle));
4613  double *rhoEPtr = rhoEData.Data<double>();
4614 
4615  int velocityHandle = exactState.GetDataIndex("velocity");
4616  pcpp::field::dataset::DataBufferType &velocityData(exactState.Field(velocityHandle));
4617  double *velocityPtr = velocityData.Data<double>();
4618 
4619  int pressureHandle = exactState.GetDataIndex("pressure");
4620  pcpp::field::dataset::DataBufferType &pressureData(exactState.Field(pressureHandle));
4621  double *pressurePtr = pressureData.Data<double>();
4622 
4623  int temperatureHandle = exactState.GetDataIndex("temperature");
4624  pcpp::field::dataset::DataBufferType &temperatureData(exactState.Field(temperatureHandle));
4625  double *temperaturePtr = temperatureData.Data<double>();
4626 
4627  // clear solution data
4628  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
4629  rhoPtr[iPoint] = 0.0;
4630  rhoEPtr[iPoint] = 0.0;
4631  pressurePtr[iPoint] = 0.0;
4632  temperaturePtr[iPoint] = 0.0;
4633  }
4634 
4635  for(int iDim = 0;iDim < numDim;iDim++){
4636  size_t dimOffset = iDim*numPointsBuffer;
4637  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
4638  rhoVPtr[iPoint+dimOffset] = 0.0;
4639  velocityPtr[iPoint+dimOffset]=0.0;
4640  }
4641  }
4642 
4643  // solve for exact solution iteratively
4644  //std::cout << "Solving for exact solution" << std::endl;
4645 
4646  // some constants, we could read these from the input file
4647  // but we know what is in there, so why bother?
4648  double mach = 2.0;
4649  double p1 = 101325;
4650  double rho1 = 1.225;
4651  double T1 = 288.15;
4652  double R = p1/rho1/T1;
4653  double gamma = 1.4;
4654  double a1 = sqrt(gamma*p1/rho1);
4655  double u1 = mach*a1;
4656  double E1 = p1/(gamma-1)/rho1 + 0.5*u1*u1;
4657 
4658 
4659  // perfect gas poiseulle exact solution
4660 
4661  // the problem dimensions, assume that flow is along the largest grid extent (2D)
4662  // in 3D, assume the flow is in the z direction
4663  double xMin = readGridExtent[0];
4664  double xMax = readGridExtent[1];
4665  double yMin = readGridExtent[2];
4666  double yMax = readGridExtent[3];
4667  double zMin = 0.;
4668  double zMax = 0.;
4669  if(numDim > 2 ) {
4670  zMin = readGridExtent[4];
4671  zMax = readGridExtent[5];
4672  }
4673 
4674  double height = 0.;
4675  double length = 0.;
4676  double width = 0.;
4677 
4678  if(numDim < 3) {
4679  if(xMax-xMin > yMax-yMin){
4680  height = yMax-yMin;
4681  length = xMax-xMin;
4682  } else {
4683  length = yMax-yMin;
4684  height = xMax-xMin;
4685  }
4686  } else {
4687  length = zMax-zMin;
4688  height = xMax-xMin;
4689  width = yMax-yMin;
4690  }
4691 
4692  double pressureRatio = 1418.55;
4693  double mu = 2.0;
4694  double pi = 3.14159265359;
4695 
4696  size_t iStart = partitionBufferInterval[0].first;
4697  size_t iEnd = partitionBufferInterval[0].second;
4698  size_t jStart = partitionBufferInterval[1].first;
4699  size_t jEnd = partitionBufferInterval[1].second;
4700  size_t kStart = 0;
4701  size_t kEnd = 0;
4702  size_t nPlane = 0;
4703 
4704  if (numDim > 2) {
4705  kStart = partitionBufferInterval[2].first;
4706  kEnd = partitionBufferInterval[2].second;
4707  nPlane = bufferSizes[0]*bufferSizes[1];
4708  }
4709 
4710  for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
4711  size_t kBufferIndex = kIndex*nPlane;
4712  size_t kPartIndex = 0;
4713  double z = 0.;
4714  if(numDim > 2 ){
4715  kPartIndex = (kIndex - kStart) + partitionInterval[2].first;
4716  z = kPartIndex*dX[2];
4717  }
4718  for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
4719  size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes[0];
4720  size_t jPartIndex = (jIndex - jStart) + partitionInterval[1].first;
4721  double y = jPartIndex*dX[1];
4722  for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
4723  size_t bufferIndex = jkBufferIndex + iIndex;
4724  size_t iPartIndex = (iIndex - iStart) + partitionInterval[0].first;
4725  double x = iPartIndex*dX[0];
4726 
4727  // x is in the flow direction
4728  // y (and z) are transverse
4729  double positionX = 0.;
4730  double positionY = 0.;
4731  double positionZ = 0.;
4732  size_t velOffSet = 0;
4733  if(direction == 0) {
4734  positionX = (xMax-x)/(xMax-xMin);
4735  positionY = y;
4736  velOffSet = 0;
4737  } else if(direction == 1){
4738  positionX = (yMax-y)/(yMax-yMin);
4739  positionY = x;
4740  velOffSet = numPointsBuffer;
4741  } else if(direction == 2){
4742  positionX = (zMax-z)/(zMax-zMin);
4743  positionY = x;
4744  positionZ = y;
4745  velOffSet = 2*numPointsBuffer;
4746  }
4747 
4748  rhoPtr[bufferIndex] = rho1;
4749  pressurePtr[bufferIndex] = p1+pressureRatio*positionX;
4750  temperaturePtr[bufferIndex] = pressurePtr[bufferIndex]/rhoPtr[bufferIndex]/R;
4751  double velocity = pressureRatio*positionY*(height-positionY)/(2*length*mu);
4752 
4753  if(numDim > 2) {
4754  // height = x
4755  // width = y
4756  // length = z
4757  double sumTerm = 0;
4758  int nMax=10;
4759  for(int n=1;n<nMax;n++){
4760  double beta = (2*n-1)*pi/height;
4761  sumTerm += sin(beta*positionY)/pow((2*n-1),3)*
4762  (sinh(beta*positionZ)+sinh(beta*(width-positionZ)))/sinh(beta*width);
4763  }
4764 
4765  sumTerm *= 4*pressureRatio/length*height*height/mu/pow(pi,3);
4766  velocity -= sumTerm;
4767  }
4768 
4769  velocityPtr[bufferIndex+velOffSet] = velocity;
4770 
4771  rhoVPtr[bufferIndex+velOffSet] = rhoPtr[bufferIndex]*velocityPtr[bufferIndex+velOffSet];
4772  rhoEPtr[bufferIndex] = pressurePtr[bufferIndex]/(gamma-1)+
4773  0.5*rhoPtr[bufferIndex]*velocity*velocity;
4774  }
4775  }
4776  }
4777  return 0;
4778 }
4779 
void TestIntegrated_PFCurvilinear3DZ(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
void TestIntegrated_VSRectilinear2DX(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
void TestIntegrated_Poiseuille2DY(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
int ApplicationDriver(ApplicationType &simulationApplication)
void TestIntegrated_ViscousShock2DX(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
int Remove(const std::string &fname)
Definition: UnixUtils.C:141
int GetDataIndex(const std::string &name) const
void TestIntegrated_PFRectilinear2DY(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
std::vector< DomainBaseType > domainvector
Definition: Simulation.H:32
void TestIntegrated_VSCurvilinear2DY(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
const std::vector< double > & DX() const
Definition: Grid.H:538
pcpp::IndexIntervalType & PartitionInterval()
Definition: Grid.H:500
void TestIntegrated_VSRectilinear2DY(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
void const size_t const size_t * gridSizes
Definition: EulerKernels.H:10
int PC2Compare(const std::string &redFileName, const std::string &blueFileName, double errTolerance, std::ostream &outStream)
Read two HDF5 files and compare the state data therein.
Definition: PC2Util.C:30
void RenewStream(std::ostringstream &outStream)
bool CheckResult(bool &localResult, pcpp::CommunicatorType &testComm)
Definition: PCPPCommUtil.C:176
int ChDir(const std::string &path)
Definition: UnixUtils.C:297
void const size_t const size_t const size_t const double const double double * y
simulation::grid::parallel_blockstructured pbsgrid_t
void const size_t const size_t const size_t const int const double const int const double * velocity
Definition: MetricKernels.H:19
void TestIntegrated_PFRectilinear3DZ(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
int GeneratePoiseuilleExact(pbsgrid_t &exactGrid, state_t &exactState, const int direction)
int OutputSingle(const std::string &fileName, const GridType &inGrid, const StateType &inState, const ConfigType &simConfig, const ConfigType &gridConfig, const ConfigType &stateConfig)
Definition: PlasCom2IO.C:12
void const size_t const size_t const size_t const double const double * x
const std::vector< size_t > & BufferSizes() const
Definition: Grid.H:459
Encapsulating class for collections of test results.
Definition: Testing.H:18
pcpp::IndexIntervalType & PartitionBufferInterval()
Definition: Grid.H:519
const std::vector< size_t > & GridSizes() const
Definition: Grid.H:469
bool FILEEXISTS(const std::string &fname)
Definition: UnixUtils.C:189
Main encapsulation of MPI.
Definition: COMM.H:62
std::string Grid(const GridType &inGrid)
Definition: Report.H:14
Testing constructs for unit testing.
simulation::state::base state_t
void TestIntegrated_VSCurvilinear3DZ(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
void TestIntegrated_Poiseuille3DZ(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
void TestIntegrated_ViscousShock3DZ(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
void const size_t const size_t * bufferSizes
Definition: MetricKernels.H:19
const std::string CWD()
Definition: UnixUtils.C:291
void UpdateResult(const std::string &name, const ValueType &result)
Updates an existing test result.
Definition: Testing.H:55
const std::vector< double > & PhysicalExtent() const
Definition: Grid.H:533
std::vector< double > & CoordinateData()
Definition: Grid.H:503
void TestIntegrated_ViscousShock2DY(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
Definition: EulerRHS.H:33
plascom2::application::domainvector DomainVector
void TestIntegrated_Poiseuille2DX(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
void TestIntegrated_PFCurvilinear2DX(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
int GenerateViscidShockExact(pbsgrid_t &exactGrid, state_t &exactState, const int direction)
Simple Block Structured Mesh object.
Definition: IndexUtil.H:21
void TestIntegrated_PFCurvilinear2DY(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
void TestIntegrated_VSRectilinear3DZ(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
void const size_t * numPointsBuffer
Definition: MetricKernels.H:19
void TestIntegrated_VSCurvilinear2DX(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
DataBufferType & Field(const std::string &name)
void TestIntegrated_PFRectilinear2DX(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)