VTK-m  2.0
CosmoToolsHaloFinder.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 // Copyright (c) 2016, Los Alamos National Security, LLC
11 // All rights reserved.
12 //
13 // Copyright 2016. Los Alamos National Security, LLC.
14 // This software was produced under U.S. Government contract DE-AC52-06NA25396
15 // for Los Alamos National Laboratory (LANL), which is operated by
16 // Los Alamos National Security, LLC for the U.S. Department of Energy.
17 // The U.S. Government has rights to use, reproduce, and distribute this
18 // software. NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC
19 // MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE
20 // USE OF THIS SOFTWARE. If software is modified to produce derivative works,
21 // such modified software should be clearly marked, so as not to confuse it
22 // with the version available from LANL.
23 //
24 // Additionally, redistribution and use in source and binary forms, with or
25 // without modification, are permitted provided that the following conditions
26 // are met:
27 //
28 // 1. Redistributions of source code must retain the above copyright notice,
29 // this list of conditions and the following disclaimer.
30 // 2. Redistributions in binary form must reproduce the above copyright notice,
31 // this list of conditions and the following disclaimer in the documentation
32 // and/or other materials provided with the distribution.
33 // 3. Neither the name of Los Alamos National Security, LLC, Los Alamos
34 // National Laboratory, LANL, the U.S. Government, nor the names of its
35 // contributors may be used to endorse or promote products derived from
36 // this software without specific prior written permission.
37 //
38 // THIS SOFTWARE IS PROVIDED BY LOS ALAMOS NATIONAL SECURITY, LLC AND
39 // CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
40 // BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
41 // FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL LOS ALAMOS
42 // NATIONAL SECURITY, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
43 // INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
44 // BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
45 // USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
46 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
47 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
48 // THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
49 //============================================================================
50 
51 #ifndef vtkm_worklet_cosmotools_cosmotools_halofinder_h
52 #define vtkm_worklet_cosmotools_cosmotools_halofinder_h
53 
55 
57 
58 namespace vtkm
59 {
60 namespace worklet
61 {
62 namespace cosmotools
63 {
64 
66 //
67 // Halo finder for all particles in domain
68 //
70 
71 template <typename T, typename StorageType>
74  vtkm::cont::ArrayHandle<T>& resultPot)
75 {
76  // Package locations for worklets
77  using CompositeLocationType =
79  CompositeLocationType location;
80  location = make_ArrayHandleCompositeVector(xLoc, yLoc, zLoc);
81 
82  vtkm::cont::ArrayHandle<vtkm::Id> leftNeighbor; // lower particle id to check for linking length
83  vtkm::cont::ArrayHandle<vtkm::Id> rightNeighbor; // upper particle id to check for linking length
85  activeMask; // mask per particle indicating active neighbor bins
86  vtkm::cont::ArrayHandle<vtkm::Id> partId; // index into all particles
87  vtkm::cont::ArrayHandle<vtkm::Id> binId; // bin id for each particle in each FOF halo
88 
89  leftNeighbor.Allocate(NUM_NEIGHBORS * nParticles);
90  rightNeighbor.Allocate(NUM_NEIGHBORS * nParticles);
91 
92  vtkm::cont::ArrayHandleConstant<bool> trueArray(true, nParticles);
93  vtkm::cont::ArrayHandleIndex indexArray(nParticles);
94 
95  // Bin all particles in domain into bins of size linking length
96  BinParticlesAll(partId, binId, leftNeighbor, rightNeighbor);
97 
98  // Mark active neighbor bins, meaning at least one particle in the bin
99  // is within linking length of the given particle indicated by mask
100  MarkActiveNeighbors<T> markActiveNeighbors(numBinsX, numBinsY, numBinsZ, NUM_NEIGHBORS, linkLen);
101  vtkm::worklet::DispatcherMapField<MarkActiveNeighbors<T>> markActiveNeighborsDispatcher(
102  markActiveNeighbors);
103  markActiveNeighborsDispatcher.Invoke(
104  indexArray, // (input) index into all particles
105  partId, // (input) particle id sorted by bin
106  binId, // (input) bin id sorted
107  partId, // (input) particle id (whole array)
108  location, // (input) location on original particle order
109  leftNeighbor, // (input) first partId for neighbor vector
110  rightNeighbor, // (input) last partId for neighbor vector
111  activeMask); // (output) mask per particle indicating valid neighbors
112 
113  // Initialize halo id of each particle to itself
114  vtkm::cont::ArrayHandle<vtkm::Id> haloIdCurrent;
116  DeviceAlgorithm::Copy(indexArray, haloIdCurrent);
117  DeviceAlgorithm::Copy(indexArray, haloIdLast);
118 
119  // rooted star is nchecked each iteration for all particles being rooted in a halo
121 
122  // Iterate over particles graft together to form halos
123  while (true)
124  {
125  // Connect each particle to another close particle to build halos
126  GraftParticles<T> graftParticles(numBinsX, numBinsY, numBinsZ, NUM_NEIGHBORS, linkLen);
127  vtkm::worklet::DispatcherMapField<GraftParticles<T>> graftParticlesDispatcher(graftParticles);
128 
129  graftParticlesDispatcher.Invoke(indexArray, // (input) index into particles
130  partId, // (input) particle id sorted by bin
131  binId, // (input) bin id sorted by bin
132  activeMask, // (input) flag indicates if neighor range is used
133  partId, // (input) particle id (whole array)
134  location, // (input) location on original particle order
135  leftNeighbor, // (input) first partId for neighbor
136  rightNeighbor, // (input) last partId for neighbor
137  haloIdCurrent); // (output)
138 #ifdef DEBUG_PRINT
139  DebugPrint("haloIdCurrent", haloIdCurrent);
140 #endif
141 
142  // Reininitialize rootedStar for each pass
143  DeviceAlgorithm::Copy(trueArray, rootedStar);
144 
145  // By comparing the haloIds from the last pass and this one
146  // determine if any particles are still migrating to halos
147  IsStar isStar;
148  vtkm::worklet::DispatcherMapField<IsStar> isStarDispatcher(isStar);
149  isStarDispatcher.Invoke(indexArray,
150  haloIdCurrent, // input (whole array)
151  haloIdLast, // input (whole array)
152  rootedStar); // output (whole array)
153 
154  // If all vertices are in rooted stars, algorithm is complete
155  bool allStars = DeviceAlgorithm::Reduce(rootedStar, true, vtkm::LogicalAnd());
156  if (allStars)
157  {
158  break;
159  }
160  else
161  // Otherwise copy current halo ids to last pass halo ids
162  {
163  PointerJump pointerJump;
164  vtkm::worklet::DispatcherMapField<PointerJump> pointerJumpDispatcher(pointerJump);
165  pointerJumpDispatcher.Invoke(indexArray, haloIdCurrent); // input (whole array)
166  DeviceAlgorithm::Copy(haloIdCurrent, haloIdLast);
167  }
168  }
169 
170  // Index into final halo id is the original particle ordering
171  // not the particles sorted by bin
172  DeviceAlgorithm::Copy(indexArray, partId);
173 #ifdef DEBUG_PRINT
174  DebugPrint("FINAL haloId", haloIdCurrent);
175  DebugPrint("FINAL partId", partId);
176 #endif
177 
178  // Call center finding on all halos using method with ReduceByKey and Scatter
179  DeviceAlgorithm::Copy(haloIdCurrent, resultHaloId);
180  MBPCenterFindingByHalo(partId, resultHaloId, resultMBP, resultPot);
181 }
182 
184 //
185 // Bin all particles in the system for halo finding
186 //
188 template <typename T, typename StorageType>
191  vtkm::cont::ArrayHandle<vtkm::Id>& leftNeighbor,
192  vtkm::cont::ArrayHandle<vtkm::Id>& rightNeighbor)
193 {
194  // Compute number of bins and ranges for each bin
195  vtkm::Vec<T, 2> result;
199  result = DeviceAlgorithm::Reduce(xLoc, xInit, vtkm::MinAndMax<T>());
200  T minX = result[0];
201  T maxX = result[1];
202  result = DeviceAlgorithm::Reduce(yLoc, yInit, vtkm::MinAndMax<T>());
203  T minY = result[0];
204  T maxY = result[1];
205  result = DeviceAlgorithm::Reduce(zLoc, zInit, vtkm::MinAndMax<T>());
206  T minZ = result[0];
207  T maxZ = result[1];
208 
209  vtkm::Id maxBins = 1048576;
210  vtkm::Id minBins = 1;
211 
212  numBinsX = static_cast<vtkm::Id>(vtkm::Floor((maxX - minX) / linkLen));
213  numBinsY = static_cast<vtkm::Id>(vtkm::Floor((maxY - minY) / linkLen));
214  numBinsZ = static_cast<vtkm::Id>(vtkm::Floor((maxZ - minZ) / linkLen));
215 
216  numBinsX = std::min(maxBins, numBinsX);
217  numBinsY = std::min(maxBins, numBinsY);
218  numBinsZ = std::min(maxBins, numBinsZ);
219 
220  numBinsX = std::max(minBins, numBinsX);
221  numBinsY = std::max(minBins, numBinsY);
222  numBinsZ = std::max(minBins, numBinsZ);
223 
224  // Compute which bin each particle is in
225  ComputeBins<T> computeBins(minX,
226  maxX, // Physical range on domain
227  minY,
228  maxY,
229  minZ,
230  maxZ,
231  numBinsX,
232  numBinsY,
233  numBinsZ); // Size of superimposed mesh
234  vtkm::worklet::DispatcherMapField<ComputeBins<T>> computeBinsDispatcher(computeBins);
235  computeBinsDispatcher.Invoke(xLoc, // input
236  yLoc, // input
237  zLoc, // input
238  binId); // output
239 
240  vtkm::cont::ArrayHandleIndex indexArray(nParticles);
241  DeviceAlgorithm::Copy(indexArray, partId);
242 
243 #ifdef DEBUG_PRINT
244  std::cout << std::endl
245  << "** BinParticlesAll (" << numBinsX << ", " << numBinsY << ", " << numBinsZ << ")"
246  << std::endl;
247  DebugPrint("xLoc", xLoc);
248  DebugPrint("yLoc", yLoc);
249  DebugPrint("zLoc", zLoc);
250  DebugPrint("partId", partId);
251  DebugPrint("binId", binId);
252  std::cout << std::endl;
253 #endif
254 
255  // Sort the particles by bin (remember that xLoc and yLoc are not sorted)
256  DeviceAlgorithm::SortByKey(binId, partId);
257 #ifdef DEBUG_PRINT
258  DebugPrint("partId", partId);
259  DebugPrint("binId", binId);
260 #endif
261 
262  // Compute indices of all left neighbor bins
263  vtkm::cont::ArrayHandleIndex countArray(nParticles);
264  ComputeNeighborBins computeNeighborBins(numBinsX, numBinsY, numBinsZ, NUM_NEIGHBORS);
265  vtkm::worklet::DispatcherMapField<ComputeNeighborBins> computeNeighborBinsDispatcher(
266  computeNeighborBins);
267  computeNeighborBinsDispatcher.Invoke(countArray, binId, leftNeighbor);
268 
269  // Compute indices of all right neighbor bins
270  ComputeBinRange computeBinRange(numBinsX);
271  vtkm::worklet::DispatcherMapField<ComputeBinRange> computeBinRangeDispatcher(computeBinRange);
272  computeBinRangeDispatcher.Invoke(leftNeighbor, rightNeighbor);
273 
274  // Convert bin range to particle range within the bins
275  DeviceAlgorithm::LowerBounds(binId, leftNeighbor, leftNeighbor);
276  DeviceAlgorithm::UpperBounds(binId, rightNeighbor, rightNeighbor);
277 }
278 
280 //
281 // Center finder for all particles given location, particle id and halo id
282 // MBP (Most Bound Particle) is particle with the minimum potential energy
283 // Method uses ReduceByKey() and Scatter()
284 //
286 template <typename T, typename StorageType>
290  vtkm::cont::ArrayHandle<T>& minPotential)
291 {
292  // Sort particles into groups according to halo id using an index into WholeArrays
293  DeviceAlgorithm::SortByKey(haloId, partId);
294 #ifdef DEBUG_PRINT
295  DebugPrint("Sorted haloId", haloId);
296  DebugPrint("Sorted partId", partId);
297 #endif
298 
299  // Find the particle in each halo with the lowest potential
300  // Compute starting and ending indices of each halo
301  vtkm::cont::ArrayHandleConstant<vtkm::Id> constArray(1, nParticles);
302  vtkm::cont::ArrayHandleIndex indexArray(nParticles);
303  vtkm::cont::ArrayHandle<vtkm::Id> uniqueHaloIds;
304  vtkm::cont::ArrayHandle<vtkm::Id> particlesPerHalo;
307  vtkm::cont::ArrayHandle<T> potential;
310 
311  // Halo ids have been sorted, reduce to find the number of particles per halo
312  DeviceAlgorithm::ReduceByKey(haloId, constArray, uniqueHaloIds, particlesPerHalo, vtkm::Add());
313 #ifdef DEBUG_PRINT
314  DebugPrint("uniqueHaloId", uniqueHaloIds);
315  DebugPrint("partPerHalo", particlesPerHalo);
316  std::cout << std::endl;
317 #endif
318 
319  // Setup the ScatterCounting worklets needed to expand the ReduceByKeyResults
320  vtkm::worklet::ScatterCounting scatter(particlesPerHalo);
321  vtkm::cont::Invoker invoke;
322 
323  // Calculate the minimum particle index per halo id and scatter
324  DeviceAlgorithm::ScanExclusive(particlesPerHalo, tempI);
325  invoke(ScatterWorklet<vtkm::Id>{}, scatter, tempI, minParticle);
326 
327  // Calculate the maximum particle index per halo id and scatter
328  DeviceAlgorithm::ScanInclusive(particlesPerHalo, tempI);
329  invoke(ScatterWorklet<vtkm::Id>{}, scatter, tempI, maxParticle);
330 
333  vtkm::cont::make_ArrayHandleTransform<IdArrayType>(maxParticle,
335 
336  DeviceAlgorithm::Copy(scaleBias, maxParticle);
337 #ifdef DEBUG_PRINT
338  DebugPrint("minParticle", minParticle);
339  DebugPrint("maxParticle", maxParticle);
340 #endif
341 
342  // Compute potentials
343  ComputePotential<T> computePotential(particleMass);
344  vtkm::worklet::DispatcherMapField<ComputePotential<T>> computePotentialDispatcher(
345  computePotential);
346 
347  computePotentialDispatcher.Invoke(indexArray,
348  partId, // input (whole array)
349  xLoc, // input (whole array)
350  yLoc, // input (whole array)
351  zLoc, // input (whole array)
352  minParticle, // input (whole array)
353  maxParticle, // input (whole array)
354  potential); // output
355 
356  // Find minimum potential for all particles in a halo and scatter
357  DeviceAlgorithm::ReduceByKey(haloId, potential, uniqueHaloIds, tempT, vtkm::Minimum());
358  invoke(ScatterWorklet<T>{}, scatter, tempT, minPotential);
359 #ifdef DEBUG_PRINT
360  DebugPrint("potential", potential);
361  DebugPrint("minPotential", minPotential);
362 #endif
363 
364  // Find the particle id matching the minimum potential (Worklet)
365  EqualsMinimumPotential<T> equalsMinimumPotential;
366  vtkm::worklet::DispatcherMapField<EqualsMinimumPotential<T>> equalsMinimumPotentialDispatcher(
367  equalsMinimumPotential);
368 
369  equalsMinimumPotentialDispatcher.Invoke(partId, potential, minPotential, mbpId);
370 
371  // Fill out entire array with center index, another reduce and scatter
373  minIndx.Allocate(nParticles);
374  DeviceAlgorithm::ReduceByKey(haloId, mbpId, uniqueHaloIds, minIndx, vtkm::Maximum());
375  invoke(ScatterWorklet<vtkm::Id>{}, scatter, minIndx, mbpId);
376 
377  // Resort particle ids and mbpId to starting order
379  DeviceAlgorithm::Copy(partId, savePartId);
380 
381  DeviceAlgorithm::SortByKey(partId, haloId);
382  DeviceAlgorithm::Copy(savePartId, partId);
383  DeviceAlgorithm::SortByKey(partId, mbpId);
384  DeviceAlgorithm::Copy(savePartId, partId);
385  DeviceAlgorithm::SortByKey(partId, minPotential);
386 
387 #ifdef DEBUG_PRINT
388  std::cout << std::endl;
389  DebugPrint("partId", partId);
390  DebugPrint("xLoc", xLoc);
391  DebugPrint("yLoc", yLoc);
392  DebugPrint("haloId", haloId);
393  DebugPrint("mbpId", mbpId);
394  DebugPrint("minPotential", minPotential);
395 #endif
396 }
397 }
398 }
399 }
400 #endif
vtkm::cont::ArrayHandle< vtkm::Id >
vtkm
Groups connected points that have the same field value.
Definition: Atomic.h:19
vtkm::MinAndMax
Binary Predicate that takes two arguments argument x, and y and returns a vtkm::Vec<T,...
Definition: BinaryOperators.h:112
vtkm::cont::ArrayHandle::Allocate
VTKM_CONT void Allocate(vtkm::Id numberOfValues, vtkm::CopyFlag preserve, vtkm::cont::Token &token) const
Allocates an array large enough to hold the given number of values.
Definition: ArrayHandle.h:465
vtkm::worklet::cosmotools::ComputeBins
Definition: ComputeBins.h:65
vtkm::worklet::cosmotools::ComputeBinRange
Definition: ComputeBinRange.h:65
vtkm::worklet::cosmotools::IsStar
Definition: IsStar.h:65
vtkm::cont::ArrayHandleCompositeVector
An ArrayHandle that combines components from other arrays.
Definition: ArrayHandleCompositeVector.h:402
vtkm::Maximum
Binary Predicate that takes two arguments argument x, and y and returns the x if x > y otherwise retu...
Definition: BinaryOperators.h:85
vtkm::cont::ArrayGetValue
VTKM_CONT T ArrayGetValue(vtkm::Id id, const vtkm::cont::ArrayHandle< T, S > &data)
Obtain a small set of values from an ArrayHandle with minimal device transfers.
Definition: ArrayGetValues.h:264
vtkm::worklet::cosmotools::EqualsMinimumPotential
Definition: EqualsMinimumPotential.h:65
vtkm::Id
vtkm::Int32 Id
Represents an ID (index into arrays).
Definition: Types.h:191
vtkm::LogicalAnd
Binary Predicate that takes two arguments argument x, and y and returns True if and only if x and y a...
Definition: BinaryPredicates.h:71
vtkm::worklet::cosmotools::ComputeNeighborBins
Definition: ComputeNeighborBins.h:65
vtkm::Floor
VTKM_EXEC_CONT vtkm::Float32 Floor(vtkm::Float32 x)
Round x to the largest integer value not greater than x.
Definition: Math.h:2204
vtkm::Add
Definition: Types.h:222
vtkm::worklet::cosmotools::ScatterWorklet
Definition: cosmotools/CosmoTools.h:139
vtkm::cont::make_ArrayHandleCompositeVector
VTKM_CONT ArrayHandleCompositeVector< ArrayTs... > make_ArrayHandleCompositeVector(const ArrayTs &... arrays)
Create a composite vector array from other arrays.
Definition: ArrayHandleCompositeVector.h:430
vtkm::worklet::cosmotools::MarkActiveNeighbors
Definition: MarkActiveNeighbors.h:66
vtkm::worklet::ScatterCounting
A scatter that maps input to some numbers of output.
Definition: ScatterCounting.h:44
vtkm::worklet::cosmotools::PointerJump
Definition: PointerJump.h:64
vtkm::worklet::DispatcherMapField
Dispatcher for worklets that inherit from WorkletMapField.
Definition: DispatcherMapField.h:25
vtkm::cont::Invoker
Allows launching any worklet without a dispatcher.
Definition: Invoker.h:41
vtkm::Vec< T, 2 >
Definition: Types.h:859
vtkm::worklet::contourtree_augmented::IdArrayType
vtkm::cont::ArrayHandle< vtkm::Id > IdArrayType
Definition: filter/scalar_topology/worklet/contourtree_augmented/Types.h:90
CosmoTools.h
vtkm::cont::ArrayHandleTransform
Implicitly transform values of one array to another with a functor.
Definition: ArrayHandleTransform.h:437
vtkm::cont::ArrayHandleConstant
An array handle with a constant value.
Definition: ArrayHandleConstant.h:63
vtkm::worklet::cosmotools::GraftParticles
Definition: GraftParticles.h:66
ArrayGetValues.h
vtkm::worklet::cosmotools::CosmoTools::BinParticlesAll
void BinParticlesAll(vtkm::cont::ArrayHandle< vtkm::Id > &partId, vtkm::cont::ArrayHandle< vtkm::Id > &binId, vtkm::cont::ArrayHandle< vtkm::Id > &leftNeighbor, vtkm::cont::ArrayHandle< vtkm::Id > &rightNeighbor)
Definition: CosmoToolsHaloFinder.h:189
vtkm::worklet::cosmotools::ScaleBiasFunctor
Definition: cosmotools/CosmoTools.h:155
vtkm::worklet::cosmotools::CosmoTools::HaloFinder
void HaloFinder(vtkm::cont::ArrayHandle< vtkm::Id > &resultHaloId, vtkm::cont::ArrayHandle< vtkm::Id > &resultMBP, vtkm::cont::ArrayHandle< T > &resultPot)
Definition: CosmoToolsHaloFinder.h:72
vtkm::worklet::cosmotools::ComputePotential
Definition: ComputePotential.h:65
vtkm::worklet::cosmotools::CosmoTools::MBPCenterFindingByHalo
void MBPCenterFindingByHalo(vtkm::cont::ArrayHandle< vtkm::Id > &partId, vtkm::cont::ArrayHandle< vtkm::Id > &haloId, vtkm::cont::ArrayHandle< vtkm::Id > &mbpId, vtkm::cont::ArrayHandle< T > &minPotential)
Definition: CosmoToolsHaloFinder.h:287
vtkm::cont::ArrayHandleIndex
An implicit array handle containing the its own indices.
Definition: ArrayHandleIndex.h:54
vtkm::Minimum
Binary Predicate that takes two arguments argument x, and y and returns the x if x < y otherwise retu...
Definition: BinaryOperators.h:99