PlasCom2  1.0
XPACC Multi-physics simluation application
Grid.C
Go to the documentation of this file.
1 #include "Grid.H"
2 #include "PCPPIntervalUtils.H"
3 #include "PCPPReport.H"
4 #include "OperatorKernels.H"
5 
6 namespace simulation {
7 
8  namespace grid {
9 
22  int parallel_blockstructured::ParallelSetup(fixtures::CommunicatorType &inComm,
23  pcpp::ParallelTopologyInfoType &cartInfo,
24  int haloDepth,std::ostream &messageStream)
25  {
26 
27  if(partitionIntervals.gridSizes.empty()){
28  messageStream << "grid::PBS::ParallelSetup: Error gridSizes empty." << std::endl;
29  return(1);
30  }
31 
32  std::vector<size_t> &gridSizes(GridSizes());
33  std::vector<bool> &isPeriodic(PeriodicDirs());
34  std::vector<size_t> &bufferSizes(BufferSizes());
35  fixtures::CommunicatorType &gridComm(Communicator());
36 
37  numDim = gridSizes.size();
38 
39  std::vector<int> &periodicDirs(cartInfo.isPeriodic);
40 
41  // Use periodic dirs indication from cartInfo only
42  if(periodicDirs.size() != numDim){
43  periodicDirs.resize(numDim,0);
44  }
45  for(int iDim = 0;iDim < numDim;iDim++){
46  isPeriodic.resize(numDim,false);
47  if(periodicDirs[iDim] > 0)
48  isPeriodic[iDim] = true;
49  }
50 
51  std::vector<int> haloDepths(2*numDim,haloDepth);
52 
53  pcpp::comm::SetupCartesianTopology(inComm,cartInfo,gridComm,messageStream);
54 
55  // Grab the local coordinates and global dimension of the Cartesian topo
56  std::vector<int> &cartCoords(gridComm.CartCoordinates());
57  std::vector<int> &cartDims(gridComm.CartDimensions());
58 
59  // Generate a report of the Cartesian setup and put it in the messageStream
60  messageStream << "grid::PBS::ParallelSetup: Cartesian setup:" << std::endl;
61  pcpp::report::CartesianSetup(messageStream,cartInfo);
62 
63  pcpp::IndexIntervalType &partitionInterval(PartitionInterval());
64  pcpp::IndexIntervalType globalInterval;
65  globalInterval.InitSimple(gridSizes);
66 
67  // === Partition the grid ===
68  if(pcpp::util::PartitionCartesianInterval(globalInterval,cartDims,cartCoords,
69  partitionInterval,messageStream)){
70  messageStream << "grid::PBS::ParallelSetup:Error: Partitioning failed." << std::endl;
71  return(1);
72  }
73 
74 
75  std::vector<size_t> extendGrid(2*numDim,0);
76 
77  std::vector<bool> haveNeighbors(2*numDim,true);
78  for(int iDim = 0;iDim < numDim; iDim++){
79  if(isPeriodic[iDim]){
80  extendGrid[2*iDim] = haloDepths[2*iDim];
81  extendGrid[2*iDim+1] = haloDepths[2*iDim+1];
82  } else {
83  if(partitionInterval[iDim].first == 0) {
84  haveNeighbors[2*iDim] = false;
85  } else {
86  extendGrid[2*iDim] = haloDepths[2*iDim];
87  }
88  if(partitionInterval[iDim].second == (gridSizes[iDim]-1)){
89  haveNeighbors[2*iDim + 1] = false;
90  } else {
91  extendGrid[2*iDim + 1] = haloDepths[2*iDim+1];
92  }
93  }
94  }
95 
96  messageStream << "grid::PBS::ParallelSetup: Configuring grid extensions (";
97  pcpp::io::DumpContents(messageStream,extendGrid,",");
98  messageStream << ")" << std::endl;
99 
100  SetDimensionExtensions(extendGrid);
101 
102  if(Finalize()){
103  messageStream << "grid::PBS::ParallelSetup:Error: finalization failed." << std::endl;
104  return(1);
105  }
106 
107  const pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
108 
109  // Set up the HALO data structure for this grid/geometry
110  simulation::grid::halo &myHalo(Halo());
111 
112  myHalo.SetGridInterval(globalInterval);
113  myHalo.SetPartitionInterval(partitionInterval);
114 
115  // Set communication directions for halo
116  myHalo.SetNeighbors(haveNeighbors);
117  // myHalo.SetPeriodicDirs(cartInfo.isPeriodic);
118 
119  // Simple create halo extents just creates uniformly sized halos everywhere
120  std::vector<pcpp::IndexIntervalType> remoteHaloExtents(myHalo.CreateRemoteHaloExtents(extendGrid));
121  std::vector<pcpp::IndexIntervalType> localHaloExtents(myHalo.CreateLocalHaloExtents(extendGrid));
122 
123  myHalo.SetLocalBufferSizes(bufferSizes);
124  myHalo.SetLocalPartitionExtent(partitionBufferInterval);
125 
126  myHalo.SetRemoteHaloExtents(remoteHaloExtents);
127  myHalo.SetLocalHaloExtents(localHaloExtents);
128  myHalo.CreateSimpleSendIndices();
129  myHalo.CreateSimpleRecvIndices();
130 
131  return(0);
132 
133  };
134 
135  int ResolveTopoName(const std::string &inName)
136  {
137  for(int iType = (int)CARTESIAN;iType < (int)NUMTOPOTYPE;iType++){
138  std::string topoName(topoNames[iType]);
139  if(inName == topoName) return(iType);
140  }
141  return((int)NUMTOPOTYPE);
142  };
143 
144  std::string ResolveTopoType(int inType)
145  {
146  if(inType >= 0 && inType < (int)NUMTOPOTYPE)
147  return(std::string(topoNames[inType]));
148  return(std::string(topoNames[(int)NUMTOPOTYPE]));
149  }
150 
151  // Take a partitionInterval, partitionBufferInterval, and partitionSubInterval and
152  // return a partitionSubBufferInterval so-to-speak i.e. subinterval wrt buffer
153  pcpp::IndexIntervalType ResolveBufferInterval(const pcpp::IndexIntervalType &partitionInterval,
154  const pcpp::IndexIntervalType &partitionBufferInterval,
155  const pcpp::IndexIntervalType &partitionSubInterval)
156  {
157  pcpp::IndexIntervalType returnValue;
158  pcpp::IndexIntervalType collisionInterval(partitionInterval.Overlap(partitionSubInterval));
159  if(collisionInterval.NNodes() == 0)
160  return(returnValue);
161  int numDim = partitionSubInterval.size();
162  returnValue = collisionInterval; // just init to correct size
163  for(int iDim = 0;iDim < numDim;iDim++){
164  collisionInterval[iDim].first -= partitionInterval[iDim].first;
165  collisionInterval[iDim].second -= partitionInterval[iDim].first;
166  returnValue[iDim].first = collisionInterval[iDim].first + partitionBufferInterval[iDim].first;
167  returnValue[iDim].second = collisionInterval[iDim].second + partitionBufferInterval[iDim].first;
168  }
169  returnValue.Sync();
170  return(returnValue);
171  };
172 
173  // Used for reading subregions from config
174  int InitSubRegionFromString(const std::string &inString,const std::vector<size_t> &gridSizes,
175  subregion &inSubRegion)
176  {
177 
178  if(inString.empty())
179  return(1);
180 
181  int numDim = gridSizes.size();
182 
183  pcpp::IndexIntervalType gridInterval;
184  gridInterval.InitSimple(gridSizes);
185 
186  std::istringstream Istr(inString);
187 
188  std::string regionName;
189  int normalDirection = 0;
190  Istr >> regionName >> normalDirection;
191  if(regionName.empty())
192  return(1);
193  inSubRegion.regionName = regionName;
194  if(std::abs(normalDirection) > numDim)
195  return(1);
196 
197  inSubRegion.normalDirection = normalDirection;
198  pcpp::IndexIntervalType &regionInterval(inSubRegion.regionInterval);
199  regionInterval.InitSimple(gridSizes); // copy grid interval and pare it down below
200 
201  for(int iDim = 0;iDim < numDim;iDim++){
202 
203  int leftBoundary = (iDim+1);
204  int rightBoundary = -leftBoundary;
205 
206  if(normalDirection == leftBoundary){
207  regionInterval[iDim].first = 0;
208  regionInterval[iDim].second = 0;
209  } else if(normalDirection == rightBoundary){
210  regionInterval[iDim].first = gridInterval[iDim].second;
211  regionInterval[iDim].second = gridInterval[iDim].second;
212  }
213 
214  // If specified, use specific indices
215  if(Istr){
216  int iStart = 0;
217  int iEnd = 0;
218  Istr >> iStart >> iEnd;
219  if((iStart != 0) && (iEnd != 0)){
220  if(iStart < 0){
221  iStart += 1;
222  regionInterval[iDim].first = gridInterval[iDim].second + iStart;
223  } else {
224  iStart -= 1;
225  regionInterval[iDim].first = iStart;
226  }
227  if(iEnd < 0){
228  iEnd += 1;
229  regionInterval[iDim].second = gridInterval[iDim].second + iEnd;
230  } else {
231  iEnd -= 1;
232  regionInterval[iDim].second = iEnd;
233  }
234 
235  // Validate
236  if(regionInterval[iDim].first > gridInterval[iDim].second)
237  return(1);
238  if(regionInterval[iDim].first > regionInterval[iDim].second)
239  return(1);
240  if(regionInterval[iDim].second > gridInterval[iDim].second)
241  return(1);
242  }
243  }
244  }
245 
246  // Sync up the created interval (always req'd when setting up by hand)
247  regionInterval.Sync();
248 
249  return(0);
250  };
251 
252  int parallel_blockstructured::ComputeUniformRectangularMetrics(std::ostream &messageStream)
253  {
254 
255  std::vector<size_t> &gridSizes(GridSizes());
256  std::vector<size_t> &bufferSizes(BufferSizes());
257  size_t numPointsBuffer = BufferSize();
258 
259  pcpp::IndexIntervalType &partInterval(PartitionInterval());
260  pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
261 
262  int numDim = Dimension(); // gridSizes.size();
263 
264  messageStream << "grid::PBS::ComputeMetrics: Computing uniform grid metrics."
265  << std::endl;
266  if(dX.size() != numDim){
267  messageStream << "grid::PBS::ComputeMetrics: Warning: grid spacings of invalid size ("
268  << dX.size() << " != " << numDim << "), attempting to generate."
269  << std::endl;
270 
271  pcpp::IndexIntervalType bufferInterval;
272  bufferInterval.InitSimple(bufferSizes);
273  const std::vector<double> &gridCoordinates(CoordinateData());
274  dX.resize(numDim);
275  if(numDim == 1){
276  int index1 = partitionBufferInterval[0].first;
277  int index2 = index1+1;
278  dX[0] = gridCoordinates[index2] - gridCoordinates[index1];
279  } else if (numDim == 2) {
280  pcpp::IndexIntervalType littleBox(partitionBufferInterval);
281  littleBox[0].second = littleBox[0].first+1;
282  littleBox[1].second = littleBox[1].first+1;
283  std::vector<size_t> littleBoxIndices;
284  bufferInterval.GetFlatIndices(littleBox,littleBoxIndices);
285  dX[0] = gridCoordinates[littleBoxIndices[1]] - gridCoordinates[littleBoxIndices[0]];
286  dX[1] = gridCoordinates[littleBoxIndices[2]+numPointsBuffer] -
287  gridCoordinates[littleBoxIndices[0]+numPointsBuffer];
288  } else {
289  pcpp::IndexIntervalType littleBox(partitionBufferInterval);
290  littleBox[0].second = littleBox[0].first+1;
291  littleBox[1].second = littleBox[1].first+1;
292  littleBox[2].second = littleBox[2].first+1;
293  std::vector<size_t> littleBoxIndices;
294  bufferInterval.GetFlatIndices(littleBox,littleBoxIndices);
295  dX[0] = gridCoordinates[littleBoxIndices[1]] - gridCoordinates[littleBoxIndices[0]];
296  dX[1] = gridCoordinates[littleBoxIndices[2]+numPointsBuffer] -
297  gridCoordinates[littleBoxIndices[0]+numPointsBuffer];
298  dX[2] = gridCoordinates[littleBoxIndices[4]+2*numPointsBuffer] -
299  gridCoordinates[littleBoxIndices[0]+2*numPointsBuffer];
300  }
301  messageStream << "Generated dX(";
302  pcpp::io::DumpContents(messageStream,dX,",");
303  messageStream << ")" << std::endl;
304  }
305 
306  gridMetric.resize(numDim,0.0);
307  gridJacobian.resize(2,1.0);
308 
309  for(int iDim = 0;iDim < numDim;iDim++){
310  if(dX[iDim] <= 0.0){
311  messageStream << "grid::PBS::ComputeMetrics: Error: invalid grid spacings (";
312  pcpp::io::DumpContents(messageStream,dX,",");
313  messageStream << "). Aborting." << std::endl;
314  return(1);
315  }
316  gridMetric[iDim] = 1/dX[iDim];
317  gridJacobian[0] *= gridMetric[iDim];
318  }
319  for(int iDim = 0;iDim < numDim;iDim++)
320  gridMetric[iDim] /= gridJacobian[0];
321 
322  gridJacobian[1] = 1.0/gridJacobian[0];
323 
324  return(0);
325  };
326 
327 
328  int parallel_blockstructured::ComputeRectilinearMetrics(std::ostream &messageStream)
329  {
330 
331  messageStream << "grid::PBS::ComputeMetrics: Computing metric for rectilinear grid."
332  << std::endl;
333  if(!stencilSet || !stencilConnectivity){
334  messageStream << "grid::PBS::ComputeMetrics: Error: Differential operators must be set "
335  << "before computing rectilinear metrics." << std::endl;
336  return(1);
337  }
338 
339  std::vector<size_t> &gridSizes(GridSizes());
340  std::vector<size_t> &bufferSizes(BufferSizes());
341 
342  pcpp::IndexIntervalType &partInterval(PartitionInterval());
343  pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
344 
345  int numDim = gridSizes.size();
346 
347  int numStencils = stencilSet->numStencils;
348  int *stencilSizes = stencilSet->stencilSizes;
349  int *stencilStarts = stencilSet->stencilStarts;
350  int *stencilOffsets = stencilSet->stencilOffsets;
351  double *stencilWeights = stencilSet->stencilWeights;
352  int numStencilValues = stencilSet->numValues;
353 
354  // Threading this routine would involve having the
355  // thread id right here. It is not certain how
356  // best to get do that. Perhaps two routines
357  // are needed, one threaded, one not. Hokey.
358  std::vector<size_t> dOpInterval;
359  partitionBufferInterval.Flatten(dOpInterval);
360  // forticate the opInterval
361  std::vector<size_t>::iterator dOpIt = dOpInterval.begin();
362  while(dOpIt != dOpInterval.end()){
363  *dOpIt += 1;
364  dOpIt++;
365  }
366 
367 
368  int numMetricComponents = numDim;
369  size_t numPointsBuffer = BufferSize();
370 
371  gridMetric.resize(numMetricComponents*numPointsBuffer,0);
372  gridJacobian.resize(2*numPointsBuffer,1.0); // J and 1/J
373 
374  // for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
375  // gridMetric[iPoint] = dX[1];
376  // gridJacobian[iPoint] = 1.0/(dX[0]*dX[1]);
377  // }
378  // for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
379  // gridMetric[numPointsBuffer+iPoint] = dX[0];
380  // gridJacobian[iPoint+numPointsBuffer] = 1.0/gridJacobian[iPoint];
381  // }
382  // return(0);
383 
384  int one = 1;
385 
386  if(numDim == 2){
387 
388  double *xData = &gridCoordinates[0];
389  double *yData = &gridCoordinates[numPointsBuffer];
390 
391 
392  // metric[1][1] = y_eta
393  int dDir = 2;
394  int metricComponent = 1;
395  size_t dirOffset = (dDir-1)*numPointsBuffer;
396  size_t metricOffset = (metricComponent-1)*numPointsBuffer;
397  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
398  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
399  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
400  stencilStarts,&numStencilValues,stencilWeights,
401  stencilOffsets,&stencilConnectivity[dirOffset],
402  yData,&gridMetric[0]);
403 
404  // metric[2][2] = x_xi
405  dDir = 1;
406  metricComponent = 2;
407  dirOffset = (dDir - 1)*numPointsBuffer;
408  metricOffset = (metricComponent - 1)*numPointsBuffer;
409  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
410  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
411  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
412  stencilStarts,&numStencilValues,stencilWeights,
413  stencilOffsets,&stencilConnectivity[dirOffset],
414  xData,&gridMetric[metricOffset]);
415 
416 
417  // 1/J = metric[1][1] * metric[2][2]
418  FC_MODULE(operators,zxy,OPERATORS,ZXY)
419  (&numDim,&numPointsBuffer,&bufferSizes[0],
420  &dOpInterval[0],&gridMetric[0],
421  &gridMetric[metricOffset],&gridJacobian[numPointsBuffer]);
422 
423 
424  // Finish off with J = 1 / (1/J)
425  double dOne = 1.0;
426  FC_MODULE(operators,yaxm1,OPERATORS,YAXM1)
427  (&numDim,&numPointsBuffer,&bufferSizes[0],
428  &dOpInterval[0],&dOne,&gridJacobian[numPointsBuffer],
429  &gridJacobian[0]);
430 
431 
432  } else if (numDim == 3) {
433 
434  double *xData = &gridCoordinates[0];
435  double *yData = &gridCoordinates[numPointsBuffer];
436  double *zData = &gridCoordinates[2*numPointsBuffer];
437 
438  // y_eta (going into J^-1 (temporarily))
439  int dDir = 2;
440  size_t dirOffset = (dDir-1)*numPointsBuffer;
441  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
442  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
443  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
444  stencilStarts,&numStencilValues,stencilWeights,
445  stencilOffsets,&stencilConnectivity[dirOffset],
446  yData,&gridJacobian[numPointsBuffer]);
447 
448  // x_xi (going into J (temporarily))
449  dDir = 1;
450  dirOffset = (dDir-1)*numPointsBuffer;
451  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
452  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
453  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
454  stencilStarts,&numStencilValues,stencilWeights,
455  stencilOffsets,&stencilConnectivity[dirOffset],
456  xData,&gridJacobian[0]);
457 
458 
459  // z_zeta (going into metric[1][1])
460  dDir = 3;
461  dirOffset = (dDir-1)*numPointsBuffer;
462  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
463  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
464  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
465  stencilStarts,&numStencilValues,stencilWeights,
466  stencilOffsets,&stencilConnectivity[dirOffset],
467  zData,&gridMetric[0]);
468 
469  // The following quantities have been calculated:
470  // * x_xi stored in J
471  // * y_eta stored in J^-1
472  // * z_zeta stored in metric[1][1]
473 
474  // Finish off metric terms using temporary storage from above
475 
476  // metric[3][3] = x_xi * y_eta = J * J^-1
477  size_t metricOffset = 2*numPointsBuffer;
478  FC_MODULE(operators,zxy,OPERATORS,ZXY)
479  (&numDim,&numPointsBuffer,&bufferSizes[0],
480  &dOpInterval[0],&gridJacobian[0],
481  &gridJacobian[numPointsBuffer],&gridMetric[metricOffset]);
482 
483  // metric[2][2] = x_xi * z_zeta = J * metric[1][1]
484  metricOffset = numPointsBuffer;
485  FC_MODULE(operators,zxy,OPERATORS,ZXY)
486  (&numDim,&numPointsBuffer,&bufferSizes[0],
487  &dOpInterval[0],&gridJacobian[0],
488  &gridMetric[0],&gridMetric[metricOffset]);
489 
490  // metric[1][1] = z_zeta * y_eta = metric[1][1] * J^-1
491  metricOffset = 0;
492  FC_MODULE(operators,yxy,OPERATORS,YXY)
493  (&numDim,&numPointsBuffer,&bufferSizes[0],
494  &dOpInterval[0],&gridJacobian[numPointsBuffer],
495  &gridMetric[metricOffset]);
496 
497  // J^-1 = x_xi * y_eta * z_zeta = J^-1 * metric[2][2]
498  metricOffset = numPointsBuffer;
499  FC_MODULE(operators,yxy,OPERATORS,YXY)
500  (&numDim,&numPointsBuffer,&bufferSizes[0],
501  &dOpInterval[0],&gridMetric[metricOffset],
502  &gridJacobian[numPointsBuffer]);
503 
504  // J = 1/J^-1
505  double dOne = 1.0;
506  FC_MODULE(operators,yaxm1,OPERATORS,YAXM1)
507  (&numDim,&numPointsBuffer,&bufferSizes[0],
508  &dOpInterval[0],&dOne,&gridJacobian[numPointsBuffer],
509  &gridJacobian[0]);
510 
511  } else if (numDim == 1){
512  return(1);
513  }
514 
515  return(0);
516  };
517 
518  // Metrics calculations are currently thread-safe but *not* themselves threaded
519  int parallel_blockstructured::ComputeMetrics(std::ostream &messageStream)
520  {
521 
522  std::vector<size_t> &gridSizes(GridSizes());
523  std::vector<size_t> &bufferSizes(BufferSizes());
524 
525  pcpp::IndexIntervalType &partInterval(PartitionInterval());
526  pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
527 
528  if(gridSizes.empty() || bufferSizes.empty() ||
529  (partInterval.NNodes() == 0) || (partitionBufferInterval.NNodes() == 0)){
530  messageStream << "grid::PBS::ComputeMetrics: Error: Grid sizes not set up."
531  << std::endl
532  << "grid::PBS::ComputeMetrics: Finalize() must be called before ComputeMetrics."
533  << std::endl;
534  return(1);
535  }
536 
537 
538  if(gridType == simulation::grid::CARTESIAN ||
539  gridType == simulation::grid::UNIRECT){
540  return(ComputeUniformRectangularMetrics(messageStream));
541  } else {
542  if(!stencilSet || !stencilConnectivity){
543  messageStream << "grid::PBS::ComputeMetrics: Error: Differential operators must be set "
544  << "before computing metrics." << std::endl;
545  return(1);
546  }
547  if(gridType == simulation::grid::RECTILINEAR)
548  return(ComputeRectilinearMetrics(messageStream));
549  }
550 
551  int numDim = gridSizes.size();
552  if(numDim == 2) {
553  return(ComputeCurvilinearMetrics2D(messageStream));
554  } else if (numDim == 3) {
555  // return(ComputeCurvilinearMetrics3D(messageStream));
556  return(CurvilinearMetrics(messageStream));
557  } else if (numDim > 1) {
558  messageStream << "grid::ComputeMetrics:Error: "
559  << "Unsupported number of dimensions."
560  << std::endl;
561  return(1);
562  }
563 
564  // Curvilinear in 1D is equivalent to rectilinear in 1D
565  return(ComputeRectilinearMetrics(messageStream));
566 
567  };
568 
569  int parallel_blockstructured::ComputeCurvilinearMetrics2D(std::ostream &messageStream)
570  {
571 
572  std::vector<size_t> &gridSizes(GridSizes());
573  std::vector<size_t> &bufferSizes(BufferSizes());
574 
575  pcpp::IndexIntervalType &partInterval(PartitionInterval());
576  pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
577 
578  if(gridSizes.empty() || bufferSizes.empty() ||
579  (partInterval.NNodes() == 0) || (partitionBufferInterval.NNodes() == 0)){
580  messageStream << "grid::PBS::ComputeCurvilinearMetrics: Error: Grid sizes not set up."
581  << std::endl
582  << "grid::PBS::ComputeCurvilinearMetrics: Finalize() must be called before ComputeMetrics."
583  << std::endl;
584  return(1);
585  }
586 
587  if(numDim != 2){
588  messageStream << "grid::PBS::ComputeMetrics:Error: Cannot create 2D curvilinear metrics for "
589  << numDim<< "-dimensional grid." << std::endl;
590  return(1);
591  }
592 
593  messageStream << "grid::PBS::ComputeMetrics: Computing 2d curvilinear metric."
594  << std::endl;
595 
596 
597 
598  int numDim = partitionIntervals.numDim;
599  std::vector<double> &gridCoordinates(CoordinateData());
600 
601 
602  int numStencils = stencilSet->numStencils;
603  int *stencilSizes = stencilSet->stencilSizes;
604  int *stencilStarts = stencilSet->stencilStarts;
605  int *stencilOffsets = stencilSet->stencilOffsets;
606  double *stencilWeights = stencilSet->stencilWeights;
607  int numStencilValues = stencilSet->numValues;
608 
609  // Threading this routine would involve having the
610  // thread id right here. It is not certain how
611  // best to get do that. Perhaps two routines
612  // are needed, one threaded, one not. Hokey.
613  std::vector<size_t> dOpInterval;
614  partitionBufferInterval.Flatten(dOpInterval);
615  // forticate the opInterval
616  std::vector<size_t>::iterator dOpIt = dOpInterval.begin();
617  while(dOpIt != dOpInterval.end()){
618  *dOpIt += 1;
619  dOpIt++;
620  }
621 
622 
623  int numMetricComponents = numDim*numDim;
624  size_t numPointsBuffer = BufferSize();
625 
626  gridMetric.resize(numMetricComponents*numPointsBuffer,0);
627  gridJacobian.resize(2*numPointsBuffer,1.0);
628 
629  int one = 1;
630 
631  double *xData = &gridCoordinates[0];
632  double *yData = &gridCoordinates[numPointsBuffer];
633 
634 
635  // metric[1][1] = y_eta
636  int dDir = 2;
637  int metricComponent = 1;
638  size_t dirOffset = (dDir-1)*numPointsBuffer;
639  size_t metricOffset = (metricComponent-1)*numPointsBuffer;
640  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
641  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
642  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
643  stencilStarts,&numStencilValues,stencilWeights,
644  stencilOffsets,&stencilConnectivity[dirOffset],
645  yData,&gridMetric[0]);
646 
647  // metric[2][2] = x_xi
648  dDir = 1;
649  metricComponent = 4;
650  dirOffset = (dDir - 1)*numPointsBuffer;
651  metricOffset = (metricComponent - 1)*numPointsBuffer;
652 
653  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
654  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
655  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
656  stencilStarts,&numStencilValues,stencilWeights,
657  stencilOffsets,&stencilConnectivity[dirOffset],
658  xData,&gridMetric[metricOffset]);
659 
660 
661  // 1/J = metric[1][1] * metric[2][2]
662  FC_MODULE(operators,zxy,OPERATORS,ZXY)
663  (&numDim,&numPointsBuffer,&bufferSizes[0],
664  &dOpInterval[0],&gridMetric[0],
665  &gridMetric[metricOffset],&gridJacobian[numPointsBuffer]);
666 
667  // if(gridType == simulation::grid::CURVILINEAR){ // calculate cross-terms
668  // metric[1][1] = gridMetric[0]
669  // metric[1][2] = gridMetric[numPointsBuffer]
670  // metric[2][1] = gridMetric[2*numPointsBuffer]
671  // metric[2][2] = gridMetric[3*numPointsBuffer]
672 
673  // deta_dy = metric[1][2] = [-]x_eta
674  double negOne = -1.0;
675  dDir = 2;
676  metricComponent = 2;
677  dirOffset = (dDir - 1)*numPointsBuffer;
678  metricOffset = (metricComponent-1)*numPointsBuffer;
679  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
680  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
681  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
682  stencilStarts,&numStencilValues,stencilWeights,
683  stencilOffsets,&stencilConnectivity[dirOffset],
684  xData,&gridMetric[metricOffset]);
685 
686  // Make the sign right
687  // metric[1][2] = -metric[1][2]
688  FC_MODULE(operators,xax,OPERATORS,XAX)
689  (&numDim,&numPointsBuffer,&bufferSizes[0],
690  &dOpInterval[0],&negOne,&gridMetric[metricOffset]);
691 
692 
693  // deta_dx = metric[2][1] = [-]y_xi
694  dDir = 1;
695  metricComponent = 3;
696  dirOffset = (dDir - 1)*numPointsBuffer;
697  metricOffset = (metricComponent-1)*numPointsBuffer;
698  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
699  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
700  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
701  stencilStarts,&numStencilValues,stencilWeights,
702  stencilOffsets,&stencilConnectivity[dirOffset],
703  yData,&gridMetric[metricOffset]);
704 
705  // Make the sign right
706  // metric[2][1] = -metric[2][1]
707  FC_MODULE(operators,xax,OPERATORS,XAX)
708  (&numDim,&numPointsBuffer,&bufferSizes[0],
709  &dOpInterval[0],&negOne,&gridMetric[metricOffset]);
710 
711  // J(temporarily) = metric[1][2]*metric[2][1]
712  FC_MODULE(operators,zxy,OPERATORS,ZXY)
713  (&numDim,&numPointsBuffer,&bufferSizes[0],
714  &dOpInterval[0],&gridMetric[numPointsBuffer],
715  &gridMetric[metricOffset],&gridJacobian[0]);
716 
717  // 1/J = y_eta*x_xi - y_xi*x_eta = 1/J - J(temporary)
718  // = metric[1][1]*metric[2][2] - metric[1][2]*metric[2][1]
719  FC_MODULE(operators,yaxpy,OPERATORS,YAXPY)
720  (&numDim,&numPointsBuffer,&bufferSizes[0],
721  &dOpInterval[0],&negOne,&gridJacobian[0],
722  &gridJacobian[numPointsBuffer]);
723 
724 
725  // Finish off with J = 1 / (1/J)
726  double dOne = 1.0;
727  FC_MODULE(operators,yaxm1,OPERATORS,YAXM1)
728  (&numDim,&numPointsBuffer,&bufferSizes[0],
729  &dOpInterval[0],&dOne,&gridJacobian[numPointsBuffer],
730  &gridJacobian[0]);
731 
732  return(0);
733 
734  };
735 
736 
737  int parallel_blockstructured::ComputeMetricIdentities(std::vector<double> &identityData,std::ostream &messageStream)
738  {
739  if(gridMetric.empty()){
740  messageStream << "PBS::ComputeMetricIdentities: Error - grid metrics empty." << std::endl;
741  return(1);
742  }
743  if(gridType < simulation::grid::CURVILINEAR){
744  identityData.resize(numDim,0.0);
745  messageStream << "PBS::ComputeMetricIdentities: Metric identities for gridType("
746  << gridType << ") are trivially zero." << std::endl;
747  return(0);
748  }
749 
750  std::vector<size_t> &gridSizes(GridSizes());
751  std::vector<size_t> &bufferSizes(BufferSizes());
752  pcpp::IndexIntervalType &partInterval(PartitionInterval());
753  pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
754 
755  int numDim = Dimension();
756  size_t numPointsBuffer = BufferSize();
757  if(numPointsBuffer == 0){
758  messageStream << "PBS::ComputeMetricIdentities: Error - no points in the grid." << std::endl;
759  return(1);
760  }
761  if(numDim < 1){
762  messageStream << "PBS::ComputeMetricIdentities: Error - grid of dimension " << numDim
763  << " not supported." << std::endl;
764  return(1);
765  }
766 
767  int numStencils = stencilSet->numStencils;
768  int *stencilSizes = stencilSet->stencilSizes;
769  int *stencilStarts = stencilSet->stencilStarts;
770  int *stencilOffsets = stencilSet->stencilOffsets;
771  double *stencilWeights = stencilSet->stencilWeights;
772  int numStencilValues = stencilSet->numValues;
773 
774  std::vector<size_t> dOpInterval;
775  partitionBufferInterval.Flatten(dOpInterval);
776 
777  // forticate the opInterval
778  std::vector<size_t>::iterator dOpIt = dOpInterval.begin();
779  while(dOpIt != dOpInterval.end()){
780  *dOpIt += 1;
781  dOpIt++;
782  }
783 
784 
785  int one = 1;
786  double a = 1.0;
787  identityData.resize(numDim*numPointsBuffer,0.0);
788  std::vector<double> identityComponent(numPointsBuffer);
789 
790  int numMetricComponents = numDim*numDim;
791  ExchangeNodalData(numMetricComponents,&gridMetric[0]);
792 
793  for(int iDim = 0;iDim < numDim;iDim++){
794  size_t iIndex = (iDim*numDim*numPointsBuffer);
795  int diffDir = iDim+1;
796  size_t dirOffset = iDim*numPointsBuffer;
797  for(int jDim = 0;jDim < numDim;jDim++){
798  size_t metricIndex = jDim*numPointsBuffer + iIndex;
799  size_t identityOffset = jDim*numPointsBuffer;
800  double *metricComponent = &gridMetric[metricIndex];
801  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
802  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
803  &diffDir,&dOpInterval[0],&numStencils,stencilSizes,
804  stencilStarts,&numStencilValues,stencilWeights,
805  stencilOffsets,&stencilConnectivity[dirOffset],
806  metricComponent,&identityComponent[0]);
807  FC_MODULE(operators,yaxpy,OPERATORS,YAXPY)
808  (&numDim,&numPointsBuffer,&bufferSizes[0],
809  &dOpInterval[0],&a,&identityComponent[0],
810  &identityData[identityOffset]);
811 
812  }
813  }
814 
815  return(0);
816  };
817 
818  int parallel_blockstructured::ExchangeNodalData(int numComponents,double *dataBuffer)
819  {
820 
821  simulation::grid::halo &myHalo(Halo());
822  fixtures::CommunicatorType &gridComm(Communicator());
823 
824  // int saveMessageId = xyzMessageId;
825  xyzMessageId = myHalo.CreateMessageBuffers(numComponents);
826 
827 
828  std::vector<int> cartNeighbors;
829  gridComm.CartNeighbors(cartNeighbors);
830 
831  myHalo.PackMessageBuffers(xyzMessageId,0,numComponents,dataBuffer);
832  myHalo.SendMessage(xyzMessageId,cartNeighbors,gridComm);
833  myHalo.ReceiveMessage(xyzMessageId,cartNeighbors,gridComm);
834  myHalo.UnPackMessageBuffers(xyzMessageId,0,numComponents,dataBuffer);
835 
836  myHalo.DestroyMessageBuffer(xyzMessageId);
837  // xyzMessageId = saveMessageId;
838 
839  return(0);
840 
841  };
842 
843 
844  int parallel_blockstructured::GradIJK(int numComponents,double *sourceBuffer,
845  double *targetBuffer,std::ostream &messageStream,bool exchange = false)
846  {
847 
848  if(!stencilSet || !stencilConnectivity){
849  messageStream << "grid::PBS::GradIJK: Error: Differential operators must be set "
850  << "before using GradIJK." << std::endl;
851  return(1);
852  }
853 
854  if(exchange){
855  if(ExchangeNodalData(numComponents,sourceBuffer)){
856  messageStream << "grid::PBS::GradIJK: Error: Data exchange failed." << std::endl;
857  return(1);
858  }
859  }
860 
861  std::vector<size_t> &gridSizes(GridSizes());
862  std::vector<size_t> &bufferSizes(BufferSizes());
863  pcpp::IndexIntervalType &partInterval(PartitionInterval());
864  pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
865 
866  int numDim = partitionIntervals.numDim;
867  int gradComponents = numComponents * numDim;
868  size_t numPointsBuffer = BufferSize();
869 
870  int numStencils = stencilSet->numStencils;
871  int *stencilSizes = stencilSet->stencilSizes;
872  int *stencilStarts = stencilSet->stencilStarts;
873  int *stencilOffsets = stencilSet->stencilOffsets;
874  double *stencilWeights = stencilSet->stencilWeights;
875  int numStencilValues = stencilSet->numValues;
876 
877  // Threading this routine would involve having the
878  // thread id right here. It is not certain how
879  // best to get do that. Perhaps two routines
880  // are needed, one threaded, one not. Hokey.
881  std::vector<size_t> dOpInterval;
882  partitionBufferInterval.Flatten(dOpInterval);
883  // forticate the opInterval
884  std::vector<size_t>::iterator dOpIt = dOpInterval.begin();
885  while(dOpIt != dOpInterval.end()){
886  *dOpIt += 1;
887  dOpIt++;
888  }
889  int one = 1.0;
890  size_t targetOffset = 0;
891  int targetComp = 0;
892  for(int iComponent = 0;iComponent < numComponents;iComponent++){
893  size_t sourceOffset = iComponent*numPointsBuffer;
894  // Put sourceBuffer(iComp)_ijk in targetBuffer
895  for(int iDir = 0;iDir < numDim;iDir++){
896  size_t dirOffset = iDir*numPointsBuffer;
897  int dDir = iDir+1;
898  targetOffset = numPointsBuffer*targetComp++;
899  // Apply operator returns sourceBuffer(iComp)_iDir
900  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
901  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
902  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
903  stencilStarts,&numStencilValues,stencilWeights,
904  stencilOffsets,&stencilConnectivity[dirOffset],
905  &sourceBuffer[sourceOffset],&targetBuffer[targetOffset]);
906  }
907  }
908  return(0);
909  };
910 
916  int parallel_blockstructured::ComputeJacobianMatrix(std::ostream &messageStream)
917  {
918 
919  int numDim = partitionIntervals.numDim;
920 
921  std::vector<size_t> &gridSizes(GridSizes());
922  std::vector<size_t> &bufferSizes(BufferSizes());
923  pcpp::IndexIntervalType &partInterval(PartitionInterval());
924  pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
925  std::vector<double> &gridCoordinates(CoordinateData());
926 
927  int numJacComponents = numDim*numDim;
928  size_t numPointsBuffer = BufferSize();
929  jacobianMatrix.resize(numJacComponents*numPointsBuffer);
930 
931  if(GradIJK(numDim,&gridCoordinates[0],&jacobianMatrix[0],
932  messageStream)){
933  messageStream << "grid::PBS::ComputeJacobianMatrix: Error, GradIJK failed."
934  << std::endl;
935  return(1);
936  }
937 
938  return(0);
939  };
940 
946  int parallel_blockstructured::ComputeJacobianMatrix2(std::ostream &messageStream)
947  {
948 
949  int numDim = partitionIntervals.numDim;
950 
951  int numComponents = numDim*numDim;
952  int numComponents2 = numDim*numComponents;
953  size_t numPointsBuffer = BufferSize();
954  jacobianMatrix2.resize(numComponents2*numPointsBuffer);
955 
956  if(GradIJK(numComponents,&jacobianMatrix[0],&jacobianMatrix2[0],
957  messageStream)){
958  messageStream << "grid::PBS::ComputeJacobianMatrix2: Error, GradIJK failed."
959  << std::endl;
960  return(1);
961  }
962 
963  return(0);
964  };
965 
966  int parallel_blockstructured::ComputeCurvilinearMetrics3D(std::ostream &messageStream){
967 
968  std::vector<size_t> &gridSizes(GridSizes());
969  std::vector<size_t> &bufferSizes(BufferSizes());
970 
971  pcpp::IndexIntervalType &partInterval(PartitionInterval());
972  pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
973 
974  if(gridSizes.empty() || bufferSizes.empty() ||
975  (partInterval.NNodes() == 0) || (partitionBufferInterval.NNodes() == 0)){
976  messageStream << "grid::PBS::ComputeCurvilinearMetrics: Error: Grid sizes not set up."
977  << std::endl
978  << "grid::PBS::ComputeCurvilinearMetrics: Finalize() must be called before ComputeMetrics."
979  << std::endl;
980  return(1);
981  }
982 
983  int numDim = bufferSizes.size();
984  if(numDim != 3){
985  messageStream << "grid::PBS::ComputeMetrics:Error: Cannot create 3D curvilinear metrics for "
986  << numDim << "-dimensional grid." << std::endl;
987  return(1);
988  }
989 
990  messageStream << "grid::PBS::ComputeMetrics: Computing 3d curvilinear metric."
991  << std::endl;
992 
993 
994  // this populates jacobianMatrix with d(xyz)/d(ijk)
995  if(ComputeJacobianMatrix(messageStream)){
996  messageStream << "grid::PBS::ComputeCurvilinearMetrics3D: Error: "
997  << "ComputeJacobianMatrix failed." << std::endl;
998  return(1);
999  }
1000 
1001  int numMetricComponents = numDim*numDim;
1002  size_t numPointsBuffer = BufferSize();
1003 
1004  gridMetric.resize(numMetricComponents*numPointsBuffer,0);
1005  gridJacobian.resize(2*numPointsBuffer,1.0);
1006 
1007 
1008  // Threading this routine would involve having the
1009  // thread id right here. It is not certain how
1010  // best to get do that. Perhaps two routines
1011  // are needed, one threaded, one not. Hokey.
1012  std::vector<size_t> dOpInterval;
1013  partitionBufferInterval.Flatten(dOpInterval);
1014  // forticate the opInterval
1015  std::vector<size_t>::iterator dOpIt = dOpInterval.begin();
1016  while(dOpIt != dOpInterval.end()){
1017  *dOpIt += 1;
1018  dOpIt++;
1019  }
1020 
1021  // Get 1/(Jacobian determinant)
1022  FC_MODULE(operators,determinant3x3,OPERATORS,DETERMINANT3X3)
1023  (&numPointsBuffer,&bufferSizes[0],
1024  &dOpInterval[0],&jacobianMatrix[0],
1025  &gridJacobian[numPointsBuffer]);
1026 
1027  // Get J = 1 / (1/J)
1028  double dOne = 1.0;
1029  FC_MODULE(operators,yaxm1,OPERATORS,YAXM1)
1030  (&numDim,&numPointsBuffer,&bufferSizes[0],
1031  &dOpInterval[0],&dOne,&gridJacobian[numPointsBuffer],
1032  &gridJacobian[0]);
1033 
1034 
1035  if(ExchangeJacobianMatrix(messageStream)){
1036  messageStream << "grid::PBS::ComputeCurvilinearMetrics3D: Error: Jacobian matrix exchange failed."
1037  << std::endl;
1038  return(1);
1039  }
1040 
1041  // this populates jacobianMatrix2 with d^2(xyz)/d(ijk)^2
1042  if(ComputeJacobianMatrix2(messageStream)){
1043  messageStream << "grid::PBS::ComputeCurvilinearMetrics3D: Error: "
1044  << "ComputeJacobianMatrix2 failed." << std::endl;
1045  return(1);
1046  }
1047 
1048  size_t componentOffset = 0;
1049  int jm1Component = 1;
1050  componentOffset = (jm1Component - 1)*numPointsBuffer;
1051  double *x_xi = &jacobianMatrix[componentOffset];
1052  jm1Component = 2;
1053  componentOffset = (jm1Component - 1)*numPointsBuffer;
1054  double *x_eta = &jacobianMatrix[componentOffset];
1055  jm1Component = 3;
1056  componentOffset = (jm1Component - 1)*numPointsBuffer;
1057  double *x_zeta = &jacobianMatrix[componentOffset];
1058  jm1Component = 4;
1059  componentOffset = (jm1Component - 1)*numPointsBuffer;
1060  double *y_xi = &jacobianMatrix[componentOffset];
1061  jm1Component = 5;
1062  componentOffset = (jm1Component - 1)*numPointsBuffer;
1063  double *y_eta = &jacobianMatrix[componentOffset];
1064  jm1Component = 6;
1065  componentOffset = (jm1Component - 1)*numPointsBuffer;
1066  double *y_zeta = &jacobianMatrix[componentOffset];
1067  jm1Component = 7;
1068  componentOffset = (jm1Component - 1)*numPointsBuffer;
1069  double *z_xi = &jacobianMatrix[componentOffset];
1070  jm1Component = 8;
1071  componentOffset = (jm1Component - 1)*numPointsBuffer;
1072  double *z_eta = &jacobianMatrix[componentOffset];
1073  jm1Component = 9;
1074  componentOffset = (jm1Component - 1)*numPointsBuffer;
1075  double *z_zeta= &jacobianMatrix[componentOffset];
1076 
1077  int j2Component = 2;
1078  componentOffset = (j2Component-1)*numPointsBuffer;
1079  double *x_xi_eta = &jacobianMatrix2[componentOffset];
1080 
1081  j2Component = 3;
1082  componentOffset = (j2Component-1)*numPointsBuffer;
1083  double *x_xi_zeta = &jacobianMatrix2[componentOffset];
1084 
1085  j2Component = 4;
1086  componentOffset = (j2Component-1)*numPointsBuffer;
1087  double *x_eta_xi = &jacobianMatrix2[componentOffset];
1088 
1089  j2Component = 6;
1090  componentOffset = (j2Component-1)*numPointsBuffer;
1091  double *x_eta_zeta = &jacobianMatrix2[componentOffset];
1092 
1093  j2Component = 7;
1094  componentOffset = (j2Component-1)*numPointsBuffer;
1095  double *x_zeta_xi = &jacobianMatrix2[componentOffset];
1096 
1097  j2Component = 8;
1098  componentOffset = (j2Component-1)*numPointsBuffer;
1099  double *x_zeta_eta = &jacobianMatrix2[componentOffset];
1100 
1101  j2Component = 11;
1102  componentOffset = (j2Component-1)*numPointsBuffer;
1103  double *y_xi_eta = &jacobianMatrix2[componentOffset];
1104 
1105  j2Component = 12;
1106  componentOffset = (j2Component-1)*numPointsBuffer;
1107  double *y_xi_zeta= &jacobianMatrix2[componentOffset];
1108 
1109  j2Component = 13;
1110  componentOffset = (j2Component-1)*numPointsBuffer;
1111  double *y_eta_xi = &jacobianMatrix2[componentOffset];
1112 
1113  j2Component = 15;
1114  componentOffset = (j2Component-1)*numPointsBuffer;
1115  double *y_eta_zeta= &jacobianMatrix2[componentOffset];
1116 
1117  j2Component = 16;
1118  componentOffset = (j2Component-1)*numPointsBuffer;
1119  double *y_zeta_xi = &jacobianMatrix2[componentOffset];
1120 
1121  j2Component = 17;
1122  componentOffset = (j2Component-1)*numPointsBuffer;
1123  double *y_zeta_eta = &jacobianMatrix2[componentOffset];
1124 
1125  j2Component = 20;
1126  componentOffset = (j2Component-1)*numPointsBuffer;
1127  double *z_xi_eta = &jacobianMatrix2[componentOffset];
1128 
1129  j2Component = 21;
1130  componentOffset = (j2Component-1)*numPointsBuffer;
1131  double *z_xi_zeta = &jacobianMatrix2[componentOffset];
1132 
1133  j2Component = 22;
1134  componentOffset = (j2Component-1)*numPointsBuffer;
1135  double *z_eta_xi = &jacobianMatrix2[componentOffset];
1136 
1137  j2Component = 24;
1138  componentOffset = (j2Component-1)*numPointsBuffer;
1139  double *z_eta_zeta = &jacobianMatrix2[componentOffset];
1140 
1141  j2Component = 25;
1142  componentOffset = (j2Component-1)*numPointsBuffer;
1143  double *z_zeta_xi = &jacobianMatrix2[componentOffset];
1144 
1145  j2Component = 26;
1146  componentOffset = (j2Component-1)*numPointsBuffer;
1147  double *z_zeta_eta = &jacobianMatrix2[componentOffset];
1148 
1149  double *xData = &gridCoordinates[0];
1150  double *yData = &gridCoordinates[numPointsBuffer];
1151  double *zData = &gridCoordinates[2*numPointsBuffer];
1152 
1153  // operator metricsum4
1156  // met[0] = y_eta_zeta*z + y_eta*z_zeta - y_zeta_eta*z - y_zeta*z_eta;
1157  FC_MODULE(operators,metricsum4,OPERATORS,METRICSUM4)
1158  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1159  y_eta,y_zeta,z_eta,z_zeta,y_eta_zeta,y_zeta_eta,zData,
1160  &gridMetric[0]);
1161  // met[1] = z_eta_zeta*x + z_eta*x_zeta - z_zeta_eta*x - z_zeta*x_eta;
1162  FC_MODULE(operators,metricsum4,OPERATORS,METRICSUM4)
1163  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1164  z_eta,z_zeta,x_eta,x_zeta,z_eta_zeta,z_zeta_eta,xData,
1165  &gridMetric[numPointsBuffer]);
1166  // met[2] = x_eta_zeta*y + x_eta*y_zeta - x_zeta_eta*y - x_zeta*y_eta;
1167  FC_MODULE(operators,metricsum4,OPERATORS,METRICSUM4)
1168  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1169  x_eta,x_zeta,y_eta,y_zeta,x_eta_zeta,x_zeta_eta,yData,
1170  &gridMetric[2*numPointsBuffer]);
1171  // met[3] = y_zeta_xi*z + y_zeta*z_xi - y_xi_zeta*z - y_xi*z_zeta;
1172  FC_MODULE(operators,metricsum4,OPERATORS,METRICSUM4)
1173  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1174  y_zeta,y_xi,z_zeta,z_xi,y_zeta_xi,y_xi_zeta,zData,
1175  &gridMetric[3*numPointsBuffer]);
1176  // met[4] = z_zeta_xi*x + z_zeta*x_xi - z_xi_zeta*x - z_xi*x_zeta;
1177  FC_MODULE(operators,metricsum4,OPERATORS,METRICSUM4)
1178  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1179  z_zeta,z_xi,x_zeta,x_xi,z_zeta_xi,z_xi_zeta,xData,
1180  &gridMetric[4*numPointsBuffer]);
1181  // met[5] = x_zeta_xi*y + x_zeta*y_xi - x_xi_zeta*y - x_xi*y_zeta;
1182  FC_MODULE(operators,metricsum4,OPERATORS,METRICSUM4)
1183  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1184  x_zeta,x_xi,y_zeta,y_xi,x_zeta_xi,x_xi_zeta,yData,
1185  &gridMetric[5*numPointsBuffer]);
1186  // met[6] = y_xi_eta*z + y_xi*z_eta - y_eta_xi*z - y_eta*z_xi;
1187  FC_MODULE(operators,metricsum4,OPERATORS,METRICSUM4)
1188  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1189  y_xi,y_eta,z_xi,z_eta,y_xi_eta,y_eta_xi,zData,
1190  &gridMetric[6*numPointsBuffer]);
1191  // met[7] = z_xi_eta*x + z_xi*x_eta - z_eta_xi*x - z_eta*x_xi;
1192  FC_MODULE(operators,metricsum4,OPERATORS,METRICSUM4)
1193  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1194  z_xi,z_eta,x_xi,x_eta,z_xi_eta,z_eta_xi,xData,
1195  &gridMetric[7*numPointsBuffer]);
1196  // met[8] = x_xi_eta*y + x_xi*y_eta - x_eta_xi*y - x_eta*y_xi;
1197  FC_MODULE(operators,metricsum4,OPERATORS,METRICSUM4)
1198  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1199  x_xi,x_eta,y_xi,y_eta,x_xi_eta,x_eta_xi,yData,
1200  &gridMetric[8*numPointsBuffer]);
1201 
1202  return(0);
1203  };
1204 
1205  int parallel_blockstructured::CurvilinearMetrics(std::ostream &messageStream){
1206 
1207  std::vector<size_t> &gridSizes(GridSizes());
1208  std::vector<size_t> &bufferSizes(BufferSizes());
1209  pcpp::IndexIntervalType bufferInterval;
1210  bufferInterval.InitSimple(bufferSizes);
1211  std::vector<size_t> fortranBufferInterval;
1212  bufferInterval.Flatten(fortranBufferInterval);
1213  std::vector<size_t>::iterator fbIt = fortranBufferInterval.begin();
1214  while(fbIt != fortranBufferInterval.end()){
1215  *fbIt = *fbIt + 1;
1216  fbIt++;
1217  }
1218 
1219  pcpp::IndexIntervalType &partInterval(PartitionInterval());
1220  pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
1221 
1222  if(gridSizes.empty() || bufferSizes.empty() ||
1223  (partInterval.NNodes() == 0) || (partitionBufferInterval.NNodes() == 0)){
1224  messageStream << "grid::PBS::ComputeCurvilinearMetrics: Error: Grid sizes not set up."
1225  << std::endl
1226  << "grid::PBS::ComputeCurvilinearMetrics: Finalize() must be called before ComputeMetrics."
1227  << std::endl;
1228  return(1);
1229  }
1230  int one = 1;
1231  double negOne = -1.0;
1232  double plusOne = 1.0;
1233 
1234  int numDim = bufferSizes.size();
1235  if(numDim != 3){
1236  messageStream << "grid::PBS::ComputeMetrics:Error: Cannot create 3D curvilinear metrics for "
1237  << numDim << "-dimensional grid." << std::endl;
1238  return(1);
1239  }
1240 
1241  messageStream << "grid::PBS::ComputeMetrics: Computing 3d curvilinear metric."
1242  << std::endl;
1243 
1244 
1245  // this populates jacobianMatrix with d(xyz)/d(ijk)
1246  if(ComputeJacobianMatrix(messageStream)){
1247  messageStream << "grid::PBS::CurvilinearMetrics: Error: "
1248  << "ComputeJacobianMatrix failed." << std::endl;
1249  return(1);
1250  }
1251 
1252  int numMetricComponents = numDim*numDim;
1253  size_t numPointsBuffer = BufferSize();
1254 
1255  gridMetric.resize(numMetricComponents*numPointsBuffer,0);
1256  gridJacobian.resize(2*numPointsBuffer,1.0);
1257 
1258 
1259  // Threading this routine would involve having the
1260  // thread id right here. It is not certain how
1261  // best to get do that. Perhaps two routines
1262  // are needed, one threaded, one not. Hokey.
1263  std::vector<size_t> dOpInterval;
1264  partitionBufferInterval.Flatten(dOpInterval);
1265  // forticate the opInterval
1266  std::vector<size_t>::iterator dOpIt = dOpInterval.begin();
1267  while(dOpIt != dOpInterval.end()){
1268  *dOpIt += 1;
1269  dOpIt++;
1270  }
1271 
1272  // Get 1/(Jacobian determinant)
1273  FC_MODULE(operators,determinant3x3,OPERATORS,DETERMINANT3X3)
1274  (&numPointsBuffer,&bufferSizes[0],
1275  &dOpInterval[0],&jacobianMatrix[0],
1276  &gridJacobian[numPointsBuffer]);
1277 
1278  // Get J = 1 / (1/J)
1279  double dOne = 1.0;
1280  FC_MODULE(operators,yaxm1,OPERATORS,YAXM1)
1281  (&numDim,&numPointsBuffer,&bufferSizes[0],
1282  &dOpInterval[0],&dOne,&gridJacobian[numPointsBuffer],
1283  &gridJacobian[0]);
1284 
1285 
1286  if(ExchangeJacobianMatrix(messageStream)){
1287  messageStream << "grid::PBS::ComputeCurvilinearMetrics3D: Error: Jacobian matrix exchange failed."
1288  << std::endl;
1289  return(1);
1290  }
1291 
1292  int numStencils = stencilSet->numStencils;
1293  int *stencilSizes = stencilSet->stencilSizes;
1294  int *stencilStarts = stencilSet->stencilStarts;
1295  int *stencilOffsets = stencilSet->stencilOffsets;
1296  double *stencilWeights = stencilSet->stencilWeights;
1297  int numStencilValues = stencilSet->numValues;
1298 
1299  size_t componentOffset = 0;
1300  int jm1Component = 1;
1301  componentOffset = (jm1Component - 1)*numPointsBuffer;
1302  double *x_xi = &jacobianMatrix[componentOffset];
1303  jm1Component = 2;
1304  componentOffset = (jm1Component - 1)*numPointsBuffer;
1305  double *x_eta = &jacobianMatrix[componentOffset];
1306  jm1Component = 3;
1307  componentOffset = (jm1Component - 1)*numPointsBuffer;
1308  double *x_zeta = &jacobianMatrix[componentOffset];
1309  jm1Component = 4;
1310  componentOffset = (jm1Component - 1)*numPointsBuffer;
1311  double *y_xi = &jacobianMatrix[componentOffset];
1312  jm1Component = 5;
1313  componentOffset = (jm1Component - 1)*numPointsBuffer;
1314  double *y_eta = &jacobianMatrix[componentOffset];
1315  jm1Component = 6;
1316  componentOffset = (jm1Component - 1)*numPointsBuffer;
1317  double *y_zeta = &jacobianMatrix[componentOffset];
1318  jm1Component = 7;
1319  componentOffset = (jm1Component - 1)*numPointsBuffer;
1320  double *z_xi = &jacobianMatrix[componentOffset];
1321  jm1Component = 8;
1322  componentOffset = (jm1Component - 1)*numPointsBuffer;
1323  double *z_eta = &jacobianMatrix[componentOffset];
1324  jm1Component = 9;
1325  componentOffset = (jm1Component - 1)*numPointsBuffer;
1326  double *z_zeta= &jacobianMatrix[componentOffset];
1327 
1328 
1329  double *xData = &gridCoordinates[0];
1330  double *yData = &gridCoordinates[numPointsBuffer];
1331  double *zData = &gridCoordinates[2*numPointsBuffer];
1332 
1333  // temporary storage
1334  jacobianMatrix2.resize(numDim*numPointsBuffer,0.0);
1335 
1336  // ! ... spatial metrics using Thomas & Lombard
1337 
1338  // ! ... xihat_x = (y_eta z)_zeta - (y_zeta z)_eta
1339  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1340  (&numDim,&numPointsBuffer,&bufferSizes[0],
1341  &fortranBufferInterval[0],y_eta,zData,&jacobianMatrix2[0]);
1342 
1343  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1344  (&numDim,&numPointsBuffer,&bufferSizes[0],
1345  &fortranBufferInterval[0],y_zeta,zData,&jacobianMatrix2[numPointsBuffer]);
1346 
1347  int dDir = 3;
1348  int metricComponent = 1;
1349  size_t dirOffset = (dDir-1)*numPointsBuffer;
1350  size_t metricOffset = (metricComponent-1)*numPointsBuffer;
1351  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1352  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1353  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1354  stencilStarts,&numStencilValues,stencilWeights,
1355  stencilOffsets,&stencilConnectivity[dirOffset],
1356  &jacobianMatrix2[0],&gridMetric[metricOffset]);
1357 
1358  dDir = 2;
1359  dirOffset = (dDir-1)*numPointsBuffer;
1360  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1361  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1362  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1363  stencilStarts,&numStencilValues,stencilWeights,
1364  stencilOffsets,&stencilConnectivity[dirOffset],
1365  &jacobianMatrix2[numPointsBuffer],&jacobianMatrix2[2*numPointsBuffer]);
1366  FC_MODULE(operators,yaxpy,OPERATORS,YAXPY)
1367  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1368  &negOne,&jacobianMatrix2[2*numPointsBuffer],&gridMetric[metricOffset]);
1369 
1370  // ! ... xihat_y = (z_eta x)_zeta - (z_zeta x)_eta
1371  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1372  (&numDim,&numPointsBuffer,&bufferSizes[0],
1373  &fortranBufferInterval[0],z_eta,xData,&jacobianMatrix2[0]);
1374  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1375  (&numDim,&numPointsBuffer,&bufferSizes[0],
1376  &fortranBufferInterval[0],z_zeta,xData,&jacobianMatrix2[numPointsBuffer]);
1377 
1378  dDir = 3;
1379  metricComponent = 2;
1380  dirOffset = (dDir-1)*numPointsBuffer;
1381  metricOffset = (metricComponent-1)*numPointsBuffer;
1382  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1383  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1384  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1385  stencilStarts,&numStencilValues,stencilWeights,
1386  stencilOffsets,&stencilConnectivity[dirOffset],
1387  &jacobianMatrix2[0],&gridMetric[metricOffset]);
1388  dDir = 2;
1389  dirOffset = (dDir-1)*numPointsBuffer;
1390  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1391  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1392  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1393  stencilStarts,&numStencilValues,stencilWeights,
1394  stencilOffsets,&stencilConnectivity[dirOffset],
1395  &jacobianMatrix2[numPointsBuffer],&jacobianMatrix2[2*numPointsBuffer]);
1396  FC_MODULE(operators,yaxpy,OPERATORS,YAXPY)
1397  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1398  &negOne,&jacobianMatrix2[2*numPointsBuffer],&gridMetric[metricOffset]);
1399 
1400  // ! ... xihat_z = (x_eta y)_zeta - (x_zeta y)_eta
1401  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1402  (&numDim,&numPointsBuffer,&bufferSizes[0],
1403  &fortranBufferInterval[0],x_eta,yData,&jacobianMatrix2[0]);
1404  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1405  (&numDim,&numPointsBuffer,&bufferSizes[0],
1406  &fortranBufferInterval[0],x_zeta,yData,&jacobianMatrix2[numPointsBuffer]);
1407 
1408  dDir = 3;
1409  metricComponent = 3;
1410  dirOffset = (dDir-1)*numPointsBuffer;
1411  metricOffset = (metricComponent-1)*numPointsBuffer;
1412  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1413  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1414  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1415  stencilStarts,&numStencilValues,stencilWeights,
1416  stencilOffsets,&stencilConnectivity[dirOffset],
1417  &jacobianMatrix2[0],&gridMetric[metricOffset]);
1418  dDir = 2;
1419  dirOffset = (dDir-1)*numPointsBuffer;
1420  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1421  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1422  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1423  stencilStarts,&numStencilValues,stencilWeights,
1424  stencilOffsets,&stencilConnectivity[dirOffset],
1425  &jacobianMatrix2[numPointsBuffer],&jacobianMatrix2[2*numPointsBuffer]);
1426  FC_MODULE(operators,yaxpy,OPERATORS,YAXPY)
1427  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1428  &negOne,&jacobianMatrix2[2*numPointsBuffer],&gridMetric[metricOffset]);
1429 
1430 
1431 
1432 
1433  // ! ... etahat_x = (y_zeta z)_xi - (y_xi z)_zeta
1434  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1435  (&numDim,&numPointsBuffer,&bufferSizes[0],
1436  &fortranBufferInterval[0],y_zeta,zData,&jacobianMatrix2[0]);
1437  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1438  (&numDim,&numPointsBuffer,&bufferSizes[0],
1439  &fortranBufferInterval[0],y_xi,zData,&jacobianMatrix2[numPointsBuffer]);
1440 
1441  dDir = 1;
1442  metricComponent = 4;
1443  dirOffset = (dDir-1)*numPointsBuffer;
1444  metricOffset = (metricComponent-1)*numPointsBuffer;
1445  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1446  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1447  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1448  stencilStarts,&numStencilValues,stencilWeights,
1449  stencilOffsets,&stencilConnectivity[dirOffset],
1450  &jacobianMatrix2[0],&gridMetric[metricOffset]);
1451  dDir = 3;
1452  dirOffset = (dDir-1)*numPointsBuffer;
1453  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1454  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1455  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1456  stencilStarts,&numStencilValues,stencilWeights,
1457  stencilOffsets,&stencilConnectivity[dirOffset],
1458  &jacobianMatrix2[numPointsBuffer],&jacobianMatrix2[2*numPointsBuffer]);
1459  FC_MODULE(operators,yaxpy,OPERATORS,YAXPY)
1460  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1461  &negOne,&jacobianMatrix2[2*numPointsBuffer],&gridMetric[metricOffset]);
1462 
1463 
1464 
1465  // ! ... etahat_y = (z_zeta x)_xi - (z_xi x)_zeta
1466  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1467  (&numDim,&numPointsBuffer,&bufferSizes[0],
1468  &fortranBufferInterval[0],z_zeta,xData,&jacobianMatrix2[0]);
1469  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1470  (&numDim,&numPointsBuffer,&bufferSizes[0],
1471  &fortranBufferInterval[0],z_xi,xData,&jacobianMatrix2[numPointsBuffer]);
1472 
1473  dDir = 1;
1474  metricComponent = 5;
1475  dirOffset = (dDir-1)*numPointsBuffer;
1476  metricOffset = (metricComponent-1)*numPointsBuffer;
1477  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1478  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1479  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1480  stencilStarts,&numStencilValues,stencilWeights,
1481  stencilOffsets,&stencilConnectivity[dirOffset],
1482  &jacobianMatrix2[0],&gridMetric[metricOffset]);
1483  dDir = 3;
1484  dirOffset = (dDir-1)*numPointsBuffer;
1485  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1486  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1487  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1488  stencilStarts,&numStencilValues,stencilWeights,
1489  stencilOffsets,&stencilConnectivity[dirOffset],
1490  &jacobianMatrix2[numPointsBuffer],&jacobianMatrix2[2*numPointsBuffer]);
1491  FC_MODULE(operators,yaxpy,OPERATORS,YAXPY)
1492  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1493  &negOne,&jacobianMatrix2[2*numPointsBuffer],&gridMetric[metricOffset]);
1494 
1495  // ! ... etahat_z = (x_zeta y)_xi - (x_xi y)_zeta
1496  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1497  (&numDim,&numPointsBuffer,&bufferSizes[0],
1498  &fortranBufferInterval[0],x_zeta,yData,&jacobianMatrix2[0]);
1499  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1500  (&numDim,&numPointsBuffer,&bufferSizes[0],
1501  &fortranBufferInterval[0],x_xi,yData,&jacobianMatrix2[numPointsBuffer]);
1502 
1503  dDir = 1;
1504  metricComponent = 6;
1505  dirOffset = (dDir-1)*numPointsBuffer;
1506  metricOffset = (metricComponent-1)*numPointsBuffer;
1507  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1508  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1509  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1510  stencilStarts,&numStencilValues,stencilWeights,
1511  stencilOffsets,&stencilConnectivity[dirOffset],
1512  &jacobianMatrix2[0],&gridMetric[metricOffset]);
1513  dDir = 3;
1514  dirOffset = (dDir-1)*numPointsBuffer;
1515  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1516  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1517  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1518  stencilStarts,&numStencilValues,stencilWeights,
1519  stencilOffsets,&stencilConnectivity[dirOffset],
1520  &jacobianMatrix2[numPointsBuffer],&jacobianMatrix2[2*numPointsBuffer]);
1521  FC_MODULE(operators,yaxpy,OPERATORS,YAXPY)
1522  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1523  &negOne,&jacobianMatrix2[2*numPointsBuffer],&gridMetric[metricOffset]);
1524 
1525 
1526  // ! ... zetahat_x = (y_xi z)_eta - (y_eta z)_xi
1527  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1528  (&numDim,&numPointsBuffer,&bufferSizes[0],
1529  &fortranBufferInterval[0],y_xi,zData,&jacobianMatrix2[0]);
1530  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1531  (&numDim,&numPointsBuffer,&bufferSizes[0],
1532  &fortranBufferInterval[0],y_eta,zData,&jacobianMatrix2[numPointsBuffer]);
1533 
1534  dDir = 2;
1535  metricComponent = 7;
1536  dirOffset = (dDir-1)*numPointsBuffer;
1537  metricOffset = (metricComponent-1)*numPointsBuffer;
1538  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1539  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1540  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1541  stencilStarts,&numStencilValues,stencilWeights,
1542  stencilOffsets,&stencilConnectivity[dirOffset],
1543  &jacobianMatrix2[0],&gridMetric[metricOffset]);
1544  dDir = 1;
1545  dirOffset = (dDir-1)*numPointsBuffer;
1546  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1547  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1548  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1549  stencilStarts,&numStencilValues,stencilWeights,
1550  stencilOffsets,&stencilConnectivity[dirOffset],
1551  &jacobianMatrix2[numPointsBuffer],&jacobianMatrix2[2*numPointsBuffer]);
1552  FC_MODULE(operators,yaxpy,OPERATORS,YAXPY)
1553  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1554  &negOne,&jacobianMatrix2[2*numPointsBuffer],&gridMetric[metricOffset]);
1555 
1556  // ! ... zetahat_y = (z_xi x)_eta - (z_eta x)_xi
1557  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1558  (&numDim,&numPointsBuffer,&bufferSizes[0],
1559  &fortranBufferInterval[0],z_xi,xData,&jacobianMatrix2[0]);
1560  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1561  (&numDim,&numPointsBuffer,&bufferSizes[0],
1562  &fortranBufferInterval[0],z_eta,xData,&jacobianMatrix2[numPointsBuffer]);
1563 
1564  dDir = 2;
1565  metricComponent = 8;
1566  dirOffset = (dDir-1)*numPointsBuffer;
1567  metricOffset = (metricComponent-1)*numPointsBuffer;
1568  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1569  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1570  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1571  stencilStarts,&numStencilValues,stencilWeights,
1572  stencilOffsets,&stencilConnectivity[dirOffset],
1573  &jacobianMatrix2[0],&gridMetric[metricOffset]);
1574  dDir = 1;
1575  dirOffset = (dDir-1)*numPointsBuffer;
1576  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1577  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1578  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1579  stencilStarts,&numStencilValues,stencilWeights,
1580  stencilOffsets,&stencilConnectivity[dirOffset],
1581  &jacobianMatrix2[numPointsBuffer],&jacobianMatrix2[2*numPointsBuffer]);
1582  FC_MODULE(operators,yaxpy,OPERATORS,YAXPY)
1583  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1584  &negOne,&jacobianMatrix2[2*numPointsBuffer],&gridMetric[metricOffset]);
1585 
1586  // ! ... zetahat_z = (x_xi y)_eta - (x_eta y)_xi
1587  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1588  (&numDim,&numPointsBuffer,&bufferSizes[0],
1589  &fortranBufferInterval[0],x_xi,yData,&jacobianMatrix2[0]);
1590  FC_MODULE(operators,zxy,OPERATORS,ZXY)
1591  (&numDim,&numPointsBuffer,&bufferSizes[0],
1592  &fortranBufferInterval[0],x_eta,yData,&jacobianMatrix2[numPointsBuffer]);
1593 
1594  dDir = 2;
1595  metricComponent = 9;
1596  dirOffset = (dDir-1)*numPointsBuffer;
1597  metricOffset = (metricComponent-1)*numPointsBuffer;
1598  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1599  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1600  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1601  stencilStarts,&numStencilValues,stencilWeights,
1602  stencilOffsets,&stencilConnectivity[dirOffset],
1603  &jacobianMatrix2[0],&gridMetric[metricOffset]);
1604  dDir = 1;
1605  dirOffset = (dDir-1)*numPointsBuffer;
1606  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
1607  (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
1608  &dDir,&dOpInterval[0],&numStencils,stencilSizes,
1609  stencilStarts,&numStencilValues,stencilWeights,
1610  stencilOffsets,&stencilConnectivity[dirOffset],
1611  &jacobianMatrix2[numPointsBuffer],&jacobianMatrix2[2*numPointsBuffer]);
1612  FC_MODULE(operators,yaxpy,OPERATORS,YAXPY)
1613  (&numDim,&numPointsBuffer,&bufferSizes[0],&dOpInterval[0],
1614  &negOne,&jacobianMatrix2[2*numPointsBuffer],&gridMetric[metricOffset]);
1615 
1616  return(0);
1617  };
1618 
1619  int parallel_blockstructured::SetupDifferentialOperator(plascom2::operators::stencilset *inStencilSet,
1620  int *inStencilConnectivity)
1621  {
1622  stencilSet = inStencilSet;
1623  stencilConnectivity = inStencilConnectivity;
1624  return(0);
1625  };
1626 
1627 
1628  int parallel_blockstructured::GenerateCoordinates(std::ostream &messageStream)
1629  {
1630 
1631  const std::vector<size_t> &gridSizes(GridSizes());
1632  const std::vector<size_t> &bufferSizes(BufferSizes());
1633  std::vector<double> &gridExtents(PhysicalExtent());
1634 
1635  pcpp::IndexIntervalType &partInterval(PartitionInterval());
1636  pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
1637 
1638  if(gridSizes.empty() || bufferSizes.empty() ||
1639  (partInterval.NNodes() == 0) || (partitionBufferInterval.NNodes() == 0)){
1640  messageStream << "grid::PBS::GenerateCoordinates: Error: Grid sizes not set up."
1641  << std::endl
1642  << "grid::PBS::GenerateCoordinates: Finalize() must be called before GenerateCoordinates."
1643  << std::endl;
1644  return(1);
1645  }
1646 
1647  int numDim = gridSizes.size();
1648 
1649  if(gridExtents.size() != 2*numDim){
1650  messageStream << "grid::PBS::GenerateCoordinates: Warning: Defaulting to physical extent [0:1] for all dimensions."
1651  << std::endl;
1652  gridExtents.resize(2*numDim,0);
1653  dX.resize(numDim,0);
1654  for(int iDim = 0;iDim < numDim;iDim++){
1655  gridExtents[2*iDim+1] = 1.0;
1656  dX[iDim] = 1.0/(double(gridSizes[iDim]-1));
1657  }
1658  }
1659 
1660  if(dX.size() != numDim){
1661  dX.resize(numDim,0);
1662  for(int iDim = 0;iDim < numDim;iDim++){
1663  dX[iDim] = (gridExtents[2*iDim+1] - gridExtents[2*iDim])/(double(gridSizes[iDim]-1));
1664  }
1665  }
1666 
1667  AllocateCoordinateData();
1668 
1669  if(scaleFactor.size() == numDim){
1670  messageStream << "grid::PBS::GenerateCoordinates: Using scalefactors (";
1671  pcpp::io::DumpContents(messageStream,scaleFactor,",");
1672  messageStream << ")" << std::endl;
1673  for(int iDim = 0;iDim < numDim;iDim++){
1674  dX[iDim] *= scaleFactor[iDim];
1675  physicalExtent[2*iDim] *= scaleFactor[iDim];
1676  physicalExtent[2*iDim+1] *= scaleFactor[iDim];
1677  }
1678  }
1679 
1680  bool isCartesian = true;
1681  for(int iDim = 0;iDim < numDim-1;iDim++){
1682  if(dX[iDim] != dX[iDim+1])
1683  isCartesian = false;
1684  }
1685 
1686  if(isCartesian) {
1687  gridType = simulation::grid::CARTESIAN;
1688  } else {
1689  gridType = simulation::grid::UNIRECT;
1690  }
1691 
1692  if(periodicDirs.size() == numDim){
1693  if(periodicLengths.size() != numDim){
1694  periodicLengths.resize(numDim,0);
1695  for(int iDim = 0;iDim < numDim;iDim++){
1696  if(periodicDirs[iDim]){
1697  periodicLengths[iDim] = (physicalExtent[2*iDim+1]-physicalExtent[2*iDim]) + dX[iDim];
1698  }
1699  }
1700  }
1701  }
1702 
1703  messageStream << "grid::PBS::GenerateCoordinates: Generating coordinates for "
1704  << (gridType == simulation::grid::CARTESIAN ? " Cartesian " : " Uniform-rectangular ")
1705  << "grid." << std::endl;
1706 
1707  if(numDim == 3){
1708 
1709  size_t bufferPlane = bufferSizes[0]*bufferSizes[1];
1710  size_t bufferLine = bufferSizes[0];
1711 
1712  size_t xStart = partInterval[0].first;
1713  size_t xEnd = partInterval[0].second;
1714  size_t yStart = partInterval[1].first;
1715  size_t yEnd = partInterval[1].second;
1716  size_t zStart = partInterval[2].first;
1717  size_t zEnd = partInterval[2].second;
1718 
1719  size_t iStart = partitionBufferInterval[0].first;
1720  size_t iEnd = partitionBufferInterval[0].second;
1721  size_t jStart = partitionBufferInterval[1].first;
1722  size_t jEnd = partitionBufferInterval[1].second;
1723  size_t kStart = partitionBufferInterval[2].first;
1724  size_t kEnd = partitionBufferInterval[2].second;
1725 
1726  size_t kIndex = kStart;
1727  size_t zIndex = zStart;
1728 
1729  for(kIndex = kStart,zIndex = zStart;kIndex <= kEnd;kIndex++,zIndex++){
1730  size_t jIndex = jStart;
1731  size_t yIndex = yStart;
1732  size_t bzIndex = kIndex * (bufferPlane);
1733  double zCoordinate = (zIndex*dX[2] + physicalExtent[4]);
1734  for(jIndex = jStart,yIndex = yStart;jIndex <= jEnd;jIndex++,yIndex++){
1735  size_t iIndex = iStart;
1736  size_t xIndex = xStart;
1737  size_t byzIndex = bzIndex + jIndex*bufferLine;
1738  double yCoordinate = (yIndex*dX[1] + physicalExtent[2]);
1739  for(iIndex = iStart,xIndex = xStart;iIndex <= iEnd;iIndex++,xIndex++){
1740  size_t bufferIndex = byzIndex + iIndex;
1741  double xCoordinate = (xIndex*dX[0] + physicalExtent[0]);
1742  gridCoordinates[bufferIndex] = xCoordinate;
1743  gridCoordinates[bufferIndex+numPointsBuffer] = yCoordinate;
1744  gridCoordinates[bufferIndex+2*numPointsBuffer] = zCoordinate;
1745  }
1746  }
1747  }
1748  } else if (numDim == 2) {
1749 
1750  size_t bufferLine = bufferSizes[0];
1751 
1752  size_t xStart = partInterval[0].first;
1753  size_t xEnd = partInterval[0].second;
1754  size_t yStart = partInterval[1].first;
1755  size_t yEnd = partInterval[1].second;
1756 
1757  size_t iStart = partitionBufferInterval[0].first;
1758  size_t iEnd = partitionBufferInterval[0].second;
1759  size_t jStart = partitionBufferInterval[1].first;
1760  size_t jEnd = partitionBufferInterval[1].second;
1761 
1762  size_t jIndex = jStart;
1763  size_t yIndex = yStart;
1764  for(jIndex = jStart,yIndex = yStart;jIndex <= jEnd;jIndex++,yIndex++){
1765  size_t iIndex = iStart;
1766  size_t xIndex = xStart;
1767  size_t byIndex = jIndex*bufferLine;
1768  double yCoordinate = yIndex*dX[1] + physicalExtent[2];
1769  for(iIndex = iStart,xIndex = xStart;iIndex <= iEnd;iIndex++,xIndex++){
1770  size_t bufferIndex = byIndex + iIndex;
1771  double xCoordinate = xIndex*dX[0] + physicalExtent[0];
1772  gridCoordinates[bufferIndex] = xCoordinate;
1773  gridCoordinates[bufferIndex+numPointsBuffer] = yCoordinate;
1774  }
1775  }
1776  } else if (numDim == 1) {
1777 
1778 
1779  size_t xStart = partInterval[0].first;
1780  size_t xEnd = partInterval[0].second;
1781 
1782  size_t iStart = partitionBufferInterval[0].first;
1783  size_t iEnd = partitionBufferInterval[0].second;
1784 
1785  size_t iIndex = iStart;
1786  size_t xIndex = xStart;
1787  for(iIndex = iStart,xIndex = xStart;iIndex <= iEnd;iIndex++,xIndex++){
1788  gridCoordinates[iIndex] = xIndex*dX[0] + physicalExtent[0];
1789  }
1790  }
1791  return(0);
1792  };
1793 
1794  size_t parallel_blockstructured::FlagMask(int *iMask,int flagValue,int flagBit) const
1795  {
1796 
1797  const std::vector<size_t> &bufferSizes(BufferSizes());
1798  const pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
1799 
1800  if(bufferSizes.empty() ||
1801  (partitionBufferInterval.NNodes() == 0) ||
1802  flagData.empty()){
1803  return(0);
1804  }
1805 
1806  pcpp::IndexIntervalType bufferInterval;
1807  bufferInterval.InitSimple(bufferSizes);
1808  std::vector<size_t> myIndices;
1809  bufferInterval.GetFlatIndices(partitionBufferInterval,myIndices);
1810  std::vector<size_t>::const_iterator indexIt = myIndices.begin();
1811  size_t retVal = 0;
1812  while(indexIt != myIndices.end()){
1813  const size_t myIndex(*indexIt++);
1814  if(flagData[myIndex] == flagValue){
1815  iMask[myIndex] |= flagBit;
1816  retVal++;
1817  }
1818  }
1819  return(retVal);
1820  };
1821 
1822  int parallel_blockstructured::GenerateGrid(int (*GridGeneratorFunction)
1823  (const std::vector<int> &,
1824  const std::vector<double> &,
1825  const std::vector<size_t> &,
1826  const std::vector<size_t> &,
1827  std::vector<double> &),
1828  std::ostream &messageStream)
1829  {
1830 
1831  const std::vector<size_t> &gridSizes(GridSizes());
1832  const std::vector<size_t> &bufferSizes(BufferSizes());
1833  std::vector<double> &gridExtents(PhysicalExtent());
1834 
1835  pcpp::IndexIntervalType &partInterval(PartitionInterval());
1836  pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
1837 
1838  const std::vector<double> &gridGenerationParameters(GridGenerationParameters());
1839  const std::vector<int> &gridGenerationOptions(GridGenerationOptions());
1840 
1841  if(gridSizes.empty() || bufferSizes.empty() ||
1842  (partInterval.NNodes() == 0) || (partitionBufferInterval.NNodes() == 0)){
1843  messageStream << "grid::PBS::GenerateGrid: Error: Grid sizes not set up."
1844  << std::endl
1845  << "grid::PBS::GenerateGrid: Finalize() must be called before GenerateGrid."
1846  << std::endl;
1847  return(1);
1848  }
1849 
1850  int numDim = gridSizes.size();
1851 
1852  AllocateCoordinateData();
1853 
1854  std::vector<double> xyz(numDim,0);
1855  std::vector<size_t> ijk(numDim,0);
1856 
1857  size_t iIndex,jIndex,kIndex,xIndex,yIndex,zIndex;
1858 
1859  if(numDim == 3){
1860 
1861 
1862  size_t bufferPlane = bufferSizes[0]*bufferSizes[1];
1863  size_t bufferLine = bufferSizes[0];
1864 
1865  size_t xStart = partInterval[0].first;
1866  size_t xEnd = partInterval[0].second;
1867  size_t yStart = partInterval[1].first;
1868  size_t yEnd = partInterval[1].second;
1869  size_t zStart = partInterval[2].first;
1870  size_t zEnd = partInterval[2].second;
1871 
1872  size_t iStart = partitionBufferInterval[0].first;
1873  size_t iEnd = partitionBufferInterval[0].second;
1874  size_t jStart = partitionBufferInterval[1].first;
1875  size_t jEnd = partitionBufferInterval[1].second;
1876  size_t kStart = partitionBufferInterval[2].first;
1877  size_t kEnd = partitionBufferInterval[2].second;
1878 
1879  for(kIndex = kStart,zIndex = zStart;kIndex <= kEnd;kIndex++,zIndex++){
1880  size_t bzIndex = kIndex * (bufferPlane);
1881  ijk[2] = zIndex;
1882  for(jIndex = jStart,yIndex = yStart;jIndex <= jEnd;jIndex++,yIndex++){
1883  ijk[1] = yIndex;
1884  size_t byzIndex = bzIndex + jIndex*bufferLine;
1885  for(iIndex = iStart,xIndex = xStart;iIndex <= iEnd;iIndex++,xIndex++){
1886  ijk[0] = xIndex;
1887  size_t bufferIndex = byzIndex + iIndex;
1888  int errorCode = GridGeneratorFunction(gridGenerationOptions,gridGenerationParameters,gridSizes,ijk,xyz);
1889  if(errorCode){
1890  messageStream << "grid::GenerateGrid: User-specified GridGeneratorFunction"
1891  << " returned error code (" << errorCode << "). Aborting." << std::endl;
1892  return(1);
1893  }
1894  gridCoordinates[bufferIndex] = xyz[0];
1895  gridCoordinates[bufferIndex+numPointsBuffer] = xyz[1];
1896  gridCoordinates[bufferIndex+2*numPointsBuffer] = xyz[2];
1897  }
1898  }
1899  }
1900  } else if (numDim == 2) {
1901 
1902  size_t bufferLine = bufferSizes[0];
1903 
1904  size_t xStart = partInterval[0].first;
1905  size_t xEnd = partInterval[0].second;
1906  size_t yStart = partInterval[1].first;
1907  size_t yEnd = partInterval[1].second;
1908 
1909  size_t iStart = partitionBufferInterval[0].first;
1910  size_t iEnd = partitionBufferInterval[0].second;
1911  size_t jStart = partitionBufferInterval[1].first;
1912  size_t jEnd = partitionBufferInterval[1].second;
1913 
1914  for(jIndex = jStart,yIndex = yStart;jIndex <= jEnd;jIndex++,yIndex++){
1915  size_t byIndex = jIndex*bufferLine;
1916  ijk[1] = yIndex;
1917  for(iIndex = iStart,xIndex = xStart;iIndex <= iEnd;iIndex++,xIndex++){
1918  ijk[0] = xIndex;
1919  size_t bufferIndex = byIndex + iIndex;
1920  int errorCode = GridGeneratorFunction(gridGenerationOptions,gridGenerationParameters,gridSizes,ijk,xyz);
1921  if(errorCode){
1922  messageStream << "grid::GenerateGrid: User-specified GridGeneratorFunction"
1923  << " returned error code (" << errorCode << "). Aborting." << std::endl;
1924  return(1);
1925  }
1926  gridCoordinates[bufferIndex] = xyz[0];
1927  gridCoordinates[bufferIndex+numPointsBuffer] = xyz[1];
1928  }
1929  }
1930  } else if (numDim == 1) {
1931 
1932 
1933  size_t xStart = partInterval[0].first;
1934  size_t xEnd = partInterval[0].second;
1935 
1936  size_t iStart = partitionBufferInterval[0].first;
1937  size_t iEnd = partitionBufferInterval[0].second;
1938  for(iIndex = iStart,xIndex = xStart;iIndex <= iEnd;iIndex++,xIndex++){
1939  ijk[0] = xIndex;
1940  int errorCode = GridGeneratorFunction(gridGenerationOptions,gridGenerationParameters,gridSizes,ijk,xyz);
1941  if(errorCode){
1942  messageStream << "grid::GenerateGrid: User-specified GridGeneratorFunction"
1943  << " returned error code (" << errorCode << "). Aborting." << std::endl;
1944  return(1);
1945  }
1946  gridCoordinates[iIndex] = xyz[0];
1947  }
1948  }
1949  return(0);
1950  };
1951 
1952  int parallel_blockstructured::GenerateCoordinates(int (*XYZFunction)
1953  (const std::vector<size_t> &ijk,std::vector<double> &xyz),
1954  std::ostream &messageStream)
1955  {
1956 
1957  const std::vector<size_t> &gridSizes(GridSizes());
1958  const std::vector<size_t> &bufferSizes(BufferSizes());
1959  std::vector<double> &gridExtents(PhysicalExtent());
1960 
1961  pcpp::IndexIntervalType &partInterval(PartitionInterval());
1962  pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
1963 
1964  if(gridSizes.empty() || bufferSizes.empty() ||
1965  (partInterval.NNodes() == 0) || (partitionBufferInterval.NNodes() == 0)){
1966  messageStream << "grid::PBS::GenerateCoordinates: Error: Grid sizes not set up."
1967  << std::endl
1968  << "grid::PBS::GenerateCoordinates: Finalize() must be called before GenerateCoordinates."
1969  << std::endl;
1970  return(1);
1971  }
1972 
1973  int numDim = gridSizes.size();
1974 
1975  AllocateCoordinateData();
1976 
1977  std::vector<double> xyz(numDim,0);
1978  std::vector<size_t> ijk(numDim,0);
1979 
1980  size_t iIndex,jIndex,kIndex,xIndex,yIndex,zIndex;
1981 
1982  if(numDim == 3){
1983 
1984 
1985  size_t bufferPlane = bufferSizes[0]*bufferSizes[1];
1986  size_t bufferLine = bufferSizes[0];
1987 
1988  size_t xStart = partInterval[0].first;
1989  size_t xEnd = partInterval[0].second;
1990  size_t yStart = partInterval[1].first;
1991  size_t yEnd = partInterval[1].second;
1992  size_t zStart = partInterval[2].first;
1993  size_t zEnd = partInterval[2].second;
1994 
1995  size_t iStart = partitionBufferInterval[0].first;
1996  size_t iEnd = partitionBufferInterval[0].second;
1997  size_t jStart = partitionBufferInterval[1].first;
1998  size_t jEnd = partitionBufferInterval[1].second;
1999  size_t kStart = partitionBufferInterval[2].first;
2000  size_t kEnd = partitionBufferInterval[2].second;
2001 
2002  for(kIndex = kStart,zIndex = zStart;kIndex <= kEnd;kIndex++,zIndex++){
2003  size_t bzIndex = kIndex * (bufferPlane);
2004  ijk[2] = zIndex;
2005  for(jIndex = jStart,yIndex = yStart;jIndex <= jEnd;jIndex++,yIndex++){
2006  ijk[1] = yIndex;
2007  size_t byzIndex = bzIndex + jIndex*bufferLine;
2008  for(iIndex = iStart,xIndex = xStart;iIndex <= iEnd;iIndex++,xIndex++){
2009  ijk[0] = xIndex;
2010  size_t bufferIndex = byzIndex + iIndex;
2011  int errorCode = XYZFunction(ijk,xyz);
2012  if(errorCode){
2013  messageStream << "grid::GenerateCoordinates: User-specified XYZ function returned error code ("
2014  << errorCode << "). Aborting." << std::endl;
2015  return(1);
2016  }
2017  gridCoordinates[bufferIndex] = xyz[0];
2018  gridCoordinates[bufferIndex+numPointsBuffer] = xyz[1];
2019  gridCoordinates[bufferIndex+2*numPointsBuffer] = xyz[2];
2020  }
2021  }
2022  }
2023  } else if (numDim == 2) {
2024 
2025  size_t bufferLine = bufferSizes[0];
2026 
2027  size_t xStart = partInterval[0].first;
2028  size_t xEnd = partInterval[0].second;
2029  size_t yStart = partInterval[1].first;
2030  size_t yEnd = partInterval[1].second;
2031 
2032  size_t iStart = partitionBufferInterval[0].first;
2033  size_t iEnd = partitionBufferInterval[0].second;
2034  size_t jStart = partitionBufferInterval[1].first;
2035  size_t jEnd = partitionBufferInterval[1].second;
2036 
2037  for(jIndex = jStart,yIndex = yStart;jIndex <= jEnd;jIndex++,yIndex++){
2038  size_t byIndex = jIndex*bufferLine;
2039  ijk[1] = yIndex;
2040  for(iIndex = iStart,xIndex = xStart;iIndex <= iEnd;iIndex++,xIndex++){
2041  ijk[0] = xIndex;
2042  size_t bufferIndex = byIndex + iIndex;
2043  int errorCode = XYZFunction(ijk,xyz);
2044  if(errorCode){
2045  messageStream << "grid::GenerateCoordinates: User-specified XYZ function returned error code ("
2046  << errorCode << "). Aborting." << std::endl;
2047  return(1);
2048  }
2049  gridCoordinates[bufferIndex] = xyz[0];
2050  gridCoordinates[bufferIndex+numPointsBuffer] = xyz[1];
2051  }
2052  }
2053  } else if (numDim == 1) {
2054 
2055 
2056  size_t xStart = partInterval[0].first;
2057  size_t xEnd = partInterval[0].second;
2058 
2059  size_t iStart = partitionBufferInterval[0].first;
2060  size_t iEnd = partitionBufferInterval[0].second;
2061  for(iIndex = iStart,xIndex = xStart;iIndex <= iEnd;iIndex++,xIndex++){
2062  ijk[0] = xIndex;
2063  int errorCode = XYZFunction(ijk,xyz);
2064  if(errorCode){
2065  messageStream << "grid::GenerateCoordinates: User-specified XYZ function returned error code ("
2066  << errorCode << "). Aborting." << std::endl;
2067  return(1);
2068  }
2069  gridCoordinates[iIndex] = xyz[0];
2070  }
2071  }
2072  return(0);
2073  };
2074 
2075  int parallel_blockstructured::ExchangeCoordinates()
2076  {
2077  // Set up the HALO data structure for this grid/geometry
2078  simulation::grid::halo &myHalo(Halo());
2079 
2080  int numDim = Dimension();
2081 
2082  // if(xyzMessageId < 0)
2083  xyzMessageId = myHalo.CreateMessageBuffers(numDim);
2084 
2085  fixtures::CommunicatorType &gridComm(Communicator());
2086 
2087  std::vector<int> cartNeighbors;
2088  gridComm.CartNeighbors(cartNeighbors);
2089 
2090  myHalo.PackMessageBuffers(xyzMessageId,0,numDim,&gridCoordinates[0]);
2091  myHalo.SendMessage(xyzMessageId,cartNeighbors,gridComm);
2092  myHalo.ReceiveMessage(xyzMessageId,cartNeighbors,gridComm);
2093  myHalo.UnPackMessageBuffers(xyzMessageId,0,numDim,&gridCoordinates[0]);
2094 
2095  myHalo.DestroyMessageBuffer(xyzMessageId);
2096 
2097  // Correct coordinates in halo for periodicity
2098  if(!periodicLengths.empty()){
2099 
2100  // std::cout << "Correcting for periodic lengths: (";
2101  // for(int iDim = 0;iDim < numDim;iDim++){
2102  // if(iDim > 0) std::cout << ",";
2103  // std::cout << periodicLengths[iDim];
2104  // }
2105  // std::cout << ")" << std::endl;
2106 
2107  size_t numPointsBuffer = BufferSize();
2108  pcpp::IndexIntervalType bufferInterval;
2109  std::vector<size_t> &bufferSizes(BufferSizes());
2110  std::vector<size_t> &gridSizes(GridSizes());
2111  bufferInterval.InitSimple(bufferSizes);
2112  pcpp::IndexIntervalType &partitionInterval(PartitionInterval());
2113  pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
2114 
2115  for(int iDim = 0;iDim < numDim;iDim++){
2116 
2117  size_t lastIndex = gridSizes[iDim] - 1;
2118  size_t dimOffset = iDim*numPointsBuffer;
2119 
2120  if(periodicDirs[iDim] && (periodicLengths[iDim] > 0)){
2121 
2122  if(partitionInterval[iDim].first == 0){ // Correct on left domain boundary
2123 
2124  pcpp::IndexIntervalType leftHaloInterval(partitionBufferInterval);
2125  leftHaloInterval[iDim].second = partitionBufferInterval[iDim].first - 1;
2126  leftHaloInterval[iDim].first = 0;
2127  std::vector<size_t> haloIndices;
2128  bufferInterval.GetFlatIndices(leftHaloInterval,haloIndices);
2129  std::vector<size_t>::iterator haloIndexIt = haloIndices.begin();
2130  while(haloIndexIt != haloIndices.end()){
2131  size_t haloIndex = *haloIndexIt++;
2132  gridCoordinates[haloIndex+dimOffset] -= periodicLengths[iDim];
2133  }
2134 
2135  } // have left domain boundary
2136 
2137  if(partitionInterval[iDim].second == lastIndex) { // Correct on right domain boundary
2138 
2139  pcpp::IndexIntervalType rightHaloInterval(partitionBufferInterval);
2140  rightHaloInterval[iDim].first = partitionBufferInterval[iDim].second + 1;
2141  rightHaloInterval[iDim].second = bufferSizes[iDim]-1;
2142  std::vector<size_t> haloIndices;
2143  bufferInterval.GetFlatIndices(rightHaloInterval,haloIndices);
2144  std::vector<size_t>::iterator haloIndexIt = haloIndices.begin();
2145  while(haloIndexIt != haloIndices.end()){
2146  size_t haloIndex = *haloIndexIt++;
2147  gridCoordinates[haloIndex+dimOffset] += periodicLengths[iDim];
2148  }
2149 
2150  } // have right domain boundary
2151 
2152  } // periodic length specified in this dimension
2153  } // loop over dimensions
2154  } // periodic lengths not empty
2155 
2156  return(0);
2157  };
2158 
2159  int parallel_blockstructured::ExchangeJacobianMatrix(std::ostream &messageStream)
2160  {
2161 
2162  int numDim = Dimension();
2163  int numComponents = numDim*numDim;
2164  size_t numPointsBuffer = BufferSize();
2165 
2166  if(jacobianMatrix.size() != (numComponents*numPointsBuffer)){
2167  messageStream << "grid::PBS::ExchangeJacobianMatrix: Error: jacobianMatrix of improper size ("
2168  << jacobianMatrix.size() << ")." << std::endl;
2169  return(1);
2170  }
2171  return(ExchangeNodalData(numComponents,&jacobianMatrix[0]));
2172  };
2173 
2174  int parallel_blockstructured::SetupThreads(std::vector<int> &numThreadPartitions)
2175  {
2176 
2177  int numThreads = 1;
2178  int threadId = 0;
2179  bool inParallel = false;
2180 
2181 #ifdef USE_OMP
2182  inParallel = omp_in_parallel();
2183  if(inParallel){
2184  threadId = omp_get_thread_num();
2185  numThreads = omp_get_num_threads();
2186  }
2187 #endif
2188 
2189  bool masterThread = (threadId == 0);
2190 
2191 #pragma omp master
2192  {
2193  // if(masterThread){
2194  SetNumThreads(numThreads);
2195  errorState = 0;
2196  }
2197 
2198  if(inParallel){
2199 #pragma omp barrier
2200  }
2201 
2202  int numDim = Dimension();
2203  std::vector<bool> partDims(numDim,true);
2204  std::vector<int>::iterator tdIt = threadDecompDirs.begin();
2205  while(tdIt != threadDecompDirs.end()){
2206  int iDim = tdIt - threadDecompDirs.begin();
2207  if(*tdIt++ == 0)
2208  partDims[iDim] = false;
2209  }
2210  pcpp::IndexIntervalType &partitionInterval(PartitionInterval());
2211  pcpp::IndexIntervalType threadPartitionInterval;
2212 
2213 
2214  if(pcpp::util::SimplePartitionInterval(partitionInterval,partDims,threadId,numThreads,
2215  threadPartitionInterval,numThreadPartitions)){
2216  errorState = 1; // errorState is shared
2217  }
2218 
2219  if(inParallel){
2220 #pragma omp barrier
2221  }
2222 
2223  if(errorState)
2224  return(1);
2225 
2226  // Get buffer intervals for thread
2227  pcpp::IndexIntervalType &partitionBufferInterval(PartitionBufferInterval());
2228  pcpp::IndexIntervalType threadPartitionBufferInterval(threadPartitionInterval);
2229 
2230  for(int iDim = 0;iDim < numDim;iDim++){
2231  threadPartitionBufferInterval[iDim].first -= partitionInterval[iDim].first;
2232  threadPartitionBufferInterval[iDim].second -= partitionInterval[iDim].first;
2233  threadPartitionBufferInterval[iDim].first += partitionBufferInterval[iDim].first;
2234  threadPartitionBufferInterval[iDim].second += partitionBufferInterval[iDim].first;
2235  }
2236  threadPartitionBufferInterval.Sync();
2237 
2238  std::vector<size_t> &bufferSizes(BufferSizes());
2239  pcpp::IndexIntervalType threadBufferInterval(threadPartitionBufferInterval);
2240  for(int iDim = 0;iDim < numDim;iDim++){
2241  if(threadBufferInterval[iDim].first == partitionBufferInterval[iDim].first)
2242  threadBufferInterval[iDim].first = 0;
2243  if(threadBufferInterval[iDim].second == partitionBufferInterval[iDim].second)
2244  threadBufferInterval[iDim].second = bufferSizes[iDim]-1;
2245  }
2246  threadBufferInterval.Sync();
2247 
2248  partitionIntervals.threadBufferIntervals[threadId] = threadBufferInterval;
2249  partitionIntervals.threadPartitionBufferIntervals[threadId] = threadPartitionBufferInterval;
2250  partitionIntervals.threadPartitionIntervals[threadId] = threadPartitionInterval;
2251 
2252  pcpp::IndexIntervalType bufferInterval;
2253  bufferInterval.InitSimple(partitionIntervals.bufferSizes);
2254 
2255  // Subregion handling - collides this thread's partition with all
2256  // subregions and records the overlap
2257  int numSubRegions = subRegions.size();
2258 
2259  std::vector<subregion>::iterator subRegionIt = subRegions.begin();
2260  while(subRegionIt != subRegions.end()){
2261 
2262  subregion &subRegion(*subRegionIt++);
2263  pcpp::IndexIntervalType &subPartitionInterval(subRegion.threadPartitionIntervals[threadId]);
2264  if(subPartitionInterval.NNodes() == 0)
2265  continue;
2266  threadPartitionInterval.Overlap(subRegion.regionInterval,subPartitionInterval);
2267  size_t numNodesOverlap = subPartitionInterval.NNodes();
2268 
2269  if(numNodesOverlap > 0){
2270  pcpp::IndexIntervalType
2271  &subPartBufferInterval(subRegion.threadPartitionBufferIntervals[threadId]);
2272  subPartBufferInterval = subPartitionInterval;
2273 
2274  for(int iDim = 0;iDim < numDim;iDim++){
2275  subPartBufferInterval[iDim].first -= partitionInterval[iDim].first;
2276  subPartBufferInterval[iDim].second -= partitionInterval[iDim].first;
2277  subPartBufferInterval[iDim].first += partitionBufferInterval[iDim].first;
2278  subPartBufferInterval[iDim].second += partitionBufferInterval[iDim].first;
2279  }
2280  subPartBufferInterval.Sync();
2281 
2282  std::vector<size_t>
2283  &subThreadBufferIndices(subRegion.threadPartitionBufferIndices[threadId]);
2284  bufferInterval.GetFlatIndices(subPartBufferInterval,subThreadBufferIndices);
2285  }
2286  }
2287  return(0);
2288  }
2289 
2290  // Should be called *inside* parallel block
2291  void parallel_blockstructured::SetNumThreads(int numThreads)
2292  {
2293 
2294  partitionIntervals.threadPartitionIntervals.resize(numThreads);
2295  partitionIntervals.threadBufferIntervals.resize(numThreads);
2296  partitionIntervals.threadPartitionBufferIntervals.resize(numThreads);
2297  std::vector<subregion>::iterator srIt = subRegions.begin();
2298  while(srIt != subRegions.end()){
2299  subregion &subRegion(*srIt++);
2300  subRegion.threadPartitionIntervals.resize(numThreads);
2301  subRegion.threadPartitionBufferIntervals.resize(numThreads);
2302  subRegion.threadPartitionBufferIndices.resize(numThreads);
2303  }
2304  // if only one thread, just copy thread-serial data into thread 0's datastructure
2305  if(numThreads == 1){
2306  partitionIntervals.threadBufferIntervals[0].InitSimple(partitionIntervals.bufferSizes);
2307  partitionIntervals.threadPartitionBufferIntervals[0] = partitionIntervals.partitionBufferInterval;
2308  std::vector<subregion>::iterator srIt = subRegions.begin();
2309  while(srIt != subRegions.end()){
2310  subregion &subRegion(*srIt++);
2311  subRegion.threadPartitionIntervals[0] = subRegion.regionPartitionInterval;
2312  subRegion.threadPartitionBufferIntervals[0] = subRegion.regionPartitionBufferInterval;
2313  subRegion.threadPartitionBufferIndices[0] = subRegion.partitionBufferIndices;
2314  }
2315  }
2316  };
2317 
2318  int parallel_blockstructured::Finalize(bool allocateCoordinateData)
2319  {
2320 
2321  // Init & validate the grid size data structures
2322  if(!partitionIntervals.gridSizes.empty()){
2323  numDim = partitionIntervals.gridSizes.size();
2324  } else {
2325  return(1);
2326  }
2327 
2328  pcpp::IndexIntervalType globalInterval;
2329  globalInterval.InitSimple(partitionIntervals.gridSizes);
2330 
2331  if(partitionIntervals.partitionInterval.empty())
2332  partitionIntervals.partitionInterval = globalInterval;
2333  else { // ensure the partition interval is a subset of the global size
2334  if(partitionIntervals.partitionInterval.size() != numDim)
2335  return(1);
2336  for(int iDim = 0;iDim < numDim;iDim++){
2337  if(partitionIntervals.partitionInterval[iDim].second > globalInterval[iDim].second)
2338  return(1);
2339  }
2340  }
2341 
2342  partitionIntervals.partitionSizes.resize(numDim);
2343  numPointsPartition = 1;
2344  for(int iDim = 0;iDim < numDim;iDim++){
2345  partitionIntervals.partitionSizes[iDim] =
2346  partitionIntervals.partitionInterval[iDim].second -
2347  partitionIntervals.partitionInterval[iDim].first + 1;
2348  if(partitionIntervals.partitionSizes[iDim] > partitionIntervals.gridSizes[iDim])
2349  return(1);
2350  numPointsPartition *= partitionIntervals.partitionSizes[iDim];
2351  }
2352  if(dimExtensions.size() != 2*numDim){
2353  dimExtensions.resize(2*numDim,0);
2354  }
2355  partitionIntervals.bufferSizes.resize(0);
2356  partitionIntervals.bufferSizes.resize(numDim,0);
2357  numPointsBuffer = 1;
2358  partitionIntervals.partitionBufferInterval.resize(numDim);
2359  for(int iDim = 0;iDim < numDim;iDim++){
2360  partitionIntervals.bufferSizes[iDim] =
2361  partitionIntervals.partitionSizes[iDim] + dimExtensions[2*iDim] +
2362  dimExtensions[2*iDim+1];
2363  partitionIntervals.partitionBufferInterval[iDim].first = dimExtensions[2*iDim];
2364  partitionIntervals.partitionBufferInterval[iDim].second =
2365  dimExtensions[2*iDim] + partitionIntervals.partitionSizes[iDim] - 1;
2366  numPointsBuffer *= partitionIntervals.bufferSizes[iDim];
2367  }
2368  partitionIntervals.partitionBufferInterval.Sync();
2369  pcpp::IndexIntervalType bufferInterval;
2370  bufferInterval.InitSimple(partitionIntervals.bufferSizes);
2371 
2372  // Subregion handling - collides this partition with all subregions and records the overlap
2373  int numSubRegions = subRegions.size();
2374  std::vector<subregion>::iterator subRegionIt = subRegions.begin();
2375  while(subRegionIt != subRegions.end()){
2376  subregion &subRegion(*subRegionIt++);
2377  partitionIntervals.partitionInterval.Overlap(subRegion.regionInterval,subRegion.regionPartitionInterval);
2378  size_t numNodesOverlap = subRegion.regionPartitionInterval.NNodes();
2379  if(numNodesOverlap > 0){
2380  subRegion.regionPartitionBufferInterval = subRegion.regionPartitionInterval;
2381  for(int iDim = 0;iDim < numDim;iDim++){
2382  size_t nPointsRegionDim =
2383  subRegion.regionPartitionInterval[iDim].second - subRegion.regionPartitionInterval[iDim].first + 1;
2384  size_t bufferOffset = subRegion.regionPartitionInterval[iDim].first -
2385  partitionIntervals.partitionInterval[iDim].first;
2386  subRegion.regionPartitionBufferInterval[iDim].first =
2387  partitionIntervals.partitionBufferInterval[iDim].first + bufferOffset;
2388  subRegion.regionPartitionBufferInterval[iDim].second =
2389  subRegion.regionPartitionBufferInterval[iDim].first + nPointsRegionDim - 1;
2390  }
2391  subRegion.regionPartitionBufferInterval.Sync();
2392  bufferInterval.GetFlatIndices(subRegion.regionPartitionBufferInterval,subRegion.partitionBufferIndices);
2393  }
2394  }
2395 
2396  if(allocateCoordinateData)
2397  AllocateCoordinateData();
2398  // SetupThreads();
2399 
2400  int numThreads = 1;
2401  bool inParallel = false;
2402 
2403 #ifdef USE_OMP
2404  inParallel = omp_in_parallel();
2405  if(inParallel){
2406  numThreads = omp_get_num_threads();
2407  }
2408 #endif
2409 
2410  // if(threadId == 0)
2411 #pragma omp master
2412  {
2413  SetNumThreads(numThreads);
2414  }
2415 
2416  if(inParallel){
2417 #pragma omp barrier
2418  }
2419 
2420  return(0);
2421  };
2422 
2423  void halo::SetLocalBufferSizes(const std::vector<size_t> &inBufferSizes)
2424  {
2425  localBufferSizes = inBufferSizes;
2426  numPointsBuffer = 1;
2427  std::vector<size_t>::iterator bsIt = localBufferSizes.begin();
2428  while(bsIt != localBufferSizes.end())
2429  numPointsBuffer *= *bsIt++;
2430  };
2431 
2432  void halo::SetRemoteHaloExtents(const std::vector<pcpp::IndexIntervalType> haloExtents)
2433  {
2434  remoteHaloExtents = haloExtents;
2435  if(!localPartitionExtent.empty()){
2436  int numDim = haloExtents.size()/2;
2437  remoteHaloBufferExtents.resize(2*numDim);
2438  for(int iDim = 0;iDim < numDim;iDim++){
2439  pcpp::IndexIntervalType &leftHaloExtent(remoteHaloExtents[2*iDim]);
2440  if(leftHaloExtent.NNodes() > 0){ // Left remote halo for iDim direction
2441  int haloSize = leftHaloExtent[iDim].second - leftHaloExtent[iDim].first + 1;
2442  pcpp::IndexIntervalType remoteHaloBufferExtent(localPartitionExtent);
2443  remoteHaloBufferExtent[iDim].first = 0;
2444  remoteHaloBufferExtent[iDim].second = haloSize - 1;
2445  remoteHaloBufferExtents[iDim*2].Copy(remoteHaloBufferExtent);
2446  }
2447  pcpp::IndexIntervalType &rightHaloExtent(remoteHaloExtents[2*iDim+1]);
2448  if(rightHaloExtent.NNodes() > 0){ // Right remote halo for iDim direction
2449  int haloSize = rightHaloExtent[iDim].second - rightHaloExtent[iDim].first + 1;
2450  pcpp::IndexIntervalType remoteHaloBufferExtent(localPartitionExtent);
2451  remoteHaloBufferExtent[iDim].second = localBufferSizes[iDim] - 1;
2452  remoteHaloBufferExtent[iDim].first = remoteHaloBufferExtent[iDim].second - haloSize + 1;
2453  remoteHaloBufferExtents[iDim*2+1].Copy(remoteHaloBufferExtent);
2454  }
2455  }
2456  }
2457  // if(!localPartitionExtent.empty()){
2458  // }
2459  };
2460 
2461  void halo::SetThreadExtent(int myThreadId,pcpp::IndexIntervalType threadInterval)
2462  {
2463  threadExtents[myThreadId] = threadInterval;
2464  if(!localPartitionExtent.empty()){
2465  int numDim = threadInterval.size();
2466  pcpp::IndexIntervalType &threadBufferExtent(threadBufferExtents[myThreadId]);
2467  threadBufferExtent = (threadInterval - globalPartitionExtent);
2468  threadBufferExtent += localPartitionExtent;
2469  for(int iDim = 0;iDim < numDim;iDim++){
2470  if(threadBufferExtent[iDim].first == localPartitionExtent[iDim].first){
2471  threadBufferExtent[iDim].first = 0;
2472  }
2473  if(threadBufferExtent[iDim].second == localPartitionExtent[iDim].second)
2474  threadBufferExtent[iDim].second = localBufferSizes[iDim] - 1;
2475  }
2476  }
2477  };
2478 
2479  void halo::SetNumThreads(int numThreadsIn)
2480  {
2481  numThreads = numThreadsIn;
2482  threadExtents.resize(numThreads,globalPartitionExtent);
2483  threadBufferExtents.resize(numThreads,localPartitionExtent);
2484  threadSendIndices.resize(numThreads);
2485  threadSendBufferIndices.resize(numThreads);
2486  threadRecvIndices.resize(numThreads);
2487  threadRecvBufferIndices.resize(numThreads);
2488  threadHaloBufferIndices.resize(numThreads);
2489  };
2490 
2491 
2492  void halo::SetLocalHaloExtents(const std::vector<pcpp::IndexIntervalType> haloExtents)
2493  {
2494  localHaloExtents = haloExtents;
2495  if(!localPartitionExtent.empty()){
2496  //assert(!bufferSizes.empty());
2497  int numDim = haloExtents.size();
2498  localHaloBufferExtents.resize(numDim);
2499  for(int iDim = 0;iDim < numDim;iDim++){
2500  pcpp::IndexIntervalType &threadSendInterval(localHaloExtents[iDim]);
2501  // Adjust to local numbering: extentWRTLocal = extentWRTGlobal - globalPartitionExtent
2502  pcpp::IndexIntervalType subThreadSendExtent(threadSendInterval - globalPartitionExtent);
2503  // Adjust to the buffer numbering: extentWRTLocalBuffer = extentWRTLocal + localPartitionInterval
2504  subThreadSendExtent += localPartitionExtent;
2505  // Obtain flat indices: bufferExtent.GetFlatIndices(subHaloExtent,sendIndices[i]);
2506  localHaloBufferExtents[iDim].Copy(subThreadSendExtent);
2507  }
2508  }
2509  };
2510 
2511  // Simple halo extent creation - takes a size in each (+/-) direction (i.e. 2*NumDim sizes)
2512  // and creates the extents that specify the remote halo regions
2513  std::vector<pcpp::IndexIntervalType> halo::CreateRemoteHaloExtents(const pcpp::IndexIntervalType &globalExtent,
2514  const pcpp::IndexIntervalType &partitionExtent,
2515  std::vector<size_t> &haloSizes)
2516  {
2517  std::vector<pcpp::IndexIntervalType> haloExtents;
2518  int numDim = haloSizes.size()/2;
2519  haloExtents.resize(numDim*2);
2520  for(int iDim = 0;iDim < numDim;iDim++){
2521  size_t numNodesDim = partitionExtent[iDim].second - partitionExtent[iDim].first + 1;
2522  bool leftHaloRequired = false;
2523  if(haveNeighbor.empty()){
2524  if(partitionExtent[iDim].first > globalExtent[iDim].first)
2525  leftHaloRequired = true;
2526  } else if (haveNeighbor[2*iDim]) {
2527  leftHaloRequired = true;
2528  }
2529  if(leftHaloRequired){
2530  size_t haloSize = haloSizes[iDim*2];
2531  pcpp::IndexIntervalType leftHaloExtent(partitionExtent);
2532  // if((partitionExtent[iDim].first - globalExtent[iDim].first) < haloSize)
2533  // haloSize = partitionExtent[iDim].first - globalExtent[iDim].first;
2534  assert(haloSize < numNodesDim);
2535  for(int hDim = 0;hDim < numDim;hDim++){
2536  if(hDim == iDim){
2537  if(partitionExtent[iDim].first == 0){
2538  leftHaloExtent[hDim].first = globalExtent[iDim].second - haloSize + 1;
2539  leftHaloExtent[hDim].second = globalExtent[iDim].second;
2540  } else {
2541  leftHaloExtent[hDim].first = partitionExtent[iDim].first - haloSize;
2542  leftHaloExtent[hDim].second = leftHaloExtent[hDim].first + haloSize - 1;
2543  }
2544  } else {
2545  leftHaloExtent[hDim] = partitionExtent[hDim];
2546  }
2547  }
2548  haloExtents[iDim*2].Copy(leftHaloExtent);
2549  }
2550 
2551  bool rightHaloRequired = false;
2552  if(haveNeighbor.empty()){
2553  if(partitionExtent[iDim].second < globalExtent[iDim].second)
2554  rightHaloRequired = true;
2555  } else if (haveNeighbor[2*iDim+1]) {
2556  rightHaloRequired = true;
2557  }
2558 
2559  if(rightHaloRequired){ // right halo will be required
2560  size_t haloSize = haloSizes[iDim*2+1];
2561  pcpp::IndexIntervalType rightHaloExtent(partitionExtent);
2562  // if((globalExtent[iDim].second - partitionExtent[iDim].second) < haloSize)
2563  // haloSize = globalExtent[iDim].second - partitionExtent[iDim].second;
2564  assert(haloSize < numNodesDim);
2565  for(int hDim = 0;hDim < numDim;hDim++){
2566  if(hDim == iDim){
2567  if(partitionExtent[iDim].second == globalExtent[iDim].second){
2568  rightHaloExtent[hDim].second = globalExtent[iDim].first + haloSize - 1;
2569  rightHaloExtent[hDim].first = globalExtent[iDim].first;
2570  } else {
2571  rightHaloExtent[hDim].first = partitionExtent[iDim].second + 1;
2572  rightHaloExtent[hDim].second = rightHaloExtent[hDim].first + haloSize - 1;
2573  }
2574  } else {
2575  rightHaloExtent[hDim] = partitionExtent[hDim];
2576  }
2577  }
2578  haloExtents[iDim*2+1].Copy(rightHaloExtent);
2579  }
2580  }
2581  return(haloExtents);
2582  };
2583 
2584  std::vector<pcpp::IndexIntervalType> halo::CreateRemoteHaloExtents(std::vector<size_t> &haloSizes)
2585  {
2586  std::vector<pcpp::IndexIntervalType> haloExtents(CreateRemoteHaloExtents(globalGridExtent,globalPartitionExtent,haloSizes));
2587  return(haloExtents);
2588  }
2589 
2590  // Simple halo extent creation - takes a size in each (+/-) direction (i.e. 2*NumDim sizes)
2591  // and creates the extents that specify the local halos (i.e. the local data upon which remote procs depend)
2592  std::vector<pcpp::IndexIntervalType> halo::CreateLocalHaloExtents(const pcpp::IndexIntervalType &globalExtent,
2593  const pcpp::IndexIntervalType &partitionExtent,
2594  std::vector<size_t> &haloSizes)
2595  {
2596  std::vector<pcpp::IndexIntervalType> haloExtents;
2597  int numDim = haloSizes.size()/2;
2598  haloExtents.resize(numDim*2);
2599  for(int iDim = 0;iDim < numDim;iDim++){
2600  size_t numNodesDim = partitionExtent[iDim].second - partitionExtent[iDim].first + 1;
2601  bool leftHaloRequired = false;
2602  if(haveNeighbor.empty()){
2603  if(partitionExtent[iDim].first > globalExtent[iDim].first) // left halo will be required
2604  leftHaloRequired = true;
2605  } else if (haveNeighbor[iDim*2]) {
2606  leftHaloRequired = true;
2607  }
2608  if(leftHaloRequired){
2609  size_t haloSize = haloSizes[iDim*2];
2610  pcpp::IndexIntervalType leftHaloExtent(partitionExtent);
2611  // if((partitionExtent[iDim].first - globalExtent[iDim].first) < haloSize)
2612  // haloSize = partitionExtent[iDim].first - globalExtent[iDim].first;
2613  assert(haloSize < numNodesDim);
2614  for(int hDim = 0;hDim < numDim;hDim++){
2615  if(hDim == iDim){
2616  leftHaloExtent[hDim].first = partitionExtent[iDim].first;
2617  leftHaloExtent[hDim].second = leftHaloExtent[hDim].first + haloSize - 1;
2618  } else {
2619  leftHaloExtent[hDim] = partitionExtent[hDim];
2620  }
2621  }
2622  haloExtents[iDim*2].Copy(leftHaloExtent);
2623  }
2624  bool rightHaloRequired = false;
2625  if(haveNeighbor.empty()){
2626  if(partitionExtent[iDim].second < globalExtent[iDim].second) // right halo will be required
2627  rightHaloRequired = true;
2628  } else if(haveNeighbor[2*iDim+1]) {
2629  rightHaloRequired = true;
2630  }
2631  if(rightHaloRequired){
2632  size_t haloSize = haloSizes[iDim*2+1];
2633  pcpp::IndexIntervalType rightHaloExtent(partitionExtent);
2634  // if((globalExtent[iDim].second - partitionExtent[iDim].second) < haloSize)
2635  // haloSize = globalExtent[iDim].second - partitionExtent[iDim].second;
2636  assert(haloSize < numNodesDim);
2637  for(int hDim = 0;hDim < numDim;hDim++){
2638  if(hDim == iDim){
2639  rightHaloExtent[hDim].first = partitionExtent[iDim].second - haloSize + 1;
2640  rightHaloExtent[hDim].second = partitionExtent[hDim].second;
2641  } else {
2642  rightHaloExtent[hDim] = partitionExtent[hDim];
2643  }
2644  }
2645  haloExtents[iDim*2+1].Copy(rightHaloExtent);
2646  }
2647  }
2648  return(haloExtents);
2649  };
2650 
2651  std::vector<pcpp::IndexIntervalType> halo::CreateLocalHaloExtents(std::vector<size_t> &haloSizes)
2652  {
2653  std::vector<pcpp::IndexIntervalType> haloExtents(CreateLocalHaloExtents(globalGridExtent,globalPartitionExtent,haloSizes));
2654  return(haloExtents);
2655  }
2656 
2657 
2659  int halo::CreateSimpleSendIndices()
2660  {
2661  haveSendData = false;
2662  size_t totalNumSend = 0;
2663  if(localHaloExtents.empty())
2664  return(0);
2665 
2666  int numSendIntervals = localHaloExtents.size();
2667  sendIndices.resize(numSendIntervals);
2668  for(int i = 0;i < numSendIntervals;i++){
2669  size_t numPointsSend = localHaloExtents[i].NNodes();
2670  if(numPointsSend > 0){
2671 
2672  // Local halo extents are wrt the global grid size. We need the flat array of source
2673  // buffer indices (i.e. wrt the local buffer size) that will be used to populate
2674  // the send buffers. If the local buffer includes "halos" or other buffer zones,
2675  // then it is a bit more complicated to obtain the source buffer indices.
2676  // The following code determines the source buffer indices.
2677  if(localPartitionExtent.empty()){ // then halos are not integrated, globalPartitionExtent
2678  // represents the entire local buffer
2679  globalPartitionExtent.GetFlatIndices(localHaloExtents[i],sendIndices[i]);
2680 
2681  } else { // the global partition interval is a sub-interval of the local buffer, deal with it
2682  //assert(!bufferSizes.empty());
2683  // Adjust to local numbering: subHaloExtent = localHaloExtents - globalPartitionExtent
2684  pcpp::IndexIntervalType subHaloExtent(localHaloExtents[i] - globalPartitionExtent);
2685  // Adjust to the buffer numbering: subHaloExtent = subHaloExtent + localPartitionInterval
2686  subHaloExtent += localPartitionExtent;
2687  // Obtain flat indices: bufferExtent.GetFlatIndices(subHaloExtent,sendIndices[i]);
2688  pcpp::IndexIntervalType bufferInterval;
2689  bufferInterval.InitSimple(localBufferSizes);
2690  bufferInterval.GetFlatIndices(subHaloExtent,sendIndices[i]);
2691  }
2692 
2693  totalNumSend += numPointsSend;
2694  }
2695  }
2696  haveSendData = (totalNumSend > 0);
2697  return(totalNumSend);
2698  }
2699 
2700 
2702  int halo::CreateThreadSendIndices(int threadId)
2703  {
2704  if(!haveSendData)
2705  return(0);
2706  if(localHaloExtents.empty())
2707  return(0);
2708  int numSendBuffers = localHaloExtents.size();
2709  threadSendIndices[threadId].resize(numSendBuffers);
2710  threadSendBufferIndices[threadId].resize(numSendBuffers);
2711  pcpp::IndexIntervalType &threadInterval(threadExtents[threadId]);
2712  for(int i = 0;i < numSendBuffers;i++){
2713  size_t numSendNodes = localHaloExtents[i].NNodes();
2714  std::vector<size_t> &threadSendInds(threadSendIndices[threadId][i]);
2715  std::vector<size_t> &threadSendBufferInds(threadSendBufferIndices[threadId][i]);
2716  if(numSendNodes > 0){
2717 
2718  //
2719  // - debugging prints
2720  // std::cout << "Num nodes in Halo = " << numSendNodes << std::endl;
2721  // std::cout << "Halo Extent: ";
2722  // localHaloExtents[i].PrettyPrint(std::cout);
2723  // std::cout << std::endl;
2724  // std::cout << "Thread Extent: ";
2725  // threadInterval.PrettyPrint(std::cout);
2726  // std::cout << std::endl;
2727  // - debugging
2728  //
2729 
2730  // collide thread Interval with the halo interval; gets the portion of the thread interval
2731  // that needs to be sent for this particular border
2732  pcpp::IndexIntervalType threadSendInterval(threadInterval.Overlap(localHaloExtents[i]));
2733  size_t numThreadSendNodes = threadSendInterval.NNodes();
2734  if(numThreadSendNodes > 0){
2735  //
2736  // debugging
2737  // std::cout << "Thread Send Extent: ";
2738  // threadSendInterval.PrettyPrint(std::cout);
2739  // std::cout << std::endl;
2740  // std::cout << "Thread collides with " << numThreadSendNodes << " nodes to send." << std::endl;
2741  // debugging
2742  //
2743 
2744  // Get the linear thread-specific source indices; these are the indices of the partition
2745  // buffer from which to source send data for this thread
2746  if(localPartitionExtent.empty()){ // then partition extent is not a sub-interval of the buffer
2747  globalPartitionExtent.GetFlatIndices(threadSendInterval,threadSendInds);
2748  } else { // partition is a subinterval of the local buffer, deal with it
2749  //assert(!bufferSizes.empty());
2750  // Adjust to local numbering: extentWRTLocal = extentWRTGlobal - globalPartitionExtent
2751  pcpp::IndexIntervalType subThreadSendExtent(threadSendInterval - globalPartitionExtent);
2752  // Adjust to the buffer numbering: extentWRTLocalBuffer = extentWRTLocal + localPartitionInterval
2753  subThreadSendExtent += localPartitionExtent;
2754  // Obtain flat indices: bufferExtent.GetFlatIndices(subHaloExtent,sendIndices[i]);
2755  pcpp::IndexIntervalType bufferInterval;
2756  bufferInterval.InitSimple(localBufferSizes);
2757  bufferInterval.GetFlatIndices(subThreadSendExtent,threadSendInds);
2758  }
2759 
2760  // Get the linear thread-specific target/pack indices; these are the indices of the send buffer
2761  // to which to source data is copied for this thread
2762  localHaloExtents[i].GetFlatIndices(threadSendInterval,threadSendBufferInds);
2763  if(threadSendInds.size() != numThreadSendNodes){
2764  exit(1);
2765  }
2766  }
2767  }
2768  }
2769  return(0);
2770  }
2771 
2772 
2773 
2774 
2775  int halo::CreateSimpleRecvIndices()
2776  {
2777  haveRecvData = false;
2778  if(remoteHaloExtents.empty()){
2779  std::cout << "halo::CreateSimpleRecvBuffers: WARNING No remote halos." << std::endl;
2780  return(0);
2781  }
2782  // DestroyRecvBuffers();
2783  size_t totalNumRecv = 0;
2784  int numRecvIntervals = remoteHaloExtents.size();
2785  // recvBuffers.resize(numRecvBuffers,NULL);
2786  recvIndices.resize(numRecvIntervals);
2787  for(int i = 0;i < numRecvIntervals;i++){
2788  std::vector<size_t> &recvIndex(recvIndices[i]);
2789  size_t numRemoteHaloNodes = remoteHaloExtents[i].NNodes();
2790  if(numRemoteHaloNodes > 0){
2791  if(localPartitionExtent.empty()){
2792  remoteHaloExtents[i].GetFlatIndices(remoteHaloExtents[i],recvIndex);
2793  } else {
2794  pcpp::IndexIntervalType localBufferExtent;
2795  localBufferExtent.InitSimple(localBufferSizes);
2796  localBufferExtent.GetFlatIndices(remoteHaloBufferExtents[i],recvIndex);
2797  }
2798  assert(recvIndex.size() == remoteHaloExtents[i].NNodes());
2799 
2800  totalNumRecv += numRemoteHaloNodes;
2801  }
2802  }
2803  haveRecvData = (totalNumRecv > 0);
2804  return(totalNumRecv);
2805  };
2806 
2807  // Destroy only works on the last buffer created, can be called
2808  // in sequence to delete all message buffers.
2809  int halo::DestroyMessageBuffer(int messageId)
2810  {
2811  if(messageId < 0) {
2812  messageId = communicationBuffers.size() - 1;
2813  } else if(messageId != communicationBuffers.size() - 1){
2814  return(1);
2815  }
2816 
2817  halo::messagebuffers &lastMessage(communicationBuffers[messageId]);
2818  int numCommunicationBorders = recvIndices.size();
2819 
2820  for(int iComm = 0;iComm < numCommunicationBorders;iComm++){
2821  std::vector<double *>::iterator mbIt = lastMessage.sendBuffers.begin();
2822  while(mbIt != lastMessage.sendBuffers.end()){
2823  delete [] *mbIt;
2824  *mbIt = NULL;
2825  mbIt++;
2826  }
2827  mbIt = lastMessage.recvBuffers.begin();
2828  while(mbIt != lastMessage.recvBuffers.end()){
2829  delete [] *mbIt;
2830  *mbIt = NULL;
2831  mbIt++;
2832  }
2833  std::vector<int *>::iterator fmbIt = lastMessage.sendFlagBuffers.begin();
2834  while(fmbIt != lastMessage.sendFlagBuffers.end()){
2835  delete [] *fmbIt;
2836  *fmbIt = NULL;
2837  fmbIt++;
2838  }
2839  fmbIt = lastMessage.recvFlagBuffers.begin();
2840  while(fmbIt != lastMessage.recvFlagBuffers.end()){
2841  delete [] *fmbIt;
2842  *fmbIt = NULL;
2843  fmbIt++;
2844  }
2845  }
2846  communicationBuffers.pop_back();
2847  return(0);
2848  };
2849 
2850  int halo::CreateMessageBuffers(int numComponents,int componentSize)
2851  {
2852 
2853  int numCommuncationBorders = recvIndices.size();
2854  if(numCommuncationBorders != sendIndices.size())
2855  return(-1);
2856  int messageId = communicationBuffers.size();
2857  halo::messagebuffers newMessageBuffers;
2858  newMessageBuffers.numComponents = numComponents;
2859  newMessageBuffers.componentSize = componentSize;
2860  if(componentSize == 8) {
2861  newMessageBuffers.sendBuffers.resize(numCommuncationBorders);
2862  newMessageBuffers.recvBuffers.resize(numCommuncationBorders);
2863  } else {
2864  newMessageBuffers.sendFlagBuffers.resize(numCommuncationBorders);
2865  newMessageBuffers.recvFlagBuffers.resize(numCommuncationBorders);
2866  }
2867  newMessageBuffers.numPointsSend.resize(numCommuncationBorders);
2868  newMessageBuffers.numPointsRecv.resize(numCommuncationBorders);
2869 
2870  for(int iComm = 0;iComm < numCommuncationBorders;iComm++){
2871  int numRecv = recvIndices[iComm].size();
2872  int numSend = sendIndices[iComm].size();
2873  newMessageBuffers.numPointsSend[iComm] = numSend;
2874  newMessageBuffers.numPointsRecv[iComm] = numRecv;
2875  if(componentSize == 8){
2876  newMessageBuffers.sendBuffers[iComm] = new double [numComponents * numSend];
2877  newMessageBuffers.recvBuffers[iComm] = new double [numComponents * numRecv];
2878  } else {
2879  newMessageBuffers.sendFlagBuffers[iComm] = new int [numComponents * numSend];
2880  newMessageBuffers.recvFlagBuffers[iComm] = new int [numComponents * numRecv];
2881  }
2882  }
2883  communicationBuffers.push_back(newMessageBuffers);
2884  return(messageId);
2885  };
2886 
2887  void halo::DestroyAll()
2888  {
2889  std::vector<messagebuffers>::iterator msgBufIt =
2890  communicationBuffers.begin();
2891  while(msgBufIt != communicationBuffers.end()){
2892  messagebuffers &msgBuffer(*msgBufIt++);
2893  std::vector<double *>::iterator bufIt = msgBuffer.sendBuffers.begin();
2894  while(bufIt != msgBuffer.sendBuffers.end()){
2895  if(*bufIt != NULL){
2896  delete [] *bufIt;
2897  *bufIt = NULL;
2898  }
2899  bufIt++;
2900  }
2901  bufIt = msgBuffer.recvBuffers.begin();
2902  while(bufIt != msgBuffer.recvBuffers.end()){
2903  if(*bufIt != NULL){
2904  delete [] *bufIt;
2905  *bufIt = NULL;
2906  }
2907  bufIt++;
2908  }
2909  std::vector<int *>::iterator flagBufIt = msgBuffer.sendFlagBuffers.begin();
2910  while(flagBufIt != msgBuffer.sendFlagBuffers.end()){
2911  if(*flagBufIt != NULL){
2912  delete [] *flagBufIt;
2913  *flagBufIt = NULL;
2914  }
2915  flagBufIt++;
2916  }
2917  flagBufIt = msgBuffer.recvFlagBuffers.begin();
2918  while(flagBufIt != msgBuffer.recvFlagBuffers.end()){
2919  if(*flagBufIt != NULL){
2920  delete [] *flagBufIt;
2921  *flagBufIt = NULL;
2922  }
2923  flagBufIt++;
2924  }
2925  }
2926  };
2927 
2937  int halo::PackMessageBuffers(const int messageId,const int componentId,
2938  const int numComponents,const double *sourceBuffer)
2939  {
2940  assert(messageId < communicationBuffers.size());
2941  halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
2942  assert((componentId+numComponents) <= messageBuffers.numComponents);
2943  int numMessages = sendIndices.size();
2944  for(int iMessage = 0;iMessage < numMessages;iMessage++){
2945  std::vector<size_t> &sendIndex(sendIndices[iMessage]);
2946  size_t numSend = sendIndex.size();
2947  for(int componentIndex = 0;componentIndex < numComponents;componentIndex++){
2948  for(size_t pointIndex = 0;pointIndex < numSend;pointIndex++){
2949  size_t sendBufferIndex = pointIndex + (componentId+componentIndex)*numSend;
2950  size_t sourceBufferIndex = componentIndex*numPointsBuffer;
2951  messageBuffers.sendBuffers[iMessage][sendBufferIndex] =
2952  sourceBuffer[sourceBufferIndex+sendIndex[pointIndex]];
2953  }
2954  }
2955  }
2956  return(0);
2957  };
2958 
2959 
2969  int halo::PackMessageBuffers(const int messageId,const int componentId,
2970  const int numComponents,const int *sourceBuffer)
2971  {
2972  assert(messageId < communicationBuffers.size());
2973  halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
2974  assert((componentId+numComponents) <= messageBuffers.numComponents);
2975  assert(messageBuffers.componentSize == 4);
2976  int numMessages = sendIndices.size();
2977  for(int iMessage = 0;iMessage < numMessages;iMessage++){
2978  std::vector<size_t> &sendIndex(sendIndices[iMessage]);
2979  size_t numSend = sendIndex.size();
2980  for(int componentIndex = 0;componentIndex < numComponents;componentIndex++){
2981  for(size_t pointIndex = 0;pointIndex < numSend;pointIndex++){
2982  size_t sendBufferIndex = pointIndex + (componentId+componentIndex)*numSend;
2983  size_t sourceBufferIndex = componentIndex*numPointsBuffer;
2984  messageBuffers.sendFlagBuffers[iMessage][sendBufferIndex] =
2985  sourceBuffer[sourceBufferIndex+sendIndex[pointIndex]];
2986  }
2987  }
2988  }
2989  return(0);
2990  };
2991 
3001  int halo::PackMessageBuffers(const int messageId,const int componentId,
3002  const int numComponents,const double *sourceBuffer,
3003  const int threadId)
3004  {
3005  assert(messageId < communicationBuffers.size());
3006  halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
3007  assert((componentId+numComponents) <= messageBuffers.numComponents);
3008  int numMessages = sendIndices.size();
3009  for(int iMessage = 0;iMessage < numMessages;iMessage++){
3010  std::vector<size_t> &threadSendIndex(threadSendIndices[threadId][iMessage]);
3011  std::vector<size_t> &threadSendBufferIndex(threadSendBufferIndices[threadId][iMessage]);
3012  std::vector<size_t> &sendIndex(sendIndices[iMessage]);
3013  size_t numThreadSend = threadSendIndex.size();
3014  size_t numSend = sendIndex.size();
3015  for(int componentIndex = 0;componentIndex < numComponents;componentIndex++){
3016  for(size_t pointIndex = 0;pointIndex < numThreadSend;pointIndex++){
3017  size_t sendBufferIndex = threadSendBufferIndex[pointIndex] + (componentId+componentIndex)*numSend;
3018  size_t sourceBufferIndex = componentIndex*numPointsBuffer + threadSendIndex[pointIndex];
3019  messageBuffers.sendBuffers[iMessage][sendBufferIndex] = sourceBuffer[sourceBufferIndex];
3020  }
3021  }
3022  }
3023  return(0);
3024  };
3025 
3035  int halo::PackMessageBuffers(const int messageId,const int componentId,
3036  const int numComponents,const int *sourceBuffer,
3037  const int threadId)
3038  {
3039  assert(messageId < communicationBuffers.size());
3040  halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
3041  assert((componentId+numComponents) <= messageBuffers.numComponents);
3042  assert(messageBuffers.componentSize == 4);
3043  int numMessages = sendIndices.size();
3044  for(int iMessage = 0;iMessage < numMessages;iMessage++){
3045  std::vector<size_t> &threadSendIndex(threadSendIndices[threadId][iMessage]);
3046  std::vector<size_t> &threadSendBufferIndex(threadSendBufferIndices[threadId][iMessage]);
3047  std::vector<size_t> &sendIndex(sendIndices[iMessage]);
3048  size_t numThreadSend = threadSendIndex.size();
3049  size_t numSend = sendIndex.size();
3050  for(int componentIndex = 0;componentIndex < numComponents;componentIndex++){
3051  for(size_t pointIndex = 0;pointIndex < numThreadSend;pointIndex++){
3052  size_t sendBufferIndex = threadSendBufferIndex[pointIndex] + (componentId+componentIndex)*numSend;
3053  size_t sourceBufferIndex = componentIndex*numPointsBuffer + threadSendIndex[pointIndex];
3054  messageBuffers.sendFlagBuffers[iMessage][sendBufferIndex] = sourceBuffer[sourceBufferIndex];
3055  }
3056  }
3057  }
3058  return(0);
3059  };
3060 
3061  int halo::SendMessage(const int messageId,const std::vector<int> &neighborRanks,
3062  fixtures::CommunicatorType &inComm)
3063  {
3064  assert(messageId < communicationBuffers.size());
3065  halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
3066  int componentSize = messageBuffers.componentSize;
3067  int numMessages = neighborRanks.size();
3068  int numActual = 0;
3069  int numComponents = messageBuffers.numComponents;
3070  if(componentSize == 8){
3071  assert(numMessages == messageBuffers.sendBuffers.size());
3072  } else {
3073  assert(numMessages == messageBuffers.sendFlagBuffers.size());
3074  }
3075  int messageTagBase = messageId*100;
3076  for(int iMessage = 0;iMessage < numMessages;iMessage++){
3077  int remoteRank = neighborRanks[iMessage];
3078  if(remoteRank >= 0){
3079  int numPointsSend = messageBuffers.numPointsSend[iMessage];
3080  int numSendItems = numPointsSend*numComponents;
3081  int messageTag = messageTagBase + iMessage;
3082  // if(remoteRank >= 0){
3083  numActual++;
3084  if(componentSize == 8) {
3085  double *sendBuffer = messageBuffers.sendBuffers[iMessage];
3086  inComm.ASendBuf(sendBuffer,numSendItems,remoteRank,messageTag);
3087  } else {
3088  int *sendBuffer = messageBuffers.sendFlagBuffers[iMessage];
3089  inComm.ASendBuf(sendBuffer,numSendItems,remoteRank,messageTag);
3090  }
3091  }
3092  }
3093  return(0);
3094  };
3095 
3096  int halo::ReceiveMessage(const int messageId,const std::vector<int> &neighborRanks,
3097  fixtures::CommunicatorType &inComm)
3098  {
3099  assert(messageId < communicationBuffers.size());
3100  halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
3101  int componentSize = messageBuffers.componentSize;
3102  int numMessages = neighborRanks.size();
3103  int numActual = 0;
3104  int numComponents = messageBuffers.numComponents;
3105  if(componentSize == 8)
3106  assert(numMessages == messageBuffers.recvBuffers.size());
3107  else
3108  assert(numMessages == messageBuffers.recvFlagBuffers.size());
3109  int messageTagBase = messageId*100;
3110  for(int iMessage = 0;iMessage < numMessages;iMessage++){
3111  int remoteRank = neighborRanks[iMessage];
3112  if(remoteRank >= 0){
3113  numActual++;
3114  int iDir = iMessage%2;
3115  int numPointsRecv = messageBuffers.numPointsSend[iMessage];
3116  int numRecvItems = numPointsRecv*numComponents;
3117  int messageTag = messageTagBase + iMessage;
3118  if(iDir == 0)
3119  messageTag++;
3120  else
3121  messageTag--;
3122  if(componentSize == 8){
3123  double *recvBuffer = messageBuffers.recvBuffers[iMessage];
3124  inComm.ARecvBuf(recvBuffer,numRecvItems,remoteRank,messageTag);
3125  } else {
3126  int *recvBuffer = messageBuffers.recvFlagBuffers[iMessage];
3127  inComm.ARecvBuf(recvBuffer,numRecvItems,remoteRank,messageTag);
3128  }
3129  }
3130  }
3131  if(numActual > 0)
3132  inComm.WaitAll();
3133  return(0);
3134  };
3135 
3145  int halo::UnPackMessageBuffers(const int messageId,const int componentId,
3146  const int numComponents,double *targetBuffer)
3147  {
3148  assert(messageId < communicationBuffers.size());
3149  halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
3150  assert((componentId+numComponents) <= messageBuffers.numComponents);
3151  int numMessages = recvIndices.size();
3152  for(int iMessage = 0;iMessage < numMessages;iMessage++){
3153  std::vector<size_t> &recvIndex(recvIndices[iMessage]);
3154  size_t numRecv = recvIndex.size();
3155  for(int componentIndex = 0;componentIndex < numComponents;componentIndex++){
3156  for(size_t pointIndex = 0;pointIndex < numRecv;pointIndex++){
3157  size_t recvBufferIndex = pointIndex + (componentId+componentIndex)*numRecv;
3158  size_t targetBufferIndex = componentIndex*numPointsBuffer + recvIndex[pointIndex];
3159  targetBuffer[targetBufferIndex] = messageBuffers.recvBuffers[iMessage][recvBufferIndex];
3160  }
3161  }
3162  }
3163  return(0);
3164  };
3165 
3175  int halo::UnPackMessageBuffers(const int messageId,const int componentId,
3176  const int numComponents,int *targetBuffer)
3177  {
3178  assert(messageId < communicationBuffers.size());
3179  halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
3180  assert((componentId+numComponents) <= messageBuffers.numComponents);
3181  assert(messageBuffers.componentSize == 4);
3182  int numMessages = recvIndices.size();
3183  for(int iMessage = 0;iMessage < numMessages;iMessage++){
3184  std::vector<size_t> &recvIndex(recvIndices[iMessage]);
3185  size_t numRecv = recvIndex.size();
3186  for(int componentIndex = 0;componentIndex < numComponents;componentIndex++){
3187  for(size_t pointIndex = 0;pointIndex < numRecv;pointIndex++){
3188  size_t recvBufferIndex = pointIndex + (componentId+componentIndex)*numRecv;
3189  size_t targetBufferIndex = componentIndex*numPointsBuffer + recvIndex[pointIndex];
3190  targetBuffer[targetBufferIndex] = messageBuffers.recvFlagBuffers[iMessage][recvBufferIndex];
3191  }
3192  }
3193  }
3194  return(0);
3195  };
3196 
3197  int halo::UnPackMessageBuffers(const int messageId,const int componentId,
3198  const int numComponents,double *targetBuffer,
3199  const int threadId)
3200  {
3201  assert(messageId < communicationBuffers.size());
3202  halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
3203  assert((componentId+numComponents) <= messageBuffers.numComponents);
3204  int numMessages = recvIndices.size();
3205  for(int iMessage = 0;iMessage < numMessages;iMessage++){
3206  std::vector<size_t> &threadRecvIndex(threadHaloBufferIndices[threadId][iMessage]);
3207  std::vector<size_t> &threadRecvBufferIndex(threadRecvBufferIndices[threadId][iMessage]);
3208  size_t recvBufferSize = recvIndices[iMessage].size();
3209  size_t numThreadRecv = threadRecvIndex.size();
3210  for(int componentIndex = 0;componentIndex < numComponents;componentIndex++){
3211  for(size_t pointIndex = 0;pointIndex < numThreadRecv;pointIndex++){
3212  size_t recvBufferIndex = threadRecvBufferIndex[pointIndex] +
3213  (componentId+componentIndex)*recvBufferSize;
3214  size_t targetBufferIndex = componentIndex*numPointsBuffer + threadRecvIndex[pointIndex];
3215  targetBuffer[targetBufferIndex] = messageBuffers.recvBuffers[iMessage][recvBufferIndex];
3216  }
3217  }
3218  }
3219  return(0);
3220  };
3221 
3222 
3223  int halo::UnPackMessageBuffers(const int messageId,const int componentId,
3224  const int numComponents,int *targetBuffer,
3225  const int threadId)
3226  {
3227  assert(messageId < communicationBuffers.size());
3228  halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
3229  assert((componentId+numComponents) <= messageBuffers.numComponents);
3230  assert(messageBuffers.componentSize == 4);
3231  int numMessages = recvIndices.size();
3232  for(int iMessage = 0;iMessage < numMessages;iMessage++){
3233  std::vector<size_t> &threadRecvIndex(threadHaloBufferIndices[threadId][iMessage]);
3234  std::vector<size_t> &threadRecvBufferIndex(threadRecvBufferIndices[threadId][iMessage]);
3235  size_t recvBufferSize = recvIndices[iMessage].size();
3236  size_t numThreadRecv = threadRecvIndex.size();
3237  for(int componentIndex = 0;componentIndex < numComponents;componentIndex++){
3238  for(size_t pointIndex = 0;pointIndex < numThreadRecv;pointIndex++){
3239  size_t recvBufferIndex = threadRecvBufferIndex[pointIndex] +
3240  (componentId+componentIndex)*recvBufferSize;
3241  size_t targetBufferIndex = componentIndex*numPointsBuffer + threadRecvIndex[pointIndex];
3242  targetBuffer[targetBufferIndex] = messageBuffers.recvFlagBuffers[iMessage][recvBufferIndex];
3243  }
3244  }
3245  }
3246  return(0);
3247  };
3248 
3249 
3252  int halo::CreateThreadRecvIndices(int threadId)
3253  {
3254  if(remoteHaloExtents.empty()){
3255  std::cout << "halo::CreateThreadRecvIndices: WARNING No remote halos." << std::endl;
3256  return(0);
3257  }
3258 
3259 
3260  // thread extent (excluding remote halos) wrt global grid size
3261  pcpp::IndexIntervalType &threadExtent(threadExtents[threadId]);
3262  // buffer extent includes remote halos
3263  pcpp::IndexIntervalType &threadBufferExtent(threadBufferExtents[threadId]);
3264 
3265 
3266  int numRecvBuffers = remoteHaloExtents.size();
3267  // if(threadRecvIndices[threadId].empty())
3268  // threadRecvIndices[threadId].resize(numRecvBuffers);
3269  threadRecvBufferIndices[threadId].resize(numRecvBuffers);
3270  threadHaloBufferIndices[threadId].resize(numRecvBuffers);
3271  // std::cout << "Thread Extent: ";
3272  // threadExtent.PrettyPrint(std::cout);
3273  // std::cout << std::endl;
3274  for(int i = 0;i < numRecvBuffers;i++){
3275  // Check to make sure thread needs to recv (i.e. it owns part of part boundary)
3276  int recvDirection = i/2;
3277  int recvOrientation = i%2;
3278  std::vector<size_t> recvIndices;
3279 
3280  pcpp::IndexIntervalType &remoteHaloExtent(remoteHaloExtents[i]);
3281  pcpp::IndexIntervalType &remoteHaloBufferExtent(remoteHaloBufferExtents[i]);
3282 
3283  std::vector<size_t> &threadRecvBufferIndex(threadRecvBufferIndices[threadId][i]);
3284  std::vector<size_t> &threadHaloBufferIndex(threadHaloBufferIndices[threadId][i]);
3285 
3286  if(!localPartitionExtent.empty()){
3287  pcpp::IndexIntervalType localBufferExtent;
3288  localBufferExtent.InitSimple(localBufferSizes);
3289  pcpp::IndexIntervalType overlapInterval(threadBufferExtent.Overlap(remoteHaloBufferExtent));
3290  if(overlapInterval.NNodes() > 0){
3291  remoteHaloBufferExtent.GetFlatIndices(overlapInterval,threadRecvBufferIndex);
3292  localBufferExtent.GetFlatIndices(overlapInterval,threadHaloBufferIndex);
3293  }
3294  } else {
3295 
3296 
3297  size_t numRemoteHaloNodes = remoteHaloExtent.NNodes();
3298  // std::cout << "numRemoteHaloNodes = " << numRemoteHaloNodes << std::endl;
3299  // std::cout << "HaloExtent(" << i << "): ";
3300  // remoteHaloExtent.PrettyPrint(std::cout);
3301  // std::cout << std::endl;
3302  if(numRemoteHaloNodes > 0){
3303 
3304  if(recvOrientation){
3305  if(threadExtent[recvDirection].second != globalPartitionExtent[recvDirection].second)
3306  continue;
3307  } else {
3308  if(threadExtent[recvDirection].first != globalPartitionExtent[recvDirection].first)
3309  continue;
3310  }
3311 
3312  pcpp::IndexIntervalType checkExtent(threadExtent);
3313  checkExtent[recvDirection] = remoteHaloExtent[recvDirection];
3314  pcpp::IndexIntervalType threadRecvExtent(checkExtent.Overlap(remoteHaloExtent));
3315  size_t totalNumComponents = 0;
3316  size_t numThreadRecvNodes = threadRecvExtent.NNodes();
3317  // std::cout << "numThreadRecvNodes = " << numThreadRecvNodes << std::endl;
3318 
3319  }
3320  }
3321  }
3322  return(0);
3323  };
3324 
3325 
3326  // int halo::Wait(CommunicatorType &inComm)
3327  // {
3328  // if(!haveRecvData)
3329  // return(0);
3330  // inComm.WaitAll();
3331  // return(0);
3332  // };
3333 
3334  }
3335 }
simulation::grid::halo::messagebuffers::recvFlagBuffers
std::vector< int * > recvFlagBuffers
Definition: Grid.H:69
simulation::grid::pbsintervals::partitionBufferInterval
pcpp::IndexIntervalType partitionBufferInterval
Partition interval wrt the local buffer dimensions.
Definition: Grid.H:374
simulation::grid::parallel_blockstructured::gridJacobian
std::vector< double > gridJacobian
Definition: Grid.H:676
simulation::grid::halo::SetPartitionInterval
void SetPartitionInterval(const pcpp::IndexIntervalType &inExtent)
Definition: Grid.H:83
simulation::grid::pbsintervals::partitionInterval
pcpp::IndexIntervalType partitionInterval
Partition interval wrt the grid dimensions.
Definition: Grid.H:372
simulation::grid::pbsintervals::partitionSizes
std::vector< size_t > partitionSizes
Number of points in each dimension for local partition.
Definition: Grid.H:370
simulation::grid::subregion::regionInterval
pcpp::IndexIntervalType regionInterval
Definition: Grid.H:39
pcpp::util::SimplePartitionInterval
int SimplePartitionInterval(const pcpp::IndexIntervalType &inInterval, std::vector< bool > partDirection, int partID, int numPart, pcpp::IndexIntervalType &outInterval, std::vector< int > &numPartitions)
Multi-dimensional interval partitioning (non-MPI)
Definition: PCPPIntervalUtils.C:151
stencilSizes
void size_t int size_t int size_t int int * stencilSizes
Definition: OperatorKernels.H:9
simulation::grid::parallel_blockstructured::gridType
int gridType
GridType.
Definition: Grid.H:654
simulation::grid::subregion::threadPartitionBufferIndices
std::vector< std::vector< size_t > > threadPartitionBufferIndices
Definition: Grid.H:46
ix::util::sizeextent::GetFlatIndices
void GetFlatIndices(const sizeextent &extent, ContainerType &indices) const
Definition: IndexUtil.H:302
simulation::grid::parallel_blockstructured::scaleFactor
std::vector< double > scaleFactor
Grid scaling factor.
Definition: Grid.H:688
simulation::grid::parallel_blockstructured::ExchangeCoordinates
int ExchangeCoordinates()
Definition: Grid.C:2075
simulation::grid::parallel_blockstructured::threadDecompDirs
std::vector< int > threadDecompDirs
Definition: Grid.H:699
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::comm::CommunicatorObject::ARecvBuf
int ARecvBuf(DataType *recvbuf, size_t nVal, unsigned int remote_rank, int tag=0)
Definition: COMM.H:163
ix::util::sizeextent::Flatten
void Flatten(ContainerType &output) const
Definition: IndexUtil.H:274
plascom2::operators::stencilset::stencilSizes
int * stencilSizes
The number of weights for each stencil.
Definition: Stencil.H:102
simulation::grid::parallel_blockstructured::ComputeCurvilinearMetrics2D
int ComputeCurvilinearMetrics2D(std::ostream &messageStream)
Definition: Grid.C:569
ix::comm::CommunicatorObject::CartNeighbors
int CartNeighbors(std::vector< int > &neighborRanks)
Definition: COMM.C:173
operators::applyoperator
subroutine applyoperator(numDim, dimSizes, numComponents, numPoints, opDir, opInterval, numStencils, stencilSizes, stencilStarts, numValues, stencilWeights, stencilOffsets, stencilID, U, dU)
applyoperator applies an operator specified as a stencil set to the provided state data ...
Definition: Operators.f90:36
operators::determinant3x3
subroutine determinant3x3(numPoints, bufferSize, bufferInterval, inMatrix, matrixDeterminant)
Computes determinant of 3x3 matrix.
Definition: Operators.f90:1449
simulation::grid::halo::messagebuffers::sendFlagBuffers
std::vector< int * > sendFlagBuffers
Definition: Grid.H:68
simulation::grid::topoNames
static const char * topoNames[]
Definition: Grid.H:30
simulation::grid::halo::CreateSimpleSendIndices
int CreateSimpleSendIndices()
Creates the flat buffer indices for each whole local halo extent.
Definition: Grid.C:2659
simulation::grid::pbsintervals::bufferSizes
std::vector< size_t > bufferSizes
Buffer size in each dimension.
Definition: Grid.H:368
simulation::grid::parallel_blockstructured::PartitionInterval
pcpp::IndexIntervalType & PartitionInterval()
Definition: Grid.H:500
simulation::grid::parallel_blockstructured::SetNumThreads
void SetNumThreads(int numThreads)
Definition: Grid.C:2291
simulation::grid::parallel_blockstructured::ComputeCurvilinearMetrics3D
int ComputeCurvilinearMetrics3D(std::ostream &messageStream)
Definition: Grid.C:966
simulation::grid::subregion::regionPartitionBufferInterval
pcpp::IndexIntervalType regionPartitionBufferInterval
Definition: Grid.H:41
simulation::grid::parallel_blockstructured::GenerateCoordinates
int GenerateCoordinates(std::ostream &)
Definition: Grid.C:1628
stencilStarts
void size_t int size_t int size_t int int int * stencilStarts
Definition: OperatorKernels.H:9
gridSizes
void const size_t const size_t * gridSizes
Definition: EulerKernels.H:10
simulation::grid::parallel_blockstructured::jacobianMatrix2
std::vector< double > jacobianMatrix2
Definition: Grid.H:678
simulation::grid::halo::messagebuffers::recvBuffers
std::vector< double * > recvBuffers
Definition: Grid.H:67
simulation::grid::parallel_blockstructured::GenerateGrid
int GenerateGrid(int(*)(const std::vector< int > &, const std::vector< double > &, const std::vector< size_t > &, const std::vector< size_t > &, std::vector< double > &), std::ostream &)
Definition: Grid.C:1822
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
plascom2::operators::stencilset::stencilOffsets
int * stencilOffsets
The offsets wrt the grid point at which the stencil is being applied.
Definition: Stencil.H:106
simulation::grid::parallel_blockstructured::CurvilinearMetrics
int CurvilinearMetrics(std::ostream &messageStream)
Definition: Grid.C:1205
simulation::grid::subregion::regionPartitionInterval
pcpp::IndexIntervalType regionPartitionInterval
Definition: Grid.H:40
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
pcpp::report::CartesianSetup
void CartesianSetup(std::ostream &outStream, const std::vector< int > &cartCoords, const std::vector< int > &cartDims)
Definition: PCPPReport.C:199
operators::yxy
subroutine yxy(numDim, numPoints, bufferSize, bufferInterval, X, Y)
YXY point-wise operator performing Y = XY (all vectors)
Definition: Operators.f90:731
simulation::grid::halo::messagebuffers::numPointsSend
std::vector< size_t > numPointsSend
Definition: Grid.H:64
simulation::grid::parallel_blockstructured::flagData
std::vector< int > flagData
Definition: Grid.H:681
simulation::grid::halo::SetNumThreads
void SetNumThreads(int numThreadsIn)
Definition: Grid.C:2479
simulation::grid::pbsintervals::gridSizes
std::vector< size_t > gridSizes
Number of points in each dimension for the grid.
Definition: Grid.H:366
simulation::grid::parallel_blockstructured::GridGenerationOptions
std::vector< int > & GridGenerationOptions()
Definition: Grid.H:649
ix::comm::CommunicatorObject::ASendBuf
int ASendBuf(DataType *sendBuf, size_t nVal, unsigned int remote_rank, int tag=-1)
Definition: COMM.H:155
simulation::grid::parallel_blockstructured::GradIJK
int GradIJK(int numComponents, double *sourceBuffer, double *targetBuffer, std::ostream &messageStream, bool exchange)
Definition: Grid.C:844
simulation::grid::parallel_blockstructured::BufferSizes
const std::vector< size_t > & BufferSizes() const
Definition: Grid.H:459
ix::util::sizeextent::NNodes
size_t NNodes() const
Definition: IndexUtil.H:254
simulation::grid::InitSubRegionFromString
int InitSubRegionFromString(const std::string &inString, const std::vector< size_t > &gridSizes, subregion &inSubRegion)
Definition: Grid.C:174
simulation::grid::parallel_blockstructured::ComputeRectilinearMetrics
int ComputeRectilinearMetrics(std::ostream &messageStream)
Definition: Grid.C:328
stencilOffsets
void size_t int size_t int size_t int int int int double int * stencilOffsets
Definition: OperatorKernels.H:9
simulation::grid::parallel_blockstructured::numPointsPartition
size_t numPointsPartition
Number of points in the partition.
Definition: Grid.H:663
simulation::grid::halo::CreateMessageBuffers
int CreateMessageBuffers(int numComponents, int componentSize=8)
Definition: Grid.C:2850
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::numPointsBuffer
size_t numPointsBuffer
Number of points in the buffer.
Definition: Grid.H:661
simulation::grid::parallel_blockstructured::ExchangeJacobianMatrix
int ExchangeJacobianMatrix(std::ostream &messageStream)
Definition: Grid.C:2159
simulation::grid::parallel_blockstructured::PartitionBufferInterval
pcpp::IndexIntervalType & PartitionBufferInterval()
Definition: Grid.H:519
simulation::grid::parallel_blockstructured::periodicDirs
std::vector< bool > periodicDirs
Definition: Grid.H:696
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
simulation::grid::subregion::threadPartitionIntervals
std::vector< pcpp::IndexIntervalType > threadPartitionIntervals
Definition: Grid.H:44
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
simulation::grid::parallel_blockstructured::FlagMask
size_t FlagMask(int *iMask, int flagValue=0, int flagBit=1) const
Definition: Grid.C:1794
simulation::grid::subregion::partitionBufferIndices
std::vector< size_t > partitionBufferIndices
Definition: Grid.H:42
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
simulation::grid::halo::SetRemoteHaloExtents
void SetRemoteHaloExtents(const std::vector< pcpp::IndexIntervalType > haloExtents)
Definition: Grid.C:2432
simulation::grid::pbsintervals::threadBufferIntervals
std::vector< pcpp::IndexIntervalType > threadBufferIntervals
Buffer interval owned by each thread (omp)
Definition: Grid.H:377
grid
Definition: Grid.f90:1
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
simulation::grid::parallel_blockstructured::jacobianMatrix
std::vector< double > jacobianMatrix
Definition: Grid.H:677
simulation::grid::parallel_blockstructured::ComputeJacobianMatrix2
int ComputeJacobianMatrix2(std::ostream &messageStream)
Calculates the gradient of the regular/naive J^-1 d^2(xyz)/d(ijk)^2.
Definition: Grid.C:946
ix::util::sizeextent::Copy
void Copy(const sizeextent &inExtent)
Definition: IndexUtil.H:228
ix::comm::CommunicatorObject
Main encapsulation of MPI.
Definition: COMM.H:62
simulation::grid::ResolveTopoName
int ResolveTopoName(const std::string &inName)
Definition: Grid.C:135
a
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 double * a
Definition: OperatorKernels.H:29
simulation::grid::parallel_blockstructured::SetupThreads
int SetupThreads(std::vector< int > &numThreadPartitions)
Definition: Grid.C:2174
plascom2::operators::stencilset
Encapsulation for a collection of operator stencils.
Definition: Stencil.H:26
simulation::grid::halo::SetThreadExtent
void SetThreadExtent(int myThreadId, pcpp::IndexIntervalType threadInterval)
Definition: Grid.C:2461
simulation::grid::parallel_blockstructured::xyzMessageId
int xyzMessageId
Definition: Grid.H:700
plascom2::operators::stencilset::numValues
int numValues
The total number of weights for all stencils (reqd for Fortran)
Definition: Stencil.H:98
simulation::grid::parallel_blockstructured::ComputeUniformRectangularMetrics
int ComputeUniformRectangularMetrics(std::ostream &messageStream)
Definition: Grid.C:252
ix::util::sizeextent::Overlap
void Overlap(const sizeextent &inextent, sizeextent &outextent) const
Definition: IndexUtil.H:324
simulation::grid::parallel_blockstructured::physicalExtent
std::vector< double > physicalExtent
Physical grid extent (when applicable)
Definition: Grid.H:684
simulation::grid::parallel_blockstructured::dX
std::vector< double > dX
Grid spacing data structure (size depends on gridTopoType)
Definition: Grid.H:686
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::parallel_blockstructured::Halo
HaloType & Halo()
Definition: Grid.H:544
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
simulation::grid::parallel_blockstructured::ExchangeNodalData
int ExchangeNodalData(int numComponents, double *dataBuffer)
Definition: Grid.C:818
bufferSizes
void const size_t const size_t * bufferSizes
Definition: MetricKernels.H:19
simulation::grid::parallel_blockstructured::subRegions
std::vector< subregion > subRegions
Sub-regions and boundaries.
Definition: Grid.H:704
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
simulation::grid::parallel_blockstructured::PeriodicDirs
std::vector< bool > & PeriodicDirs()
Definition: Grid.H:576
simulation::grid::parallel_blockstructured::PhysicalExtent
const std::vector< double > & PhysicalExtent() const
Definition: Grid.H:533
simulation::grid::pbsintervals::threadPartitionIntervals
std::vector< pcpp::IndexIntervalType > threadPartitionIntervals
Partition intervals owned by thread (omp)
Definition: Grid.H:381
operators::zxy
subroutine zxy(numDim, numPoints, bufferSize, bufferInterval, X, Y, Z)
ZXY point-wise operator performing Z = XY (all vectors)
Definition: Operators.f90:679
simulation::grid::parallel_blockstructured::CoordinateData
std::vector< double > & CoordinateData()
Definition: Grid.H:503
simulation::grid::parallel_blockstructured::partitionIntervals
pbsintervals partitionIntervals
The relevant partition intervals.
Definition: Grid.H:665
simulation::grid::parallel_blockstructured::dimExtensions
std::vector< size_t > dimExtensions
Number of "halo/ghost" points in each (+/-) of the 2*numDim directions.
Definition: Grid.H:667
plascom2::operators::stencilset::stencilWeights
double * stencilWeights
The stencil weights.
Definition: Stencil.H:108
stencilWeights
void size_t int size_t int size_t int int int int double * stencilWeights
Definition: OperatorKernels.H:9
numStencils
void size_t int size_t int size_t int * numStencils
Definition: OperatorKernels.H:9
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
simulation::grid::parallel_blockstructured::AllocateCoordinateData
void AllocateCoordinateData()
Definition: Grid.H:524
simulation::grid::pbsintervals::threadPartitionBufferIntervals
std::vector< pcpp::IndexIntervalType > threadPartitionBufferIntervals
Partition buffer interval owned by each thread (omp)
Definition: Grid.H:379
simulation::grid::halo::messagebuffers::numPointsRecv
std::vector< size_t > numPointsRecv
Definition: Grid.H:65
yaxpy
subroutine yaxpy(N1, N2, N3, localInterval, a, X, Y)
Definition: RKUtil.f90:113
operators::yaxm1
subroutine yaxm1(numDim, numPoints, bufferSize, bufferInterval, a, X, Y)
YAXM1 point-wise operator performing Y = a/X (scalar a)
Definition: Operators.f90:1396
simulation::grid::ResolveBufferInterval
pcpp::IndexIntervalType ResolveBufferInterval(const pcpp::IndexIntervalType &partitionInterval, const pcpp::IndexIntervalType &partitionBufferInterval, const pcpp::IndexIntervalType &partitionSubInterval)
Definition: Grid.C:153
simulation::grid::ResolveTopoType
std::string ResolveTopoType(int inType)
Definition: Grid.C:144
simulation::grid::parallel_blockstructured::Dimension
int Dimension() const
Definition: Grid.H:428
ix::util::sizeextent
Simple Block Structured Mesh object.
Definition: IndexUtil.H:21
simulation::grid::halo::SetGridInterval
void SetGridInterval(const pcpp::IndexIntervalType &inExtent)
Definition: Grid.H:82
simulation::grid::parallel_blockstructured::stencilConnectivity
int * stencilConnectivity
Definition: Grid.H:691
simulation::grid::subregion::threadPartitionBufferIntervals
std::vector< pcpp::IndexIntervalType > threadPartitionBufferIntervals
Definition: Grid.H:45
simulation::grid::parallel_base::Communicator
fixtures::CommunicatorType & Communicator() const
Definition: Grid.H:350
numComponents
void size_t int * numComponents
Definition: OperatorKernels.H:9
operators::metricsum4
subroutine metricsum4(numDim, numPoints, bufferSize, bufferInterval, buf1, buf2, buf3, buf4, buf5, buf6, buf7, metricSum)
Computes buf1*buf4 - buf2*buf3 + buf7*(buf5 - buf6)
Definition: Operators.f90:1536
simulation::grid::parallel_blockstructured::errorState
int errorState
Definition: Grid.H:706
simulation::grid::parallel_blockstructured::GridGenerationParameters
std::vector< double > & GridGenerationParameters()
Definition: Grid.H:643
operators::xax
subroutine xax(numDim, numPoints, bufferSize, bufferInterval, a, X)
XAX point-wise operator performing X = aX (scalar a)
Definition: Operators.f90:1119
simulation::grid::parallel_blockstructured::gridCoordinates
std::vector< double > gridCoordinates
The coordinate buffer.
Definition: Grid.H:671
simulation::grid::parallel_blockstructured::periodicLengths
std::vector< double > periodicLengths
Definition: Grid.H:697
simulation::grid::parallel_blockstructured::Finalize
int Finalize(bool allocateCoordinateData=false)
Definition: Grid.C:2318
simulation::grid::parallel_blockstructured::gridMetric
std::vector< double > gridMetric
Grid metrics (sizes depend on gridTopoType)
Definition: Grid.H:675
ix::comm::CommunicatorObject::CartDimensions
std::vector< int > & CartDimensions()
Definition: COMM.H:118
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
simulation::grid::halo::DestroyMessageBuffer
int DestroyMessageBuffer(int messageId)
Definition: Grid.C:2809
simulation::grid::halo::SetLocalBufferSizes
void SetLocalBufferSizes(const std::vector< size_t > &inBufferSizes)
Definition: Grid.C:2423
FC_MODULE
void FC_MODULE(euler, flux1d, EULER, FLUX1D)(const int *numDim
simulation::grid::halo::CreateThreadRecvIndices
int CreateThreadRecvIndices(int threadId)
Create receive buffers for state data from neighboring processors.
Definition: Grid.C:3252
ix::comm::CommunicatorObject::CartCoordinates
std::vector< int > & CartCoordinates()
Definition: COMM.H:117
simulation::grid::parallel_blockstructured::stencilSet
plascom2::operators::stencilset * stencilSet
Differential operator.
Definition: Grid.H:690
simulation::grid::halo::CreateThreadSendIndices
int CreateThreadSendIndices(int threadId)
Creates send buffers for each whole local halo extent.
Definition: Grid.C:2702