vtk-m/vtkm/worklet/cosmotools/CosmoToolsCenterFinder.h
Allison Vacanti 3eb91af961 Use ArrayGetValue where possible in worklets.
Replace code such as `myArray.GetPortalControl().Get(0)` with
`vtkm::cont::ArrayGetValue(0, myArray)`, which will just retrieve
the single value without transferring the entire contents of
`myArray` to the host. This should prevent memory thrashing.

Also did some general style clean-ups and fixed loops that called
`GetPortalControl` in the loop body.
2019-08-22 16:42:17 -04:00

427 lines
20 KiB
C++

//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
// Copyright (c) 2016, Los Alamos National Security, LLC
// All rights reserved.
//
// Copyright 2016. Los Alamos National Security, LLC.
// This software was produced under U.S. Government contract DE-AC52-06NA25396
// for Los Alamos National Laboratory (LANL), which is operated by
// Los Alamos National Security, LLC for the U.S. Department of Energy.
// The U.S. Government has rights to use, reproduce, and distribute this
// software. NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC
// MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE
// USE OF THIS SOFTWARE. If software is modified to produce derivative works,
// such modified software should be clearly marked, so as not to confuse it
// with the version available from LANL.
//
// Additionally, redistribution and use in source and binary forms, with or
// without modification, are permitted provided that the following conditions
// are met:
//
// 1. Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// 2. Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// 3. Neither the name of Los Alamos National Security, LLC, Los Alamos
// National Laboratory, LANL, the U.S. Government, nor the names of its
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY LOS ALAMOS NATIONAL SECURITY, LLC AND
// CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
// BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL LOS ALAMOS
// NATIONAL SECURITY, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
// USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
// THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//============================================================================
#ifndef vtkm_worklet_cosmotools_cosmotools_centerfinder_h
#define vtkm_worklet_cosmotools_cosmotools_centerfinder_h
#include <vtkm/worklet/cosmotools/CosmoTools.h>
#include <vtkm/cont/ArrayGetValues.h>
namespace vtkm
{
namespace worklet
{
namespace cosmotools
{
///////////////////////////////////////////////////////////////////////////////
//
// Center finder for particles in FOF halo using estimations but with exact final answer
// MBP (Most Bound Particle) is particle with the minimum potential energy
//
///////////////////////////////////////////////////////////////////////////////
template <typename T, typename StorageType>
vtkm::Id CosmoTools<T, StorageType>::MBPCenterFinderMxN(T* mxnPotential)
{
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;
// Bin all particles in the halo into bins of size linking length
BinParticlesHalo(partId, binId, uniqueBins, partPerBin, particleOffset, binX, binY, binZ);
#ifdef DEBUG_PRINT
DebugPrint("uniqueBins", uniqueBins);
DebugPrint("partPerBin", partPerBin);
#endif
// Compute the estimated potential per bin using 27 contiguous bins
vtkm::cont::ArrayHandle<T> partPotential;
MBPCenterFindingByKey(binId, partId, partPotential);
// Reduce by key to get the estimated minimum potential per bin within 27 neighbors
vtkm::cont::ArrayHandle<vtkm::Id> tempId;
vtkm::cont::ArrayHandle<T> minPotential;
DeviceAlgorithm::ReduceByKey(binId, partPotential, tempId, minPotential, vtkm::Minimum());
// Reduce by key to get the estimated maximum potential per bin within 27 neighbors
vtkm::cont::ArrayHandle<T> maxPotential;
DeviceAlgorithm::ReduceByKey(binId, partPotential, tempId, maxPotential, vtkm::Maximum());
#ifdef DEBUG_PRINT
DebugPrint("minPotential", minPotential);
DebugPrint("maxPotential", maxPotential);
#endif
// Compute potentials estimate for a bin using all other bins
// Particles in the other bins are located at the closest point to this bin
vtkm::cont::ArrayHandleIndex uniqueIndex(uniqueBins.GetNumberOfValues());
vtkm::cont::ArrayHandle<T> bestEstPotential;
vtkm::cont::ArrayHandle<T> worstEstPotential;
// Initialize each bin potential with the nxn for that bin
DeviceAlgorithm::Copy(minPotential, bestEstPotential);
DeviceAlgorithm::Copy(maxPotential, worstEstPotential);
// Estimate only across the uniqueBins that contain particles
ComputePotentialBin<T> computePotentialBin(uniqueBins.GetNumberOfValues(), particleMass, linkLen);
vtkm::worklet::DispatcherMapField<ComputePotentialBin<T>> computePotentialBinDispatcher(
computePotentialBin);
computePotentialBinDispatcher.Invoke(uniqueIndex, // input
partPerBin, // input (whole array)
binX, // input (whole array)
binY, // input (whole array)
binZ, // input (whole array)
bestEstPotential, // input/output
worstEstPotential); // input/output
#ifdef DEBUG_PRINT
DebugPrint("bestEstPotential", bestEstPotential);
DebugPrint("worstEstPotential", worstEstPotential);
std::cout << "Number of bestEstPotential " << bestEstPotential.GetNumberOfValues() << std::endl;
std::cout << "Number of worstEstPotential " << worstEstPotential.GetNumberOfValues() << std::endl;
#endif
// Sort everything by the best estimated potential per bin
vtkm::cont::ArrayHandle<T> tempBest;
DeviceAlgorithm::Copy(bestEstPotential, tempBest);
DeviceAlgorithm::SortByKey(tempBest, worstEstPotential);
// Use the worst estimate for the first selected bin to compare to best of all others
// Any bin that passes is a candidate for having the MBP
T cutoffPotential = vtkm::cont::ArrayGetValue(0, worstEstPotential);
worstEstPotential.ReleaseResources();
tempBest.ReleaseResources();
vtkm::cont::ArrayHandle<vtkm::Id> candidate;
DeviceAlgorithm::Copy(vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, nParticles), candidate);
SetCandidateParticles<T> setCandidateParticles(cutoffPotential);
vtkm::worklet::DispatcherMapField<SetCandidateParticles<T>> setCandidateParticlesDispatcher(
setCandidateParticles);
setCandidateParticlesDispatcher.Invoke(bestEstPotential, // input
particleOffset, // input
partPerBin, // input
candidate); // output (whole array)
// Copy the M candidate particles to a new array
vtkm::cont::ArrayHandle<vtkm::Id> mparticles;
DeviceAlgorithm::CopyIf(partId, candidate, mparticles);
// Compute potentials only on the candidate particles
vtkm::cont::ArrayHandle<T> mpotential;
ComputePotentialOnCandidates<T> computePotentialOnCandidates(nParticles, particleMass);
vtkm::worklet::DispatcherMapField<ComputePotentialOnCandidates<T>>
computePotentialOnCandidatesDispatcher(computePotentialOnCandidates);
computePotentialOnCandidatesDispatcher.Invoke(mparticles,
xLoc, // input (whole array)
yLoc, // input (whole array)
zLoc, // input (whole array)
mpotential); // output
// Of the M candidate particles which has the minimum potential
DeviceAlgorithm::SortByKey(mpotential, mparticles);
#ifdef DEBUG_PRINT
DebugPrint("mparticles", mparticles);
DebugPrint("mpotential", mpotential);
#endif
// Return the found MBP particle and its potential
vtkm::Id mxnMBP = vtkm::cont::ArrayGetValue(0, mparticles);
*mxnPotential = vtkm::cont::ArrayGetValue(0, mpotential);
return mxnMBP;
}
///////////////////////////////////////////////////////////////////////////////
//
// Bin particles in one halo for quick MBP finding
//
///////////////////////////////////////////////////////////////////////////////
template <typename T, typename StorageType>
void CosmoTools<T, 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)
{
// Compute number of bins and ranges for each bin
vtkm::Vec<T, 2> xRange(vtkm::cont::ArrayGetValue(0, xLoc));
vtkm::Vec<T, 2> yRange(vtkm::cont::ArrayGetValue(0, yLoc));
vtkm::Vec<T, 2> zRange(vtkm::cont::ArrayGetValue(0, zLoc));
xRange = DeviceAlgorithm::Reduce(xLoc, xRange, vtkm::MinAndMax<T>());
T minX = xRange[0];
T maxX = xRange[1];
yRange = DeviceAlgorithm::Reduce(yLoc, yRange, vtkm::MinAndMax<T>());
T minY = yRange[0];
T maxY = yRange[1];
zRange = DeviceAlgorithm::Reduce(zLoc, zRange, vtkm::MinAndMax<T>());
T minZ = zRange[0];
T maxZ = zRange[1];
numBinsX = static_cast<vtkm::Id>(vtkm::Floor((maxX - minX) / linkLen));
numBinsY = static_cast<vtkm::Id>(vtkm::Floor((maxY - minY) / linkLen));
numBinsZ = static_cast<vtkm::Id>(vtkm::Floor((maxZ - minZ) / linkLen));
vtkm::Id maxBins = 1048576;
numBinsX = std::min(maxBins, numBinsX);
numBinsY = std::min(maxBins, numBinsY);
numBinsZ = std::min(maxBins, numBinsZ);
vtkm::Id minBins = 1;
numBinsX = std::max(minBins, numBinsX);
numBinsY = std::max(minBins, numBinsY);
numBinsZ = std::max(minBins, numBinsZ);
#ifdef DEBUG_PRINT
std::cout << std::endl
<< "** BinParticlesHalo (" << numBinsX << ", " << numBinsY << ", " << numBinsZ << ") ("
<< minX << ", " << minY << ", " << minZ << ") (" << maxX << ", " << maxY << ", " << maxZ
<< ")" << std::endl;
#endif
// Compute which bin each particle is in
ComputeBins<T> computeBins(minX,
maxX, // Physical range on domain
minY,
maxY,
minZ,
maxZ,
numBinsX,
numBinsY,
numBinsZ); // Size of superimposed mesh
vtkm::worklet::DispatcherMapField<ComputeBins<T>> computeBinsDispatcher(computeBins);
computeBinsDispatcher.Invoke(xLoc, // input
yLoc, // input
zLoc, // input
binId); // output
vtkm::cont::ArrayHandleIndex indexArray(nParticles);
DeviceAlgorithm::Copy(indexArray, partId);
#ifdef DEBUG_PRINT
DebugPrint("xLoc", xLoc);
DebugPrint("yLoc", yLoc);
DebugPrint("zLoc", zLoc);
DebugPrint("partId", partId);
DebugPrint("binId", binId);
#endif
// Sort the particles by bin
DeviceAlgorithm::SortByKey(binId, partId);
// Count the number of particles per bin
vtkm::cont::ArrayHandleConstant<vtkm::Id> constArray(1, nParticles);
DeviceAlgorithm::ReduceByKey(binId, constArray, uniqueBins, partPerBin, vtkm::Add());
#ifdef DEBUG_PRINT
DebugPrint("sorted binId", binId);
DebugPrint("sorted partId", partId);
DebugPrint("uniqueBins", uniqueBins);
DebugPrint("partPerBin", partPerBin);
#endif
// Calculate the bin indices
ComputeBinIndices<T> computeBinIndices(numBinsX, numBinsY, numBinsZ);
vtkm::worklet::DispatcherMapField<ComputeBinIndices<T>> computeBinIndicesDispatcher(
computeBinIndices);
computeBinIndicesDispatcher.Invoke(uniqueBins, // input
binX, // input
binY, // input
binZ); // input
DeviceAlgorithm::ScanExclusive(partPerBin, particleOffset);
}
///////////////////////////////////////////////////////////////////////////////
//
// Center finder for all particles given location, particle id and key id
// Assumed that key and particles are already sorted
// MBP (Most Bound Particle) is particle with the minimum potential energy
// Method uses ScanInclusiveByKey() and ArrayHandleReverse
//
///////////////////////////////////////////////////////////////////////////////
template <typename T, typename StorageType>
void CosmoTools<T, StorageType>::MBPCenterFindingByKey(vtkm::cont::ArrayHandle<vtkm::Id>& keyId,
vtkm::cont::ArrayHandle<vtkm::Id>& partId,
vtkm::cont::ArrayHandle<T>& minPotential)
{
// Compute starting and ending indices of each key (bin or halo)
vtkm::cont::ArrayHandleIndex indexArray(nParticles);
vtkm::cont::ArrayHandle<T> potential;
vtkm::cont::ArrayHandleReverse<vtkm::cont::ArrayHandle<vtkm::Id>> keyReverse(keyId);
vtkm::cont::ArrayHandleReverse<vtkm::cont::ArrayHandle<T>> minPotReverse(minPotential);
// Compute indices of all left neighbor bins per bin not per particle
vtkm::cont::ArrayHandle<vtkm::Id> leftNeighbor;
vtkm::cont::ArrayHandle<vtkm::Id> rightNeighbor;
leftNeighbor.Allocate(NUM_NEIGHBORS * nParticles);
rightNeighbor.Allocate(NUM_NEIGHBORS * nParticles);
vtkm::cont::ArrayHandleIndex countArray(nParticles);
ComputeNeighborBins computeNeighborBins(numBinsX, numBinsY, numBinsZ, NUM_NEIGHBORS);
vtkm::worklet::DispatcherMapField<ComputeNeighborBins> computeNeighborBinsDispatcher(
computeNeighborBins);
computeNeighborBinsDispatcher.Invoke(countArray, keyId, leftNeighbor);
// Compute indices of all right neighbor bins
ComputeBinRange computeBinRange(numBinsX);
vtkm::worklet::DispatcherMapField<ComputeBinRange> computeBinRangeDispatcher(computeBinRange);
computeBinRangeDispatcher.Invoke(leftNeighbor, rightNeighbor);
// Convert bin range to particle range within the bins
DeviceAlgorithm::LowerBounds(keyId, leftNeighbor, leftNeighbor);
DeviceAlgorithm::UpperBounds(keyId, rightNeighbor, rightNeighbor);
#ifdef DEBUG_PRINT
DebugPrint("leftNeighbor", leftNeighbor);
DebugPrint("rightNeighbor", rightNeighbor);
#endif
// Initialize halo id of each particle to itself
// Compute potentials on particles in 27 neighbors to find minimum
ComputePotentialNeighbors<T> computePotentialNeighbors(
numBinsX, numBinsY, numBinsZ, NUM_NEIGHBORS, particleMass);
vtkm::worklet::DispatcherMapField<ComputePotentialNeighbors<T>>
computePotentialNeighborsDispatcher(computePotentialNeighbors);
computePotentialNeighborsDispatcher.Invoke(indexArray,
keyId, // input (whole array)
partId, // input (whole array)
xLoc, // input (whole array)
yLoc, // input (whole array)
zLoc, // input (whole array)
leftNeighbor, // input (whole array)
rightNeighbor, // input (whole array)
potential); // output
// Find minimum potential for all particles in a halo
DeviceAlgorithm::ScanInclusiveByKey(keyId, potential, minPotential, vtkm::Minimum());
DeviceAlgorithm::ScanInclusiveByKey(keyReverse, minPotReverse, minPotReverse, vtkm::Minimum());
#ifdef DEBUG_PRINT
DebugPrint("potential", potential);
DebugPrint("minPotential", minPotential);
#endif
// Find the particle id matching the minimum potential
vtkm::cont::ArrayHandle<vtkm::Id> centerId;
EqualsMinimumPotential<T> equalsMinimumPotential;
vtkm::worklet::DispatcherMapField<EqualsMinimumPotential<T>> equalsMinimumPotentialDispatcher(
equalsMinimumPotential);
equalsMinimumPotentialDispatcher.Invoke(partId, potential, minPotential, centerId);
}
///////////////////////////////////////////////////////////////////////////////
//
// Center finder for particles in a single halo given location and particle id
// MBP (Most Bound Particle) is particle with the minimum potential energy
// Method uses ScanInclusiveByKey() and ArrayHandleReverse
//
///////////////////////////////////////////////////////////////////////////////
template <typename T, typename StorageType>
vtkm::Id CosmoTools<T, StorageType>::MBPCenterFinderNxN(T* nxnPotential)
{
vtkm::cont::ArrayHandle<T> potential;
vtkm::cont::ArrayHandle<T> minPotential;
vtkm::cont::ArrayHandleReverse<vtkm::cont::ArrayHandle<T>> minPotReverse(minPotential);
vtkm::cont::ArrayHandleIndex particleIndex(nParticles);
// Compute potentials (Worklet)
ComputePotentialNxN<T> computePotentialHalo(nParticles, particleMass);
vtkm::worklet::DispatcherMapField<ComputePotentialNxN<T>> computePotentialHaloDispatcher(
computePotentialHalo);
computePotentialHaloDispatcher.Invoke(particleIndex, // input
xLoc, // input (whole array)
yLoc, // input (whole array)
zLoc, // input (whole array)
potential); // output
// Find minimum potential for all particles in a halo
DeviceAlgorithm::ScanInclusive(potential, minPotential, vtkm::Minimum());
DeviceAlgorithm::ScanInclusive(minPotReverse, minPotReverse, vtkm::Minimum());
// Find the particle id matching the minimum potential
vtkm::cont::ArrayHandle<vtkm::Id> centerId;
EqualsMinimumPotential<T> equalsMinimumPotential;
vtkm::worklet::DispatcherMapField<EqualsMinimumPotential<T>> equalsMinimumPotentialDispatcher(
equalsMinimumPotential);
equalsMinimumPotentialDispatcher.Invoke(particleIndex, potential, minPotential, centerId);
// Fill out entire array with center index
vtkm::cont::ArrayHandleReverse<vtkm::cont::ArrayHandle<vtkm::Id>> centerIdReverse(centerId);
DeviceAlgorithm::ScanInclusive(centerId, centerId, vtkm::Maximum());
DeviceAlgorithm::ScanInclusive(centerIdReverse, centerIdReverse, vtkm::Maximum());
vtkm::Id nxnMBP = vtkm::cont::ArrayGetValue(0, centerId);
*nxnPotential = vtkm::cont::ArrayGetValue(nxnMBP, potential);
return nxnMBP;
}
}
}
}
#endif