PlasCom2  1.0
XPACC Multi-physics simluation application
TestGrid.C
Go to the documentation of this file.
1 #include "Testing.H"
2 #include "Simulation.H"
3 #include "OperatorKernels.H"
4 #include "Stencil.H"
5 #include "PCPPCommUtil.H"
6 #include "PCPPReport.H"
7 #include "PCPPIntervalUtils.H"
8 #include "PCPPHDF5.H"
9 #include "Report.H"
10 #include "PC2IO.H"
11 #include "GridGen.H"
12 
13 // Check result is used to check test results across
14 // all processors in parallel tests.
15 using pcpp::comm::CheckResult;
16 
17 
18 void TestGrid_SubRegion(ix::test::results &serialUnitResults)
19 {
20 
21  std::ostringstream messageStream;
22  simulation::grid::subregion subRegion;
23 
24  std::vector<size_t> gridSizes(3,100);
25  pcpp::IndexIntervalType gridInterval;
26  gridInterval.InitSimple(gridSizes);
27 
28  pcpp::IndexIntervalType &subRegionInterval(subRegion.regionInterval);
29 
30  bool subregionTest = true;
31 
32  std::string initString1("testRegion 1");
33  if(simulation::grid::InitSubRegionFromString(initString1,gridSizes,subRegion))
34  subregionTest = false;
35 
36  if(subregionTest){
37  if(subRegion.regionName != "testRegion")
38  subregionTest = false;
39  if(subRegion.normalDirection != 1)
40  subregionTest = false;
41  if(subRegionInterval[0].first != 0 &&
42  subRegionInterval[0].second != 0)
43  subregionTest = false;
44  if(subRegionInterval[1] != gridInterval[1])
45  subregionTest = false;
46  if(subRegionInterval[2] != gridInterval[2])
47  subregionTest = false;
48  }
49 
50  std::string initString0("testRegion 0");
51  if(simulation::grid::InitSubRegionFromString(initString0,gridSizes,subRegion))
52  subregionTest = false;
53 
54  if(subregionTest){
55  if(subRegion.regionName != "testRegion")
56  subregionTest = false;
57  if(subRegion.normalDirection != 0)
58  subregionTest = false;
59  if(subRegionInterval != gridInterval)
60  subregionTest = false;
61  }
62 
63  std::string initString_1("testRegion -1");
64  if(simulation::grid::InitSubRegionFromString(initString_1,gridSizes,subRegion))
65  subregionTest = false;
66 
67  if(subregionTest){
68  if(subRegion.regionName != "testRegion")
69  subregionTest = false;
70  if(subRegion.normalDirection != -1)
71  subregionTest = false;
72  if(subRegionInterval[0].first != gridInterval[0].second)
73  subregionTest = false;
74  if(subRegionInterval[0].second != gridInterval[0].second)
75  subregionTest = false;
76  if(subRegionInterval[1] != gridInterval[1])
77  subregionTest = false;
78  if(subRegionInterval[2] != gridInterval[2])
79  subregionTest = false;
80  }
81 
82 
83  std::string initString12("testRegion 1 1 4");
84  if(simulation::grid::InitSubRegionFromString(initString12,gridSizes,subRegion))
85  subregionTest = false;
86 
87 
88  if(subregionTest){
89  if(subRegion.regionName != "testRegion")
90  subregionTest = false;
91  if(subRegion.normalDirection != 1)
92  subregionTest = false;
93  if(subRegionInterval[0].first != 0 &&
94  subRegionInterval[0].second != 3)
95  subregionTest = false;
96  if(subRegionInterval[1] != gridInterval[1])
97  subregionTest = false;
98  if(subRegionInterval[2] != gridInterval[2])
99  subregionTest = false;
100  }
101 
102  std::string initString_12("testRegion -1 -4 -1");
103  if(simulation::grid::InitSubRegionFromString(initString_12,gridSizes,subRegion))
104  subregionTest = false;
105 
106  if(subregionTest){
107  if(subRegion.regionName != "testRegion")
108  subregionTest = false;
109  if(subRegion.normalDirection != -1)
110  subregionTest = false;
111  if(subRegionInterval[0].first != (gridInterval[0].second-3))
112  subregionTest = false;
113  if(subRegionInterval[0].second != gridInterval[0].second)
114  subregionTest = false;
115  if(subRegionInterval[1] != gridInterval[1])
116  subregionTest = false;
117  if(subRegionInterval[2] != gridInterval[2])
118  subregionTest = false;
119  }
120 
121 
122  std::string initString13("testRegion 1 1 4 10 15");
123  if(simulation::grid::InitSubRegionFromString(initString13,gridSizes,subRegion))
124  subregionTest = false;
125 
126 
127  if(subregionTest){
128  if(subRegion.regionName != "testRegion")
129  subregionTest = false;
130  if(subRegion.normalDirection != 1)
131  subregionTest = false;
132  if(subRegionInterval[0].first != 0 &&
133  subRegionInterval[0].second != 3)
134  subregionTest = false;
135  if(subRegionInterval[1].first != 9)
136  subregionTest = false;
137  if(subRegionInterval[1].second != 14)
138  subregionTest = false;
139  if(subRegionInterval[2] != gridInterval[2])
140  subregionTest = false;
141  }
142 
143  std::string initString14("testRegion 1 1 4 10 15 20 23");
144  if(simulation::grid::InitSubRegionFromString(initString14,gridSizes,subRegion))
145  subregionTest = false;
146 
147 
148  if(subregionTest){
149  if(subRegion.regionName != "testRegion")
150  subregionTest = false;
151  if(subRegion.normalDirection != 1)
152  subregionTest = false;
153  if(subRegionInterval[0].first != 0 &&
154  subRegionInterval[0].second != 3)
155  subregionTest = false;
156  if(subRegionInterval[1].first != 9)
157  subregionTest = false;
158  if(subRegionInterval[1].second != 14)
159  subregionTest = false;
160  if(subRegionInterval[2].first != 19)
161  subregionTest = false;
162  if(subRegionInterval[2].second != 22)
163  subregionTest = false;
164  }
165 
166  std::string initString_14("testRegion -1 -4 -1 1 10 2 20");
167  if(simulation::grid::InitSubRegionFromString(initString_14,gridSizes,subRegion))
168  subregionTest = false;
169 
170 
171  if(subregionTest){
172  if(subRegion.regionName != "testRegion")
173  subregionTest = false;
174  if(subRegion.normalDirection != -1)
175  subregionTest = false;
176  if(subRegionInterval[0].first != (gridInterval[0].second-3) &&
177  subRegionInterval[0].second != gridInterval[0].second)
178  subregionTest = false;
179  if(subRegionInterval[1].first != 0)
180  subregionTest = false;
181  if(subRegionInterval[1].second != 9)
182  subregionTest = false;
183  if(subRegionInterval[2].first != 1)
184  subregionTest = false;
185  if(subRegionInterval[2].second != 19)
186  subregionTest = false;
187  }
188 
189  std::string initString_24("testRegion -2 1 4 -10 -1 2 20");
190  if(simulation::grid::InitSubRegionFromString(initString_24,gridSizes,subRegion))
191  subregionTest = false;
192 
193 
194  if(subregionTest){
195  if(subRegion.regionName != "testRegion")
196  subregionTest = false;
197  if(subRegion.normalDirection != -2)
198  subregionTest = false;
199  if(subRegionInterval[0].first != 0 &&
200  subRegionInterval[0].second != 3)
201  subregionTest = false;
202  if(subRegionInterval[1].first != (gridInterval[1].second - 9))
203  subregionTest = false;
204  if(subRegionInterval[1].second != gridInterval[1].second)
205  subregionTest = false;
206  if(subRegionInterval[2].first != 1)
207  subregionTest = false;
208  if(subRegionInterval[2].second != 19)
209  subregionTest = false;
210  }
211 
212  std::string initString02("testRegion 0 1 -1 1 10 1 -1");
213  if(simulation::grid::InitSubRegionFromString(initString02,gridSizes,subRegion))
214  subregionTest = false;
215  if(subRegion.normalDirection != 0)
216  subregionTest = false;
217  if(subRegionInterval[0] != gridInterval[0])
218  subregionTest = false;
219  if(subRegionInterval[1].first != 0 &&
220  subRegionInterval[1].second != 9)
221  subregionTest = false;
222  if(subRegionInterval[2] != gridInterval[2])
223  subregionTest = false;
224 
225  serialUnitResults.UpdateResult("Grid:SubRegion:InitSubRegionFromString",subregionTest);
226 }
227 
228 
229 void TestGrid_CartesianMetric(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
230 {
231  std::ostringstream messageStream;
232 
233  fixtures::ConfigurationType testConfig;
234 
235  typedef simulation::grid::parallel_blockstructured pbsgrid_t;
236  typedef plascom2::operators::sbp::base operator_t;
237  typedef simulation::grid::halo halo_t;
238  typedef pcpp::ParallelGlobalType global_t;
239  typedef pcpp::CommunicatorType comm_t;
240  typedef pcpp::IndexIntervalType interval_t;
241 
242  global_t myGlobal;
243  myGlobal.Init("TestGrid_CartesianMetric",testComm);
244  myGlobal.SetVerbLevel(3);
245  myGlobal.Profiling(true);
246 
247  std::vector<bool> testResults(2,true);
248  std::vector<size_t> allSizes(3,0);
249 
250  // Test metrics for grids of this size
251  allSizes[0] = 101;
252  allSizes[1] = 21;
253  allSizes[2] = 11;
254 
255  for(int numDim = 2;numDim < 4;numDim++){
256 
257  bool sectionResult = true;
258 
259  pbsgrid_t testGrid;
260  int haloDepth = 2;
261  pcpp::ParallelTopologyInfoType cartInfo;
262 
263  // Make a non-periodic cartesian grid & topology
264  cartInfo.numDimensions = numDim;
265  cartInfo.isPeriodic.resize(numDim,0);
266 
267  // Set grid sizes
268  size_t numGlobalNodes = 1;
269  size_t numGlobalCells = 1;
270  std::vector<size_t> gridSizes(numDim,0);
271  for(int iDim = 0;iDim < numDim;iDim++){
272  gridSizes[iDim] = allSizes[iDim];
273  numGlobalNodes *= gridSizes[iDim];
274  numGlobalCells *= (gridSizes[iDim]-1);
275  }
276  testGrid.SetGridSizes(gridSizes);
277 
278  messageStream << "Setting up grid with sizes: (";
279  pcpp::io::DumpContents(messageStream,gridSizes,",");
280  messageStream << ")" << std::endl;
281  myGlobal.StdOut(messageStream.str());
282  pcpp::io::RenewStream(messageStream);
283 
284 
285  // This call does everything to set up the grid and required halos
286  messageStream << "====== Grid Parallel Setup ========= " << std::endl;
287  if(testGrid.ParallelSetup(testComm,cartInfo,haloDepth,messageStream)){
288  myGlobal.ErrOut("Grid::ParallelSetup failed with messages:\n");
289  myGlobal.ErrOut(messageStream.str());
290  pcpp::io::RenewStream(messageStream);
291  sectionResult = false;
292  }
293  sectionResult = CheckResult(sectionResult,testComm);
294  if(!sectionResult){
295  testResults[numDim-2] = false;
296  continue;
297  }
298  messageStream << "==================================== " << std::endl;
299 
300  // Generate a report of the Cartesian grid and topo setup
301  messageStream << "------- Grid & Topo Report --------- " << std::endl;
302  pcpp::report::CartesianSetup(messageStream,cartInfo);
303  messageStream << simulation::report::Grid(testGrid) << std::endl
304  << "------------------------------------ " << std::endl;
305  myGlobal.StdOut(messageStream.str());
306  pcpp::io::RenewStream(messageStream);
307 
308  // Set up grid type and generate coordinates
309  std::vector<double> physicalGridExtent(2*numDim,0);
310  std::vector<double> dX(numDim,0);
311  double expectedJacobian = 1.0;
312  std::vector<double> expectedMetric(numDim,0);
313  for(int iDim = 0;iDim < numDim;iDim++){
314  physicalGridExtent[iDim*2+1] = 1.0;
315  dX[iDim] = 1.0/(double(gridSizes[iDim]-1));
316  expectedMetric[iDim] = 1.0/dX[iDim];
317  expectedJacobian *= expectedMetric[iDim];
318  }
319  for(int iDim = 0;iDim < numDim;iDim++){
320  expectedMetric[iDim] *= 1.0/expectedJacobian;
321  }
322 
323  testGrid.SetPhysicalExtent(physicalGridExtent);
324 
325  if(testGrid.GenerateCoordinates(messageStream)){
326  myGlobal.StdOut("Coordinate generation failed with messages:\n");
327  myGlobal.StdOut(messageStream.str());
328  myGlobal.StdOut("\nSkipping rest of this test.\n");
329  pcpp::io::RenewStream(messageStream);
330  sectionResult = false;
331  }
332  sectionResult = CheckResult(sectionResult,testComm);
333  if(!sectionResult){
334  testResults[numDim-2] = false;
335  continue;
336  }
337 
338  myGlobal.StdOut(messageStream.str());
339  pcpp::io::RenewStream(messageStream);
340 
341  if(testGrid.ComputeMetrics(messageStream)){
342  myGlobal.StdOut("Metrics generation failed with messages:\n");
343  myGlobal.StdOut(messageStream.str());
344  myGlobal.StdOut("\nSkipping rest of this test.\n");
345  pcpp::io::RenewStream(messageStream);
346  sectionResult = false;
347  }
348  sectionResult = CheckResult(sectionResult,testComm);
349  if(!sectionResult){
350  testResults[numDim-2] = false;
351  continue;
352  }
353 
354  myGlobal.StdOut(messageStream.str());
355  pcpp::io::RenewStream(messageStream);
356 
357  std::vector<double> &gridMetric(testGrid.Metric());
358  std::vector<double> &gridJacobian(testGrid.Jacobian());
359 
360  if(gridMetric != expectedMetric){
361  messageStream << "Grid metric not as expected. (";
362  pcpp::io::DumpContents(messageStream,gridMetric,",");
363  messageStream << ") != (";
364  pcpp::io::DumpContents(messageStream,expectedMetric,",");
365  messageStream << "). Test failed for " << numDim << " dimensions."
366  << std::endl;
367  myGlobal.StdOut(messageStream.str());
368  pcpp::io::RenewStream(messageStream);
369  sectionResult = false;
370  }
371  sectionResult = CheckResult(sectionResult,testComm);
372  if(!sectionResult){
373  testResults[numDim-2] = false;
374  continue;
375  }
376 
377  if(gridJacobian.size() != 2){
378  myGlobal.StdOut("Grid Jacobian not of proper size.");
379  sectionResult = false;
380  }
381  sectionResult = CheckResult(sectionResult,testComm);
382  if(!sectionResult){
383  testResults[numDim-2] = false;
384  continue;
385  }
386 
387  if(gridJacobian[0] != expectedJacobian){
388  messageStream << "Grid Jacobian does not have expected value ("
389  << gridJacobian[0] << " != " << expectedJacobian
390  << "). Test failed." << std::endl;
391  myGlobal.StdOut(messageStream.str());
392  sectionResult = false;
393  }
394  if(gridJacobian[1] != 1.0/expectedJacobian){
395  messageStream << "Grid Jacobian does not have expected value ("
396  << gridJacobian[1] << " != " << 1.0/expectedJacobian
397  << "). Test failed." << std::endl;
398  myGlobal.StdOut(messageStream.str());
399  sectionResult = false;
400  }
401 
402  sectionResult = CheckResult(sectionResult,testComm);
403  if(!sectionResult){
404  testResults[numDim-2] = false;
405  continue;
406  }
407 
408  }
409 
410  parallelUnitResults.UpdateResult("Grid:PBS:CartesianMetric2D",testResults[0]);
411  parallelUnitResults.UpdateResult("Grid:PBS:CartesianMetric3D",testResults[1]);
412 
413 };
414 
415 
416 template<int numDim>
417 int GenerateRectilinearMetric1(std::vector<size_t> &ijk,std::vector<double> &xyz,double &jacm1)
418 {
419  double scale[] = {101,51,21};
420  double dxDXi = 2.0*(ijk[0]+1)/(scale[0]*scale[0]);
421  double dyDEta = 4.0*(ijk[1]+1)/(scale[1]*scale[1]);
422  xyz[0] = dyDEta;
423  xyz[1] = dxDXi;
424  jacm1 = xyz[0]*xyz[1];
425  if (numDim == 3){
426  double dzDZeta = 6.0*(ijk[2]+1)/(scale[2]*scale[2]);
427  xyz[0] = dyDEta*dzDZeta;
428  xyz[1] = dxDXi*dzDZeta;
429  xyz[2] = dxDXi*dyDEta;
430  jacm1 = xyz[0]*xyz[1]*xyz[2];
431  }
432  return(0);
433 };
434 
435 void TestGrid_RectilinearMetric(ix::test::results &parallelUnitResults,
436  pcpp::CommunicatorType &testComm)
437 {
438  std::ostringstream messageStream;
439 
440  fixtures::ConfigurationType testConfig;
441 
442  typedef simulation::grid::parallel_blockstructured pbsgrid_t;
443  typedef plascom2::operators::sbp::base operator_t;
444  typedef simulation::grid::halo halo_t;
445  typedef pcpp::ParallelGlobalType global_t;
446  typedef pcpp::CommunicatorType comm_t;
447  typedef pcpp::IndexIntervalType interval_t;
448 
449  global_t myGlobal;
450  myGlobal.Init("TestGrid_RectilinearMetric",testComm);
451  myGlobal.SetVerbLevel(3);
452  myGlobal.Profiling(true);
453 
454  std::vector<bool> testResults(2,true);
455  std::vector<size_t> allSizes(3,0);
456  std::vector<double> scale(3,0);
457 
458  // Test metrics for grids of this size
459  allSizes[0] = 101;
460  allSizes[1] = 51;
461  allSizes[2] = 21;
462 
463  for(int numDim = 2;numDim < 4;numDim++){
464 
465  bool sectionResults = true;
466 
467  pbsgrid_t testGrid;
468  int interiorSBPOrder = 4;
469  int haloDepth = interiorSBPOrder+1;
470  pcpp::ParallelTopologyInfoType cartInfo;
471 
472 
473  // Make a non-periodic cartesian grid & topology
474  cartInfo.numDimensions = numDim;
475  cartInfo.isPeriodic.resize(numDim,0);
476 
477  // Set grid sizes
478  size_t numGlobalNodes = 1;
479  size_t numGlobalCells = 1;
480  std::vector<size_t> gridSizes(numDim,0);
481  std::vector<double> gridScale(numDim,1.0);
482  for(int iDim = 0;iDim < numDim;iDim++){
483  gridSizes[iDim] = allSizes[iDim];
484  gridScale[iDim] = double(allSizes[iDim]);
485  numGlobalNodes *= gridSizes[iDim];
486  numGlobalCells *= (gridSizes[iDim]-1);
487  }
488  testGrid.SetGridSizes(gridSizes);
489 
490 
491  messageStream << "Setting up grid with sizes: (";
492  pcpp::io::DumpContents(messageStream,gridSizes,",");
493  messageStream << ")" << std::endl;
494  myGlobal.StdOut(messageStream.str());
495  pcpp::io::RenewStream(messageStream);
496 
497 
498  // This call does everything to set up the grid and required halos
499  messageStream << "====== Grid Parallel Setup ========= " << std::endl;
500  if(testGrid.ParallelSetup(testComm,cartInfo,haloDepth,messageStream)){
501  myGlobal.ErrOut("Grid::ParallelSetup failed with messages:\n");
502  myGlobal.ErrOut(messageStream.str());
503  pcpp::io::RenewStream(messageStream);
504  sectionResults = false;
505  }
506  sectionResults = CheckResult(sectionResults,testComm);
507  if(!sectionResults){
508  testResults[numDim-2] = false;
509  continue;
510  }
511  messageStream << "==================================== " << std::endl;
512 
513  interval_t &partitionBufferInterval(testGrid.PartitionBufferInterval());
514  interval_t &partitionInterval(testGrid.PartitionInterval());
515  std::vector<size_t> partitionSizes(partitionInterval.Sizes());
516  std::vector<size_t> partitionStarts(partitionInterval.Starts());
517  std::vector<size_t> partitionBufferStarts(partitionBufferInterval.Starts());
518 
519  // Generate a report of the Cartesian grid and topo setup
520  messageStream << "------- Grid & Topo Report --------- " << std::endl;
521  pcpp::report::CartesianSetup(messageStream,cartInfo);
522  messageStream << simulation::report::Grid(testGrid) << std::endl;
523  messageStream << "------------------------------------ " << std::endl;
524  myGlobal.StdOut(messageStream.str());
525  pcpp::io::RenewStream(messageStream);
526 
527  // Set up grid coordinates and expected metric
528  size_t numPointsBuffer = testGrid.BufferSize();
529 
530  std::vector<double> expectedMetric(numDim*numPointsBuffer,0);
531  std::vector<double> expectedJacobian(2*numPointsBuffer,0);
532  std::vector<double> pointMetric(numDim,0);
533  std::vector<size_t> pointIndices(numDim,0);
534 
535  testGrid.SetType(simulation::grid::RECTILINEAR);
536 
537  std::ostringstream topoTypeStream;
538  topoTypeStream << numDim << "DSMesh";
539 
540  std::vector<size_t> &bufferSizes(testGrid.BufferSizes());
541  std::vector<size_t> rGridSizes(testGrid.GridSizes());
542  std::reverse(rGridSizes.begin(),rGridSizes.end());
543 
544  if(numDim == 2) {
545 
546  // testGrid.GenerateCoordinates(GenerateRectilinearGrid1<2>,messageStream);
547  testGrid.GenerateCoordinates(gridgen::RectilinearGrid1<2,101,51>,messageStream);
548  testGrid.ExchangeCoordinates();
549 
550  size_t nPlane = bufferSizes[0];
551  for(size_t iY = 0;iY < partitionSizes[1];iY++){
552  pointIndices[1] = iY + partitionStarts[1];
553  size_t bufferY = iY + partitionBufferStarts[1];
554  for(size_t iX = 0;iX < partitionSizes[0];iX++){
555  pointIndices[0] = iX + partitionStarts[0];
556  size_t bufferX = iX + partitionBufferStarts[0];
557  size_t bufferIndex = bufferY*nPlane + bufferX;
558  pointMetric[0] = 0;
559  pointMetric[1] = 0;
560  GenerateRectilinearMetric1<2>(pointIndices,pointMetric,
561  expectedJacobian[bufferIndex+numPointsBuffer]);
562  expectedMetric[bufferIndex] = pointMetric[0];
563  expectedMetric[bufferIndex+numPointsBuffer] = pointMetric[1];
564  }
565  }
566 
567 
568  // ====== write the grid to file, just so we can look at it ====
569  if(ix::sys::FILEEXISTS("rectilinear2d.h5")){
570  ix::sys::Remove("rectilinear2d.h5");
571  ix::sys::Remove("rectilinear2d.xdmf");
572  }
573  pcpp::io::hdf5::base hdf5File("rectilinear2d.h5",testComm);
574  pcpp::io::hdf5::WriteGrid(testGrid,"rectilinear2d","",hdf5File);
575  hdf5File.Close();
576  std::ofstream xdmfStream;
577  xdmfStream.open("rectilinear2d.xdmf");
578  pcpp::io::xdmf::OpenFileTag(xdmfStream);
579  pcpp::io::xdmf::OpenGridTag("rectilinear2d","Uniform",0,xdmfStream);
580  pcpp::io::xdmf::WriteGridSection(topoTypeStream.str(),"rectilinear2d.h5","rectilinear2d",
581  rGridSizes,xdmfStream);
582  pcpp::io::xdmf::CloseGridTag(xdmfStream);
583  pcpp::io::xdmf::CloseFileTag(xdmfStream);
584  xdmfStream.close();
585  // ====== ---------------------- ====
586 
587  } else if(numDim == 3) {
588 
589  testGrid.GenerateCoordinates(gridgen::RectilinearGrid1<3,101,51,21>,messageStream);
590  testGrid.ExchangeCoordinates();
591 
592  size_t nPlane = bufferSizes[0]*bufferSizes[1];
593  for(size_t iZ = 0;iZ < partitionSizes[2];iZ++){
594  pointIndices[2] = iZ + partitionStarts[2];
595  size_t bufferZ = iZ + partitionBufferStarts[2];
596  for(size_t iY = 0;iY < partitionSizes[1];iY++){
597  pointIndices[1] = iY + partitionStarts[1];
598  size_t bufferY = iY + partitionBufferStarts[1];
599  for(size_t iX = 0;iX < partitionSizes[0];iX++){
600  pointIndices[0] = iX + partitionStarts[0];
601  size_t bufferX = iX + partitionBufferStarts[0];
602  size_t bufferIndex = bufferZ*nPlane + bufferY*bufferSizes[0] + bufferX;
603  pointMetric[0] = 0;
604  pointMetric[1] = 0;
605  GenerateRectilinearMetric1<3>(pointIndices,pointMetric,
606  expectedJacobian[bufferIndex+numPointsBuffer]);
607  expectedMetric[bufferIndex] = pointMetric[0];
608  expectedMetric[bufferIndex+numPointsBuffer] = pointMetric[1];
609  expectedMetric[bufferIndex+2*numPointsBuffer] = pointMetric[2];
610  }
611  }
612  }
613 
614  // ====== write the grid to file, just so we can look at it ====
615  if(ix::sys::FILEEXISTS("rectilinear3d.h5")){
616  ix::sys::Remove("rectilinear3d.h5");
617  ix::sys::Remove("rectilinear3d.xdmf");
618  }
619  pcpp::io::hdf5::base hdf5File("rectilinear3d.h5",testComm);
620  pcpp::io::hdf5::WriteGrid(testGrid,"rectilinear3d","",hdf5File);
621  hdf5File.Close();
622  std::ofstream xdmfStream;
623  xdmfStream.open("rectilinear3d.xdmf");
624  pcpp::io::xdmf::OpenFileTag(xdmfStream);
625  pcpp::io::xdmf::OpenGridTag("rectilinear3d","Uniform",0,xdmfStream);
626  pcpp::io::xdmf::WriteGridSection(topoTypeStream.str(),"rectilinear3d.h5","rectilinear3d",
627  rGridSizes,xdmfStream);
628  pcpp::io::xdmf::CloseGridTag(xdmfStream);
629  pcpp::io::xdmf::CloseFileTag(xdmfStream);
630  xdmfStream.close();
631  // ====== ---------------------- ====
632  }
633 
634 
635 
636  myGlobal.StdOut(messageStream.str());
637  pcpp::io::RenewStream(messageStream);
638 
639  // Set up differential operator for grid
640  operator_t sbpOperator;
641  plascom2::operators::sbp::Initialize(sbpOperator,interiorSBPOrder);
642  int connBoundaryDepth = (sbpOperator.numStencils-2)/2;
643  // Forticate the stencil starts
644  for(int iStencil = 0;iStencil < sbpOperator.numStencils;iStencil++)
645  sbpOperator.stencilStarts[iStencil]++;
646 
647  std::vector<int> stencilConnectivity(numDim*numPointsBuffer,0);
648  std::vector<size_t> opInterval;
649  partitionBufferInterval.Flatten(opInterval);
650  std::vector<size_t>::iterator dOpIt = opInterval.begin();
651  while(dOpIt != opInterval.end()){
652  *dOpIt += 1;
653  dOpIt++;
654  }
655 
656  if(plascom2::operators::sbp::CreateStencilConnectivity(numDim,&bufferSizes[0],&opInterval[0],
657  connBoundaryDepth,&stencilConnectivity[0],true)){
658  messageStream << "CreateStencilConnectivity failed. Aborting the test." << std::endl;
659  myGlobal.StdOut(messageStream.str());
660  pcpp::io::RenewStream(messageStream);
661  sectionResults = false;
662  }
663  sectionResults = CheckResult(sectionResults,testComm);
664  if(!sectionResults){
665  testResults[numDim-2] = false;
666  continue;
667  }
668 
669  testGrid.SetupDifferentialOperator(&sbpOperator,&stencilConnectivity[0]);
670 
671  if(testGrid.ComputeMetrics(messageStream)){
672  myGlobal.StdOut("Metrics generation failed with messages:\n");
673  myGlobal.StdOut(messageStream.str());
674  myGlobal.StdOut("\nSkipping rest of this test.\n");
675  pcpp::io::RenewStream(messageStream);
676  sectionResults = false;
677  }
678  sectionResults = CheckResult(sectionResults,testComm);
679  if(!sectionResults){
680  testResults[numDim-2] = false;
681  continue;
682  }
683 
684  myGlobal.StdOut(messageStream.str());
685  pcpp::io::RenewStream(messageStream);
686 
687  std::vector<double> &gridMetric(testGrid.Metric());
688  std::vector<double> &gridJacobian(testGrid.Jacobian());
689 
690 
691  if(numDim == 2){
692 
693  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
694 
695  size_t iPoint2 = iPoint + numPointsBuffer;
696  double error11 = std::abs(gridMetric[iPoint] - expectedMetric[iPoint]);
697  double error22 = std::abs(gridMetric[iPoint2] - expectedMetric[iPoint2]);
698 
699  if(error11 > 1e-12){
700  sectionResults = false;
701  messageStream << "gridMetric11[" << iPoint << "] = "
702  << gridMetric[iPoint] << "," << expectedMetric[iPoint]
703  << std::endl;
704  }
705  if(error22 > 1e-12){
706  sectionResults = false;
707  messageStream << "gridMetric22[" << iPoint << "] = "
708  << gridMetric[iPoint2] << "," << expectedMetric[iPoint2]
709  << std::endl;
710  }
711  }
712  sectionResults = CheckResult(sectionResults,testComm);
713  } else if (numDim == 3){
714 
715  for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
716 
717  size_t iPoint2 = iPoint + numPointsBuffer;
718  size_t iPoint3 = iPoint + 2*numPointsBuffer;
719  double error11 = std::abs(gridMetric[iPoint] - expectedMetric[iPoint]);
720  double error22 = std::abs(gridMetric[iPoint2] - expectedMetric[iPoint2]);
721  double error33 = std::abs(gridMetric[iPoint3] - expectedMetric[iPoint3]);
722 
723  if(error11 > 1e-12){
724  sectionResults = false;
725  messageStream << "gridMetric11[" << iPoint << "] = "
726  << gridMetric[iPoint] << "," << expectedMetric[iPoint]
727  << std::endl;
728  }
729  if(error22 > 1e-12){
730  sectionResults = false;
731  messageStream << "gridMetric22[" << iPoint << "] = "
732  << gridMetric[iPoint2] << "," << expectedMetric[iPoint2]
733  << std::endl;
734  }
735  if(error33 > 1e-12){
736  sectionResults = false;
737  messageStream << "gridMetric33[" << iPoint << "] = "
738  << gridMetric[iPoint3] << "," << expectedMetric[iPoint3]
739  << std::endl;
740  }
741  }
742  sectionResults = CheckResult(sectionResults,testComm);
743  }
744  if(!sectionResults){
745  testResults[numDim-2] = false;
746  messageStream << "Test failed in " << numDim << " dimensions." << std::endl;
747  myGlobal.StdOut(messageStream.str());
748  pcpp::io::RenewStream(messageStream);
749  continue;
750  }
751 
752  // if(gridJacobian.size() != 1){
753  // myGlobal.StdOut("Grid Jacobian not of proper size.");
754  // sectionResults = false;
755  // continue;
756  // }
757  // if(gridJacobian[0] != expectedJacobian){
758  // messageStream << "Grid Jacobian does not have expected value ("
759  // << gridJacobian[0] << " != " << expectedJacobian
760  // << "). Test failed." << std::endl;
761  // myGlobal.StdOut(messageStream.str());
762  // sectionResults = false;
763  // continue;
764  // }
765 
766  }
767 
768  if(testResults[0])
769  myGlobal.StdOut("2D Test passed.\n");
770  if(testResults[1])
771  myGlobal.StdOut("3D Test passed.\n");
772 
773  parallelUnitResults.UpdateResult("Grid:PBS:Rectilinear2D",testResults[0]);
774  parallelUnitResults.UpdateResult("Grid:PBS:Rectilinear3D",testResults[1]);
775 
776 };
777 
778 #define CLSIZE1 101
779 #define CLSIZE2 101
780 #define CLSIZE3 1001
781 
782 
783 template<int numDim>
784 int GenerateCurvilinearMetric1(std::vector<size_t> &ijk,
785  std::vector<double> &met,
786  std::vector<double> &jac,
787  double &jacm1)
788 {
789 
790  double scale[] = {CLSIZE1,CLSIZE2,CLSIZE3};
791  double r0 = 0.25;
792  double rmax = 1.25;
793  double rt = 2.0;
794  double dr = (rmax - r0)/(double(scale[0]-1));
795  double d0 = (2.0*M_PI)/(double(scale[1]));
796  double dp = (2.0*M_PI)/(double(scale[2]));
797  double r = r0 + ijk[0]*dr;
798  double a0 = ijk[1]*d0;
799 
800 
801  if(numDim == 2){
802  double dxDXi = dr*std::cos(a0);
803  double dxDEta = -1.0*d0*r*std::sin(a0);
804  double dyDXi = dr*std::sin(a0);
805  double dyDEta = d0*r*std::cos(a0);
806 
807  met[0] = dyDEta;
808  met[1] = -1.0*dxDEta;
809  met[2] = -1.0*dyDXi;
810  met[3] = dxDXi;
811 
812  jac[0] = dxDXi;
813  jac[1] = dxDEta;
814  jac[2] = dyDXi;
815  jac[3] = dyDEta;
816 
817  jacm1 = met[0]*met[3] - met[1]*met[2];
818 
819  if(std::abs(jacm1 - r*dr*d0) > 1e-12){
820  std::cout << "Warning: Errors detected in expected 2d Jacobian." << std::endl;
821  }
822  } else if (numDim == 3){
823 
824  double phi = ijk[2]*dp;
825 
826  double x = (rt + r*std::cos(a0))*std::sin(phi);
827 
828  double x_xi = dr*std::cos(a0)*std::sin(phi);
829  double x_xi_eta = -1.0*d0*dr*std::sin(a0)*std::sin(phi);
830  double x_xi_zeta = dr*std::cos(a0)*dp*std::cos(phi);
831  jac[0] = x_xi;
832 
833  double x_eta = -1.0*d0*r*std::sin(a0)*std::sin(phi);
834  double x_eta_xi = -1.0*d0*dr*std::sin(a0)*std::sin(phi);
835  double x_eta_zeta = -1.0*dp*d0*r*std::sin(a0)*std::cos(phi);
836  jac[1] = x_eta;
837 
838  double x_zeta = (rt + r*std::cos(a0))*dp*std::cos(phi);
839  double x_zeta_xi = dr*std::cos(a0)*dp*std::cos(phi);
840  double x_zeta_eta = -1.0*r*d0*std::sin(a0)*dp*std::cos(phi);
841  jac[2] = x_zeta;
842 
843  double y = (rt + r*std::cos(a0))*std::cos(phi);
844 
845  double y_xi = dr*std::cos(a0)*std::cos(phi);
846  double y_xi_eta = -1.0*dr*d0*std::sin(a0)*std::cos(phi);
847  double y_xi_zeta = -1.0*dr*dp*std::cos(a0)*std::sin(phi);
848  jac[3] = y_xi;
849 
850  double y_eta = -1.0*d0*r*std::sin(a0)*std::cos(phi);
851  double y_eta_xi = -1.0*d0*dr*std::sin(a0)*std::cos(phi);
852  double y_eta_zeta = d0*dp*r*std::sin(a0)*std::sin(phi);
853  jac[4] = y_eta;
854 
855  double y_zeta = -1.0*dp*(rt+r*std::cos(a0))*std::sin(phi);
856  double y_zeta_xi = -1.0*dp*dr*std::cos(a0)*std::sin(phi);
857  double y_zeta_eta = dp*r*d0*std::sin(a0)*std::sin(phi);
858  jac[5] = y_zeta;
859 
860  double z = r*std::sin(a0);
861 
862  double z_xi = dr*std::sin(a0);
863  double z_xi_eta = dr*d0*std::cos(a0);
864  double z_xi_zeta = 0.0;
865  jac[6] = z_xi;
866 
867  double z_eta = d0*r*std::cos(a0);
868  double z_eta_xi = d0*dr*std::cos(a0);
869  double z_eta_zeta = 0.0;
870  jac[7] = z_eta;
871 
872  double z_zeta = 0.0;
873  double z_zeta_xi = 0.0;
874  double z_zeta_eta = 0.0;
875  jac[8] = z_zeta;
876 
877  // Expected grid metric tensor
878  met[0] = y_eta_zeta*z + y_eta*z_zeta - y_zeta_eta*z - y_zeta*z_eta;
879  met[1] = z_eta_zeta*x + z_eta*x_zeta - z_zeta_eta*x - z_zeta*x_eta;
880  met[2] = x_eta_zeta*y + x_eta*y_zeta - x_zeta_eta*y - x_zeta*y_eta;
881  met[3] = y_zeta_xi*z + y_zeta*z_xi - y_xi_zeta*z - y_xi*z_zeta;
882  met[4] = z_zeta_xi*x + z_zeta*x_xi - z_xi_zeta*x - z_xi*x_zeta;
883  met[5] = x_zeta_xi*y + x_zeta*y_xi - x_xi_zeta*y - x_xi*y_zeta;
884  met[6] = y_xi_eta*z + y_xi*z_eta - y_eta_xi*z - y_eta*z_xi;
885  met[7] = z_xi_eta*x + z_xi*x_eta - z_eta_xi*x - z_eta*x_xi;
886  met[8] = x_xi_eta*y + x_xi*y_eta - x_eta_xi*y - x_eta*y_xi;
887 
888  jacm1 =
889  x_xi*(z_zeta*y_eta - z_eta*y_zeta) -
890  y_xi*(z_zeta*x_eta - z_eta*x_zeta) +
891  z_xi*(y_zeta*x_eta - y_eta*x_zeta);
892 
893  double expectedJac = r*(rt+r*std::cos(a0))*dr*dp*d0;
894 
895  double jacErr = std::abs(jacm1 - expectedJac)/expectedJac;
896 
897  if(jacErr > 1e-3){
898  std::cout << "Warning: Errors detected in expected 3d Jacobian ("
899  << jacm1 << "," << expectedJac << ")." << std::endl;
900  }
901 
902  }
903  if(jacm1 < 0){
904  std::cout << "Warning: Expected Jacobian is negative." << std::endl;
905  }
906  return(0);
907 };
908 
909 template<int numDim>
910 int GenerateCurvilinearMetric2(std::vector<size_t> &ijk,
911  std::vector<double> &xyz,double &jacm1)
912 {
913 
914  double scale[] = {CLSIZE1,CLSIZE2,CLSIZE3};
915  double I = (ijk[0] + 1)/scale[0];
916  double J = (ijk[1] + 1)/scale[1];
917  // xyz[0] = I*I + J;
918  // xyz[1] = 2.0*J*J + I;
919 
920  double dxDXi = 2.0*I/scale[0];
921  double dxDEta = 1.0/scale[1];
922  double dyDXi = 1.0/scale[0];
923  double dyDEta = 4.0*J/scale[1];
924 
925  xyz[0] = dyDEta;
926  xyz[1] = -1.0*dxDEta;
927  xyz[2] = -1.0*dyDXi;
928  xyz[3] = dxDXi;
929 
930  jacm1 = xyz[0]*xyz[3] - xyz[1]*xyz[2];
931 
932  if(jacm1 < 0){
933  std::cout << "Warning: Negative Jacobian." << std::endl;
934  }
935  if (numDim == 3){
936  xyz[0] = 0;
937  xyz[1] = 0;
938  xyz[2] = 0;
939  xyz[3] = 0;
940  xyz[4] = 0;
941  xyz[5] = 0;
942  xyz[6] = 0;
943  xyz[7] = 0;
944  xyz[8] = 0;
945  }
946  return(0);
947 };
948 
949 void TestGrid_CurvilinearMetric(ix::test::results &parallelUnitResults,
950  pcpp::CommunicatorType &testComm)
951 {
952 
953  std::ostringstream messageStream;
954 
955  fixtures::ConfigurationType testConfig;
956 
957  typedef simulation::grid::parallel_blockstructured pbsgrid_t;
958  typedef simulation::state::base state_t;
959  typedef simulation::grid::halo halo_t;
960  typedef plascom2::operators::sbp::base operator_t;
961  typedef pcpp::ParallelGlobalType global_t;
962  typedef pcpp::CommunicatorType comm_t;
963  typedef pcpp::IndexIntervalType interval_t;
964 
965  global_t myGlobal;
966  myGlobal.Init("TestGrid_CurvilinearMetric",testComm);
967  myGlobal.SetVerbLevel(3);
968  myGlobal.Profiling(true);
969 
970  std::vector<bool> testResults(2,true);
971  std::vector<size_t> allSizes(3,0);
972  std::vector<double> scale(3,0);
973 
974  // Test metrics for grids of this size
975  allSizes[0] = CLSIZE1;
976  allSizes[1] = CLSIZE2;
977  allSizes[2] = CLSIZE3;
978 
979  for(int numDim = 2;numDim <= 3;numDim++){
980 
981  bool sectionResults = true;
982 
983  pbsgrid_t testGrid;
984  int interiorSBPOrder = 4;
985  int haloDepth = interiorSBPOrder+1;
986  pcpp::ParallelTopologyInfoType cartInfo;
987 
988 
989  // Make a non-periodic cartesian grid & topology
990  cartInfo.numDimensions = numDim;
991  cartInfo.isPeriodic.resize(numDim,1);
992  cartInfo.isPeriodic[0] = 0;
993 
994 
995  // Set grid sizes
996  size_t numGlobalNodes = 1;
997  size_t numGlobalCells = 1;
998  std::vector<size_t> gridSizes(numDim,0);
999  for(int iDim = 0;iDim < numDim;iDim++){
1000  gridSizes[iDim] = allSizes[iDim];
1001  numGlobalNodes *= gridSizes[iDim];
1002  numGlobalCells *= (gridSizes[iDim]-1);
1003  }
1004  testGrid.SetGridSizes(gridSizes);
1005 
1006 
1007  messageStream << "Setting up grid with sizes: (";
1008  pcpp::io::DumpContents(messageStream,gridSizes,",");
1009  messageStream << ")" << std::endl;
1010  myGlobal.StdOut(messageStream.str());
1011  pcpp::io::RenewStream(messageStream);
1012 
1013 
1014  // This call does everything to set up the grid and required halos
1015  messageStream << "====== Grid Parallel Setup ========= " << std::endl;
1016  if(testGrid.ParallelSetup(testComm,cartInfo,haloDepth,messageStream)){
1017  myGlobal.ErrOut("Grid::ParallelSetup failed with messages:\n");
1018  myGlobal.ErrOut(messageStream.str());
1019  pcpp::io::RenewStream(messageStream);
1020  sectionResults = false;
1021  }
1022  sectionResults = CheckResult(sectionResults,testComm);
1023  if(!sectionResults){
1024  testResults[numDim-2] = false;
1025  continue;
1026  }
1027  messageStream << "==================================== " << std::endl;
1028 
1029  // Generate a report of the Cartesian grid and topo setup
1030  messageStream << "------- Grid & Topo Report --------- " << std::endl;
1031  pcpp::report::CartesianSetup(messageStream,cartInfo);
1032  messageStream << simulation::report::Grid(testGrid) << std::endl;
1033  messageStream << "------------------------------------ " << std::endl;
1034  myGlobal.StdOut(messageStream.str());
1035  pcpp::io::RenewStream(messageStream);
1036 
1037 
1038  testComm.Barrier();
1039 
1040  interval_t &partitionBufferInterval(testGrid.PartitionBufferInterval());
1041  interval_t &partitionInterval(testGrid.PartitionInterval());
1042  std::vector<size_t> partitionSizes(partitionInterval.Sizes());
1043  std::vector<size_t> partitionStarts(partitionInterval.Starts());
1044  std::vector<size_t> partitionBufferStarts(partitionBufferInterval.Starts());
1045  std::vector<size_t> &bufferSizes(testGrid.BufferSizes());
1046  interval_t bufferInterval;
1047  bufferInterval.InitSimple(bufferSizes);
1048  std::vector<size_t> partitionBufferIndices;
1049  bufferInterval.GetFlatIndices(partitionBufferInterval,partitionBufferIndices);
1050 
1051  // Set up grid coordinates and expected metric
1052  size_t numPointsBuffer = testGrid.BufferSize();
1053 
1054  int numMetricComponents = numDim*numDim;
1055  std::vector<double> expectedMetric(numMetricComponents*numPointsBuffer,0);
1056  std::vector<double> metricIdentities;
1057  std::vector<double> expectedJacobian(2*numPointsBuffer,0);
1058  std::vector<double> expectedJacobianMatrix(numMetricComponents*numPointsBuffer,0);
1059  std::vector<double> pointMetric(numMetricComponents,0);
1060  std::vector<double> pointMatrix(numMetricComponents,0);
1061  std::vector<size_t> pointIndices(numDim,0);
1062  testGrid.SetType(simulation::grid::CURVILINEAR);
1063 
1064  state_t testState;
1065  testState.AddField("expectedMetric",'n',numMetricComponents,8,"");
1066  testState.AddField("metricIdentities",'n',numDim,8,"");
1067  testState.AddField("expectedJacobian",'n',2,8,"");
1068  testState.AddField("expectedJacobianMatrix",'n',numMetricComponents,8,"");
1069  testState.AddField("gridMetric",'n',numMetricComponents,8,"");
1070  testState.AddField("metricError",'n',numMetricComponents,8,"");
1071  testState.AddField("gridJacobian",'n',2,8,"");
1072  testState.AddField("jacobianMatrix",'n',numMetricComponents,8,"");
1073  testState.AddField("sconn",'n',numDim,4,"");
1074  testState.AddField("xyz",'n',numDim,8,"");
1075  testState.AddField("ijk",'n',numDim,4,"");
1076  testState.AddField("bufferIndex1",'n',1,4,"");
1077  testState.AddField("jacobianError",'n',2,8,"");
1078  testState.AddField("jacobianMatrixError",'n',numMetricComponents,8,"");
1079  testState.Create(numPointsBuffer,0);
1080 
1081  testState.SetFieldBuffer("expectedMetric",&expectedMetric[0]);
1082  testState.SetFieldBuffer("expectedJacobian",&expectedJacobian[0]);
1083  testState.SetFieldBuffer("expectedJacobianMatrix",&expectedJacobianMatrix[0]);
1084 
1085  int *ijkFieldPtr = testState.GetFieldBuffer<int>("ijk");
1086  int *bufferIndex1 = testState.GetFieldBuffer<int>("bufferIndex1");
1087 
1088  std::ostringstream topoTypeStream;
1089  topoTypeStream << numDim << "DSMesh";
1090 
1091  std::vector<size_t> rGridSizes(testGrid.GridSizes());
1092  std::reverse(rGridSizes.begin(),rGridSizes.end());
1093 
1094  testComm.Barrier();
1095  messageStream << "Running " << numDim << "-dimensional test." << std::endl;
1096  myGlobal.StdOut(messageStream.str());
1097  pcpp::io::RenewStream(messageStream);
1098 
1099  if(numDim == 2) {
1100 
1101  testComm.Barrier();
1102  myGlobal.StdOut("Generating 2D grid.\n");
1103 
1104  testGrid.GenerateCoordinates(gridgen::CurvilinearGrid1<2,CLSIZE1,CLSIZE2>,messageStream);
1105 
1106  testComm.Barrier();
1107  myGlobal.StdOut("Exchanging coordinates.\n");
1108 
1109  testGrid.ExchangeCoordinates();
1110 
1111  testComm.Barrier();
1112  myGlobal.StdOut("Generating expected grid metrics.\n");
1113 
1114  size_t nPlane = bufferSizes[0];
1115  for(size_t iY = 0;iY < partitionSizes[1];iY++){
1116  pointIndices[1] = iY + partitionStarts[1];
1117  size_t bufferY = iY + partitionBufferStarts[1];
1118  for(size_t iX = 0;iX < partitionSizes[0];iX++){
1119  pointIndices[0] = iX + partitionStarts[0];
1120  size_t bufferX = iX + partitionBufferStarts[0];
1121  size_t bufferIndex = bufferY*nPlane + bufferX;
1122  ijkFieldPtr[bufferIndex] = pointIndices[0];
1123  ijkFieldPtr[bufferIndex+numPointsBuffer] = pointIndices[1];
1124  bufferIndex1[bufferIndex] = bufferIndex;
1125  for(int iComp = 0;iComp < numMetricComponents;iComp++){
1126  pointMetric[iComp] = 0;
1127  }
1128  GenerateCurvilinearMetric1<2>(pointIndices,pointMetric,pointMatrix,
1129  expectedJacobian[bufferIndex+numPointsBuffer]);
1130  expectedJacobian[bufferIndex] = 1.0/expectedJacobian[bufferIndex+numPointsBuffer];
1131  for(int iComp = 0;iComp < numMetricComponents;iComp++){
1132  expectedMetric[bufferIndex+iComp*numPointsBuffer] = pointMetric[iComp];
1133  expectedJacobianMatrix[bufferIndex+iComp*numPointsBuffer] = pointMatrix[iComp];
1134  }
1135  }
1136  }
1137 
1138  } else if(numDim == 3) {
1139 
1140  testComm.Barrier();
1141  myGlobal.StdOut("Generating 3D grid.\n");
1142 
1143  testGrid.GenerateCoordinates(gridgen::CurvilinearGrid1<3,CLSIZE1,CLSIZE2,CLSIZE3>,messageStream);
1144 
1145  testComm.Barrier();
1146  myGlobal.StdOut("Exchanging coordinates.\n");
1147 
1148  testGrid.ExchangeCoordinates();
1149 
1150  testComm.Barrier();
1151  myGlobal.StdOut("Generating expected grid metrics.\n");
1152 
1153  size_t nPlane = bufferSizes[0]*bufferSizes[1];
1154  for(size_t iZ = 0;iZ < partitionSizes[2];iZ++){
1155  pointIndices[2] = iZ + partitionStarts[2];
1156  size_t bufferZ = iZ + partitionBufferStarts[2];
1157  for(size_t iY = 0;iY < partitionSizes[1];iY++){
1158  pointIndices[1] = iY + partitionStarts[1];
1159  size_t bufferY = iY + partitionBufferStarts[1];
1160  for(size_t iX = 0;iX < partitionSizes[0];iX++){
1161  pointIndices[0] = iX + partitionStarts[0];
1162  size_t bufferX = iX + partitionBufferStarts[0];
1163  size_t bufferIndex = bufferZ*nPlane + bufferY*bufferSizes[0] + bufferX;
1164  ijkFieldPtr[bufferIndex] = pointIndices[0];
1165  ijkFieldPtr[bufferIndex+numPointsBuffer] = pointIndices[1];
1166  ijkFieldPtr[bufferIndex+2*numPointsBuffer] = pointIndices[2];
1167  bufferIndex1[bufferIndex] = bufferIndex;
1168  for(int iComp = 0;iComp < numMetricComponents;iComp++){
1169  pointMetric[iComp] = 0;
1170  }
1171  GenerateCurvilinearMetric1<3>(pointIndices,pointMetric,pointMatrix,
1172  expectedJacobian[bufferIndex+numPointsBuffer]);
1173  expectedJacobian[bufferIndex] = 1.0/expectedJacobian[bufferIndex+numPointsBuffer];
1174  for(int iComp = 0;iComp < numMetricComponents;iComp++){
1175  expectedMetric[bufferIndex+iComp*numPointsBuffer] = pointMetric[iComp];
1176  expectedJacobianMatrix[bufferIndex+iComp*numPointsBuffer] = pointMatrix[iComp];
1177  }
1178  }
1179  }
1180  }
1181  }
1182  std::vector<double> &gridCoordinates(testGrid.CoordinateData());
1183  testState.SetFieldBuffer("xyz",&gridCoordinates[0]);
1184 
1185 
1186  testComm.Barrier();
1187  myGlobal.StdOut(messageStream.str());
1188  pcpp::io::RenewStream(messageStream);
1189 
1190  myGlobal.StdOut("Setting up differential operators for grid.\n");
1191 
1192  // Set up differential operator for grid
1193  operator_t sbpOperator;
1194  plascom2::operators::sbp::Initialize(sbpOperator,interiorSBPOrder);
1195  int connBoundaryDepth = (sbpOperator.numStencils-2)/2;
1196  // Forticate the stencil starts
1197  for(int iStencil = 0;iStencil < sbpOperator.numStencils;iStencil++)
1198  sbpOperator.stencilStarts[iStencil]++;
1199 
1200  std::vector<int> stencilConnectivity(numDim*numPointsBuffer,0);
1201  testState.SetFieldBuffer("sconn",&stencilConnectivity[0]);
1202 
1203  std::vector<size_t> opInterval;
1204  partitionBufferInterval.Flatten(opInterval);
1205  std::vector<size_t>::iterator dOpIt = opInterval.begin();
1206  while(dOpIt != opInterval.end()){
1207  *dOpIt += 1;
1208  dOpIt++;
1209  }
1210 
1211  if(plascom2::operators::sbp::CreateStencilConnectivity(numDim,&bufferSizes[0],&opInterval[0],
1212  connBoundaryDepth,&cartInfo.isPeriodic[0],
1213  &stencilConnectivity[0],true)){
1214  messageStream << "CreateStencilConnectivity failed. Aborting the test." << std::endl;
1215  myGlobal.StdOut(messageStream.str());
1216  pcpp::io::RenewStream(messageStream);
1217  sectionResults = false;
1218  }
1219  sectionResults = CheckResult(sectionResults,testComm);
1220  if(!sectionResults){
1221  testResults[numDim-2] = false;
1222  continue;
1223  }
1224  testGrid.SetupDifferentialOperator(&sbpOperator,&stencilConnectivity[0]);
1225 
1226  testComm.Barrier();
1227  myGlobal.StdOut("Computing Jacobian matrix.\n");
1228  // Compute the inverse Jacobian matrix
1229  if(testGrid.ComputeJacobianMatrix(messageStream)){
1230  myGlobal.StdOut("Jacobian matrix generation failed with messages:\n");
1231  myGlobal.StdOut(messageStream.str());
1232  myGlobal.StdOut("Skipping the rest of this test.\n");
1233  pcpp::io::RenewStream(messageStream);
1234  sectionResults = false;
1235  }
1236  sectionResults = CheckResult(sectionResults,testComm);
1237  if(!sectionResults){
1238  testResults[numDim-2] = false;
1239  continue;
1240  }
1241 
1242  testComm.Barrier();
1243  myGlobal.StdOut(messageStream.str());
1244  pcpp::io::RenewStream(messageStream);
1245 
1246  myGlobal.StdOut("Computing metrics.\n");
1247  // Compute the actual metrics and Jacobian
1248  if(testGrid.ComputeMetrics(messageStream)){
1249  myGlobal.StdOut("Metrics generation failed with messages:\n");
1250  myGlobal.StdOut(messageStream.str());
1251  myGlobal.StdOut("\nSkipping rest of this test.\n");
1252  pcpp::io::RenewStream(messageStream);
1253  sectionResults = false;
1254  }
1255  sectionResults = CheckResult(sectionResults,testComm);
1256  if(!sectionResults){
1257  testResults[numDim-2] = false;
1258  continue;
1259  }
1260  testComm.Barrier();
1261  myGlobal.StdOut(messageStream.str());
1262  pcpp::io::RenewStream(messageStream);
1263 
1264  myGlobal.StdOut("Computing metric identities.\n");
1265  if(testGrid.ComputeMetricIdentities(metricIdentities,messageStream)){
1266  myGlobal.StdOut("Compute metric identities failed with messages:\n");
1267  myGlobal.StdOut(messageStream.str());
1268  myGlobal.StdOut("\nSkipping rest of this test.\n");
1269  pcpp::io::RenewStream(messageStream);
1270  sectionResults = false;
1271  }
1272  sectionResults = CheckResult(sectionResults,testComm);
1273  if(!sectionResults){
1274  testResults[numDim-2] = false;
1275  continue;
1276  }
1277  testComm.Barrier();
1278  myGlobal.StdOut(messageStream.str());
1279  pcpp::io::RenewStream(messageStream);
1280 
1281  myGlobal.StdOut("Checking Jacobian matrix.\n");
1282  // Preliminary - check the components of J^-1. If these
1283  // aren't right, the metrics have no chance.
1284  std::vector<double> &jm1Matrix(testGrid.JacobianMatrix());
1285  testState.SetFieldBuffer("jacobianMatrix",&jm1Matrix[0]);
1286 
1287  double *jm1MatrixError = testState.GetFieldBuffer<double>("jacobianMatrixError");
1288  std::vector<size_t>::iterator indexIt = partitionBufferIndices.begin();
1289  while(indexIt != partitionBufferIndices.end()){
1290  size_t iPoint = *indexIt++;
1291  for(int iComponent = 0;iComponent < numMetricComponents;iComponent++){
1292  size_t iMetricPoint = iPoint + iComponent*numPointsBuffer;
1293  jm1MatrixError[iMetricPoint] = std::abs(jm1Matrix[iMetricPoint] -
1294  expectedJacobianMatrix[iMetricPoint]);
1295  if(jm1MatrixError[iMetricPoint] > 1e-7){
1296  sectionResults = false;
1297  }
1298  }
1299  }
1300  testComm.Barrier();
1301  sectionResults = CheckResult(sectionResults,testComm);
1302  if(!sectionResults){
1303  testResults[numDim-2] = false;
1304  messageStream << numDim << "D Jacobian matrix test failed." << std::endl;
1305  myGlobal.StdOut(messageStream.str());
1306  pcpp::io::RenewStream(messageStream);
1307  }
1308  sectionResults = true;
1309 
1310  myGlobal.StdOut("Checking Jacobian determinant.\n");
1311  std::vector<double> &gridJacobian(testGrid.Jacobian());
1312  testState.SetFieldBuffer("gridJacobian",&gridJacobian[0]);
1313  bool jacSection = true;
1314  if(gridJacobian.size() != 2*numPointsBuffer){
1315  myGlobal.StdOut("Grid Jacobian not of proper size.");
1316  jacSection = false;
1317  }
1318  jacSection = CheckResult(jacSection,testComm);
1319  if(!jacSection){
1320  testResults[numDim-2] = false;
1321  messageStream << numDim << "D Jacobian size test failed." << std::endl;
1322  myGlobal.StdOut(messageStream.str());
1323  pcpp::io::RenewStream(messageStream);
1324  // continue;
1325  }
1326  if(jacSection){
1327  double *jacobianError = testState.GetFieldBuffer<double>("jacobianError");
1328  indexIt = partitionBufferIndices.begin();
1329  while(indexIt != partitionBufferIndices.end()){
1330  size_t iPoint = *indexIt++;
1331  for(int iComponent = 0;iComponent < 2;iComponent++){
1332  size_t iMetricPoint = iPoint + iComponent*numPointsBuffer;
1333  jacobianError[iMetricPoint] = std::abs(gridJacobian[iMetricPoint] -
1334  expectedJacobian[iMetricPoint]);
1335  jacobianError[iMetricPoint] /= expectedJacobian[iMetricPoint]; // positive-definite != 0
1336  if(jacobianError[iMetricPoint] > 1e-6){
1337  jacSection = false;
1338  }
1339  }
1340  }
1341  }
1342  testComm.Barrier();
1343  jacSection = CheckResult(sectionResults,testComm);
1344  if(!jacSection){
1345  testResults[numDim-2] = false;
1346  messageStream << numDim << "D Jacobian accuracy test failed." << std::endl;
1347  myGlobal.StdOut(messageStream.str());
1348  pcpp::io::RenewStream(messageStream);
1349  }
1350 
1351  // Now check the actual metrics
1352  testComm.Barrier();
1353  myGlobal.StdOut("Checking the metrics.\n");
1354  std::vector<double> &gridMetric(testGrid.Metric());
1355  testState.SetFieldBuffer("gridMetric",&gridMetric[0]);
1356 
1357  double *metricError = testState.GetFieldBuffer<double>("metricError");
1358  std::vector<bool> metricTest(numMetricComponents,true);
1359 
1360  indexIt = partitionBufferIndices.begin();
1361  while(indexIt != partitionBufferIndices.end()){
1362  size_t iPoint = *indexIt++;
1363  for(int iComponent = 0;iComponent < numMetricComponents;iComponent++){
1364  size_t iMetricPoint = iPoint + iComponent*numPointsBuffer;
1365  metricError[iMetricPoint] = std::abs(gridMetric[iMetricPoint] -
1366  expectedMetric[iMetricPoint]);
1367  if(metricError[iMetricPoint] > 1e-6){
1368  metricTest[iComponent] = false;
1369  sectionResults = false;
1370  }
1371  }
1372  }
1373  for(int iComp = 0;iComp < numMetricComponents;iComp++){
1374  if(!metricTest[iComp]){
1375  messageStream << "Metric component " << iComp+1 << " failed test."
1376  << std::endl;
1377  myGlobal.StdOut(messageStream.str());
1378  pcpp::io::RenewStream(messageStream);
1379  }
1380  }
1381  testComm.Barrier();
1382  sectionResults = CheckResult(sectionResults,testComm);
1383  if(!sectionResults){
1384  testResults[numDim-2] = false;
1385  messageStream << numDim << "D metric accuracy test failed." << std::endl;
1386  myGlobal.StdOut(messageStream.str());
1387  pcpp::io::RenewStream(messageStream);
1388  // continue;
1389  }
1390 
1391  // Check the metric identities
1392  testComm.Barrier();
1393  myGlobal.StdOut("Checking the metric identities.\n");
1394  testState.SetFieldBuffer("metricIdentities",&metricIdentities[0]);
1395 
1396  std::vector<bool> metricIdentityTest(numDim,true);
1397  std::vector<double> maxIdentityError(numDim,0.0);
1398  double identityTolerance = 5e-11;
1399  indexIt = partitionBufferIndices.begin();
1400  while(indexIt != partitionBufferIndices.end()){
1401  size_t iPoint = *indexIt++;
1402  for(int iComponent = 0;iComponent < numDim;iComponent++){
1403  size_t idPoint = iPoint + iComponent*numPointsBuffer;
1404  double identityMagnitude = std::abs(metricIdentities[idPoint]);
1405  if(identityMagnitude > identityTolerance){
1406  metricIdentityTest[iComponent] = false;
1407  sectionResults = false;
1408  }
1409  if(identityMagnitude > maxIdentityError[iComponent])
1410  maxIdentityError[iComponent] = identityMagnitude;
1411  }
1412  }
1413  for(int iComp = 0;iComp < numDim;iComp++){
1414  messageStream << "Max identity error for component " << iComp+1 << " = "
1415  << maxIdentityError[iComp] << std::endl;
1416  if(!metricIdentityTest[iComp]){
1417  messageStream << "Metric identity component " << iComp+1 << " failed test."
1418  << std::endl;
1419  myGlobal.StdOut(messageStream.str());
1420  pcpp::io::RenewStream(messageStream);
1421  }
1422  }
1423  testComm.Barrier();
1424  sectionResults = CheckResult(sectionResults,testComm);
1425  if(!sectionResults){
1426  testResults[numDim-2] = false;
1427  messageStream << numDim << "D metric test failed." << std::endl;
1428  myGlobal.StdOut(messageStream.str());
1429  pcpp::io::RenewStream(messageStream);
1430  // continue;
1431  }
1432  messageStream << std::endl;
1433  myGlobal.StdOut(messageStream.str());
1434  pcpp::io::RenewStream(messageStream);
1435 
1436  myGlobal.StdOut("Writing test results to HDF5 file.\n");
1437  // Write the testing data for visualization
1438  std::ostringstream Ostr;
1439  Ostr << "curvilinear" << numDim << "d";
1440  const std::string domainName(Ostr.str());
1441  const std::string hdfFileName(domainName+".h5");
1442  const std::string xdmfFileName(domainName+".xdmf");
1443  const std::string geometryName("curvilinear");
1444  const std::string gridName("test1");
1445  if(ix::sys::FILEEXISTS(hdfFileName)){
1446  ix::sys::Remove(hdfFileName);
1447  ix::sys::Remove(xdmfFileName);
1448  }
1449  fixtures::ConfigurationType testConfig;
1450  state_t paramState;
1451  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
1452  gridName,testGrid,testState,paramState,
1453  testConfig,messageStream)){
1454 
1455  myGlobal.StdOut("Output to HDF5 failed:\n");
1456  // myGlobal.StdOut(messageStream.str());
1457  myGlobal.StdOut("Continuing....\n");
1458 
1459  }
1460  pcpp::io::RenewStream(messageStream);
1461  }
1462  testComm.Barrier();
1463  myGlobal.StdOut("Tests done.\n");
1464  parallelUnitResults.UpdateResult("Grid:PBS:Curvilinear2D",testResults[0]);
1465  parallelUnitResults.UpdateResult("Grid:PBS:Curvilinear3D",testResults[1]);
1466 
1467 };
1468 
1469 
1470 void TestGrid_CurvilinearVGWavy(ix::test::results &parallelUnitResults,
1471  pcpp::CommunicatorType &testComm)
1472 {
1473 
1474  std::ostringstream messageStream;
1475 
1476  fixtures::ConfigurationType testConfig;
1477 
1478  typedef simulation::grid::parallel_blockstructured pbsgrid_t;
1479  typedef simulation::state::base state_t;
1480  typedef simulation::grid::halo halo_t;
1481  typedef plascom2::operators::sbp::base operator_t;
1482  typedef pcpp::ParallelGlobalType global_t;
1483  typedef pcpp::CommunicatorType comm_t;
1484  typedef pcpp::IndexIntervalType interval_t;
1485 
1486  global_t myGlobal;
1487  myGlobal.Init("TestGrid_CurvilinearMetricVG",testComm);
1488  myGlobal.SetVerbLevel(3);
1489  myGlobal.Profiling(true);
1490 
1491  std::vector<bool> testResults(2,true);
1492  std::vector<size_t> allSizes(3,0);
1493  std::vector<double> scale(3,0);
1494 
1495  // Test metrics for grids of this size
1496  allSizes[0] = 21;
1497  allSizes[1] = 21;
1498  allSizes[2] = 21;
1499 
1500  int numDim = 3;
1501 
1502  bool sectionResults = true;
1503 
1504  pbsgrid_t testGrid;
1505  int interiorSBPOrder = 4;
1506  int haloDepth = interiorSBPOrder+1;
1507  pcpp::ParallelTopologyInfoType cartInfo;
1508 
1509 
1510  // Make a non-periodic cartesian grid & topology
1511  cartInfo.numDimensions = numDim;
1512  cartInfo.isPeriodic.resize(numDim,1);
1513 
1514 
1515  // Set grid sizes
1516  size_t numGlobalNodes = 1;
1517  size_t numGlobalCells = 1;
1518  std::vector<size_t> gridSizes(numDim,0);
1519  for(int iDim = 0;iDim < numDim;iDim++){
1520  gridSizes[iDim] = allSizes[iDim];
1521  numGlobalNodes *= gridSizes[iDim];
1522  numGlobalCells *= (gridSizes[iDim]-1);
1523  }
1524  testGrid.SetGridSizes(gridSizes);
1525 
1526 
1527  messageStream << "Setting up grid with sizes: (";
1528  pcpp::io::DumpContents(messageStream,gridSizes,",");
1529  messageStream << ")" << std::endl;
1530  myGlobal.StdOut(messageStream.str());
1531  pcpp::io::RenewStream(messageStream);
1532 
1533 
1534  // This call does everything to set up the grid and required halos
1535  messageStream << "====== Grid Parallel Setup ========= " << std::endl;
1536  if(testGrid.ParallelSetup(testComm,cartInfo,haloDepth,messageStream)){
1537  myGlobal.ErrOut("Grid::ParallelSetup failed with messages:\n");
1538  myGlobal.ErrOut(messageStream.str());
1539  pcpp::io::RenewStream(messageStream);
1540  sectionResults = false;
1541  }
1542  sectionResults = CheckResult(sectionResults,testComm);
1543  if(!sectionResults){
1544  testResults[numDim-2] = false;
1545  }
1546  messageStream << "==================================== " << std::endl;
1547 
1548  // Generate a report of the Cartesian grid and topo setup
1549  messageStream << "------- Grid & Topo Report --------- " << std::endl;
1550  pcpp::report::CartesianSetup(messageStream,cartInfo);
1551  messageStream << simulation::report::Grid(testGrid) << std::endl;
1552  messageStream << "------------------------------------ " << std::endl;
1553  myGlobal.StdOut(messageStream.str());
1554  pcpp::io::RenewStream(messageStream);
1555 
1556 
1557  testComm.Barrier();
1558 
1559  interval_t &partitionBufferInterval(testGrid.PartitionBufferInterval());
1560  interval_t &partitionInterval(testGrid.PartitionInterval());
1561  std::vector<size_t> partitionSizes(partitionInterval.Sizes());
1562  std::vector<size_t> partitionStarts(partitionInterval.Starts());
1563  std::vector<size_t> partitionBufferStarts(partitionBufferInterval.Starts());
1564  std::vector<size_t> &bufferSizes(testGrid.BufferSizes());
1565  interval_t bufferInterval;
1566  bufferInterval.InitSimple(bufferSizes);
1567  std::vector<size_t> partitionBufferIndices;
1568  bufferInterval.GetFlatIndices(partitionBufferInterval,partitionBufferIndices);
1569 
1570  // Set up grid coordinates and expected metric
1571  size_t numPointsBuffer = testGrid.BufferSize();
1572 
1573  int numMetricComponents = numDim*numDim;
1574  std::vector<double> expectedMetric(numMetricComponents*numPointsBuffer,0);
1575  std::vector<double> metricIdentities;
1576  std::vector<double> expectedJacobian(2*numPointsBuffer,0);
1577  std::vector<double> expectedJacobianMatrix(numMetricComponents*numPointsBuffer,0);
1578  std::vector<double> pointMetric(numMetricComponents,0);
1579  std::vector<double> pointMatrix(numMetricComponents,0);
1580  std::vector<size_t> pointIndices(numDim,0);
1581  testGrid.SetType(simulation::grid::CURVILINEAR);
1582 
1583  state_t testState;
1584  testState.AddField("expectedMetric",'n',numMetricComponents,8,"");
1585  testState.AddField("metricIdentities",'n',numDim,8,"");
1586  testState.AddField("expectedJacobian",'n',2,8,"");
1587  testState.AddField("expectedJacobianMatrix",'n',numMetricComponents,8,"");
1588  testState.AddField("gridMetric",'n',numMetricComponents,8,"");
1589  testState.AddField("metricError",'n',numMetricComponents,8,"");
1590  testState.AddField("gridJacobian",'n',2,8,"");
1591  testState.AddField("jacobianMatrix",'n',numMetricComponents,8,"");
1592  testState.AddField("sconn",'n',numDim,4,"");
1593  testState.AddField("xyz",'n',numDim,8,"");
1594  testState.AddField("ijk",'n',numDim,4,"");
1595  testState.AddField("bufferIndex1",'n',1,4,"");
1596  testState.AddField("jacobianError",'n',2,8,"");
1597  testState.AddField("jacobianMatrixError",'n',numMetricComponents,8,"");
1598  testState.Create(numPointsBuffer,0);
1599 
1600  testState.SetFieldBuffer("expectedMetric",&expectedMetric[0]);
1601  testState.SetFieldBuffer("expectedJacobian",&expectedJacobian[0]);
1602  testState.SetFieldBuffer("expectedJacobianMatrix",&expectedJacobianMatrix[0]);
1603 
1604  int *ijkFieldPtr = testState.GetFieldBuffer<int>("ijk");
1605  int *bufferIndex1 = testState.GetFieldBuffer<int>("bufferIndex1");
1606 
1607  std::ostringstream topoTypeStream;
1608  topoTypeStream << numDim << "DSMesh";
1609 
1610  std::vector<size_t> rGridSizes(testGrid.GridSizes());
1611  std::reverse(rGridSizes.begin(),rGridSizes.end());
1612 
1613  testComm.Barrier();
1614  myGlobal.StdOut(messageStream.str());
1615  pcpp::io::RenewStream(messageStream);
1616  myGlobal.StdOut("Generating 3D grid.\n");
1617  double xA = 1.0;
1618  double xOmega = 1.0/4.0;
1619  double xL = 4.0;
1620  testGrid.GenerateCoordinates(gridgen::VGWavy<21,21,21,4>,messageStream);
1621  std::vector<double> periodicLengths(numDim,xL);
1622  testGrid.SetPeriodicLengths(periodicLengths);
1623  testComm.Barrier();
1624  myGlobal.StdOut("Exchanging coordinates.\n");
1625  testGrid.ExchangeCoordinates();
1626  testComm.Barrier();
1627 
1628  std::vector<double> &gridCoordinates(testGrid.CoordinateData());
1629  testState.SetFieldBuffer("xyz",&gridCoordinates[0]);
1630 
1631 
1632  testComm.Barrier();
1633  myGlobal.StdOut(messageStream.str());
1634  pcpp::io::RenewStream(messageStream);
1635 
1636  myGlobal.StdOut("Setting up differential operators for grid.\n");
1637 
1638  // Set up differential operator for grid
1639  operator_t sbpOperator;
1640  plascom2::operators::sbp::Initialize(sbpOperator,interiorSBPOrder);
1641  int connBoundaryDepth = (sbpOperator.numStencils-2)/2;
1642  // Forticate the stencil starts
1643  for(int iStencil = 0;iStencil < sbpOperator.numStencils;iStencil++)
1644  sbpOperator.stencilStarts[iStencil]++;
1645 
1646  std::vector<int> stencilConnectivity(numDim*numPointsBuffer,0);
1647  testState.SetFieldBuffer("sconn",&stencilConnectivity[0]);
1648 
1649  std::vector<size_t> opInterval;
1650  partitionBufferInterval.Flatten(opInterval);
1651  std::vector<size_t>::iterator dOpIt = opInterval.begin();
1652  while(dOpIt != opInterval.end()){
1653  *dOpIt += 1;
1654  dOpIt++;
1655  }
1656 
1657  if(plascom2::operators::sbp::CreateStencilConnectivity(numDim,&bufferSizes[0],&opInterval[0],
1658  connBoundaryDepth,&cartInfo.isPeriodic[0],
1659  &stencilConnectivity[0],true)){
1660  messageStream << "CreateStencilConnectivity failed. Aborting the test." << std::endl;
1661  myGlobal.StdOut(messageStream.str());
1662  pcpp::io::RenewStream(messageStream);
1663  sectionResults = false;
1664  }
1665  sectionResults = CheckResult(sectionResults,testComm);
1666  if(!sectionResults){
1667  testResults[numDim-2] = false;
1668  }
1669  testGrid.SetupDifferentialOperator(&sbpOperator,&stencilConnectivity[0]);
1670 
1671  testComm.Barrier();
1672  myGlobal.StdOut("Computing Jacobian matrix.\n");
1673  // Compute the inverse Jacobian matrix
1674  if(testGrid.ComputeJacobianMatrix(messageStream)){
1675  myGlobal.StdOut("Jacobian matrix generation failed with messages:\n");
1676  myGlobal.StdOut(messageStream.str());
1677  myGlobal.StdOut("Skipping the rest of this test.\n");
1678  pcpp::io::RenewStream(messageStream);
1679  sectionResults = false;
1680  }
1681  sectionResults = CheckResult(sectionResults,testComm);
1682  if(!sectionResults){
1683  testResults[numDim-2] = false;
1684  }
1685 
1686  testComm.Barrier();
1687  myGlobal.StdOut(messageStream.str());
1688  pcpp::io::RenewStream(messageStream);
1689 
1690  myGlobal.StdOut("Computing metrics.\n");
1691  // Compute the actual metrics and Jacobian
1692  if(testGrid.ComputeMetrics(messageStream)){
1693  myGlobal.StdOut("Metrics generation failed with messages:\n");
1694  myGlobal.StdOut(messageStream.str());
1695  myGlobal.StdOut("\nSkipping rest of this test.\n");
1696  pcpp::io::RenewStream(messageStream);
1697  sectionResults = false;
1698  }
1699  sectionResults = CheckResult(sectionResults,testComm);
1700  if(!sectionResults){
1701  testResults[numDim-2] = false;
1702  }
1703  testComm.Barrier();
1704  myGlobal.StdOut(messageStream.str());
1705  pcpp::io::RenewStream(messageStream);
1706 
1707  myGlobal.StdOut("Computing metric identities.\n");
1708  if(testGrid.ComputeMetricIdentities(metricIdentities,messageStream)){
1709  myGlobal.StdOut("Compute metric identities failed with messages:\n");
1710  myGlobal.StdOut(messageStream.str());
1711  myGlobal.StdOut("\nSkipping rest of this test.\n");
1712  pcpp::io::RenewStream(messageStream);
1713  sectionResults = false;
1714  }
1715  sectionResults = CheckResult(sectionResults,testComm);
1716  if(!sectionResults){
1717  testResults[numDim-2] = false;
1718  }
1719  testComm.Barrier();
1720  myGlobal.StdOut(messageStream.str());
1721  pcpp::io::RenewStream(messageStream);
1722 
1723  // Check the metric identities
1724  testComm.Barrier();
1725  myGlobal.StdOut("Checking the metric identities.\n");
1726  testState.SetFieldBuffer("metricIdentities",&metricIdentities[0]);
1727 
1728  std::vector<bool> metricIdentityTest(numDim,true);
1729  std::vector<double> maxIdentityError(numDim,0.0);
1730  double identityTolerance = 5e-11;
1731  std::vector<size_t>::iterator indexIt = partitionBufferIndices.begin();
1732  while(indexIt != partitionBufferIndices.end()){
1733  size_t iPoint = *indexIt++;
1734  for(int iComponent = 0;iComponent < numDim;iComponent++){
1735  size_t idPoint = iPoint + iComponent*numPointsBuffer;
1736  double identityMagnitude = std::abs(metricIdentities[idPoint]);
1737  if(identityMagnitude > identityTolerance){
1738  metricIdentityTest[iComponent] = false;
1739  sectionResults = false;
1740  }
1741  if(identityMagnitude > maxIdentityError[iComponent])
1742  maxIdentityError[iComponent] = identityMagnitude;
1743  }
1744  }
1745  for(int iComp = 0;iComp < numDim;iComp++){
1746  messageStream << "Max identity error for component " << iComp+1 << " = "
1747  << maxIdentityError[iComp] << std::endl;
1748  if(!metricIdentityTest[iComp]){
1749  messageStream << "Metric identity component " << iComp+1 << " failed test."
1750  << std::endl;
1751  myGlobal.StdOut(messageStream.str());
1752  pcpp::io::RenewStream(messageStream);
1753  }
1754  }
1755  testComm.Barrier();
1756  sectionResults = CheckResult(sectionResults,testComm);
1757  if(!sectionResults){
1758  testResults[numDim-2] = false;
1759  messageStream << numDim << "D metric test failed." << std::endl;
1760  myGlobal.StdOut(messageStream.str());
1761  pcpp::io::RenewStream(messageStream);
1762  // continue;
1763  }
1764  messageStream << std::endl;
1765  myGlobal.StdOut(messageStream.str());
1766  pcpp::io::RenewStream(messageStream);
1767 
1768  myGlobal.StdOut("Writing test results to HDF5 file.\n");
1769  // Write the testing data for visualization
1770  std::ostringstream Ostr;
1771  Ostr << "vgwavy3d";
1772  const std::string domainName(Ostr.str());
1773  const std::string hdfFileName(domainName+".h5");
1774  const std::string xdmfFileName(domainName+".xdmf");
1775  const std::string geometryName("vgwavy");
1776  const std::string gridName("test1");
1777  if(ix::sys::FILEEXISTS(hdfFileName)){
1778  ix::sys::Remove(hdfFileName);
1779  ix::sys::Remove(xdmfFileName);
1780  }
1781 
1782  // fixtures::ConfigurationType testConfig;
1783  state_t paramState;
1784  if(plascom2::io::hdf5::OutputSingle(hdfFileName,domainName,geometryName,
1785  gridName,testGrid,testState,paramState,
1786  testConfig,messageStream)){
1787 
1788  myGlobal.StdOut("Output to HDF5 failed:\n");
1789  // myGlobal.StdOut(messageStream.str());
1790  myGlobal.StdOut("Continuing....\n");
1791 
1792  }
1793  pcpp::io::RenewStream(messageStream);
1794  testComm.Barrier();
1795  myGlobal.StdOut("Tests done.\n");
1796  parallelUnitResults.UpdateResult("Grid:PBS:MetricCancellationVGWavy",testResults[1]);
1797 };
1798 
1799 
1800 void TestGrid_PBS_IntegratedHalo(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
1801 {
1802 
1803  std::ostringstream messageStream;
1804 
1805  fixtures::ConfigurationType testConfig;
1806 
1807  // typedef simulation::grid::parallel_uniform_blockstructured grid_t;
1808  typedef simulation::grid::parallel_blockstructured pbsgrid_t;
1809  typedef plascom2::operators::sbp::base operator_t;
1810  typedef simulation::grid::halo halo_t;
1811  typedef pcpp::ParallelGlobalType global_t;
1812  typedef pcpp::CommunicatorType comm_t;
1813  typedef pcpp::IndexIntervalType interval_t;
1814 
1815  global_t myGlobal;
1816  myGlobal.Init("TestIntegratedHalo",testComm);
1817  myGlobal.SetVerbLevel(3);
1818  myGlobal.Profiling(true);
1819 
1820  int numProc = testComm.Size();
1821 
1822  // grid_t testGrid;
1823  pbsgrid_t pbsGrid;
1824  // halo_t testHalo;
1825  // comm_t &gridComm(testGrid.Communicator());
1826  comm_t &pbsGridComm(pbsGrid.Communicator());
1827 
1828  int numDim = 3;
1829  std::vector<size_t> gridSize(numDim,0);
1830  std::vector<double> dX(numDim,0);
1831 
1832  size_t numGlobalNodes = 1;
1833  size_t numGlobalCells = 1;
1834  for(int iDim = 0;iDim < numDim;iDim++){
1835  // 100x50x25
1836  gridSize[iDim] = static_cast<size_t>(std::pow(2.0,(numDim-iDim-1)%numDim)*25);
1837  dX[iDim] = 1.0/static_cast<double>(gridSize[iDim]-1);
1838  numGlobalNodes *= gridSize[iDim];
1839  numGlobalCells *= (gridSize[iDim]-1);
1840  }
1841  // testGrid.SetGlobalSize(gridSize);
1842  pbsGrid.SetGridSizes(gridSize);
1843 
1844  // testGrid.SetGridSpacing(dX);
1845 
1846  operator_t sbpOperator;
1847  plascom2::operators::sbp::Initialize(sbpOperator,4);
1848  int boundaryDepth = (sbpOperator.numStencils-1)/2;
1849  std::vector<int> haloDepths(2*numDim,boundaryDepth);
1850  std::vector<int> boundaryDepths(2*numDim,boundaryDepth);
1851  // testGrid.SetHaloDepth(haloDepths);
1852  // testGrid.SetBoundaryDepth(boundaryDepths);
1853 
1854  // Set up parallel Cartesian topology
1855  // references to some grid-specific data structures
1856  pcpp::ParallelTopologyInfoType cartInfo;
1857  cartInfo.numDimensions = numDim;
1858  pcpp::ParallelTopologyInfoType pbsCartInfo;
1859  pbsCartInfo.numDimensions = numDim;
1860  pbsCartInfo.isPeriodic.resize(numDim,1);
1861  // pcpp::comm::SetupCartesianTopology(testComm,cartInfo,gridComm,messageStream);
1862  pcpp::comm::SetupCartesianTopology(testComm,pbsCartInfo,pbsGridComm,messageStream);
1863 
1864  // Grab the local coordinates and global dimension of the Cartesian topo
1865  std::vector<int> &cartCoords(pbsGridComm.CartCoordinates());
1866  std::vector<int> &cartDims(pbsGridComm.CartDimensions());
1867 
1868  // Grab the +/- neighbors in each Cartesian direction
1869  std::vector<int> cartNeighbors;
1870  pbsGridComm.CartNeighbors(cartNeighbors);
1871 
1872  // Generate a report of the Cartesian setup and put it to stdout
1873  pcpp::report::CartesianSetup(messageStream,pbsCartInfo);
1874  myGlobal.StdOut(messageStream.str(),2);
1875  pcpp::io::RenewStream(messageStream);
1876 
1877  interval_t globalInterval;
1878  globalInterval.InitSimple(gridSize);
1879 
1880  // interval_t &partitionInterval(pbsGrid.PartitionInterval());
1881  interval_t &pbsPartInterval(pbsGrid.PartitionInterval());
1882 
1883  // === Partition the grid ===
1884  // if(pcpp::util::PartitionCartesianInterval(globalInterval,cartDims,cartCoords,partitionInterval,messageStream)){
1885  // myGlobal.ErrOut("Partitioning failed.\n");
1886  // return;
1887  // }
1888  if(pcpp::util::PartitionCartesianInterval(globalInterval,cartDims,cartCoords,pbsPartInterval,messageStream)){
1889  myGlobal.ErrOut("Partitioning failed.\n");
1890  return;
1891  }
1892  std::vector<size_t> extendGrid(2*numDim,0);
1893  std::vector<bool> isPeriodic(numDim,true);
1894  size_t numPointsPart = pbsPartInterval.NNodes();
1895  for(int iDim = 0;iDim < numDim; iDim++){
1896  if(isPeriodic[iDim]){
1897  extendGrid[2*iDim] = haloDepths[2*iDim];
1898  extendGrid[2*iDim+1] = haloDepths[2*iDim+1];
1899  } else {
1900  if(pbsPartInterval[iDim].first == 0)
1901  extendGrid[2*iDim] = haloDepths[2*iDim];
1902  if(pbsPartInterval[iDim].second == (gridSize[iDim]-1))
1903  extendGrid[2*iDim + 1] = haloDepths[2*iDim+1];
1904  }
1905  }
1906  pbsGrid.SetDimensionExtensions(extendGrid);
1907  pbsGrid.Finalize();
1908 
1909  const std::vector<size_t> &bufferSizes(pbsGrid.BufferSizes());
1910  const pcpp::IndexIntervalType &localPartitionInterval(pbsGrid.PartitionBufferInterval());
1911 
1912  size_t numPointsBuffer = 1;
1913  std::vector<size_t>::const_iterator bsIt = bufferSizes.begin();
1914  while(bsIt != bufferSizes.end())
1915  numPointsBuffer *= *bsIt++;
1916 
1917  bool bufferSizeTest = true;
1918  if(numPointsBuffer != pbsGrid.BufferSize())
1919  bufferSizeTest = false;
1920  parallelUnitResults.UpdateResult("Grid:PBS:BufferSize",bufferSizeTest);
1921 
1922  std::cout << "Number of Points in Partition: " << numPointsPart << std::endl;
1923  std::cout << "Buffer Size: (";
1924  pcpp::io::DumpContents(std::cout,bufferSizes,",");
1925  std::cout << ")" << std::endl;
1926  bool testResult = true;
1927  for(int iDim = 0;iDim < numDim;iDim++){
1928  if(bufferSizes[iDim] != (pbsPartInterval.NNodes(iDim) + extendGrid[2*iDim] + extendGrid[2*iDim + 1]))
1929  testResult = false;
1930  }
1931  parallelUnitResults.UpdateResult("Grid:PBS:GridSizing",testResult);
1932 
1933  // Set up the HALO data structure
1934  halo_t testHalo(globalInterval,pbsPartInterval);
1935 
1936  // This forces the halo to have neighbors even on
1937  // partitions with domain boundaries.
1938  std::vector<bool> haveNeighbors(2*numDim,true);
1939  testHalo.SetNeighbors(haveNeighbors);
1940 
1941  std::vector<pcpp::IndexIntervalType> remoteHaloExtents(testHalo.CreateRemoteHaloExtents(extendGrid));
1942  std::vector<pcpp::IndexIntervalType> localHaloExtents(testHalo.CreateLocalHaloExtents(extendGrid));
1943 
1944  testHalo.SetLocalBufferSizes(bufferSizes);
1945  testHalo.SetLocalPartitionExtent(localPartitionInterval);
1946 
1947  testHalo.SetRemoteHaloExtents(remoteHaloExtents);
1948  testHalo.SetLocalHaloExtents(localHaloExtents);
1949 
1950  int myRank = testComm.Rank();
1951  int numRemoteHalos = remoteHaloExtents.size();
1952  int numLocalHalos = localHaloExtents.size();
1953  for(int iProc = 0;iProc < numProc;iProc++){
1954  if(iProc == myRank) {
1955  std::cout << "Processor " << iProc << " Summary " << std::endl
1956  << "+++++++++++++++++++++++++++++++++++" << std::endl;
1957  std::cout << "Partition Interval: ";
1958  pbsPartInterval.PrettyPrint(std::cout);
1959  std::cout << std::endl;
1960  std::cout << "Number of Remote Halos: " << numRemoteHalos << std::endl;
1961  std::cout << "Number of Local Halos: " << numLocalHalos << std::endl;
1962  for(int iHalo = 0;iHalo < numRemoteHalos;iHalo++){
1963  std::cout << "Remote Halo(" << iHalo << "):";
1964  remoteHaloExtents[iHalo].PrettyPrint(std::cout);
1965  std::cout << std::endl;
1966  std::cout << "Local Halo(" << iHalo << "):";
1967  localHaloExtents[iHalo].PrettyPrint(std::cout);
1968  std::cout << std::endl;
1969  }
1970  }
1971  testComm.Barrier();
1972  }
1973 
1974  const std::vector<pcpp::IndexIntervalType> &remoteHaloBufferExtents(testHalo.RemoteHaloBufferExtents());
1975  const std::vector<pcpp::IndexIntervalType> &localHaloBufferExtents(testHalo.LocalHaloBufferExtents());
1976 
1977  bool haloBufferExtentTest = true;
1978  for(int iProc = 0;iProc < numProc;iProc++){
1979  if(iProc == myRank) {
1980  std::cout << "Processor " << iProc << " Buffer Summary " << std::endl
1981  << "+++++++++++++++++++++++++++++++++++" << std::endl;
1982  std::cout << "Partition Buffer Interval: ";
1983  localPartitionInterval.PrettyPrint(std::cout);
1984  std::cout << std::endl;
1985 
1986  for(int iHalo = 0;iHalo < numRemoteHalos;iHalo++){
1987 
1988  int ijkDirection = iHalo/2;
1989  int leftRight = iHalo%2;
1990 
1991  const pcpp::IndexIntervalType &remoteHaloBufferExtent(remoteHaloBufferExtents[iHalo]);
1992  const pcpp::IndexIntervalType &localHaloBufferExtent(localHaloBufferExtents[iHalo]);
1993 
1994  std::cout << "Remote Halo Buffer(" << iHalo << "):";
1995  remoteHaloBufferExtent.PrettyPrint(std::cout);
1996  std::cout << std::endl;
1997  std::cout << "Local Halo Buffer(" << iHalo << "):";
1998  localHaloBufferExtents[iHalo].PrettyPrint(std::cout);
1999  std::cout << std::endl;
2000 
2001  for(int iDim = 0;iDim < numDim;iDim++){
2002  if(iDim == ijkDirection){
2003  if(leftRight == 0) { // left halo
2004  if(remoteHaloBufferExtent[iDim].first != 0)
2005  haloBufferExtentTest = false;
2006  if(remoteHaloBufferExtent[iDim].second != (extendGrid[iHalo] - 1))
2007  haloBufferExtentTest = false;
2008  if(localHaloBufferExtent[iDim].first != localPartitionInterval[iDim].first)
2009  haloBufferExtentTest = false;
2010  if(localHaloBufferExtent[iDim].second != (2*extendGrid[iHalo] - 1))
2011  haloBufferExtentTest = false;
2012  } else { // right halo
2013  if(remoteHaloBufferExtent[iDim].second != (bufferSizes[iDim]-1))
2014  haloBufferExtentTest = false;
2015  if(remoteHaloBufferExtent[iDim].first != (bufferSizes[iDim]-extendGrid[iHalo]))
2016  haloBufferExtentTest = false;
2017  if(localHaloBufferExtent[iDim].second != localPartitionInterval[iDim].second)
2018  haloBufferExtentTest = false;
2019  if(localHaloBufferExtent[iDim].first != (localPartitionInterval[iDim].second - extendGrid[iHalo] + 1))
2020  haloBufferExtentTest = false;
2021  }
2022  } else { // some other direction
2023  if(remoteHaloBufferExtent[iDim] != localPartitionInterval[iDim])
2024  haloBufferExtentTest = false;
2025  if(localHaloBufferExtent[iDim] != localPartitionInterval[iDim])
2026  haloBufferExtentTest = false;
2027  }
2028  }
2029  }
2030  }
2031  testComm.Barrier();
2032  }
2033 
2034  haloBufferExtentTest = CheckResult(haloBufferExtentTest,testComm);
2035 
2036  parallelUnitResults.UpdateResult("Grid:PBS:HaloExtents",haloBufferExtentTest);
2037 
2038  testHalo.CreateSimpleSendIndices();
2039  const std::vector<std::vector<size_t> > &sendIndices(testHalo.SendIndices());
2040  bool sendIndicesTest = true;
2041  pcpp::IndexIntervalType bufferInterval;
2042  bufferInterval.InitSimple(bufferSizes);
2043 
2044  if(sendIndices.size() != (2*numDim))
2045  sendIndicesTest = false;
2046  for(int iDim = 0;iDim < numDim;iDim++){
2047 
2048 
2049  int leftDir = 2*iDim;
2050  int rightDir = 2*iDim + 1;
2051  size_t numPlane = 1;
2052 
2053  for(int sDim = 1;sDim <= numDim-1;sDim++){
2054  int orthoDim = (iDim+sDim)%numDim;
2055  numPlane *= pbsPartInterval.NNodes(orthoDim);
2056  }
2057 
2058  int leftPointsExpected = numPlane*extendGrid[2*iDim];
2059  int rightPointsExpected = numPlane*extendGrid[2*iDim+1];
2060 
2061  const std::vector<size_t> &leftSendIndices(sendIndices[leftDir]);
2062  const std::vector<size_t> &rightSendIndices(sendIndices[rightDir]);
2063 
2064  int numPointsSendLeft = leftSendIndices.size();
2065  int numPointsSendRight = leftSendIndices.size();
2066 
2067  if(numPointsSendLeft != leftPointsExpected)
2068  sendIndicesTest = false;
2069  if(numPointsSendRight != rightPointsExpected)
2070  sendIndicesTest = false;
2071 
2072  std::vector<size_t> expectedLeftIndices;
2073  bufferInterval.GetFlatIndices(localHaloBufferExtents[leftDir],expectedLeftIndices);
2074  if(leftSendIndices != expectedLeftIndices)
2075  sendIndicesTest = false;
2076 
2077  std::vector<size_t> expectedRightIndices;
2078  bufferInterval.GetFlatIndices(localHaloBufferExtents[rightDir],expectedRightIndices);
2079  if(rightSendIndices != expectedRightIndices)
2080  sendIndicesTest = false;
2081  }
2082 
2083  sendIndicesTest = CheckResult(sendIndicesTest,testComm);
2084  parallelUnitResults.UpdateResult("Grid:PBS:HaloSendIndices",sendIndicesTest);
2085 
2086  testHalo.CreateSimpleRecvIndices();
2087  const std::vector<std::vector<size_t> > &recvIndices(testHalo.RecvIndices());
2088  bool recvIndicesTest = true;
2089 
2090  if(recvIndices.size() != (2*numDim))
2091  recvIndicesTest = false;
2092 
2093  for(int iDim = 0;iDim < numDim;iDim++){
2094 
2095 
2096  int leftDir = 2*iDim;
2097  int rightDir = 2*iDim + 1;
2098 
2099  size_t numPlane = 1;
2100  for(int sDim = 1;sDim <= numDim-1;sDim++){
2101  int orthoDim = (iDim+sDim)%numDim;
2102  numPlane *= pbsPartInterval.NNodes(orthoDim);
2103  }
2104 
2105  int leftPointsExpected = numPlane*extendGrid[2*iDim];
2106  int rightPointsExpected = numPlane*extendGrid[2*iDim+1];
2107 
2108  const std::vector<size_t> &leftRecvIndices(recvIndices[leftDir]);
2109  const std::vector<size_t> &rightRecvIndices(recvIndices[rightDir]);
2110 
2111  int numPointsRecvLeft = leftRecvIndices.size();
2112  int numPointsRecvRight = leftRecvIndices.size();
2113 
2114  if(numPointsRecvLeft != leftPointsExpected)
2115  recvIndicesTest = false;
2116  if(numPointsRecvRight != rightPointsExpected)
2117  recvIndicesTest = false;
2118 
2119  std::vector<size_t> expectedLeftIndices;
2120  bufferInterval.GetFlatIndices(remoteHaloBufferExtents[leftDir],expectedLeftIndices);
2121  if(leftRecvIndices != expectedLeftIndices)
2122  recvIndicesTest = false;
2123 
2124  std::vector<size_t> expectedRightIndices;
2125  bufferInterval.GetFlatIndices(remoteHaloBufferExtents[rightDir],expectedRightIndices);
2126  if(rightRecvIndices != expectedRightIndices)
2127  recvIndicesTest = false;
2128  }
2129 
2130  recvIndicesTest = CheckResult(recvIndicesTest,testComm);
2131  parallelUnitResults.UpdateResult("Grid:PBS:HaloRecvIndices",recvIndicesTest);
2132 
2133  // Create simple message buffers - assuming the grid coordinates will
2134  // be used to test
2135  int numComponents = numDim;
2136  int messageId = testHalo.CreateMessageBuffers(numComponents);
2137  bool createMessageBuffersTest = true;
2138  if(messageId != 0)
2139  createMessageBuffersTest = false;
2140  std::vector<halo_t::messagebuffers> &communicationBuffers(testHalo.CommunicationBuffers());
2141  if(communicationBuffers.size() != 1)
2142  createMessageBuffersTest = false;
2143  halo_t::messagebuffers &messageBuffers(communicationBuffers[0]);
2144  if(messageBuffers.numComponents != numComponents)
2145  createMessageBuffersTest = false;
2146  if(messageBuffers.sendBuffers.size() != 2*numDim)
2147  createMessageBuffersTest = false;
2148  if(messageBuffers.recvBuffers.size() != 2*numDim)
2149  createMessageBuffersTest = false;
2150 
2151  createMessageBuffersTest = CheckResult(createMessageBuffersTest,testComm);
2152  parallelUnitResults.UpdateResult("Grid:PBS:CreateMessageBuffers",createMessageBuffersTest);
2153 
2154  bool packSendBuffersTest = true;
2155  std::vector<double> gridCoordinates(numDim*numPointsBuffer,0.0);
2156  // pbsPartInterval: interval wrt global grid size
2157  // localPartitionInterval: interval wrt local buffer size
2158  // 3D-specific test!!!
2159  if(numDim == 3) {
2160  size_t partStartX = pbsPartInterval[0].first;
2161  size_t partEndX = pbsPartInterval[0].second;
2162  size_t partStartY = pbsPartInterval[1].first;
2163  size_t partEndY = pbsPartInterval[1].second;
2164  size_t partStartZ = pbsPartInterval[2].first;
2165  size_t partEndZ = pbsPartInterval[2].second;
2166  size_t bufferStartX = localPartitionInterval[0].first;
2167  size_t bufferEndX = localPartitionInterval[0].second;
2168  size_t bufferStartY = localPartitionInterval[1].first;
2169  size_t bufferEndY = localPartitionInterval[1].second;
2170  size_t bufferStartZ = localPartitionInterval[2].first;
2171  size_t bufferEndZ = localPartitionInterval[2].second;
2172  size_t buffPlane = bufferSizes[0]*bufferSizes[1];
2173  size_t buffRow = bufferSizes[0];
2174  size_t iPartZ,iPartY,iPartX,iBufX, iBufY, iBufZ;
2175  for(iPartZ = partStartZ, iBufZ = bufferStartZ;
2176  iPartZ <= partEndZ;iPartZ++,iBufZ++){
2177  size_t bufZ = iBufZ*buffPlane;
2178  for(iPartY = partStartY, iBufY = bufferStartY;
2179  iPartY <= partEndY;iPartY++,iBufY++){
2180  size_t bufYZ = bufZ + iBufY*buffRow;
2181  for(iPartX = partStartX, iBufX = bufferStartX;
2182  iPartX <= partEndX;iPartX++,iBufX++){
2183  size_t bufIndex = bufYZ + iBufX;
2184  gridCoordinates[bufIndex] = iPartX;
2185  gridCoordinates[bufIndex+numPointsBuffer] = iPartY;
2186  gridCoordinates[bufIndex+2*numPointsBuffer] = iPartZ;
2187  }
2188  }
2189  }
2190  testHalo.PackMessageBuffers(messageId,0,3,&gridCoordinates[0]);
2191  }
2192 
2193  for(int iDim = 0;iDim < numDim;iDim++){
2194 
2195  int leftDir = 2*iDim;
2196  int rightDir = 2*iDim+1;
2197 
2198  size_t numPlane = 1;
2199  for(int sDim = 1;sDim <= numDim-1;sDim++){
2200  int orthoDim = (iDim+sDim)%numDim;
2201  numPlane *= pbsPartInterval.NNodes(orthoDim);
2202  }
2203 
2204  int leftPointsExpected = numPlane*extendGrid[leftDir];
2205  int rightPointsExpected = numPlane*extendGrid[rightDir];
2206 
2207  double *leftSendBuffer = messageBuffers.sendBuffers[leftDir];
2208  double *rightSendBuffer = messageBuffers.sendBuffers[rightDir];
2209 
2210  int numLeftSendPoints = messageBuffers.numPointsSend[leftDir];
2211  int numRightSendPoints = messageBuffers.numPointsSend[rightDir];
2212 
2213  if(leftPointsExpected != numLeftSendPoints)
2214  packSendBuffersTest = false;
2215  if(numRightSendPoints != rightPointsExpected)
2216  packSendBuffersTest = false;
2217 
2218  // Make sure packed data is right
2219  if(numDim == 3){
2220  size_t iStart = localHaloExtents[leftDir][0].first;
2221  size_t iEnd = localHaloExtents[leftDir][0].second;
2222  size_t jStart = localHaloExtents[leftDir][1].first;
2223  size_t jEnd = localHaloExtents[leftDir][1].second;
2224  size_t kStart = localHaloExtents[leftDir][2].first;
2225  size_t kEnd = localHaloExtents[leftDir][2].second;
2226  size_t iPoint = 0;
2227  for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
2228  for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
2229  for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
2230  if(leftSendBuffer[iPoint] != iIndex)
2231  packSendBuffersTest = false;
2232  if(leftSendBuffer[iPoint + numLeftSendPoints] != jIndex)
2233  packSendBuffersTest = false;
2234  if(leftSendBuffer[iPoint + 2*numLeftSendPoints] != kIndex)
2235  packSendBuffersTest = false;
2236  iPoint++;
2237  }
2238  }
2239  }
2240 
2241  iStart = localHaloExtents[rightDir][0].first;
2242  iEnd = localHaloExtents[rightDir][0].second;
2243  jStart = localHaloExtents[rightDir][1].first;
2244  jEnd = localHaloExtents[rightDir][1].second;
2245  kStart = localHaloExtents[rightDir][2].first;
2246  kEnd = localHaloExtents[rightDir][2].second;
2247  iPoint = 0;
2248  for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
2249  for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
2250  for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
2251  if(rightSendBuffer[iPoint] != iIndex)
2252  packSendBuffersTest = false;
2253  if(rightSendBuffer[iPoint + numRightSendPoints] != jIndex)
2254  packSendBuffersTest = false;
2255  if(rightSendBuffer[iPoint + 2*numRightSendPoints] != kIndex)
2256  packSendBuffersTest = false;
2257  iPoint++;
2258  }
2259  }
2260  }
2261  }
2262  }
2263 
2264 
2265  // Send Test w/Recv Buffer check
2266  testHalo.SendMessage(0,cartNeighbors,pbsGridComm);
2267  testHalo.ReceiveMessage(0,cartNeighbors,pbsGridComm);
2268  testHalo.UnPackMessageBuffers(0,0,3,&gridCoordinates[0]);
2269 
2270  bool recvDataTest = true;
2271  for(int iDim = 0;iDim < numDim;iDim++){
2272 
2273  int leftDir = 2*iDim;
2274  int rightDir = 2*iDim+1;
2275 
2276  size_t numPlane = 1;
2277  for(int sDim = 1;sDim <= numDim-1;sDim++){
2278  int orthoDim = (iDim+sDim)%numDim;
2279  numPlane *= pbsPartInterval.NNodes(orthoDim);
2280  }
2281 
2282  int leftPointsExpected = numPlane*extendGrid[leftDir];
2283  int rightPointsExpected = numPlane*extendGrid[rightDir];
2284 
2285  double *leftRecvBuffer = messageBuffers.recvBuffers[leftDir];
2286  double *rightRecvBuffer = messageBuffers.recvBuffers[rightDir];
2287 
2288  int numLeftRecvPoints = messageBuffers.numPointsRecv[leftDir];
2289  int numRightRecvPoints = messageBuffers.numPointsRecv[rightDir];
2290 
2291  if(leftPointsExpected != numLeftRecvPoints)
2292  recvDataTest = false;
2293  if(numRightRecvPoints != rightPointsExpected)
2294  recvDataTest = false;
2295 
2296  // Make sure recvd data is right
2297  if(numDim == 3){
2298  size_t iStart = remoteHaloExtents[leftDir][0].first;
2299  size_t iEnd = remoteHaloExtents[leftDir][0].second;
2300  size_t jStart = remoteHaloExtents[leftDir][1].first;
2301  size_t jEnd = remoteHaloExtents[leftDir][1].second;
2302  size_t kStart = remoteHaloExtents[leftDir][2].first;
2303  size_t kEnd = remoteHaloExtents[leftDir][2].second;
2304 
2305  size_t ibStart = remoteHaloBufferExtents[leftDir][0].first;
2306  size_t ibEnd = remoteHaloBufferExtents[leftDir][0].second;
2307  size_t jbStart = remoteHaloBufferExtents[leftDir][1].first;
2308  size_t jbEnd = remoteHaloBufferExtents[leftDir][1].second;
2309  size_t kbStart = remoteHaloBufferExtents[leftDir][2].first;
2310  size_t kbEnd = remoteHaloBufferExtents[leftDir][2].second;
2311 
2312  size_t iPoint = 0;
2313  size_t kIndex, kbIndex, jIndex, jbIndex, iIndex,ibIndex;
2314  for(kIndex = kStart,kbIndex = kbStart;kIndex <= kEnd;kIndex++,kbIndex++){
2315  for(jIndex = jStart, jbIndex = jbStart;jIndex <= jEnd;jIndex++,jbIndex++){
2316  for(iIndex = iStart,ibIndex = ibStart;iIndex <= iEnd;iIndex++,ibIndex++){
2317  size_t bufferIndex = (kbIndex * bufferSizes[0]*bufferSizes[1]) +
2318  (jbIndex * bufferSizes[0]) + ibIndex;
2319 
2320  if(leftRecvBuffer[iPoint] != iIndex)
2321  recvDataTest = false;
2322  if(leftRecvBuffer[iPoint + numLeftRecvPoints] != jIndex)
2323  recvDataTest = false;
2324  if(leftRecvBuffer[iPoint + 2*numLeftRecvPoints] != kIndex)
2325  recvDataTest = false;
2326 
2327  if(gridCoordinates[bufferIndex] != iIndex)
2328  recvDataTest = false;
2329  if(gridCoordinates[bufferIndex + numPointsBuffer] != jIndex)
2330  recvDataTest = false;
2331  if(gridCoordinates[bufferIndex + 2*numPointsBuffer] != kIndex)
2332  recvDataTest = false;
2333 
2334  iPoint++;
2335  }
2336  }
2337  }
2338 
2339  iStart = remoteHaloExtents[rightDir][0].first;
2340  iEnd = remoteHaloExtents[rightDir][0].second;
2341  jStart = remoteHaloExtents[rightDir][1].first;
2342  jEnd = remoteHaloExtents[rightDir][1].second;
2343  kStart = remoteHaloExtents[rightDir][2].first;
2344  kEnd = remoteHaloExtents[rightDir][2].second;
2345 
2346  ibStart = remoteHaloBufferExtents[rightDir][0].first;
2347  ibEnd = remoteHaloBufferExtents[rightDir][0].second;
2348  jbStart = remoteHaloBufferExtents[rightDir][1].first;
2349  jbEnd = remoteHaloBufferExtents[rightDir][1].second;
2350  kbStart = remoteHaloBufferExtents[rightDir][2].first;
2351  kbEnd = remoteHaloBufferExtents[rightDir][2].second;
2352 
2353  iPoint = 0;
2354 
2355  for(kIndex = kStart,kbIndex = kbStart;kIndex <= kEnd;kIndex++,kbIndex++){
2356  for(jIndex = jStart, jbIndex = jbStart;jIndex <= jEnd;jIndex++,jbIndex++){
2357  for(iIndex = iStart,ibIndex = ibStart;iIndex <= iEnd;iIndex++,ibIndex++){
2358  size_t bufferIndex = (kbIndex * bufferSizes[0]*bufferSizes[1]) +
2359  (jbIndex * bufferSizes[0]) + ibIndex;
2360 
2361  if(rightRecvBuffer[iPoint] != iIndex)
2362  recvDataTest = false;
2363  if(rightRecvBuffer[iPoint + numRightRecvPoints] != jIndex)
2364  recvDataTest = false;
2365  if(rightRecvBuffer[iPoint + 2*numRightRecvPoints] != kIndex)
2366  recvDataTest = false;
2367 
2368  if(gridCoordinates[bufferIndex] != iIndex)
2369  recvDataTest = false;
2370  if(gridCoordinates[bufferIndex + numPointsBuffer] != jIndex)
2371  recvDataTest = false;
2372  if(gridCoordinates[bufferIndex + 2*numPointsBuffer] != kIndex)
2373  recvDataTest = false;
2374 
2375  iPoint++;
2376  }
2377  }
2378  }
2379  }
2380  }
2381 
2382  recvDataTest = CheckResult(recvDataTest,testComm);
2383  parallelUnitResults.UpdateResult("Grid:PBS:RecvHaloData",recvDataTest);
2384 
2385 
2386  myGlobal.Finalize();
2387 
2388 }
gridMetric
void const size_t const size_t const size_t const int const int const double * gridMetric
Definition: EulerKernels.H:10
pcpp::io::xdmf::OpenFileTag
void OpenFileTag(std::ostream &outStream)
Definition: PCPPHDF5.C:12
CLSIZE2
#define CLSIZE2
Definition: TestGrid.C:779
simulation::grid::parallel_blockstructured::Jacobian
std::vector< double > & Jacobian()
Definition: Grid.H:624
simulation::grid::halo::LocalHaloBufferExtents
const std::vector< pcpp::IndexIntervalType > & LocalHaloBufferExtents() const
Definition: Grid.H:155
pcpp::io::xdmf::CloseFileTag
void CloseFileTag(std::ostream &outStream)
Definition: PCPPHDF5.C:20
plascom2::state_t
simulation::state::base state_t
Definition: PlasCom2.H:17
simulation::grid::subregion::regionInterval
pcpp::IndexIntervalType regionInterval
Definition: Grid.H:39
ix::util::sizeextent::GetFlatIndices
void GetFlatIndices(const sizeextent &extent, ContainerType &indices) const
Definition: IndexUtil.H:302
simulation::grid::parallel_blockstructured::Metric
std::vector< double > & Metric()
Definition: Grid.H:622
simulation::grid::parallel_blockstructured::ExchangeCoordinates
int ExchangeCoordinates()
Definition: Grid.C:2075
simulation::grid::halo::CreateSimpleRecvIndices
int CreateSimpleRecvIndices()
Definition: Grid.C:2775
simulation::grid::halo::CreateLocalHaloExtents
std::vector< pcpp::IndexIntervalType > CreateLocalHaloExtents(const pcpp::IndexIntervalType &globalExtent, const pcpp::IndexIntervalType &partitionExtent, std::vector< size_t > &haloSizes)
Definition: Grid.C:2592
ix::sys::Remove
int Remove(const std::string &fname)
Definition: UnixUtils.C:141
ix::global::GlobalObj::ErrOut
virtual int ErrOut(const std::string &outstr)
Definition: Global.H:456
GenerateCurvilinearMetric1
int GenerateCurvilinearMetric1(std::vector< size_t > &ijk, std::vector< double > &met, std::vector< double > &jac, double &jacm1)
Definition: TestGrid.C:784
TestGrid_CurvilinearVGWavy
void TestGrid_CurvilinearVGWavy(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
Definition: TestGrid.C:1470
simulation::grid::halo::RecvIndices
const std::vector< std::vector< size_t > > & RecvIndices() const
Definition: Grid.H:159
simulation::grid::halo::CreateSimpleSendIndices
int CreateSimpleSendIndices()
Creates the flat buffer indices for each whole local halo extent.
Definition: Grid.C:2659
operator_t
plascom2::operators::sbp::base operator_t
Definition: TestMaxwellSolver.C:1568
ix::global::GlobalObj::Profiling
virtual bool Profiling()
Get profiling state.
Definition: Global.H:229
GenerateCurvilinearMetric2
int GenerateCurvilinearMetric2(std::vector< size_t > &ijk, std::vector< double > &xyz, double &jacm1)
Definition: TestGrid.C:910
simulation::grid::parallel_blockstructured::PartitionInterval
pcpp::IndexIntervalType & PartitionInterval()
Definition: Grid.H:500
simulation::grid::halo::CommunicationBuffers
std::vector< messagebuffers > & CommunicationBuffers()
Definition: Grid.H:157
pcpp::io::xdmf::CloseGridTag
void CloseGridTag(std::ostream &outStream)
Definition: PCPPHDF5.C:32
simulation::grid::parallel_blockstructured::GenerateCoordinates
int GenerateCoordinates(std::ostream &)
Definition: Grid.C:1628
gridSizes
void const size_t const size_t * gridSizes
Definition: EulerKernels.H:10
ix::util::RenewStream
void RenewStream(std::ostringstream &outStream)
Definition: primitive_utilities.H:638
pcpp::comm::CheckResult
bool CheckResult(bool &localResult, pcpp::CommunicatorType &testComm)
Definition: PCPPCommUtil.C:176
simulation::grid::halo::messagebuffers::recvBuffers
std::vector< double * > recvBuffers
Definition: Grid.H:67
y
void const size_t const size_t const size_t const double const double double * y
Definition: OperatorKernels.H:33
simulation::grid::parallel_blockstructured::SetDimensionExtensions
void SetDimensionExtensions(const std::vector< size_t > &inExtensions)
Definition: Grid.H:492
simulation::grid::halo::SendMessage
int SendMessage(const int messageId, const std::vector< int > &neighborRanks, fixtures::CommunicatorType &inComm)
Definition: Grid.C:3061
plascom2::operators::stencilset::stencilStarts
int * stencilStarts
The starting index into the stencilWeight and stencilOffset arrays for each stencil.
Definition: Stencil.H:104
simulation::grid::halo::UnPackMessageBuffers
int UnPackMessageBuffers(const int messageId, const int componentId, const int numComponents, double *targetBuffer)
Unpack generic recv buffers for specified message.
Definition: Grid.C:3145
comm_t
pcpp::CommunicatorType comm_t
Definition: TestRK4Advancer.C:17
pcpp::report::CartesianSetup
void CartesianSetup(std::ostream &outStream, const std::vector< int > &cartCoords, const std::vector< int > &cartDims)
Definition: PCPPReport.C:199
pcpp::field::dataset::Create
virtual size_t Create(size_t number_of_nodes=0, size_t number_of_cells=0)
Definition: PCPPFieldData.H:578
opInterval
void const size_t const size_t const size_t * opInterval
Definition: EulerKernels.H:10
GenerateRectilinearMetric1
int GenerateRectilinearMetric1(std::vector< size_t > &ijk, std::vector< double > &xyz, double &jacm1)
Definition: TestGrid.C:417
simulation::grid::halo::messagebuffers::numPointsSend
std::vector< size_t > numPointsSend
Definition: Grid.H:64
interval_t
pcpp::IndexIntervalType interval_t
Definition: TestMaxwellSolver.C:1567
simulation::grid::parallel_blockstructured::SetGridSizes
void SetGridSizes(const std::vector< size_t > &inSize)
Definition: Grid.H:452
global_t
pcpp::ParallelGlobalType global_t
Definition: TestHDF5.C:16
TestGrid_CartesianMetric
void TestGrid_CartesianMetric(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
Definition: TestGrid.C:229
ix::global::ParallelGlobalObj::Init
virtual int Init(const std::string &name, CommunicatorType &incomm)
Definition: Global.H:578
plascom2::io::hdf5::OutputSingle
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
x
void const size_t const size_t const size_t const double const double * x
Definition: OperatorKernels.H:33
pcpp::io::hdf5::WriteGrid
int WriteGrid(const GridType &inGrid, const std::string &gridName, const std::string &parentPath, base &hdfFile)
Definition: PCPPHDF5.H:1131
simulation::grid::parallel_blockstructured::BufferSizes
const std::vector< size_t > & BufferSizes() const
Definition: Grid.H:459
simulation::grid::halo::SendIndices
const std::vector< std::vector< size_t > > & SendIndices() const
Definition: Grid.H:158
simulation::grid::halo::RemoteHaloBufferExtents
const std::vector< pcpp::IndexIntervalType > & RemoteHaloBufferExtents() const
Definition: Grid.H:154
pcpp::io::xdmf::OpenGridTag
void OpenGridTag(const std::string &gridName, const std::string &gridType, double inTime, std::ostream &outStream)
Definition: PCPPHDF5.C:26
simulation::grid::InitSubRegionFromString
int InitSubRegionFromString(const std::string &inString, const std::vector< size_t > &gridSizes, subregion &inSubRegion)
Definition: Grid.C:174
ix::test::results
Encapsulating class for collections of test results.
Definition: Testing.H:18
simulation::grid::parallel_blockstructured::SetPeriodicLengths
void SetPeriodicLengths(const std::vector< double > &inLengths)
Definition: Grid.H:588
simulation::grid::halo::CreateMessageBuffers
int CreateMessageBuffers(int numComponents, int componentSize=8)
Definition: Grid.C:2850
CLSIZE1
#define CLSIZE1
Definition: TestGrid.C:778
pcpp::comm::paralleltopologyinfo::isPeriodic
std::vector< int > isPeriodic
Definition: PCPPCommUtil.H:24
simulation::grid::parallel_blockstructured::SetupDifferentialOperator
int SetupDifferentialOperator(plascom2::operators::stencilset *inStencilSet, int *inStencilConnectivity)
Definition: Grid.C:1619
simulation::grid::parallel_blockstructured::ComputeJacobianMatrix
int ComputeJacobianMatrix(std::ostream &messageStream)
Calculates the regular/naive Jacobian matrix d(xyz)/d(ijk)
Definition: Grid.C:916
simulation::grid::parallel_blockstructured::PartitionBufferInterval
pcpp::IndexIntervalType & PartitionBufferInterval()
Definition: Grid.H:519
simulation::grid::parallel_blockstructured::GridSizes
const std::vector< size_t > & GridSizes() const
Definition: Grid.H:469
simulation::grid::halo::PackMessageBuffers
int PackMessageBuffers(const int messageId, const int componentId, const int numComponents, const double *sourceBuffer)
Pack generic send buffers for specified message.
Definition: Grid.C:2937
simulation::grid::halo::SetLocalPartitionExtent
void SetLocalPartitionExtent(const pcpp::IndexIntervalType &inExtent)
Definition: Grid.H:90
TestGrid_SubRegion
void TestGrid_SubRegion(ix::test::results &serialUnitResults)
Definition: TestGrid.C:18
simulation::grid::parallel_blockstructured::ComputeMetrics
int ComputeMetrics(std::ostream &messageStream)
Definition: Grid.C:519
simulation::grid::subregion::regionName
std::string regionName
Definition: Grid.H:37
ix::global::GlobalObj::StdOut
virtual int StdOut(const std::string &outstr, unsigned char inlev=1)
Definition: Global.H:439
simulation::grid::halo::CreateRemoteHaloExtents
std::vector< pcpp::IndexIntervalType > CreateRemoteHaloExtents(const pcpp::IndexIntervalType &globalExtent, const pcpp::IndexIntervalType &partitionExtent, std::vector< size_t > &haloSizes)
Definition: Grid.C:2513
CLSIZE3
#define CLSIZE3
Definition: TestGrid.C:780
simulation::grid::halo::SetRemoteHaloExtents
void SetRemoteHaloExtents(const std::vector< pcpp::IndexIntervalType > haloExtents)
Definition: Grid.C:2432
TestGrid_CurvilinearMetric
void TestGrid_CurvilinearMetric(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
Definition: TestGrid.C:949
simulation::grid::halo::SetNeighbors
void SetNeighbors(const std::vector< bool > &inNeighbors)
Definition: Grid.H:145
simulation::grid::halo::messagebuffers::sendBuffers
std::vector< double * > sendBuffers
Definition: Grid.H:66
ix::sys::FILEEXISTS
bool FILEEXISTS(const std::string &fname)
Definition: UnixUtils.C:189
ix::comm::CommunicatorObject
Main encapsulation of MPI.
Definition: COMM.H:62
simulation::report::Grid
std::string Grid(const GridType &inGrid)
Definition: Report.H:14
plascom2::operators::sbp::CreateStencilConnectivity
int CreateStencilConnectivity(int numDim, size_t *dimSizes, size_t *opInterval, int boundaryDepth, int *stencilID, bool fortranInterval=false)
Creates simple stencil connectivity assuming all domain boundaries are physical boundaries.
Definition: Stencil.C:840
Testing.H
Testing constructs for unit testing.
plascom2::operators::stencilset
Encapsulation for a collection of operator stencils.
Definition: Stencil.H:26
pcpp::field::dataset::GetFieldBuffer
BufferDataType * GetFieldBuffer(const std::string &fieldName)
Definition: PCPPFieldData.H:746
simulation::grid::parallel_blockstructured::SetType
void SetType(int inType)
Definition: Grid.H:596
gridJacobian
void const size_t const size_t const size_t const int const double * gridJacobian
Definition: MetricKernels.H:9
TestGrid_PBS_IntegratedHalo
void TestGrid_PBS_IntegratedHalo(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
Definition: TestGrid.C:1800
apptools::DumpContents
void DumpContents(std::ostream &Ostr, const ContainerType &c, std::string del="\)
Dump container contents.
Definition: AppTools.H:117
simulation::grid::parallel_blockstructured::ComputeMetricIdentities
int ComputeMetricIdentities(std::vector< double > &identityData, std::ostream &messageStream)
Definition: Grid.C:737
simulation::grid::parallel_blockstructured::ParallelSetup
int ParallelSetup(fixtures::CommunicatorType &inComm, pcpp::ParallelTopologyInfoType &cartInfo, int haloDepth, std::ostream &messageStream)
Perform cartesian communicator initialization for a grid.
Definition: Grid.C:22
simulation::grid::halo::SetLocalHaloExtents
void SetLocalHaloExtents(const std::vector< pcpp::IndexIntervalType > haloExtents)
Definition: Grid.C:2492
simulation::grid::halo::ReceiveMessage
int ReceiveMessage(const int messageId, const std::vector< int > &neighborRanks, fixtures::CommunicatorType &inComm)
Definition: Grid.C:3096
bufferSizes
void const size_t const size_t * bufferSizes
Definition: MetricKernels.H:19
bufferInterval
void size_t int size_t int size_t int int int int double int int double double *void size_t int size_t int int int int int double int size_t size_t size_t double double *void size_t int size_t int size_t size_t int double int double double *void size_t size_t size_t * bufferInterval
Definition: OperatorKernels.H:29
pcpp::util::PartitionCartesianInterval
int PartitionCartesianInterval(const pcpp::IndexIntervalType &globalExtent, const std::vector< int > &cartDims, const std::vector< int > &cartCoords, pcpp::IndexIntervalType &partExtent, std::ostream &messageStream)
Get local sub-interval of an n-dimensional integer interval.
Definition: PCPPIntervalUtils.C:372
ix::test::results::UpdateResult
void UpdateResult(const std::string &name, const ValueType &result)
Updates an existing test result.
Definition: Testing.H:55
simulation::grid::parallel_blockstructured::CoordinateData
std::vector< double > & CoordinateData()
Definition: Grid.H:503
ix::util::sizeextent::PrettyPrint
std::ostream & PrettyPrint(std::ostream &outStream) const
Definition: IndexUtil.C:6
pcpp::comm::SetupCartesianTopology
int SetupCartesianTopology(pcpp::CommunicatorType &parallelCommunicator, pcpp::ParallelTopologyInfoType &topologyInfo, pcpp::CommunicatorType &topoCommunicator, std::ostream &messageStream)
Sets up a communicator with Cartesian topology.
Definition: PCPPCommUtil.C:48
plascom2::operators::stencilset::numStencils
int numStencils
The number of stencils (e.g. interior + boundary)
Definition: Stencil.H:93
pcpp::io::xdmf::WriteGridSection
int WriteGridSection(const std::string &topoType, const std::string &geomType, const std::string &fileName, const std::string &gridPath, const std::vector< size_t > &gridSize, std::ostream &outStream)
Definition: PCPPHDF5.C:38
ix::global::ParallelGlobalObj::Finalize
virtual int Finalize()
Finalizes the global object, and it&#39;s profiler object.
Definition: Global.H:622
ix::global::GlobalObj::SetVerbLevel
void SetVerbLevel(unsigned char l)
Definition: Global.H:393
testfixtures::viscid::pbsgrid_t
simulation::grid::parallel_blockstructured pbsgrid_t
Definition: ViscidTestFixtures.H:10
simulation::grid::parallel_blockstructured::JacobianMatrix
std::vector< double > & JacobianMatrix()
Definition: Grid.H:626
plascom2::operators::sbp::Initialize
int Initialize(base &stencilSet, int interiorOrder)
Initialize the sbp::base stencilset with the SBP operator of given order.
Definition: Stencil.C:360
simulation::grid::halo::messagebuffers::numPointsRecv
std::vector< size_t > numPointsRecv
Definition: Grid.H:65
plascom2::halo_t
simulation::grid::halo halo_t
Definition: PlasCom2.H:18
ix::util::sizeextent
Simple Block Structured Mesh object.
Definition: IndexUtil.H:21
pcpp::field::dataset::AddField
void AddField(const std::string &name, char loc, unsigned int ncomp, unsigned int dsize, const std::string &unit)
Definition: PCPPFieldData.H:869
simulation::grid::parallel_base::Communicator
fixtures::CommunicatorType & Communicator() const
Definition: Grid.H:350
numComponents
void size_t int * numComponents
Definition: OperatorKernels.H:9
TestGrid_RectilinearMetric
void TestGrid_RectilinearMetric(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
Definition: TestGrid.C:435
pcpp::field::dataset::SetFieldBuffer
void SetFieldBuffer(const std::string &name, void *buf)
Definition: PCPPFieldData.H:716
simulation::grid::parallel_blockstructured::SetPhysicalExtent
void SetPhysicalExtent(const std::vector< double > &inExtent)
Definition: Grid.H:535
simulation::grid::parallel_blockstructured::Finalize
int Finalize(bool allocateCoordinateData=false)
Definition: Grid.C:2318
ix::util::sizeextent::InitSimple
void InitSimple(const ContainerType &inSize)
Definition: IndexUtil.H:169
simulation::grid::parallel_blockstructured::BufferSize
size_t BufferSize() const
Definition: Grid.H:430
numPointsBuffer
void const size_t * numPointsBuffer
Definition: MetricKernels.H:19
simulation::grid::halo::SetLocalBufferSizes
void SetLocalBufferSizes(const std::vector< size_t > &inBufferSizes)
Definition: Grid.C:2423
simulation::grid::parallel_blockstructured
Definition: Grid.H:387