VTK-m  2.2
BoundsMap.h
Go to the documentation of this file.
1 //============================================================================
2 // Copyright (c) Kitware, Inc.
3 // All rights reserved.
4 // See LICENSE.txt for details.
5 //
6 // This software is distributed WITHOUT ANY WARRANTY; without even
7 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
8 // PURPOSE. See the above copyright notice for more information.
9 //============================================================================
10 
11 #ifndef vtk_m_filter_flow_internal_BoundsMap_h
12 #define vtk_m_filter_flow_internal_BoundsMap_h
13 
14 #include <vtkm/Bounds.h>
16 #include <vtkm/cont/DataSet.h>
19 #include <vtkm/cont/Field.h>
21 
22 #include <vtkm/thirdparty/diy/diy.h>
23 
24 #ifdef VTKM_ENABLE_MPI
25 #include <mpi.h>
26 #include <vtkm/thirdparty/diy/mpi-cast.h>
27 #endif
28 
29 #include <algorithm>
30 #include <iostream>
31 #include <set>
32 #include <vector>
33 
34 namespace vtkm
35 {
36 namespace filter
37 {
38 namespace flow
39 {
40 namespace internal
41 {
42 
43 class VTKM_ALWAYS_EXPORT BoundsMap
44 {
45 public:
46  BoundsMap() {}
47 
48  BoundsMap(const vtkm::cont::DataSet& dataSet) { this->Init({ dataSet }); }
49 
50  BoundsMap(const vtkm::cont::DataSet& dataSet, const vtkm::Id& blockId)
51  {
52  this->Init({ dataSet }, { blockId });
53  }
54 
55  BoundsMap(const std::vector<vtkm::cont::DataSet>& dataSets) { this->Init(dataSets); }
56 
57  BoundsMap(const vtkm::cont::PartitionedDataSet& pds) { this->Init(pds.GetPartitions()); }
58 
59  BoundsMap(const vtkm::cont::PartitionedDataSet& pds, const std::vector<vtkm::Id>& blockIds)
60  {
61  this->Init(pds.GetPartitions(), blockIds);
62  }
63 
64  vtkm::Bounds GetGlobalBounds() const { return this->GlobalBounds; }
65 
66  vtkm::Bounds GetBlockBounds(vtkm::Id idx) const
67  {
68  VTKM_ASSERT(idx >= 0 && static_cast<std::size_t>(idx) < this->BlockBounds.size());
69 
70  return this->BlockBounds[static_cast<std::size_t>(idx)];
71  }
72 
73  vtkm::Id GetLocalBlockId(vtkm::Id idx) const
74  {
75  VTKM_ASSERT(idx >= 0 && idx < this->LocalNumBlocks);
76  return this->LocalIDs[static_cast<std::size_t>(idx)];
77  }
78 
79  std::vector<int> FindRank(vtkm::Id blockId) const
80  {
81  auto it = this->BlockToRankMap.find(blockId);
82  if (it == this->BlockToRankMap.end())
83  return {};
84 
85  return it->second;
86  }
87 
88  std::vector<vtkm::Id> FindBlocks(const vtkm::Vec3f& p) const { return this->FindBlocks(p, -1); }
89 
90  std::vector<vtkm::Id> FindBlocks(const vtkm::Vec3f& p,
91  const std::vector<vtkm::Id>& ignoreBlocks) const
92  {
93  vtkm::Id ignoreID = (ignoreBlocks.empty() ? -1 : ignoreBlocks[0]);
94  return FindBlocks(p, ignoreID);
95  }
96 
97  std::vector<vtkm::Id> FindBlocks(const vtkm::Vec3f& p, vtkm::Id ignoreBlock) const
98  {
99  std::vector<vtkm::Id> blockIDs;
100  if (this->GlobalBounds.Contains(p))
101  {
102  vtkm::Id blockId = 0;
103  for (auto& it : this->BlockBounds)
104  {
105  if (blockId != ignoreBlock && it.Contains(p))
106  blockIDs.emplace_back(blockId);
107  blockId++;
108  }
109  }
110 
111  return blockIDs;
112  }
113 
114  vtkm::Id GetTotalNumBlocks() const { return this->TotalNumBlocks; }
115  vtkm::Id GetLocalNumBlocks() const { return this->LocalNumBlocks; }
116 
117 private:
118  void Init(const std::vector<vtkm::cont::DataSet>& dataSets, const std::vector<vtkm::Id>& blockIds)
119  {
120  if (dataSets.size() != blockIds.size())
121  throw vtkm::cont::ErrorFilterExecution("Number of datasets and block ids must match");
122 
123  this->LocalIDs = blockIds;
124  this->LocalNumBlocks = dataSets.size();
125 
126  vtkmdiy::mpi::communicator comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
127 
128  //1. Get the min/max blockId
129  vtkm::Id locMinId = 0, locMaxId = 0;
130  if (!this->LocalIDs.empty())
131  {
132  locMinId = *std::min_element(this->LocalIDs.begin(), this->LocalIDs.end());
133  locMaxId = *std::max_element(this->LocalIDs.begin(), this->LocalIDs.end());
134  }
135 
136  vtkm::Id globalMinId = 0, globalMaxId = 0;
137 
138  vtkmdiy::mpi::all_reduce(comm, locMinId, globalMinId, vtkmdiy::mpi::minimum<vtkm::Id>{});
139  vtkmdiy::mpi::all_reduce(comm, locMaxId, globalMaxId, vtkmdiy::mpi::maximum<vtkm::Id>{});
140  if (globalMinId != 0 || (globalMaxId - globalMinId) < 1)
141  throw vtkm::cont::ErrorFilterExecution("Invalid block ids");
142 
143  //2. Find out how many blocks everyone has.
144  std::vector<vtkm::Id> locBlockCounts(comm.size(), 0), globalBlockCounts(comm.size(), 0);
145  locBlockCounts[comm.rank()] = this->LocalIDs.size();
146  vtkmdiy::mpi::all_reduce(comm, locBlockCounts, globalBlockCounts, std::plus<vtkm::Id>{});
147 
148  //note: there might be duplicates...
149  vtkm::Id globalNumBlocks =
150  std::accumulate(globalBlockCounts.begin(), globalBlockCounts.end(), vtkm::Id{ 0 });
151 
152  //3. given the counts per rank, calc offset for this rank.
153  vtkm::Id offset = 0;
154  for (int i = 0; i < comm.rank(); i++)
155  offset += globalBlockCounts[i];
156 
157  //4. calc the blocks on each rank.
158  std::vector<vtkm::Id> localBlockIds(globalNumBlocks, 0);
159  vtkm::Id idx = 0;
160  for (const auto& bid : this->LocalIDs)
161  localBlockIds[offset + idx++] = bid;
162 
163  //use an MPI_Alltoallv instead.
164  std::vector<vtkm::Id> globalBlockIds(globalNumBlocks, 0);
165  vtkmdiy::mpi::all_reduce(comm, localBlockIds, globalBlockIds, std::plus<vtkm::Id>{});
166 
167 
168  //5. create a rank -> blockId map.
169  // rankToBlockIds[rank] = {this->LocalIDs on rank}.
170  std::vector<std::vector<vtkm::Id>> rankToBlockIds(comm.size());
171 
172  offset = 0;
173  for (int rank = 0; rank < comm.size(); rank++)
174  {
175  vtkm::Id numBIds = globalBlockCounts[rank];
176  rankToBlockIds[rank].resize(numBIds);
177  for (vtkm::Id i = 0; i < numBIds; i++)
178  rankToBlockIds[rank][i] = globalBlockIds[offset + i];
179 
180  offset += numBIds;
181  }
182 
183  //6. there might be duplicates, so count number of UNIQUE blocks.
184  std::set<vtkm::Id> globalUniqueBlockIds;
185  globalUniqueBlockIds.insert(globalBlockIds.begin(), globalBlockIds.end());
186  this->TotalNumBlocks = globalUniqueBlockIds.size();
187 
188  //Build a vector of : blockIdsToRank[blockId] = {ranks that have this blockId}
189  std::vector<std::vector<vtkm::Id>> blockIdsToRank(this->TotalNumBlocks);
190  for (int rank = 0; rank < comm.size(); rank++)
191  {
192  for (const auto& bid : rankToBlockIds[rank])
193  {
194  blockIdsToRank[bid].push_back(rank);
195  this->BlockToRankMap[bid].push_back(rank);
196  }
197  }
198 
199  this->Build(dataSets);
200  }
201 
202  void Init(const std::vector<vtkm::cont::DataSet>& dataSets)
203  {
204  this->LocalNumBlocks = dataSets.size();
205 
206  vtkm::cont::AssignerPartitionedDataSet assigner(this->LocalNumBlocks);
207  this->TotalNumBlocks = assigner.nblocks();
208  std::vector<int> ids;
209 
210  vtkmdiy::mpi::communicator Comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
211  assigner.local_gids(Comm.rank(), ids);
212  for (const auto& i : ids)
213  this->LocalIDs.emplace_back(static_cast<vtkm::Id>(i));
214 
215  for (vtkm::Id id = 0; id < this->TotalNumBlocks; id++)
216  this->BlockToRankMap[id] = { assigner.rank(static_cast<int>(id)) };
217  this->Build(dataSets);
218  }
219 
220  void Build(const std::vector<vtkm::cont::DataSet>& dataSets)
221  {
222  std::vector<vtkm::Float64> vals(static_cast<std::size_t>(this->TotalNumBlocks * 6), 0);
223  std::vector<vtkm::Float64> vals2(vals.size());
224 
225  std::vector<vtkm::Float64> localMins((this->TotalNumBlocks * 3),
226  std::numeric_limits<vtkm::Float64>::max());
227  std::vector<vtkm::Float64> localMaxs((this->TotalNumBlocks * 3),
228  -std::numeric_limits<vtkm::Float64>::max());
229 
230  for (std::size_t i = 0; i < this->LocalIDs.size(); i++)
231  {
232  const vtkm::cont::DataSet& ds = dataSets[static_cast<std::size_t>(i)];
234 
235  vtkm::Id localID = this->LocalIDs[i];
236  localMins[localID * 3 + 0] = bounds.X.Min;
237  localMins[localID * 3 + 1] = bounds.Y.Min;
238  localMins[localID * 3 + 2] = bounds.Z.Min;
239  localMaxs[localID * 3 + 0] = bounds.X.Max;
240  localMaxs[localID * 3 + 1] = bounds.Y.Max;
241  localMaxs[localID * 3 + 2] = bounds.Z.Max;
242  }
243 
244  std::vector<vtkm::Float64> globalMins, globalMaxs;
245 
246 #ifdef VTKM_ENABLE_MPI
247  globalMins.resize(this->TotalNumBlocks * 3);
248  globalMaxs.resize(this->TotalNumBlocks * 3);
249 
250  vtkmdiy::mpi::communicator comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
251 
252  vtkmdiy::mpi::all_reduce(comm, localMins, globalMins, vtkmdiy::mpi::minimum<vtkm::Float64>{});
253  vtkmdiy::mpi::all_reduce(comm, localMaxs, globalMaxs, vtkmdiy::mpi::maximum<vtkm::Float64>{});
254 #else
255  globalMins = localMins;
256  globalMaxs = localMaxs;
257 #endif
258 
259  this->BlockBounds.resize(static_cast<std::size_t>(this->TotalNumBlocks));
260  this->GlobalBounds = vtkm::Bounds();
261 
262  std::size_t idx = 0;
263  for (auto& block : this->BlockBounds)
264  {
265  block = vtkm::Bounds(globalMins[idx + 0],
266  globalMaxs[idx + 0],
267  globalMins[idx + 1],
268  globalMaxs[idx + 1],
269  globalMins[idx + 2],
270  globalMaxs[idx + 2]);
271  this->GlobalBounds.Include(block);
272  idx += 3;
273  }
274  }
275 
276  vtkm::Id LocalNumBlocks = 0;
277  std::vector<vtkm::Id> LocalIDs;
278  std::map<vtkm::Id, std::vector<vtkm::Int32>> BlockToRankMap;
279  vtkm::Id TotalNumBlocks = 0;
280  std::vector<vtkm::Bounds> BlockBounds;
281  vtkm::Bounds GlobalBounds;
282 };
283 
284 }
285 }
286 }
287 } // namespace vtkm::filter::flow::internal
288 
289 #endif //vtk_m_filter_flow_internal_BoundsMap_h
vtkm::cont::PartitionedDataSet::GetPartitions
const std::vector< vtkm::cont::DataSet > & GetPartitions() const
Get an STL vector of all DataSet objects stored in PartitionedDataSet.
AssignerPartitionedDataSet.h
vtkm
Groups connected points that have the same field value.
Definition: Atomic.h:19
VTKM_ASSERT
#define VTKM_ASSERT(condition)
Definition: Assert.h:43
vtkm::cont::AssignerPartitionedDataSet
Assigner for PartitionedDataSet partitions.
Definition: AssignerPartitionedDataSet.h:49
vtkm::cont::ErrorFilterExecution
This class is primarily intended to filters to throw in the control environment to indicate an execut...
Definition: ErrorFilterExecution.h:27
vtkm::cont::DataSet
Contains and manages the geometric data structures that VTK-m operates on.
Definition: DataSet.h:57
EnvironmentTracker.h
ErrorFilterExecution.h
Bounds.h
Field.h
vtkm::Id
vtkm::Int64 Id
Base type to use to index arrays.
Definition: Types.h:227
vtkm::Bounds
Represent an axis-aligned 3D bounds in space.
Definition: Bounds.h:29
vtkm::cont::EnvironmentTracker::GetCommunicator
static const vtkmdiy::mpi::communicator & GetCommunicator()
vtkm::Vec< vtkm::FloatDefault, 3 >
vtkm::Range::Min
vtkm::Float64 Min
The minumum value of the range (inclusive).
Definition: Range.h:34
PartitionedDataSet.h
vtkm::cont::DataSet::GetCoordinateSystem
vtkm::cont::CoordinateSystem GetCoordinateSystem(vtkm::Id index=0) const
vtkm::Bounds::X
vtkm::Range X
The range of values in the X direction.
Definition: Bounds.h:33
vtkm::cont::CoordinateSystem::GetBounds
vtkm::Bounds GetBounds() const
Definition: CoordinateSystem.h:128
VTKM_ALWAYS_EXPORT
#define VTKM_ALWAYS_EXPORT
Definition: ExportMacros.h:89
DataSet.h
vtkm::cont::PartitionedDataSet
Comprises a set of vtkm::cont::DataSet objects.
Definition: PartitionedDataSet.h:26