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_centerfinder_h
52 #define vtkm_worklet_cosmotools_cosmotools_centerfinder_h
53 
54 #include <vtkm/worklet/cosmotools/CosmoTools.h>
55 
56 #include <vtkm/cont/ArrayGetValues.h>
57 
58 namespace vtkm
59 {
60 namespace worklet
61 {
62 namespace cosmotools
63 {
64 
65 ///////////////////////////////////////////////////////////////////////////////
66 //
67 // Center finder for particles in FOF halo using estimations but with exact final answer
68 // MBP (Most Bound Particle) is particle with the minimum potential energy
69 //
70 ///////////////////////////////////////////////////////////////////////////////
71 
72 template <typename T, typename StorageType>
MBPCenterFinderMxN(T * mxnPotential)73 vtkm::Id CosmoTools<T, StorageType>::MBPCenterFinderMxN(T* mxnPotential)
74 {
75   vtkm::cont::ArrayHandle<vtkm::Id> partId;
76   vtkm::cont::ArrayHandle<vtkm::Id> binId;
77 
78   vtkm::cont::ArrayHandle<vtkm::Id> uniqueBins;
79   vtkm::cont::ArrayHandle<vtkm::Id> partPerBin;
80   vtkm::cont::ArrayHandle<vtkm::Id> particleOffset;
81 
82   vtkm::cont::ArrayHandle<vtkm::Id> binX;
83   vtkm::cont::ArrayHandle<vtkm::Id> binY;
84   vtkm::cont::ArrayHandle<vtkm::Id> binZ;
85 
86   // Bin all particles in the halo into bins of size linking length
87   BinParticlesHalo(partId, binId, uniqueBins, partPerBin, particleOffset, binX, binY, binZ);
88 #ifdef DEBUG_PRINT
89   DebugPrint("uniqueBins", uniqueBins);
90   DebugPrint("partPerBin", partPerBin);
91 #endif
92 
93   // Compute the estimated potential per bin using 27 contiguous bins
94   vtkm::cont::ArrayHandle<T> partPotential;
95   MBPCenterFindingByKey(binId, partId, partPotential);
96 
97   // Reduce by key to get the estimated minimum potential per bin within 27 neighbors
98   vtkm::cont::ArrayHandle<vtkm::Id> tempId;
99   vtkm::cont::ArrayHandle<T> minPotential;
100   DeviceAlgorithm::ReduceByKey(binId, partPotential, tempId, minPotential, vtkm::Minimum());
101 
102   // Reduce by key to get the estimated maximum potential per bin within 27 neighbors
103   vtkm::cont::ArrayHandle<T> maxPotential;
104   DeviceAlgorithm::ReduceByKey(binId, partPotential, tempId, maxPotential, vtkm::Maximum());
105 #ifdef DEBUG_PRINT
106   DebugPrint("minPotential", minPotential);
107   DebugPrint("maxPotential", maxPotential);
108 #endif
109 
110   // Compute potentials estimate for a bin using all other bins
111   // Particles in the other bins are located at the closest point to this bin
112   vtkm::cont::ArrayHandleIndex uniqueIndex(uniqueBins.GetNumberOfValues());
113   vtkm::cont::ArrayHandle<T> bestEstPotential;
114   vtkm::cont::ArrayHandle<T> worstEstPotential;
115 
116   // Initialize each bin potential with the nxn for that bin
117   DeviceAlgorithm::Copy(minPotential, bestEstPotential);
118   DeviceAlgorithm::Copy(maxPotential, worstEstPotential);
119 
120   // Estimate only across the uniqueBins that contain particles
121   ComputePotentialBin<T> computePotentialBin(uniqueBins.GetNumberOfValues(), particleMass, linkLen);
122   vtkm::worklet::DispatcherMapField<ComputePotentialBin<T>> computePotentialBinDispatcher(
123     computePotentialBin);
124 
125   computePotentialBinDispatcher.Invoke(uniqueIndex,        // input
126                                        partPerBin,         // input (whole array)
127                                        binX,               // input (whole array)
128                                        binY,               // input (whole array)
129                                        binZ,               // input (whole array)
130                                        bestEstPotential,   // input/output
131                                        worstEstPotential); // input/output
132 #ifdef DEBUG_PRINT
133   DebugPrint("bestEstPotential", bestEstPotential);
134   DebugPrint("worstEstPotential", worstEstPotential);
135   std::cout << "Number of bestEstPotential " << bestEstPotential.GetNumberOfValues() << std::endl;
136   std::cout << "Number of worstEstPotential " << worstEstPotential.GetNumberOfValues() << std::endl;
137 #endif
138 
139   // Sort everything by the best estimated potential per bin
140   vtkm::cont::ArrayHandle<T> tempBest;
141   DeviceAlgorithm::Copy(bestEstPotential, tempBest);
142   DeviceAlgorithm::SortByKey(tempBest, worstEstPotential);
143 
144   // Use the worst estimate for the first selected bin to compare to best of all others
145   // Any bin that passes is a candidate for having the MBP
146   T cutoffPotential = vtkm::cont::ArrayGetValue(0, worstEstPotential);
147   worstEstPotential.ReleaseResources();
148   tempBest.ReleaseResources();
149 
150   vtkm::cont::ArrayHandle<vtkm::Id> candidate;
151   DeviceAlgorithm::Copy(vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, nParticles), candidate);
152 
153   SetCandidateParticles<T> setCandidateParticles(cutoffPotential);
154   vtkm::worklet::DispatcherMapField<SetCandidateParticles<T>> setCandidateParticlesDispatcher(
155     setCandidateParticles);
156   setCandidateParticlesDispatcher.Invoke(bestEstPotential, // input
157                                          particleOffset,   // input
158                                          partPerBin,       // input
159                                          candidate);       // output (whole array)
160 
161   // Copy the M candidate particles to a new array
162   vtkm::cont::ArrayHandle<vtkm::Id> mparticles;
163   DeviceAlgorithm::CopyIf(partId, candidate, mparticles);
164 
165   // Compute potentials only on the candidate particles
166   vtkm::cont::ArrayHandle<T> mpotential;
167   ComputePotentialOnCandidates<T> computePotentialOnCandidates(nParticles, particleMass);
168   vtkm::worklet::DispatcherMapField<ComputePotentialOnCandidates<T>>
169     computePotentialOnCandidatesDispatcher(computePotentialOnCandidates);
170 
171   computePotentialOnCandidatesDispatcher.Invoke(mparticles,
172                                                 xLoc,        // input (whole array)
173                                                 yLoc,        // input (whole array)
174                                                 zLoc,        // input (whole array)
175                                                 mpotential); // output
176 
177   // Of the M candidate particles which has the minimum potential
178   DeviceAlgorithm::SortByKey(mpotential, mparticles);
179 #ifdef DEBUG_PRINT
180   DebugPrint("mparticles", mparticles);
181   DebugPrint("mpotential", mpotential);
182 #endif
183 
184   // Return the found MBP particle and its potential
185   vtkm::Id mxnMBP = vtkm::cont::ArrayGetValue(0, mparticles);
186   *mxnPotential = vtkm::cont::ArrayGetValue(0, mpotential);
187 
188   return mxnMBP;
189 }
190 
191 ///////////////////////////////////////////////////////////////////////////////
192 //
193 // Bin particles in one halo for quick MBP finding
194 //
195 ///////////////////////////////////////////////////////////////////////////////
196 template <typename T, typename StorageType>
BinParticlesHalo(vtkm::cont::ArrayHandle<vtkm::Id> & partId,vtkm::cont::ArrayHandle<vtkm::Id> & binId,vtkm::cont::ArrayHandle<vtkm::Id> & uniqueBins,vtkm::cont::ArrayHandle<vtkm::Id> & partPerBin,vtkm::cont::ArrayHandle<vtkm::Id> & particleOffset,vtkm::cont::ArrayHandle<vtkm::Id> & binX,vtkm::cont::ArrayHandle<vtkm::Id> & binY,vtkm::cont::ArrayHandle<vtkm::Id> & binZ)197 void CosmoTools<T, StorageType>::BinParticlesHalo(vtkm::cont::ArrayHandle<vtkm::Id>& partId,
198                                                   vtkm::cont::ArrayHandle<vtkm::Id>& binId,
199                                                   vtkm::cont::ArrayHandle<vtkm::Id>& uniqueBins,
200                                                   vtkm::cont::ArrayHandle<vtkm::Id>& partPerBin,
201                                                   vtkm::cont::ArrayHandle<vtkm::Id>& particleOffset,
202                                                   vtkm::cont::ArrayHandle<vtkm::Id>& binX,
203                                                   vtkm::cont::ArrayHandle<vtkm::Id>& binY,
204                                                   vtkm::cont::ArrayHandle<vtkm::Id>& binZ)
205 {
206   // Compute number of bins and ranges for each bin
207   vtkm::Vec<T, 2> xRange(vtkm::cont::ArrayGetValue(0, xLoc));
208   vtkm::Vec<T, 2> yRange(vtkm::cont::ArrayGetValue(0, yLoc));
209   vtkm::Vec<T, 2> zRange(vtkm::cont::ArrayGetValue(0, zLoc));
210   xRange = DeviceAlgorithm::Reduce(xLoc, xRange, vtkm::MinAndMax<T>());
211   T minX = xRange[0];
212   T maxX = xRange[1];
213   yRange = DeviceAlgorithm::Reduce(yLoc, yRange, vtkm::MinAndMax<T>());
214   T minY = yRange[0];
215   T maxY = yRange[1];
216   zRange = DeviceAlgorithm::Reduce(zLoc, zRange, vtkm::MinAndMax<T>());
217   T minZ = zRange[0];
218   T maxZ = zRange[1];
219 
220   numBinsX = static_cast<vtkm::Id>(vtkm::Floor((maxX - minX) / linkLen));
221   numBinsY = static_cast<vtkm::Id>(vtkm::Floor((maxY - minY) / linkLen));
222   numBinsZ = static_cast<vtkm::Id>(vtkm::Floor((maxZ - minZ) / linkLen));
223 
224   vtkm::Id maxBins = 1048576;
225   numBinsX = std::min(maxBins, numBinsX);
226   numBinsY = std::min(maxBins, numBinsY);
227   numBinsZ = std::min(maxBins, numBinsZ);
228 
229   vtkm::Id minBins = 1;
230   numBinsX = std::max(minBins, numBinsX);
231   numBinsY = std::max(minBins, numBinsY);
232   numBinsZ = std::max(minBins, numBinsZ);
233 
234 #ifdef DEBUG_PRINT
235   std::cout << std::endl
236             << "** BinParticlesHalo (" << numBinsX << ", " << numBinsY << ", " << numBinsZ << ") ("
237             << minX << ", " << minY << ", " << minZ << ") (" << maxX << ", " << maxY << ", " << maxZ
238             << ")" << std::endl;
239 #endif
240 
241   // Compute which bin each particle is in
242   ComputeBins<T> computeBins(minX,
243                              maxX, // Physical range on domain
244                              minY,
245                              maxY,
246                              minZ,
247                              maxZ,
248                              numBinsX,
249                              numBinsY,
250                              numBinsZ); // Size of superimposed mesh
251   vtkm::worklet::DispatcherMapField<ComputeBins<T>> computeBinsDispatcher(computeBins);
252   computeBinsDispatcher.Invoke(xLoc,   // input
253                                yLoc,   // input
254                                zLoc,   // input
255                                binId); // output
256 
257   vtkm::cont::ArrayHandleIndex indexArray(nParticles);
258   DeviceAlgorithm::Copy(indexArray, partId);
259 
260 #ifdef DEBUG_PRINT
261   DebugPrint("xLoc", xLoc);
262   DebugPrint("yLoc", yLoc);
263   DebugPrint("zLoc", zLoc);
264   DebugPrint("partId", partId);
265   DebugPrint("binId", binId);
266 #endif
267 
268   // Sort the particles by bin
269   DeviceAlgorithm::SortByKey(binId, partId);
270 
271   // Count the number of particles per bin
272   vtkm::cont::ArrayHandleConstant<vtkm::Id> constArray(1, nParticles);
273   DeviceAlgorithm::ReduceByKey(binId, constArray, uniqueBins, partPerBin, vtkm::Add());
274 #ifdef DEBUG_PRINT
275   DebugPrint("sorted binId", binId);
276   DebugPrint("sorted partId", partId);
277   DebugPrint("uniqueBins", uniqueBins);
278   DebugPrint("partPerBin", partPerBin);
279 #endif
280 
281   // Calculate the bin indices
282   ComputeBinIndices<T> computeBinIndices(numBinsX, numBinsY, numBinsZ);
283   vtkm::worklet::DispatcherMapField<ComputeBinIndices<T>> computeBinIndicesDispatcher(
284     computeBinIndices);
285 
286   computeBinIndicesDispatcher.Invoke(uniqueBins, // input
287                                      binX,       // input
288                                      binY,       // input
289                                      binZ);      // input
290 
291   DeviceAlgorithm::ScanExclusive(partPerBin, particleOffset);
292 }
293 
294 ///////////////////////////////////////////////////////////////////////////////
295 //
296 // Center finder for all particles given location, particle id and key id
297 // Assumed that key and particles are already sorted
298 // MBP (Most Bound Particle) is particle with the minimum potential energy
299 // Method uses ScanInclusiveByKey() and ArrayHandleReverse
300 //
301 ///////////////////////////////////////////////////////////////////////////////
302 template <typename T, typename StorageType>
MBPCenterFindingByKey(vtkm::cont::ArrayHandle<vtkm::Id> & keyId,vtkm::cont::ArrayHandle<vtkm::Id> & partId,vtkm::cont::ArrayHandle<T> & minPotential)303 void CosmoTools<T, StorageType>::MBPCenterFindingByKey(vtkm::cont::ArrayHandle<vtkm::Id>& keyId,
304                                                        vtkm::cont::ArrayHandle<vtkm::Id>& partId,
305                                                        vtkm::cont::ArrayHandle<T>& minPotential)
306 {
307   // Compute starting and ending indices of each key (bin or halo)
308   vtkm::cont::ArrayHandleIndex indexArray(nParticles);
309   vtkm::cont::ArrayHandle<T> potential;
310 
311   vtkm::cont::ArrayHandleReverse<vtkm::cont::ArrayHandle<vtkm::Id>> keyReverse(keyId);
312   vtkm::cont::ArrayHandleReverse<vtkm::cont::ArrayHandle<T>> minPotReverse(minPotential);
313 
314   // Compute indices of all left neighbor bins per bin not per particle
315   vtkm::cont::ArrayHandle<vtkm::Id> leftNeighbor;
316   vtkm::cont::ArrayHandle<vtkm::Id> rightNeighbor;
317   leftNeighbor.Allocate(NUM_NEIGHBORS * nParticles);
318   rightNeighbor.Allocate(NUM_NEIGHBORS * nParticles);
319 
320   vtkm::cont::ArrayHandleIndex countArray(nParticles);
321   ComputeNeighborBins computeNeighborBins(numBinsX, numBinsY, numBinsZ, NUM_NEIGHBORS);
322   vtkm::worklet::DispatcherMapField<ComputeNeighborBins> computeNeighborBinsDispatcher(
323     computeNeighborBins);
324   computeNeighborBinsDispatcher.Invoke(countArray, keyId, leftNeighbor);
325 
326   // Compute indices of all right neighbor bins
327   ComputeBinRange computeBinRange(numBinsX);
328   vtkm::worklet::DispatcherMapField<ComputeBinRange> computeBinRangeDispatcher(computeBinRange);
329   computeBinRangeDispatcher.Invoke(leftNeighbor, rightNeighbor);
330 
331   // Convert bin range to particle range within the bins
332   DeviceAlgorithm::LowerBounds(keyId, leftNeighbor, leftNeighbor);
333   DeviceAlgorithm::UpperBounds(keyId, rightNeighbor, rightNeighbor);
334 #ifdef DEBUG_PRINT
335   DebugPrint("leftNeighbor", leftNeighbor);
336   DebugPrint("rightNeighbor", rightNeighbor);
337 #endif
338 
339   // Initialize halo id of each particle to itself
340   // Compute potentials on particles in 27 neighbors to find minimum
341   ComputePotentialNeighbors<T> computePotentialNeighbors(
342     numBinsX, numBinsY, numBinsZ, NUM_NEIGHBORS, particleMass);
343   vtkm::worklet::DispatcherMapField<ComputePotentialNeighbors<T>>
344     computePotentialNeighborsDispatcher(computePotentialNeighbors);
345 
346   computePotentialNeighborsDispatcher.Invoke(indexArray,
347                                              keyId,         // input (whole array)
348                                              partId,        // input (whole array)
349                                              xLoc,          // input (whole array)
350                                              yLoc,          // input (whole array)
351                                              zLoc,          // input (whole array)
352                                              leftNeighbor,  // input (whole array)
353                                              rightNeighbor, // input (whole array)
354                                              potential);    // output
355 
356   // Find minimum potential for all particles in a halo
357   DeviceAlgorithm::ScanInclusiveByKey(keyId, potential, minPotential, vtkm::Minimum());
358   DeviceAlgorithm::ScanInclusiveByKey(keyReverse, minPotReverse, minPotReverse, vtkm::Minimum());
359 #ifdef DEBUG_PRINT
360   DebugPrint("potential", potential);
361   DebugPrint("minPotential", minPotential);
362 #endif
363 
364   // Find the particle id matching the minimum potential
365   vtkm::cont::ArrayHandle<vtkm::Id> centerId;
366   EqualsMinimumPotential<T> equalsMinimumPotential;
367   vtkm::worklet::DispatcherMapField<EqualsMinimumPotential<T>> equalsMinimumPotentialDispatcher(
368     equalsMinimumPotential);
369 
370   equalsMinimumPotentialDispatcher.Invoke(partId, potential, minPotential, centerId);
371 }
372 
373 ///////////////////////////////////////////////////////////////////////////////
374 //
375 // Center finder for particles in a single halo given location and particle id
376 // MBP (Most Bound Particle) is particle with the minimum potential energy
377 // Method uses ScanInclusiveByKey() and ArrayHandleReverse
378 //
379 ///////////////////////////////////////////////////////////////////////////////
380 template <typename T, typename StorageType>
MBPCenterFinderNxN(T * nxnPotential)381 vtkm::Id CosmoTools<T, StorageType>::MBPCenterFinderNxN(T* nxnPotential)
382 {
383   vtkm::cont::ArrayHandle<T> potential;
384   vtkm::cont::ArrayHandle<T> minPotential;
385 
386   vtkm::cont::ArrayHandleReverse<vtkm::cont::ArrayHandle<T>> minPotReverse(minPotential);
387 
388   vtkm::cont::ArrayHandleIndex particleIndex(nParticles);
389 
390   // Compute potentials (Worklet)
391   ComputePotentialNxN<T> computePotentialHalo(nParticles, particleMass);
392   vtkm::worklet::DispatcherMapField<ComputePotentialNxN<T>> computePotentialHaloDispatcher(
393     computePotentialHalo);
394 
395   computePotentialHaloDispatcher.Invoke(particleIndex, // input
396                                         xLoc,          // input (whole array)
397                                         yLoc,          // input (whole array)
398                                         zLoc,          // input (whole array)
399                                         potential);    // output
400 
401   // Find minimum potential for all particles in a halo
402   DeviceAlgorithm::ScanInclusive(potential, minPotential, vtkm::Minimum());
403   DeviceAlgorithm::ScanInclusive(minPotReverse, minPotReverse, vtkm::Minimum());
404 
405   // Find the particle id matching the minimum potential
406   vtkm::cont::ArrayHandle<vtkm::Id> centerId;
407   EqualsMinimumPotential<T> equalsMinimumPotential;
408   vtkm::worklet::DispatcherMapField<EqualsMinimumPotential<T>> equalsMinimumPotentialDispatcher(
409     equalsMinimumPotential);
410 
411   equalsMinimumPotentialDispatcher.Invoke(particleIndex, potential, minPotential, centerId);
412 
413   // Fill out entire array with center index
414   vtkm::cont::ArrayHandleReverse<vtkm::cont::ArrayHandle<vtkm::Id>> centerIdReverse(centerId);
415   DeviceAlgorithm::ScanInclusive(centerId, centerId, vtkm::Maximum());
416   DeviceAlgorithm::ScanInclusive(centerIdReverse, centerIdReverse, vtkm::Maximum());
417 
418   vtkm::Id nxnMBP = vtkm::cont::ArrayGetValue(0, centerId);
419   *nxnPotential = vtkm::cont::ArrayGetValue(nxnMBP, potential);
420 
421   return nxnMBP;
422 }
423 }
424 }
425 }
426 #endif
427