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