Block-Structured AMR Software Framework
AMReX_ParticleDataAdaptorI.H
Go to the documentation of this file.
1 #include "Profiler.h"
2 #include "Error.h"
3 #include "SVTKUtils.h"
4 #include "MeshMetadata.h"
5 // svtk includes
6 #include <svtkPolyData.h>
7 #include <svtkDataSetAttributes.h>
8 #include <svtkCellData.h>
9 #include <svtkPointData.h>
10 #include <svtkMultiBlockDataSet.h>
11 
12 
13 
14 // local includes
16 #include <AMReX_InSituUtils.H>
17 
18 namespace amrex
19 {
20 //-----------------------------------------------------------------------------
21 template <typename ParticleType, int NArrayReal, int NArrayInt>
22 ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>*
23 ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::New()
24 {
25  auto result = new ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>;
26  result->InitializeObjectBase();
27  return result;
28 }
29 // senseiNewMacro(ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>);
30 
31 //-----------------------------------------------------------------------------
32 template <typename ParticleType, int NArrayReal, int NArrayInt>
33 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::SetDataSource(
35  const std::map<std::string, std::vector<int>> & rStructs,
36  const std::map<std::string, int> & iStructs,
37  const std::map<std::string, std::vector<int>> & rArrays,
38  const std::map<std::string, int> & iArrays)
39 {
40  // set the data source on the particles
41  this->ReleaseData();
42  this->m_particles = particles;
43 
44  // set the array names
45  int ret = this->SetArrayNames(rStructs, iStructs, rArrays, iArrays);
46  if(ret)
47  {
48  SENSEI_ERROR("problem with array names in SetDataSource");
49  return ret;
50  }
51  return ret;
52 }
53 
54 //-----------------------------------------------------------------------------
55 template <typename ParticleType, int NArrayReal, int NArrayInt>
56 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::SetArrayNames(
57  const std::map<std::string, std::vector<int>> & rStructs,
58  const std::map<std::string, int> & iStructs,
59  const std::map<std::string, std::vector<int>> & rArrays,
60  const std::map<std::string, int> & iArrays)
61 {
62  if(rStructs.size() <= ParticleType::NReal)
63  {
64  // check that no indices from the rStructs map exceed allowable value
65  for(auto s : rStructs)
66  {
67  for(auto i : s.second)
68  {
69  if(i >= ParticleType::NReal)
70  {
71  SENSEI_ERROR("rStruct index exceeds internal storage size");
72  return -1;
73  }
74  }
75  }
76  m_realStructs = rStructs;
77  }
78  else
79  {
80  SENSEI_ERROR("rStructs array size exceeds internal storage size");
81  return -1;
82  }
83 
84  if(iStructs.size() <= ParticleType::NInt)
85  {
86  // check that no indices from the iStructs map exceed allowable value
87  for(auto s : iStructs)
88  {
89  if(s.second >= ParticleType::NInt)
90  {
91  SENSEI_ERROR("iStructs index exceeds internal storage size");
92  return -1;
93  }
94  }
95  m_intStructs = iStructs;
96  }
97  else
98  {
99  SENSEI_ERROR("iStructs array size exceeds internal storage size");
100  return -1;
101  }
102 
103 
104  if(rArrays.size() <= NArrayReal)
105  {
106  // check that no indices from the rArrays map exceed allowable value
107  for(auto s : rArrays)
108  {
109  for(auto i : s.second)
110  {
111  if(i >= NArrayReal)
112  {
113  SENSEI_ERROR("rArrays index exceeds internal storage size");
114  return -1;
115  }
116  }
117  }
118  m_realArrays = rArrays;
119  }
120  else
121  {
122  SENSEI_ERROR("rArrays array size exceeds internal storage size");
123  return -1;
124  }
125  if(iArrays.size() <= NArrayInt)
126  {
127  // check that no indices from the iArrays map exceed allowable value
128  for(auto s : iArrays)
129  {
130  if(s.second >= NArrayInt)
131  {
132  SENSEI_ERROR("iArray index exceeds internal storage size");
133  return -1;
134  }
135  }
136  m_intArrays = iArrays;
137  }
138  else
139  {
140  SENSEI_ERROR("iArrays array size exceeds internal storage size");
141  return -1;
142  }
143  return 0;
144 }
145 //-----------------------------------------------------------------------------
146 template <typename ParticleType, int NArrayReal, int NArrayInt>
147 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::GetNumberOfMeshes(unsigned int &numMeshes)
148 {
149  numMeshes = 1;
150  return 0;
151 }
152 
153 //-----------------------------------------------------------------------------
154 #if SENSEI_VERSION_MAJOR < 3
155 template <typename ParticleType, int NArrayReal, int NArrayInt>
156 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::GetMeshName(
157  unsigned int id, std::string &meshName)
158 {
159  meshName = m_particlesName;
160  return 0;
161 }
162 #endif
163 
164 //-----------------------------------------------------------------------------
165 #if SENSEI_VERSION_MAJOR < 3
166 template <typename ParticleType, int NArrayReal, int NArrayInt>
167 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::GetMeshHasGhostNodes(
168  const std::string &meshName,
169  int &nLayers)
170 {
171  nLayers = 0;
172  return 0;
173 }
174 #endif
175 
176 //-----------------------------------------------------------------------------
177 #if SENSEI_VERSION_MAJOR < 3
178 template <typename ParticleType, int NArrayReal, int NArrayInt>
179 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::GetMeshHasGhostCells(
180  const std::string &meshName,
181  int &nLayers)
182 {
183  nLayers = 0;
184  return 0;
185 }
186 #endif
187 
188 //-----------------------------------------------------------------------------
189 #if SENSEI_VERSION_MAJOR < 3
190 template <typename ParticleType, int NArrayReal, int NArrayInt>
191 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::GetNumberOfArrays(
192  const std::string &meshName,
193  int association,
194  unsigned int &numberOfArrays)
195 {
196  numberOfArrays = 0;
197  if(association == svtkDataObject::POINT)
198  {
199  numberOfArrays = m_realStructs.size()
200  + m_intStructs.size()
201  + m_realArrays.size()
202  + m_intArrays.size();
203  }
204  return 0;
205 }
206 #endif
207 //-----------------------------------------------------------------------------
208 #if SENSEI_VERSION_MAJOR < 3
209 template <typename ParticleType, int NArrayReal, int NArrayInt>
210 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::GetArrayName(
211  const std::string &meshName,
212  int association,
213  unsigned int index,
214  std::string &arrayName)
215 {
216  if(association == svtkDataObject::POINT)
217  {
218  if(index < m_realStructs.size())
219  {
220  auto a = m_realStructs.begin() + index;
221  arrayName = *(a).first;
222  return 0;
223  }
224  if(index < m_intStructs.size())
225  {
226  int ind = index - m_realStructs.size();
227  auto a = m_intStructs.begin() + ind;
228  arrayName = *(a).first;
229  return 0;
230  }
231  if(index < m_realArrays.size())
232  {
233  int ind = index - m_realStructs.size() - m_intStructs.size();
234  auto a = m_realArrays.begin() + ind;
235  arrayName = *(a).first;
236  return 0;
237  }
238  if(index < m_intArrays.size())
239  {
240  int ind = index - m_realStructs.size() - m_intStructs.size() - m_realArrays.size();
241  auto a = m_intArrays.begin() + ind;
242  arrayName = *(a).first;
243  return 0;
244  }
245  }
246 
247  return -1;
248 }
249 #endif
250 
251 //-----------------------------------------------------------------------------
252 template <typename ParticleType, int NArrayReal, int NArrayInt>
253 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::GetMesh(
254  const std::string &meshName,
255  bool structureOnly,
256  svtkDataObject *&mesh)
257 {
258  mesh = nullptr;
259  int nprocs = 1;
260  int rank = 0;
261  MPI_Comm_size(this->GetCommunicator(), &nprocs);
262  MPI_Comm_rank(this->GetCommunicator(), &rank);
263 
264  if (meshName != m_particlesName)
265  {
266  SENSEI_ERROR("No mesh named \"" << meshName << "\"")
267  return -1;
268  }
269  svtkMultiBlockDataSet* mb = svtkMultiBlockDataSet::New();
270 
271  if (structureOnly)
272  {
273  mesh = mb;
274  return 0;
275  }
276 
277  mb->SetNumberOfBlocks(nprocs);
278  svtkPolyData *pd = BuildParticles();
279  mb->SetBlock(rank, pd);
280  pd->Delete();
281  mesh = mb;
282 
283  return 0;
284 }
285 
286 //-----------------------------------------------------------------------------
287 template <typename ParticleType, int NArrayReal, int NArrayInt>
288 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::AddGhostNodesArray(
289  svtkDataObject*,
290  const std::string &meshName)
291 {
292  if (meshName != m_particlesName)
293  {
294  SENSEI_ERROR("no mesh named \"" << meshName << "\"")
295  return -1;
296  }
297  return 0;
298 }
299 
300 //-----------------------------------------------------------------------------
301 template <typename ParticleType, int NArrayReal, int NArrayInt>
302 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::AddGhostCellsArray(
303  svtkDataObject*,
304  const std::string &meshName)
305 {
306  if (meshName != m_particlesName)
307  {
308  SENSEI_ERROR("no mesh named \"" << meshName << "\"")
309  return -1;
310  }
311  return 0;
312 }
313 
314 //-----------------------------------------------------------------------------
315 template <typename ParticleType, int NArrayReal, int NArrayInt>
316 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::AddArray(
317  svtkDataObject* mesh,
318  const std::string &meshName,
319  int association,
320  const std::string &arrayName)
321 {
322  if (meshName != m_particlesName)
323  {
324  SENSEI_ERROR("no mesh named \"" << meshName << "\"");
325  return -1;
326  }
327 
328  if (association != svtkDataObject::POINT)
329  {
330  SENSEI_ERROR("Invalid association " << association);
331  return -1;
332  }
333 
334  if(m_realStructs.find(arrayName) != m_realStructs.end())
335  {
336  return this->AddParticlesAOSRealArray(arrayName, mesh);
337  }
338 
339  if(m_intStructs.find(arrayName) != m_intStructs.end())
340  {
341  return this->AddParticlesAOSIntArray(arrayName, mesh);
342  }
343 
344  if(m_realArrays.find(arrayName) != m_realArrays.end())
345  {
346  return this->AddParticlesSOARealArray(arrayName, mesh);
347  }
348 
349  if(m_intArrays.find(arrayName) != m_intArrays.end())
350  {
351  return this->AddParticlesSOAIntArray(arrayName, mesh);
352  }
353 
354  SENSEI_ERROR("Invalid array name " << arrayName);
355  return -1;
356 }
357 
358 //-----------------------------------------------------------------------------
359 template <typename ParticleType, int NArrayReal, int NArrayInt>
360 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::ReleaseData()
361 {
362  this->m_particles = nullptr;
363  return 0;
364 }
365 
366 //-----------------------------------------------------------------------------
367 #if SENSEI_VERSION_MAJOR >= 3
368 template <typename ParticleType, int NArrayReal, int NArrayInt>
369 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::GetMeshMetadata(
370  unsigned int id,
371  sensei::MeshMetadataPtr &metadata)
372 {
373  sensei::TimeEvent<64> event("AmrMeshDataAdaptor::GetMeshMetadata");
374 
375  int nprocs = 1;
376  int rank = 0;
377  MPI_Comm_size(this->GetCommunicator(), &nprocs);
378  MPI_Comm_rank(this->GetCommunicator(), &rank);
379  if (id != 0 && id != 1)
380  {
381  SENSEI_ERROR("invalid mesh id " << id)
382  return -1;
383  }
384 
385 
386  // AMR data is always expected to be a global view
387  metadata->GlobalView = false; // tells if the information describes data
388  // on this rank or all ranks. Passed into
389  // Set methods, Get methods generate the
390  // desired view.
391 
392  // name of mesh (all)
393  metadata->MeshName = m_particlesName;
394 
395  // container mesh type (all)
396  metadata->MeshType = SVTK_MULTIBLOCK_DATA_SET;
397 
398  // block mesh type (all)
399  metadata->BlockType = SVTK_POLY_DATA;
400 
401  // global number of blocks (all)
402  metadata->NumBlocks = nprocs;
403 
404  // number of blocks on each rank (all)
405  metadata->NumBlocksLocal = {1};
406 
407  // global cell index space extent (Cartesian, AMR, optional)
408  // std::array<int,6> Extent;
409 
410  // global bounding box (all, optional)
411  // std::array<double,6> Bounds;
412 
413  // type enum of point data (unstructured, optional)
414 #ifdef AMREX_SINGLE_PRECISION_PARTICLES
415  metadata->CoordinateType = SVTK_FLOAT;
416 #else
417  metadata->CoordinateType = SVTK_DOUBLE;
418 #endif
419 
420  // total number of points in all blocks (all, optional)
421  // long NumPoints;
422 
423  // total number of cells in all blocks (all, optional)
424  // long NumCells;
425 
426  // total cell array size in all blocks (all, optional)
427  // long CellArraySize;
428 
429  // number of arrays (all)
430  metadata->NumArrays = m_realStructs.size()
431  + m_intStructs.size()
432  + m_realArrays.size()
433  + m_intArrays.size();
434 
435  // number of ghost cell layers (all)
436  metadata->NumGhostCells = 0;
437 
438  // number of ghost node layers (all)
439  metadata->NumGhostNodes = 0;
440 
441  // number of AMR levels (AMR)
442  // metadata->NumLevels;
443 
444  // non zero if the mesh does not change in time (all)
445  metadata->StaticMesh = 0;
446 
447  // name of each data array (all)
448  metadata->ArrayName = {};
449  for(auto s : m_realStructs)
450  {
451  metadata->ArrayName.push_back(s.first);
452  }
453  for(auto s : m_intStructs)
454  {
455  metadata->ArrayName.push_back(s.first);
456  }
457  for(auto s : m_realArrays)
458  {
459  metadata->ArrayName.push_back(s.first);
460  }
461  for(auto s : m_intArrays)
462  {
463  metadata->ArrayName.push_back(s.first);
464  }
465 
466  // centering of each data array (all)
467  metadata->ArrayCentering = {};
468  for(auto s : m_realStructs)
469  {
470  metadata->ArrayCentering.push_back(svtkDataObject::POINT);
471  }
472  for(auto s : m_intStructs)
473  {
474  metadata->ArrayCentering.push_back(svtkDataObject::POINT);
475  }
476  for(auto s : m_realArrays)
477  {
478  metadata->ArrayCentering.push_back(svtkDataObject::POINT);
479  }
480  for(auto s : m_intArrays)
481  {
482  metadata->ArrayCentering.push_back(svtkDataObject::POINT);
483  }
484 
485  // number of components of each array (all)
486  metadata->ArrayComponents = {};
487  for(auto s : m_realStructs)
488  {
489  metadata->ArrayComponents.push_back(s.second.size());
490  }
491  for(auto s : m_intStructs)
492  {
493  metadata->ArrayComponents.push_back(1);
494  }
495  for(auto s : m_realArrays)
496  {
497  metadata->ArrayComponents.push_back(s.second.size());
498  }
499  for(auto s : m_intArrays)
500  {
501  metadata->ArrayComponents.push_back(1);
502  }
503 
504  // type enum of each data array (all)
505  metadata->ArrayType = {};
506  for(auto s : m_realStructs)
507  {
508 #ifdef AMREX_SINGLE_PRECISION_PARTICLES
509  metadata->ArrayType.push_back(SVTK_FLOAT);
510 #else
511  metadata->ArrayType.push_back(SVTK_DOUBLE);
512 #endif
513  }
514  for(auto s : m_intStructs)
515  {
516  metadata->ArrayType.push_back(SVTK_INT);
517  }
518  for(auto s : m_realArrays)
519  {
520 #ifdef AMREX_SINGLE_PRECISION_PARTICLES
521  metadata->ArrayType.push_back(SVTK_FLOAT);
522 #else
523  metadata->ArrayType.push_back(SVTK_DOUBLE);
524 #endif
525  }
526  for(auto s : m_intArrays)
527  {
528  metadata->ArrayType.push_back(SVTK_INT);
529  }
530 
531  // global min,max of each array (all, optional)
532  // metadata->ArrayRange = {};
533 
534  // rank where each block resides (all, optional)
535  metadata->BlockOwner = {rank};
536 
537  // global id of each block (all, optional)
538  metadata->BlockIds = metadata->BlockOwner;
539 
540  // number of points for each block (all, optional)
541  auto nptsOnProc = this->m_particles->TotalNumberOfParticles(true, true);
542  metadata->BlockNumPoints = {nptsOnProc};
543 
544  // number of cells for each block (all, optional)
545  metadata->BlockNumCells = {nptsOnProc};
546 
547  // cell array size for each block (unstructured, optional)
548  metadata->BlockCellArraySize = {nptsOnProc};
549 
550  // note: for AMR BlockExtents and BlockBounds are always global
551  // index space extent of each block [i0,i1, j0,j1, k0,k1] (Cartesian, AMR, optional)
552  // std::vector<std::array<int,6>> BlockExtents;
553 
554  // bounds of each block [x0,x1, y0,y1, z0,z1] (all, optional)
555  // std::vector<std::array<double,6>> BlockBounds;
556  if (metadata->Flags.BlockBoundsSet())
557  {
558  std::array<double,6> bounds({
560  std::numeric_limits<double>::lowest(),
561 #if (AMREX_SPACEDIM == 1)
562  0.0, 0.0, 0.0, 0.0
563 #elif (AMREX_SPACEDIM == 2)
565  std::numeric_limits<double>::lowest(),
566  0.0, 0.0
567 #elif (AMREX_SPACEDIM == 3)
569  std::numeric_limits<double>::lowest(),
571  std::numeric_limits<double>::lowest()
572 #endif
573  });
574 
575  // loop over levels in the particle container
576  const auto& particles = this->m_particles->GetParticles();
577  for (int lev = 0; lev < particles.size(); lev++)
578  {
579  // loop over ParticleTiles on level
580  auto& pmap = particles[lev];
581  for (const auto& kv : pmap)
582  {
583  // loop over particles in ParticleTile
584  auto& particle_tile = kv.second;
585  auto ptd = particle_tile.getConstParticleTileData();
586 
587  // visit only the "real" particles, skip the "neighbor" particles.
588  long long numReal = particle_tile.numRealParticles();
589  for (long long i = 0; i < numReal; ++i)
590  {
591  const auto &part = ptd[i];
592  if (part.id() > 0)
593  {
594 #if (AMREX_SPACEDIM == 1)
595  bounds[0] = bounds[0] > part.pos(0) ? part.pos(0) : bounds[0];
596  bounds[1] = bounds[1] < part.pos(0) ? part.pos(0) : bounds[1];
597 #elif (AMREX_SPACEDIM == 2)
598  bounds[0] = bounds[0] > part.pos(0) ? part.pos(0) : bounds[0];
599  bounds[1] = bounds[1] < part.pos(0) ? part.pos(0) : bounds[1];
600  bounds[2] = bounds[2] > part.pos(1) ? part.pos(1) : bounds[2];
601  bounds[3] = bounds[3] < part.pos(1) ? part.pos(1) : bounds[3];
602 #elif (AMREX_SPACEDIM == 3)
603  bounds[0] = bounds[0] > part.pos(0) ? part.pos(0) : bounds[0];
604  bounds[1] = bounds[1] < part.pos(0) ? part.pos(0) : bounds[1];
605  bounds[2] = bounds[2] > part.pos(1) ? part.pos(1) : bounds[2];
606  bounds[3] = bounds[3] < part.pos(1) ? part.pos(1) : bounds[3];
607  bounds[4] = bounds[4] > part.pos(2) ? part.pos(2) : bounds[4];
608  bounds[5] = bounds[5] < part.pos(2) ? part.pos(2) : bounds[5];
609 #endif
610  }
611  }
612  }
613  }
614 
615  metadata->BlockBounds = {bounds};
616  }
617 
618  // min max of each array on each block.
619  // indexed by block then array. (all, optional)
620  // std::vector<std::vector<std::array<double,2>>> BlockArrayRange;
621  if (metadata->Flags.BlockArrayRangeSet())
622  {
623  SENSEI_ERROR("BlockArrayRange requested but not yet implemented.")
624  }
625 
626  // refinement ratio in i,j, and k directions for each level (AMR)
627  // std::vector<std::array<int,3>> RefRatio;
628 
629  // number of blocks in each level (AMR)
630  // std::vector<int> BlocksPerLevel;
631 
632  // AMR level of each block (AMR)
633  // std::vector<int> BlockLevel;
634 
635  // flag indicating presence of a periodic boundary in the i,j,k direction (all)
636  // std::array<int,3> PeriodicBoundary;
637 
638 #if defined(SENSEI_DEBUG)
639  metadata->Validate(this->GetCommunicator(), sensei::MeshMetadataFlags());
640  metadata->ToStream(std::cerr);
641 #endif
642 
643  return 0;
644 }
645 #endif
646 
647 //-----------------------------------------------------------------------------
648 template <typename ParticleType, int NArrayReal, int NArrayInt>
649 svtkPolyData* ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::BuildParticles()
650 {
651  // return particle data pd
652  svtkPolyData* pd = svtkPolyData::New();
653 
654  const auto& particles = this->m_particles->GetParticles();
655  long long numParticles = this->m_particles->TotalNumberOfParticles(true, true);
656 
657  // allocate vertex storage for particles
658 #ifdef AMREX_SINGLE_PRECISION_PARTICLES
659  svtkNew<svtkFloatArray> coords;
660 #else
661  svtkNew<svtkDoubleArray> coords;
662 #endif
663  coords->SetName("coords");
664  coords->SetNumberOfComponents(3);
665  coords->SetNumberOfTuples(numParticles);
666 #ifdef AMREX_SINGLE_PRECISION_PARTICLES
667  float *pCoords = coords->GetPointer(0);
668 #else
669  double *pCoords = coords->GetPointer(0);
670 #endif
671 
672  // use this to index into the SVTK array as we copy level by level and tile by
673  // tile
674  long long ptId = 0;
675 
676  // allocate connectivity array for particles
677  svtkNew<svtkCellArray> vertex;
678  vertex->AllocateExact(numParticles, 1);
679 
680  // points->SetNumberOfPoints(numParticles);
681 
682  // loop over levels in the particle container
683  for (int lev = 0; lev < particles.size(); lev++)
684  {
685  using MyParIter = ParIter_impl<ParticleType, NArrayReal, NArrayInt>;
686  for (MyParIter pti(*this->m_particles, lev); pti.isValid(); ++pti)
687  {
688  auto ptd = pti.GetParticleTile().getParticleTileData();
689  auto numReal = pti.numParticles();
690 
691  for (long long i = 0; i < numReal; ++i)
692  {
693  const auto& part = make_particle<ParticleType>{}(ptd, i);
694  if (part.id() > 0)
695  {
696  // add a vertex type cell
697  vertex->InsertNextCell(1);
698  vertex->InsertCellPoint(ptId);
699  // copy the particle coordinates
700 #if (AMREX_SPACEDIM == 1)
701  pCoords[0] = part.pos(0);
702  pCoords[1] = amrex_particle_real(0);
703  pCoords[2] = amrex_particle_real(0);
704 #elif (AMREX_SPACEDIM == 2)
705  pCoords[0] = part.pos(0);
706  pCoords[1] = part.pos(1);
707  pCoords[2] = amrex_particle_real(0);
708 #elif (AMREX_SPACEDIM == 3)
709  pCoords[0] = part.pos(0);
710  pCoords[1] = part.pos(1);
711  pCoords[2] = part.pos(2);
712 #endif
713  pCoords += 3;
714  ++ptId;
715  }
716  }
717  }
718  }
719 
720  // pass the particle coordinates into SVTK's point data structure.
721  svtkNew<svtkPoints> points;
722  points->SetData(coords);
723 
724  // add point and vertex data to output mesh
725  pd->SetPoints(points);
726  pd->SetVerts(vertex);
727 
728  return pd;
729 }
730 
731 //-----------------------------------------------------------------------------
732 template <typename ParticleType, int NArrayReal, int NArrayInt>
733 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::AddParticlesIDArray(
734  svtkDataObject* mesh)
735 {
736  auto svtk_particles = dynamic_cast<svtkPolyData*>(mesh);
737  const auto& particles = this->m_particles->GetParticles();
738  auto nptsOnProc = this->m_particles->TotalNumberOfParticles(true, true);
739 
740  // allocate a SVTK array for the data
741  svtkNew<svtkIntArray> idArray;
742  idArray->SetName("id");
743  idArray->SetNumberOfComponents(1);
744  idArray->SetNumberOfValues(nptsOnProc);
745 
746  // get pointer underlying to data
747  int *partIds = idArray->GetPointer(0);
748 
749  // loop over particles and append their cpu value to the list
750  using MyParIter = ParIter_impl<ParticleType, NArrayReal, NArrayInt>;
751  long ptId = 0;
752  for (int level = 0; level < particles.size(); ++level)
753  {
754  for (MyParIter pti(*this->m_particles, level); pti.isValid(); ++pti)
755  {
756  auto ptd = pti.GetParticleTile().getParticleTileData();
757  auto numReal = pti.numParticles();
758  for (long long i = 0; i < numReal; ++i)
759  {
760  const auto &part = make_particle<ParticleType>{}(ptd, i);
761  if (part.id() > 0)
762  {
763  partIds[i] = part.id();
764  }
765  }
766  partIds += numReal;
767  }
768  }
769 
770  // the association for this array is svtkDataObject::POINT
771  svtk_particles->GetPointData()->AddArray(idArray);
772 
773  return 0;
774 }
775 
776 //-----------------------------------------------------------------------------
777 template <typename ParticleType, int NArrayReal, int NArrayInt>
778 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::AddParticlesCPUArray(
779  svtkDataObject* mesh)
780 {
781  auto svtk_particles = dynamic_cast<svtkPolyData*>(mesh);
782  const auto& particles = this->m_particles->GetParticles();
783  auto nptsOnProc = this->m_particles->TotalNumberOfParticles(true, true);
784 
785  // allocate a SVTK array for the data
786  svtkNew<svtkIntArray> cpuArray;
787  cpuArray->SetName("cpu");
788  cpuArray->SetNumberOfComponents(1);
789  cpuArray->SetNumberOfValues(nptsOnProc);
790 
791  // get pointer to underlying data
792  int* partCpu = cpuArray->GetPointer(0);
793 
794  // loop over particles and append their cpu value to the list
795  using MyParIter = ParIter_impl<ParticleType, NArrayReal, NArrayInt>;
796  for (int level = 0; level < particles.size(); ++level)
797  {
798  for (MyParIter pti(*this->m_particles, level); pti.isValid(); ++pti)
799  {
800  auto ptd = pti.GetParticleTile().getParticleTileData();
801  auto numReal = pti.numParticles();
802  for (long long i = 0; i < numReal; ++i)
803  {
804  const auto &part = make_particle<ParticleType>{}(ptd, i);
805  if (part.id() > 0)
806  {
807  partCpu[i] = part.cpu();
808  }
809  }
810  partCpu += numReal;
811  }
812  }
813 
814  // the association for this array is svtkDataObject::POINT
815  svtk_particles->GetPointData()->AddArray(cpuArray);
816 
817  return 0;
818 }
819 
820 //-----------------------------------------------------------------------------
821 template <typename ParticleType, int NArrayReal, int NArrayInt>
822 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::AddParticlesSOARealArray(
823  const std::string &arrayName,
824  svtkDataObject* mesh)
825 {
826  const long nParticles = this->m_particles->TotalNumberOfParticles(true, true);
827 
828  // check that the name of the arrays is listed in the m_realArrays
829  RealDataMapType::iterator ait = m_realArrays.find(arrayName);
830  if (ait == m_realArrays.end())
831  {
832  SENSEI_ERROR("No real SOA named \"" << arrayName << "\"");
833  return -1;
834  }
835 
836  // get the indices of the structs on the particles
837  const std::vector<int> &indices = ait->second;
838  int nComps = indices.size();
839 
840  // check that the indices are within the templated storage spec
841  for(auto i : indices)
842  {
843  if(i >= NArrayReal)
844  {
845  SENSEI_ERROR("Index out of bounds for real SOA named \"" << arrayName << "\"");
846  return -1;
847  }
848  }
849 
850  // allocate the svtkArray
851 #ifdef AMREX_SINGLE_PRECISION_PARTICLES
852  svtkNew<svtkFloatArray> data;
853 #else
854  svtkNew<svtkDoubleArray> data;
855 #endif
856  data->SetName(arrayName.c_str());
857  data->SetNumberOfComponents(nComps);
858  data->SetNumberOfTuples(nParticles);
859 
860 #ifdef AMREX_SINGLE_PRECISION_PARTICLES
861  float* pData = data->GetPointer(0);
862 #else
863  double* pData = data->GetPointer(0);
864 #endif
865 
866  using MyParIter = ParIter_impl<ParticleType, NArrayReal, NArrayInt>;
867  for (int level = 0; level < this->m_particles->numLevels(); ++level)
868  {
869  for (MyParIter pti(*this->m_particles, level); pti.isValid(); ++pti)
870  {
871  auto& particle_attributes = pti.GetStructOfArrays();
872  auto ptd = pti.GetParticleTile().getParticleTileData();
873 
874  auto numReal = pti.numParticles();
875  // shuffle from the AMReX component order
876  for (int j = 0; j < nComps; ++j)
877  {
878  int compInd = indices[j];
879  // component to copy
880  const auto &realData = particle_attributes.GetRealData(compInd);
881 
882  for (long long i = 0; i < numReal; ++i)
883  {
884  const auto &part = ptd[i];
885  if (part.id() > 0)
886  {
887  pData[i*nComps + j] = realData[i];
888  }
889  }
890  }
891  pData += numReal * nComps;
892  }
893  }
894 
895  // get the block of the mesh
896  int rank = 0;
897  MPI_Comm_rank(this->GetCommunicator(), &rank);
898 
899  auto blocks = dynamic_cast<svtkMultiBlockDataSet*>(mesh);
900 
901  auto block = dynamic_cast<svtkPolyData*>(blocks->GetBlock(rank));
902  block->GetPointData()->AddArray(data);
903 
904  return 0;
905 }
906 
907 //-----------------------------------------------------------------------------
908 template <typename ParticleType, int NArrayReal, int NArrayInt>
909 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::AddParticlesSOAIntArray(
910  const std::string &arrayName,
911  svtkDataObject* mesh)
912 {
913  // get the particles from the particle container
914  auto nptsOnProc = this->m_particles->TotalNumberOfParticles(true, true);
915 
916  // check that the name of the arrays is listed in the m_intArrays
917  IntDataMapType::iterator ait = m_intArrays.find(arrayName);
918  if (ait == m_intArrays.end())
919  {
920  SENSEI_ERROR("No int SOA named \"" << arrayName << "\"");
921  return -1;
922  }
923 
924  // get the index of the structs on the particles
925  int index = ait->second;
926 
927  // check that the indices are within the templated storage spec
928  if(index >= NArrayInt)
929  {
930  SENSEI_ERROR("Index out of bounds for int SOA named \"" << arrayName << "\"");
931  return -1;
932  }
933 
934  svtkNew<svtkIntArray> data;
935  data->SetName(arrayName.c_str());
936  data->SetNumberOfComponents(1);
937  data->SetNumberOfValues(nptsOnProc);
938  int* pData = data->GetPointer(0);
939 
940  // fill array
941  using MyParIter = ParIter_impl<ParticleType, NArrayReal, NArrayInt>;
942  for (int level = 0; level< this->m_particles->numLevels(); level++)
943  {
944  for (MyParIter pti(*this->m_particles, level); pti.isValid(); ++pti)
945  {
946  auto& particle_attributes = pti.GetStructOfArrays();
947  auto ptd = pti.GetParticleTile().getParticleTileData();
948 
949  auto numReal = pti.numParticles();
950  // shuffle from the AMReX component order
951  // component to copy
952  const auto &intData = particle_attributes.GetIntData(index);
953 
954  for (long long i = 0; i < numReal; ++i)
955  {
956  const auto &part = ptd[i];
957  if (part.id() > 0)
958  {
959  pData[i] = intData[i];
960  }
961  }
962  pData += numReal;
963  }
964  }
965 
966  // get the block of the mesh
967  int rank = 0;
968  MPI_Comm_rank(this->GetCommunicator(), &rank);
969 
970  auto blocks = dynamic_cast<svtkMultiBlockDataSet*>(mesh);
971 
972  auto block = dynamic_cast<svtkPolyData*>(blocks->GetBlock(rank));
973  block->GetPointData()->AddArray(data);
974 
975  return 0;
976 }
977 
978 //-----------------------------------------------------------------------------
979 template <typename ParticleType, int NArrayReal, int NArrayInt>
980 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::AddParticlesAOSRealArray(
981  const std::string &arrayName,
982  svtkDataObject* mesh)
983 {
984  // get the particles from the particle container
985  const auto& particles = this->m_particles->GetParticles();
986  auto nptsOnProc = this->m_particles->TotalNumberOfParticles(true, true);
987 
988  // check that the name of the struct data is listed in the m_realStructs
989  RealDataMapType::iterator ait = m_realStructs.find(arrayName);
990  if (ait == m_realStructs.end())
991  {
992  SENSEI_ERROR("No real AOS named \"" << arrayName << "\"");
993  return -1;
994  }
995 
996  // get the indices of the struct on the particles
997  std::vector<int> indices = ait->second;
998  int nComps = indices.size();
999 
1000  // check that the indices are within the templated storage spec
1001  for (auto i : indices)
1002  {
1003  if (i >= ParticleType::NReal)
1004  {
1005  SENSEI_ERROR("Index out of bounds for real AOS named \"" << arrayName << "\"");
1006  return -1;
1007  }
1008  }
1009 
1010  // allocate the svtk array
1011 #ifdef AMREX_SINGLE_PRECISION_PARTICLES
1012  svtkNew<svtkFloatArray> data;
1013 #else
1014  svtkNew<svtkDoubleArray> data;
1015 #endif
1016 
1017  data->SetName(arrayName.c_str());
1018  data->SetNumberOfComponents(nComps);
1019  data->SetNumberOfTuples(nptsOnProc);
1020 
1021 #ifdef AMREX_SINGLE_PRECISION_PARTICLES
1022  float *pData = data->GetPointer(0);
1023 #else
1024  double *pData = data->GetPointer(0);
1025 #endif
1026 
1027  if constexpr(!ParticleType::is_soa_particle) {
1028 
1029  // copy the data from each level
1030  using MyParIter = ParIter_impl<ParticleType, NArrayReal, NArrayInt>;
1031  for (int level = 0; level<particles.size();level++)
1032  {
1033  for (MyParIter pti(*this->m_particles, level); pti.isValid(); ++pti)
1034  {
1035  auto& aos = pti.GetArrayOfStructs();
1036 
1037  auto numReal = pti.numParticles();
1038  // shuffle from the AMReX component order
1039  for (int j = 0; j < nComps; ++j)
1040  {
1041  for (long long i = 0; i < numReal; ++i)
1042  {
1043  const auto &part = aos[i];
1044  if (part.id() > 0)
1045  {
1046  pData[i*nComps + j] = part.rdata(indices[j]);
1047  }
1048  }
1049  }
1050  pData += numReal * nComps;
1051  }
1052  }
1053 
1054  // get the block of the mesh
1055  int rank = 0;
1056  MPI_Comm_rank(this->GetCommunicator(), &rank);
1057 
1058  auto blocks = dynamic_cast<svtkMultiBlockDataSet*>(mesh);
1059 
1060  auto block = dynamic_cast<svtkPolyData*>(blocks->GetBlock(rank));
1061  block->GetPointData()->AddArray(data);
1062 
1063  }
1064 
1065  return 0;
1066 }
1067 
1068 //-----------------------------------------------------------------------------
1069 template <typename ParticleType, int NArrayReal, int NArrayInt>
1070 int ParticleDataAdaptor<ParticleType, NArrayReal, NArrayInt>::AddParticlesAOSIntArray(
1071  const std::string &arrayName,
1072  svtkDataObject* mesh)
1073 {
1074  // get the particles from the particle container
1075  const auto& particles = this->m_particles->GetParticles();
1076 
1077  auto nptsOnProc = this->m_particles->TotalNumberOfParticles(true, true);
1078 
1079  // check that the name of the struct data is listed in the m_intStructs
1080  IntDataMapType::iterator ait = m_intStructs.find(arrayName);
1081  if (ait == m_intStructs.end())
1082  {
1083  SENSEI_ERROR("No int AOS named \"" << arrayName << "\"");
1084  return -1;
1085  }
1086 
1087  // get the index of the struct on the particles
1088  int index = ait->second;
1089 
1090  // check that the index are within the templated storage spec
1091  if(index >= ParticleType::NInt)
1092  {
1093  SENSEI_ERROR("Index out of bounds for int AOS named \"" << arrayName << "\"");
1094  return -1;
1095  }
1096 
1097  if constexpr(!ParticleType::is_soa_particle) {
1098 
1099  // allocate svtkArray
1100  svtkNew<svtkIntArray> data;
1101  data->SetName(arrayName.c_str());
1102  data->SetNumberOfComponents(1);
1103  data->SetNumberOfValues(nptsOnProc);
1104  int* pData = data->GetPointer(0);
1105 
1106  using MyParIter = ParIter_impl<ParticleType, NArrayReal, NArrayInt>;
1107  for (int level = 0; level<particles.size(); ++level)
1108  {
1109  for (MyParIter pti(*this->m_particles, level); pti.isValid(); ++pti)
1110  {
1111  const auto& aos = pti.GetArrayOfStructs();
1112 
1113  long long numReal = pti.numParticles();
1114  for (long long i = 0; i < numReal; ++i)
1115  {
1116  const auto &part = aos[i];
1117  if (part.id() > 0)
1118  {
1119  pData[i] = part.idata(index);
1120  }
1121  }
1122  pData += numReal;
1123  }
1124  }
1125 
1126  // get the block of the mesh
1127  int rank = 0;
1128  MPI_Comm_rank(this->GetCommunicator(), &rank);
1129 
1130  auto blocks = dynamic_cast<svtkMultiBlockDataSet*>(mesh);
1131 
1132  auto block = dynamic_cast<svtkPolyData*>(blocks->GetBlock(rank));
1133  block->GetPointData()->AddArray(data);
1134 
1135  }
1136 
1137  return 0;
1138 }
1139 
1140 }
double amrex_particle_real
Definition: AMReX_REAL.H:64
if(!(yy_init))
Definition: amrex_iparser.lex.nolint.H:935
A distributed container for Particles sorted onto the levels, grids, and tiles of a block-structured ...
Definition: AMReX_ParticleContainer.H:145
const Vector< ParticleLevel > & GetParticles() const
Return the underlying Vector (over AMR levels) of ParticleLevels. Const version.
Definition: AMReX_ParticleContainer.H:1004
@ max
Definition: AMReX_ParallelReduce.H:17
Definition: AMReX_Amr.cpp:49
Definition: AMReX_MakeParticle.H:16