Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-10-05 01:49:02 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Eliminate need for SpatialDecomposition in contour tree filters.

Browse Source 
This commit address issue #713 and removes the need to pass information
about the spatial decomposition (block origins, block sizes, block index
and blocks per dimension) to the contour tree filter. Block origin
information is added to CellSetStructured (as GlobalPointSize) and the
distributed contour tree filter will get this information from
CellSetStructured instead of expecting it as parameters to the
constructor. Information about blocks per dimension and block indices
are computed from the information in CellSetStructure albeit requiring
global communication across all ranks. To avoid this communication cost,
the caller of the filter may explicitly specify this information via the
SetBlockIndices() method.
This commit is contained in:
Gunther H. Weber 2022-07-13 11:27:56 -07:00
parent 30dda88e59
commit cd3fe3d47d
22 changed files with 913 additions and 538 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
examples/contour_tree_augmented/ContourTreeApp.cxx
View File

@ -548,14 +548,8 @@ int main(int argc, char* argv[])
static_cast<vtkm::Id>(dims[1]),
static_cast<vtkm::Id>(0));
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockIndices;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockOrigins;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockSizes;
localBlockIndices.Allocate(blocksPerRank);
localBlockOrigins.Allocate(blocksPerRank);
localBlockSizes.Allocate(blocksPerRank);
auto localBlockIndicesPortal = localBlockIndices.WritePortal();
auto localBlockOriginsPortal = localBlockOrigins.WritePortal();
auto localBlockSizesPortal = localBlockSizes.WritePortal();
{
vtkm::Id lastDimSize =
@ -604,30 +598,29 @@ int main(int argc, char* argv[])
vtkm::Vec<ValueType, 2> origin(0, blockIndex * blockSize);
vtkm::Vec<ValueType, 2> spacing(1, 1);
ds = dsb.Create(vdims, origin, spacing);
vtkm::cont::CellSetStructured<2> cs;
cs.SetPointDimensions(vdims);
cs.SetGlobalPointDimensions(vtkm::Id2{ globalSize[0], globalSize[1] });
cs.SetGlobalPointIndexStart(vtkm::Id2{ 0, blockIndex * blockSize });
ds.SetCellSet(cs);
localBlockIndicesPortal.Set(localBlockIndex, vtkm::Id3(blockIndex, 0, 0));
localBlockOriginsPortal.Set(localBlockIndex,
vtkm::Id3((blockStart / blockSliceSize), 0, 0));
localBlockSizesPortal.Set(localBlockIndex,
vtkm::Id3(currBlockSize, static_cast<vtkm::Id>(dims[0]), 0));
localBlockIndicesPortal.Set(localBlockIndex, vtkm::Id3(0, blockIndex, 0));
}
// 3D data
else
{
vtkm::Id3 vdims;
vdims[0] = static_cast<vtkm::Id>(dims[1]);
vdims[1] = static_cast<vtkm::Id>(dims[0]);
vdims[0] = static_cast<vtkm::Id>(dims[0]);
vdims[1] = static_cast<vtkm::Id>(dims[1]);
vdims[2] = static_cast<vtkm::Id>(currBlockSize);
vtkm::Vec<ValueType, 3> origin(0, 0, (blockIndex * blockSize));
vtkm::Vec<ValueType, 3> origin(0, 0, blockIndex * blockSize);
vtkm::Vec<ValueType, 3> spacing(1, 1, 1);
ds = dsb.Create(vdims, origin, spacing);
localBlockIndicesPortal.Set(localBlockIndex, vtkm::Id3(0, 0, blockIndex));
localBlockOriginsPortal.Set(localBlockIndex,
vtkm::Id3(0, 0, (blockStart / blockSliceSize)));
localBlockSizesPortal.Set(
localBlockIndex,
vtkm::Id3(static_cast<vtkm::Id>(dims[0]), static_cast<vtkm::Id>(dims[1]), currBlockSize));
vtkm::cont::CellSetStructured<3> cs;
cs.SetPointDimensions(vdims);
cs.SetGlobalPointDimensions(globalSize);
cs.SetGlobalPointIndexStart(vtkm::Id3{ 0, 0, blockIndex * blockSize });
ds.SetCellSet(cs);
}
std::vector<vtkm::Float32> subValues((values.begin() + blockStart),
@ -648,8 +641,7 @@ int main(int argc, char* argv[])
computeRegularStructure);
#ifdef WITH_MPI
filter.SetSpatialDecomposition(
blocksPerDim, globalSize, localBlockIndices, localBlockOrigins, localBlockSizes);
filter.SetBlockIndices(blocksPerDim, localBlockIndices);
#endif
filter.SetActiveField("values");

60
examples/contour_tree_distributed/ContourTreeApp.cxx
View File

@ -401,14 +401,8 @@ int main(int argc, char* argv[])
vtkm::Id3 globalSize;
vtkm::Id3 blocksPerDim;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockIndices;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockOrigins;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockSizes;
localBlockIndices.Allocate(blocksPerRank);
localBlockOrigins.Allocate(blocksPerRank);
localBlockSizes.Allocate(blocksPerRank);
auto localBlockIndicesPortal = localBlockIndices.WritePortal();
auto localBlockOriginsPortal = localBlockOrigins.WritePortal();
auto localBlockSizesPortal = localBlockSizes.WritePortal();
// Read the pre-split data files
if (preSplitFiles)
@ -595,6 +589,11 @@ int main(int argc, char* argv[])
static_cast<ValueType>(offset[1]) };
const vtkm::Vec<ValueType, 2> v_spacing{ 1, 1 };
ds = dsb.Create(v_dims, v_origin, v_spacing);
vtkm::cont::CellSetStructured<2> cs;
cs.SetPointDimensions(v_dims);
cs.SetGlobalPointDimensions(vtkm::Id2{ globalSize[0], globalSize[1] });
cs.SetGlobalPointIndexStart(vtkm::Id2{ offset[0], offset[1] });
ds.SetCellSet(cs);
}
else
{
@ -607,6 +606,11 @@ int main(int argc, char* argv[])
static_cast<ValueType>(offset[2]) };
vtkm::Vec<ValueType, 3> v_spacing(1, 1, 1);
ds = dsb.Create(v_dims, v_origin, v_spacing);
vtkm::cont::CellSetStructured<3> cs;
cs.SetPointDimensions(v_dims);
cs.SetGlobalPointDimensions(globalSize);
cs.SetGlobalPointIndexStart(vtkm::Id3{ offset[0], offset[1], offset[2] });
ds.SetCellSet(cs);
}
ds.AddPointField("values", values);
// and add to partition
@ -617,14 +621,6 @@ int main(int argc, char* argv[])
vtkm::Id3{ static_cast<vtkm::Id>(blockIndex[0]),
static_cast<vtkm::Id>(blockIndex[1]),
static_cast<vtkm::Id>(nDims == 3 ? blockIndex[2] : 0) });
localBlockOriginsPortal.Set(blockNo,
vtkm::Id3{ static_cast<vtkm::Id>(offset[0]),
static_cast<vtkm::Id>(offset[1]),
static_cast<vtkm::Id>(nDims == 3 ? offset[2] : 0) });
localBlockSizesPortal.Set(blockNo,
vtkm::Id3{ static_cast<vtkm::Id>(dims[0]),
static_cast<vtkm::Id>(dims[1]),
static_cast<vtkm::Id>(nDims == 3 ? dims[2] : 0) });
if (blockNo == 0)
{
@ -753,7 +749,6 @@ int main(int argc, char* argv[])
: vtkm::Id3(static_cast<vtkm::Id>(dims[0]),
static_cast<vtkm::Id>(dims[1]),
static_cast<vtkm::Id>(1));
std::cout << blocksPerDim << " " << globalSize << std::endl;
{
vtkm::Id lastDimSize =
(nDims == 2) ? static_cast<vtkm::Id>(dims[1]) : static_cast<vtkm::Id>(dims[2]);
@ -801,12 +796,12 @@ int main(int argc, char* argv[])
vtkm::Vec<ValueType, 2> origin(0, blockIndex * blockSize);
vtkm::Vec<ValueType, 2> spacing(1, 1);
ds = dsb.Create(vdims, origin, spacing);
vtkm::cont::CellSetStructured<2> cs;
cs.SetPointDimensions(vdims);
cs.SetGlobalPointDimensions(vtkm::Id2{ globalSize[0], globalSize[1] });
cs.SetGlobalPointIndexStart(vtkm::Id2{ 0, (blockStart / blockSliceSize) });
ds.SetCellSet(cs);
localBlockIndicesPortal.Set(localBlockIndex, vtkm::Id3(0, blockIndex, 0));
localBlockOriginsPortal.Set(localBlockIndex,
vtkm::Id3(0, (blockStart / blockSliceSize), 0));
localBlockSizesPortal.Set(localBlockIndex,
vtkm::Id3(static_cast<vtkm::Id>(dims[0]), currBlockSize, 0));
}
// 3D data
else
@ -818,14 +813,12 @@ int main(int argc, char* argv[])
vtkm::Vec<ValueType, 3> origin(0, 0, (blockIndex * blockSize));
vtkm::Vec<ValueType, 3> spacing(1, 1, 1);
ds = dsb.Create(vdims, origin, spacing);
vtkm::cont::CellSetStructured<3> cs;
cs.SetPointDimensions(vdims);
cs.SetGlobalPointDimensions(globalSize);
cs.SetGlobalPointIndexStart(vtkm::Id3(0, 0, blockStart / blockSliceSize));
ds.SetCellSet(cs);
localBlockIndicesPortal.Set(localBlockIndex, vtkm::Id3(0, 0, blockIndex));
localBlockOriginsPortal.Set(localBlockIndex,
vtkm::Id3(0, 0, (blockStart / blockSliceSize)));
localBlockSizesPortal.Set(localBlockIndex,
vtkm::Id3(static_cast<vtkm::Id>(dims[0]),
static_cast<vtkm::Id>(dims[1]),
currBlockSize));
}
std::vector<vtkm::Float32> subValues((values.begin() + blockStart),
@ -863,13 +856,9 @@ int main(int argc, char* argv[])
prevTime = currTime;
// Convert the mesh of values into contour tree, pairs of vertex ids
vtkm::filter::ContourTreeUniformDistributed filter(blocksPerDim,
globalSize,
localBlockIndices,
localBlockOrigins,
localBlockSizes,
timingsLogLevel,
treeLogLevel);
vtkm::filter::scalar_topology::ContourTreeUniformDistributed filter(timingsLogLevel,
treeLogLevel);
filter.SetBlockIndices(blocksPerDim, localBlockIndices);
filter.SetUseBoundaryExtremaOnly(useBoundaryExtremaOnly);
filter.SetUseMarchingCubes(useMarchingCubes);
filter.SetAugmentHierarchicalTree(augmentHierarchicalTree);
@ -895,8 +884,7 @@ int main(int argc, char* argv[])
{
if (computeHierarchicalVolumetricBranchDecomposition)
{
vtkm::filter::scalar_topology::DistributedBranchDecompositionFilter bd_filter(
blocksPerDim, globalSize, localBlockIndices, localBlockOrigins, localBlockSizes);
vtkm::filter::scalar_topology::DistributedBranchDecompositionFilter bd_filter;
auto bd_result = bd_filter.Execute(result);
for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no)

20
vtkm/cont/CellSetStructured.h
View File

@ -51,6 +51,11 @@ public:
this->Structure.SetPointDimensions(dimensions);
}
void SetGlobalPointDimensions(SchedulingRangeType dimensions)
{
this->Structure.SetGlobalPointDimensions(dimensions);
}
void SetGlobalPointIndexStart(SchedulingRangeType start)
{
this->Structure.SetGlobalPointIndexStart(start);
@ -58,8 +63,18 @@ public:
SchedulingRangeType GetPointDimensions() const { return this->Structure.GetPointDimensions(); }
SchedulingRangeType GetGlobalPointDimensions() const
{
return this->Structure.GetGlobalPointDimensions();
}
SchedulingRangeType GetCellDimensions() const { return this->Structure.GetCellDimensions(); }
SchedulingRangeType GetGlobalCellDimensions() const
{
return this->Structure.GetGlobalCellDimensions();
}
SchedulingRangeType GetGlobalPointIndexStart() const
{
return this->Structure.GetGlobalPointIndexStart();
@ -225,17 +240,20 @@ public:
static VTKM_CONT void save(BinaryBuffer& bb, const Type& cs)
{
vtkmdiy::save(bb, cs.GetPointDimensions());
vtkmdiy::save(bb, cs.GetGlobalPointDimensions());
vtkmdiy::save(bb, cs.GetGlobalPointIndexStart());
}
static VTKM_CONT void load(BinaryBuffer& bb, Type& cs)
{
typename Type::SchedulingRangeType dims, start;
typename Type::SchedulingRangeType dims, gdims, start;
vtkmdiy::load(bb, dims);
vtkmdiy::load(bb, gdims);
vtkmdiy::load(bb, start);
cs = Type{};
cs.SetPointDimensions(dims);
cs.SetGlobalPointDimensions(gdims);
cs.SetGlobalPointIndexStart(start);
}
};

20
vtkm/cont/PartitionedDataSet.cxx
View File

@ -15,6 +15,13 @@
#include <vtkm/cont/ErrorBadValue.h>
#include <vtkm/cont/Field.h>
#include <vtkm/cont/PartitionedDataSet.h>
#include <vtkm/internal/Configure.h>
// clang-format off
VTKM_THIRDPARTY_PRE_INCLUDE
#include <vtkm/thirdparty/diy/diy.h>
VTKM_THIRDPARTY_POST_INCLUDE
// clang-format on
namespace vtkm
{
@ -73,6 +80,19 @@ vtkm::Id PartitionedDataSet::GetNumberOfPartitions() const
return static_cast<vtkm::Id>(this->Partitions.size());
}
VTKM_CONT
vtkm::Id PartitionedDataSet::GetGlobalNumberOfPartitions() const
{
#ifdef VTKM_ENABLE_MPI
auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
vtkm::Id globalSize = 0;
vtkmdiy::mpi::all_reduce(comm, GetNumberOfPartitions(), globalSize, std::plus<vtkm::Id>{});
return globalSize;
#else
return GetNumberOfPartitions();
#endif
}
VTKM_CONT
const vtkm::cont::DataSet& PartitionedDataSet::GetPartition(vtkm::Id blockId) const
{

28
vtkm/cont/PartitionedDataSet.h
View File

@ -32,16 +32,16 @@ public:
using reference = typename StorageVec::reference;
using const_reference = typename StorageVec::const_reference;
/// create a new PartitionedDataSet containng a single DataSet @a ds
/// Create a new PartitionedDataSet containng a single DataSet @a ds.
VTKM_CONT
PartitionedDataSet(const vtkm::cont::DataSet& ds);
/// create a new PartitionedDataSet with the existing one @a src
/// Create a new PartitionedDataSet with the existing one @a src.
VTKM_CONT
PartitionedDataSet(const vtkm::cont::PartitionedDataSet& src);
/// create a new PartitionedDataSet with a DataSet vector @a partitions.
/// Create a new PartitionedDataSet with a DataSet vector @a partitions.
VTKM_CONT
explicit PartitionedDataSet(const std::vector<vtkm::cont::DataSet>& partitions);
/// create a new PartitionedDataSet with the capacity set to be @a size
/// Create a new PartitionedDataSet with the capacity set to be @a size.
VTKM_CONT
explicit PartitionedDataSet(vtkm::Id size);
@ -53,33 +53,41 @@ public:
VTKM_CONT
~PartitionedDataSet();
/// get the field @a field_name from partition @a partition_index
/// Get the field @a field_name from partition @a partition_index.
VTKM_CONT
vtkm::cont::Field GetField(const std::string& field_name, int partition_index) const;
/// Get number of DataSet objects stored in this PartitionedDataSet.
VTKM_CONT
vtkm::Id GetNumberOfPartitions() const;
/// Get number of partations across all MPI ranks.
/// @warning This method requires global communication (MPI_Allreduce) if MPI is enabled.
VTKM_CONT
vtkm::Id GetGlobalNumberOfPartitions() const;
/// Get the DataSet @a partId.
VTKM_CONT
const vtkm::cont::DataSet& GetPartition(vtkm::Id partId) const;
/// Get an STL vector of all DataSet objects stored in PartitionedDataSet.
VTKM_CONT
const std::vector<vtkm::cont::DataSet>& GetPartitions() const;
/// add DataSet @a ds to the end of the contained DataSet vector
/// Add DataSet @a ds to the end of the contained DataSet vector.
VTKM_CONT
void AppendPartition(const vtkm::cont::DataSet& ds);
/// add DataSet @a ds to position @a index of the contained DataSet vector
/// Add DataSet @a ds to position @a index of the contained DataSet vector.
VTKM_CONT
void InsertPartition(vtkm::Id index, const vtkm::cont::DataSet& ds);
/// replace the @a index positioned element of the contained DataSet vector
/// with @a ds
/// Replace the @a index positioned element of the contained DataSet vector
/// with @a ds.
VTKM_CONT
void ReplacePartition(vtkm::Id index, const vtkm::cont::DataSet& ds);
/// append the DataSet vector "partitions" to the end of the contained one
/// Append the DataSet vector @a partitions to the end of the contained one.
VTKM_CONT
void AppendPartitions(const std::vector<vtkm::cont::DataSet>& partitions);

1
vtkm/filter/scalar_topology/CMakeLists.txt
View File

@ -16,6 +16,7 @@ set(scalar_topology_headers
set(scalar_topology_sources
internal/BranchDecompositionBlock.cxx
internal/ComputeBlockIndices.cxx
internal/ComputeDistributedBranchDecompositionFunctor.cxx
ContourTreeUniform.cxx
ContourTreeUniformAugmented.cxx

117
vtkm/filter/scalar_topology/ContourTreeUniformAugmented.cxx
View File

@ -51,6 +51,7 @@
//==============================================================================
#include <vtkm/filter/scalar_topology/ContourTreeUniformAugmented.h>
#include <vtkm/filter/scalar_topology/internal/ComputeBlockIndices.h>
#include <vtkm/filter/scalar_topology/worklet/ContourTreeUniformAugmented.h>
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/meshtypes/ContourTreeMesh.h>
@ -83,12 +84,9 @@ ContourTreeAugmented::ContourTreeAugmented(bool useMarchingCubes,
this->SetOutputFieldName("resultData");
}
void ContourTreeAugmented::SetSpatialDecomposition(
void ContourTreeAugmented::SetBlockIndices(
vtkm::Id3 blocksPerDim,
vtkm::Id3 globalSize,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockIndices,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockOrigins,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockSizes)
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockIndices)
{
if (this->MultiBlockTreeHelper)
{
@ -96,7 +94,7 @@ void ContourTreeAugmented::SetSpatialDecomposition(
}
this->MultiBlockTreeHelper =
std::make_unique<vtkm::worklet::contourtree_distributed::MultiBlockContourTreeHelper>(
blocksPerDim, globalSize, localBlockIndices, localBlockOrigins, localBlockSizes);
blocksPerDim, localBlockIndices);
}
const vtkm::worklet::contourtree_augmented::ContourTree& ContourTreeAugmented::GetContourTree()
@ -234,7 +232,7 @@ VTKM_CONT void ContourTreeAugmented::PreExecute(const vtkm::cont::PartitionedDat
{
if (this->MultiBlockTreeHelper)
{
if (this->MultiBlockTreeHelper->GetGlobalNumberOfBlocks(input) !=
if (input.GetGlobalNumberOfPartitions() !=
this->MultiBlockTreeHelper->GetGlobalNumberOfBlocks())
{
throw vtkm::cont::ErrorFilterExecution(
@ -246,6 +244,12 @@ VTKM_CONT void ContourTreeAugmented::PreExecute(const vtkm::cont::PartitionedDat
"Global number of block in MultiBlock dataset does not match the SpatialDecomposition");
}
}
else
{
// No block indices set -> compute information automatically later
this->MultiBlockTreeHelper =
std::make_unique<vtkm::worklet::contourtree_distributed::MultiBlockContourTreeHelper>(input);
}
}
//-----------------------------------------------------------------------------
@ -268,10 +272,6 @@ VTKM_CONT void ContourTreeAugmented::DoPostExecute(const vtkm::cont::Partitioned
unsigned int compRegularStruct =
(this->ComputeRegularStructure > 0) ? this->ComputeRegularStructure : 2;
auto localBlocksOriginPortal =
this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.LocalBlockOrigins.ReadPortal();
auto localBlocksSizesPortal =
this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.LocalBlockSizes.ReadPortal();
for (std::size_t bi = 0; bi < static_cast<std::size_t>(input.GetNumberOfPartitions()); bi++)
{
// create the local contour tree mesh
@ -279,20 +279,25 @@ VTKM_CONT void ContourTreeAugmented::DoPostExecute(const vtkm::cont::Partitioned
auto currBlock = input.GetPartition(static_cast<vtkm::Id>(bi));
auto currField =
currBlock.GetField(this->GetActiveFieldName(), this->GetActiveFieldAssociation());
vtkm::Id3 pointDimensions, globalPointDimensions, globalPointIndexStart;
currBlock.GetCellSet().CastAndCallForTypes<VTKM_DEFAULT_CELL_SET_LIST_STRUCTURED>(
vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(),
pointDimensions,
globalPointDimensions,
globalPointIndexStart);
//const vtkm::cont::ArrayHandle<T,StorageType> &fieldData = currField.GetData().Cast<vtkm::cont::ArrayHandle<T,StorageType> >();
vtkm::cont::ArrayHandle<T> fieldData;
vtkm::cont::ArrayCopy(currField.GetData(), fieldData);
auto currContourTreeMesh = vtkm::worklet::contourtree_distributed::MultiBlockContourTreeHelper::
ComputeLocalContourTreeMesh<T>(
this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.LocalBlockOrigins.ReadPortal()
.Get(static_cast<vtkm::Id>(bi)),
this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.LocalBlockSizes.ReadPortal().Get(
static_cast<vtkm::Id>(bi)),
this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.GlobalSize,
fieldData,
MultiBlockTreeHelper->LocalContourTrees[bi],
MultiBlockTreeHelper->LocalSortOrders[bi],
compRegularStruct);
ComputeLocalContourTreeMesh<T>(globalPointIndexStart,
pointDimensions,
globalPointDimensions,
fieldData,
MultiBlockTreeHelper->LocalContourTrees[bi],
MultiBlockTreeHelper->LocalSortOrders[bi],
compRegularStruct);
localContourTreeMeshes[bi] = currContourTreeMesh;
// create the local data block structure
localDataBlocks[bi] = new vtkm::worklet::contourtree_distributed::ContourTreeBlockData<T>();
@ -303,10 +308,9 @@ VTKM_CONT void ContourTreeAugmented::DoPostExecute(const vtkm::cont::Partitioned
localDataBlocks[bi]->NeighborConnectivity = currContourTreeMesh->NeighborConnectivity;
localDataBlocks[bi]->NeighborOffsets = currContourTreeMesh->NeighborOffsets;
localDataBlocks[bi]->MaxNeighbors = currContourTreeMesh->MaxNeighbors;
localDataBlocks[bi]->BlockOrigin = localBlocksOriginPortal.Get(static_cast<vtkm::Id>(bi));
localDataBlocks[bi]->BlockSize = localBlocksSizesPortal.Get(static_cast<vtkm::Id>(bi));
localDataBlocks[bi]->GlobalSize =
this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.GlobalSize;
localDataBlocks[bi]->BlockOrigin = globalPointIndexStart;
localDataBlocks[bi]->BlockSize = pointDimensions;
localDataBlocks[bi]->GlobalSize = globalPointDimensions;
// We need to augment at least with the boundary vertices when running in parallel
localDataBlocks[bi]->ComputeRegularStructure = compRegularStruct;
}
@ -320,33 +324,32 @@ VTKM_CONT void ContourTreeAugmented::DoPostExecute(const vtkm::cont::Partitioned
// Compute the gids for our local blocks
using RegularDecomposer = vtkmdiy::RegularDecomposer<vtkmdiy::DiscreteBounds>;
const vtkm::filter::scalar_topology::internal::SpatialDecomposition& spatialDecomp =
this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition;
const auto numDims = spatialDecomp.NumberOfDimensions();
// ... division vector
RegularDecomposer::DivisionsVector diyDivisions(numDims);
for (vtkm::IdComponent d = 0;
d < static_cast<vtkm::IdComponent>(spatialDecomp.NumberOfDimensions());
++d)
RegularDecomposer::DivisionsVector diyDivisions;
std::vector<int> vtkmdiyLocalBlockGids;
vtkmdiy::DiscreteBounds diyBounds(0);
if (this->MultiBlockTreeHelper->BlocksPerDimension[0] == -1)
{
diyDivisions[d] = static_cast<int>(spatialDecomp.BlocksPerDimension[d]);
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
"BlocksPerDimension not set. Computing block indices "
"from information in CellSetStructured.");
diyBounds = vtkm::filter::scalar_topology::internal::ComputeBlockIndices(
input, diyDivisions, vtkmdiyLocalBlockGids);
}
// ... coordinates of local blocks
auto localBlockIndicesPortal = spatialDecomp.LocalBlockIndices.ReadPortal();
std::vector<vtkm::Id> vtkmdiyLocalBlockGids(static_cast<size_t>(input.GetNumberOfPartitions()));
for (vtkm::Id bi = 0; bi < input.GetNumberOfPartitions(); bi++)
else
{
RegularDecomposer::DivisionsVector diyCoords(static_cast<size_t>(numDims));
auto currentCoords = localBlockIndicesPortal.Get(bi);
for (vtkm::IdComponent d = 0; d < numDims; ++d)
{
diyCoords[d] = static_cast<int>(currentCoords[d]);
}
vtkmdiyLocalBlockGids[static_cast<size_t>(bi)] =
RegularDecomposer::coords_to_gid(diyCoords, diyDivisions);
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
"BlocksPerDimension set. Using information provided by caller.");
diyBounds = vtkm::filter::scalar_topology::internal::ComputeBlockIndices(
input,
this->MultiBlockTreeHelper->BlocksPerDimension,
this->MultiBlockTreeHelper->LocalBlockIndices,
diyDivisions,
vtkmdiyLocalBlockGids);
}
int numDims = diyBounds.min.dimension();
int globalNumberOfBlocks =
std::accumulate(diyDivisions.cbegin(), diyDivisions.cend(), 1, std::multiplies<int>{});
// Add my local blocks to the vtkmdiy master.
for (std::size_t bi = 0; bi < static_cast<std::size_t>(input.GetNumberOfPartitions()); bi++)
@ -361,16 +364,15 @@ VTKM_CONT void ContourTreeAugmented::DoPostExecute(const vtkm::cont::Partitioned
RegularDecomposer::BoolVector wrap(3, false);
RegularDecomposer::CoordinateVector ghosts(3, 1);
RegularDecomposer decomposer(static_cast<int>(numDims),
spatialDecomp.GetVTKmDIYBounds(),
static_cast<int>(spatialDecomp.GetGlobalNumberOfBlocks()),
diyBounds,
globalNumberOfBlocks,
shareFace,
wrap,
ghosts,
diyDivisions);
// Define which blocks live on which rank so that vtkmdiy can manage them
vtkmdiy::DynamicAssigner assigner(
comm, static_cast<int>(size), static_cast<int>(spatialDecomp.GetGlobalNumberOfBlocks()));
vtkmdiy::DynamicAssigner assigner(comm, static_cast<int>(size), globalNumberOfBlocks);
for (vtkm::Id bi = 0; bi < input.GetNumberOfPartitions(); bi++)
{
assigner.set_rank(static_cast<int>(rank),
@ -394,6 +396,13 @@ VTKM_CONT void ContourTreeAugmented::DoPostExecute(const vtkm::cont::Partitioned
if (rank == 0)
{
vtkm::Id3 dummy1, globalPointDimensions, dummy2;
vtkm::cont::DataSet firstDS = input.GetPartition(0);
firstDS.GetCellSet().CastAndCallForTypes<VTKM_DEFAULT_CELL_SET_LIST_STRUCTURED>(
vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(),
dummy1,
globalPointDimensions,
dummy2);
// Now run the contour tree algorithm on the last block to compute the final tree
vtkm::Id currNumIterations;
vtkm::worklet::contourtree_augmented::ContourTree currContourTree;
@ -412,12 +421,12 @@ VTKM_CONT void ContourTreeAugmented::DoPostExecute(const vtkm::cont::Partitioned
contourTreeMeshOut.MaxNeighbors = localDataBlocks[0]->MaxNeighbors;
// Construct the mesh boundary exectuion object needed for boundary augmentation
vtkm::Id3 minIdx(0, 0, 0);
vtkm::Id3 maxIdx = this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.GlobalSize;
vtkm::Id3 maxIdx = globalPointDimensions;
maxIdx[0] = maxIdx[0] - 1;
maxIdx[1] = maxIdx[1] - 1;
maxIdx[2] = maxIdx[2] > 0 ? (maxIdx[2] - 1) : 0;
auto meshBoundaryExecObj = contourTreeMeshOut.GetMeshBoundaryExecutionObject(
this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.GlobalSize, minIdx, maxIdx);
auto meshBoundaryExecObj =
contourTreeMeshOut.GetMeshBoundaryExecutionObject(globalPointDimensions, minIdx, maxIdx);
// Run the worklet to compute the final contour tree
worklet.Run(
contourTreeMeshOut.SortedValues, // Unused param. Provide something to keep API happy

21
vtkm/filter/scalar_topology/ContourTreeUniformAugmented.h
View File

@ -114,19 +114,24 @@ public:
///
/// Note: Only used when running on a multi-block dataset.
/// @param[in] blocksPerDim Number of data blocks used in each data dimension
/// @param[in] globalSize Global extends of the input mesh (i.e., number of mesh points in each dimension)
/// @param[in] localBlockIndices Array with the (x,y,z) index of each local data block with
/// with respect to blocksPerDim
/// @param[in] localBlockOrigins Array with the (x,y,z) origin (with regard to mesh index) of each
/// local data block
/// @param[in] localBlockSizes Array with the sizes (i.e., extends in number of mesh points) of each
/// local data block
VTKM_CONT
void SetBlockIndices(vtkm::Id3 blocksPerDim,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockIndices);
VTKM_CONT
VTKM_DEPRECATED(1.9,
"Set PointSize, GlobalPointOrigin, and GlobalPointSize in CellSetStructured and "
"optionally use SetBlockIndices.")
void SetSpatialDecomposition(vtkm::Id3 blocksPerDim,
vtkm::Id3 globalSize,
vtkm::Id3,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockIndices,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockOrigins,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockSizes);
const vtkm::cont::ArrayHandle<vtkm::Id3>&,
const vtkm::cont::ArrayHandle<vtkm::Id3>&)
{
SetBlockIndices(blocksPerDim, localBlockIndices);
}
//@{
/// Get the contour tree computed by the filter

234
vtkm/filter/scalar_topology/ContourTreeUniformDistributed.cxx
View File

@ -51,6 +51,7 @@
//==============================================================================
// vtkm includes
#include <vtkm/cont/CellSetStructured.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/Timer.h>
@ -61,7 +62,7 @@
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/meshtypes/mesh_boundary/MeshBoundaryContourTreeMesh.h>
// distributed contour tree includes
#include <vtkm/filter/scalar_topology/internal/SpatialDecomposition.h>
#include <vtkm/filter/scalar_topology/internal/ComputeBlockIndices.h>
#include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/BoundaryTree.h>
#include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/BoundaryTreeMaker.h>
#include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/CombineHyperSweepBlockFunctor.h>
@ -165,10 +166,10 @@ void SaveHierarchicalTreeDot(
//-----------------------------------------------------------------------------
ContourTreeUniformDistributed::ContourTreeUniformDistributed(
vtkm::Id3 blocksPerDim,
vtkm::Id3 globalSize,
vtkm::Id3, // globalSize, -> Now in CellSetStructured
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockIndices,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockOrigins,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockSizes,
const vtkm::cont::ArrayHandle<vtkm::Id3>&, // localBlockOrigins, -> Use from CellSetStructured
const vtkm::cont::ArrayHandle<vtkm::Id3>&, // localBlockSizes, -> Use from CellSetStructured
vtkm::cont::LogLevel timingsLogLevel,
vtkm::cont::LogLevel treeLogLevel)
: UseBoundaryExtremaOnly(true)
@ -177,19 +178,33 @@ ContourTreeUniformDistributed::ContourTreeUniformDistributed(
, SaveDotFiles(false)
, TimingsLogLevel(timingsLogLevel)
, TreeLogLevel(treeLogLevel)
, MultiBlockSpatialDecomposition(blocksPerDim,
globalSize,
localBlockIndices,
localBlockOrigins,
localBlockSizes)
, LocalMeshes(static_cast<std::size_t>(localBlockSizes.GetNumberOfValues()))
, LocalContourTrees(static_cast<std::size_t>(localBlockSizes.GetNumberOfValues()))
, LocalBoundaryTrees(static_cast<std::size_t>(localBlockSizes.GetNumberOfValues()))
, LocalInteriorForests(static_cast<std::size_t>(localBlockSizes.GetNumberOfValues()))
, BlocksPerDimension(blocksPerDim)
, LocalBlockIndices(localBlockIndices)
, LocalMeshes()
, LocalContourTrees()
, LocalBoundaryTrees()
, LocalInteriorForests()
{
this->SetOutputFieldName("resultData");
}
ContourTreeUniformDistributed::ContourTreeUniformDistributed(vtkm::cont::LogLevel timingsLogLevel,
vtkm::cont::LogLevel treeLogLevel)
: UseBoundaryExtremaOnly(true)
, UseMarchingCubes(false)
, AugmentHierarchicalTree(false)
, SaveDotFiles(false)
, TimingsLogLevel(timingsLogLevel)
, TreeLogLevel(treeLogLevel)
, BlocksPerDimension(vtkm::Id3{ -1, -1, -1 })
, LocalBlockIndices()
, LocalMeshes()
, LocalContourTrees()
, LocalBoundaryTrees()
, LocalInteriorForests()
{
this->SetOutputFieldName("resultData");
}
//-----------------------------------------------------------------------------
// Functions used in PrepareForExecution() to compute the local contour
@ -237,7 +252,7 @@ void ContourTreeUniformDistributed::ComputeLocalTree(
template <typename T, typename StorageType, typename MeshType, typename MeshBoundaryExecType>
void ContourTreeUniformDistributed::ComputeLocalTreeImpl(
const vtkm::Id blockIndex,
const vtkm::cont::DataSet&, // input,
const vtkm::cont::DataSet& ds, // input,
const vtkm::cont::ArrayHandle<T, StorageType>& field,
MeshType& mesh,
MeshBoundaryExecType& meshBoundaryExecObject)
@ -278,10 +293,14 @@ void ContourTreeUniformDistributed::ComputeLocalTreeImpl(
// Create an IdRelabeler since we are using a DataSetMesh type here, we don't need
// the IdRelabeler for the BRACT construction when we are using a ContourTreeMesh.
vtkm::Id3 pointDimensions, globalPointDimensions, globalPointIndexStart;
ds.GetCellSet().CastAndCallForTypes<VTKM_DEFAULT_CELL_SET_LIST_STRUCTURED>(
vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(),
pointDimensions,
globalPointDimensions,
globalPointIndexStart);
auto localToGlobalIdRelabeler = vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler(
this->MultiBlockSpatialDecomposition.LocalBlockOrigins.ReadPortal().Get(blockIndex),
this->MultiBlockSpatialDecomposition.LocalBlockSizes.ReadPortal().Get(blockIndex),
this->MultiBlockSpatialDecomposition.GlobalSize);
globalPointIndexStart, pointDimensions, globalPointDimensions);
// Initialize the BoundaryTreeMaker
auto boundaryTreeMaker =
vtkm::worklet::contourtree_distributed::BoundaryTreeMaker<MeshType, MeshBoundaryExecType>(
@ -322,9 +341,9 @@ void ContourTreeUniformDistributed::ComputeLocalTreeImpl(
"Before Fan In",
mesh,
field,
this->MultiBlockSpatialDecomposition.LocalBlockOrigins.ReadPortal().Get(blockIndex),
this->MultiBlockSpatialDecomposition.LocalBlockSizes.ReadPortal().Get(blockIndex),
this->MultiBlockSpatialDecomposition.GlobalSize);
globalPointIndexStart,
pointDimensions,
globalPointDimensions);
bractFile << bractString << std::endl;
}
@ -461,18 +480,40 @@ vtkm::cont::DataSet ContourTreeUniformDistributed::DoExecute(const vtkm::cont::D
VTKM_CONT void ContourTreeUniformDistributed::PreExecute(
const vtkm::cont::PartitionedDataSet& input)
{
if (vtkm::filter::scalar_topology::internal::SpatialDecomposition::GetGlobalNumberOfBlocks(
input) != this->MultiBlockSpatialDecomposition.GetGlobalNumberOfBlocks())
// TODO/FIXME: The following may be too expensive for a "sanity" check as it
// requires global communication
auto globalNumberOfPartitions = input.GetGlobalNumberOfPartitions();
if (globalNumberOfPartitions < 2)
{
throw vtkm::cont::ErrorFilterExecution(
"Global number of blocks in MultiBlock dataset does not match the SpatialDecomposition");
throw vtkm::cont::ErrorFilterExecution("ContourTreeUniformDistributed filter expects a "
"PartitionedDataSet with at least two partitions.");
}
if (this->MultiBlockSpatialDecomposition.GetLocalNumberOfBlocks() !=
input.GetNumberOfPartitions())
if (this->BlocksPerDimension[0] != -1)
{
throw vtkm::cont::ErrorFilterExecution(
"Local number of blocks in MultiBlock dataset does not match the SpatialDecomposition");
if (this->BlocksPerDimension[1] < 1 || this->BlocksPerDimension[2] < 1)
{
throw vtkm::cont::ErrorFilterExecution("Invalid input BlocksPerDimension.");
}
if (globalNumberOfPartitions !=
this->BlocksPerDimension[0] * this->BlocksPerDimension[1] * this->BlocksPerDimension[2])
{
throw vtkm::cont::ErrorFilterExecution("Global number of blocks in data set does not match "
"expected value based on BlocksPerDimension");
}
if (this->LocalBlockIndices.GetNumberOfValues() != input.GetNumberOfPartitions())
{
throw vtkm::cont::ErrorFilterExecution("Local number of partitions in data set does not "
"match number of specified blocks indices.");
}
}
// Allocate vectors
this->LocalMeshes.resize(static_cast<std::size_t>(input.GetGlobalNumberOfPartitions()));
this->LocalContourTrees.resize(static_cast<std::size_t>(input.GetGlobalNumberOfPartitions()));
this->LocalBoundaryTrees.resize(static_cast<std::size_t>(input.GetGlobalNumberOfPartitions()));
this->LocalInteriorForests.resize(static_cast<std::size_t>(input.GetGlobalNumberOfPartitions()));
}
vtkm::cont::PartitionedDataSet ContourTreeUniformDistributed::DoExecutePartitions(
@ -546,15 +587,6 @@ VTKM_CONT void ContourTreeUniformDistributed::PostExecute(
vtkm::cont::Timer timer;
timer.Start();
// We are running in parallel and need to merge the contour tree in PostExecute
// TODO/FIXME: This filter should only be used in a parallel setting with more
// than one block. Is there a better way to enforce this? thrown an exception
// instead of an empty return? What is the appropriate exception?
if (this->MultiBlockSpatialDecomposition.GetGlobalNumberOfBlocks() == 1)
{
return;
}
auto field = // TODO/FIXME: Correct for more than one block per rank?
input.GetPartition(0).GetField(this->GetActiveFieldName(), this->GetActiveFieldAssociation());
@ -602,12 +634,15 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric(
inputContourTreeMaster.foreach (
[&](DistributedContourTreeBlockData* currInBlock, const vtkmdiy::Master::ProxyWithLink&) {
vtkm::Id blockNo = currInBlock->LocalBlockNo;
const vtkm::cont::DataSet& currDS = hierarchicalTreeOutputDataSet[blockNo];
// The block size and origin may be modified during the FanIn so we need to use the
// size and origin from the original decomposition instead of looking it up in the currInBlock
auto currBlockSize =
this->MultiBlockSpatialDecomposition.LocalBlockSizes.ReadPortal().Get(blockNo);
auto currBlockOrigin =
this->MultiBlockSpatialDecomposition.LocalBlockOrigins.ReadPortal().Get(blockNo);
vtkm::Id3 pointDimensions, globalPointDimensions, globalPointIndexStart;
currDS.GetCellSet().CastAndCallForTypes<VTKM_DEFAULT_CELL_SET_LIST_STRUCTURED>(
vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(),
pointDimensions,
globalPointDimensions,
globalPointIndexStart);
// NOTE: Use dummy link to make DIY happy. The dummy link is never used, since all
// communication is via RegularDecomposer, which sets up its own links. No need
@ -618,9 +653,9 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric(
currInBlock->GlobalBlockId,
new HyperSweepBlock(blockNo,
currInBlock->GlobalBlockId,
currBlockOrigin,
currBlockSize,
this->MultiBlockSpatialDecomposition.GlobalSize,
globalPointIndexStart,
pointDimensions,
globalPointDimensions,
*currInBlock->HierarchicalAugmenter.AugmentedTree),
new vtkmdiy::Link());
});
@ -774,6 +809,7 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric(
});
}
//-----------------------------------------------------------------------------
template <typename FieldType>
VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
@ -809,30 +845,32 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
// 1.1.2 Compute the gids for our local blocks
using RegularDecomposer = vtkmdiy::RegularDecomposer<vtkmdiy::DiscreteBounds>;
const auto& spatialDecomp = this->MultiBlockSpatialDecomposition;
const auto numDims = spatialDecomp.NumberOfDimensions();
// ... compute division vector for global domain
RegularDecomposer::DivisionsVector diyDivisions(numDims);
for (vtkm::IdComponent d = 0; d < static_cast<vtkm::IdComponent>(numDims); ++d)
RegularDecomposer::DivisionsVector diyDivisions;
std::vector<int> vtkmdiyLocalBlockGids;
vtkmdiy::DiscreteBounds diyBounds(0);
if (this->BlocksPerDimension[0] == -1)
{
diyDivisions[d] = static_cast<int>(spatialDecomp.BlocksPerDimension[d]);
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
"BlocksPerDimension not set. Computing block indices "
"from information in CellSetStructured.");
diyBounds = vtkm::filter::scalar_topology::internal::ComputeBlockIndices(
input, diyDivisions, vtkmdiyLocalBlockGids);
}
// ... compute coordinates of local blocks
auto localBlockIndicesPortal = spatialDecomp.LocalBlockIndices.ReadPortal();
std::vector<int> vtkmdiyLocalBlockGids(static_cast<size_t>(input.GetNumberOfPartitions()));
for (vtkm::Id bi = 0; bi < input.GetNumberOfPartitions(); bi++)
else
{
RegularDecomposer::DivisionsVector diyCoords(static_cast<size_t>(numDims));
auto currentCoords = localBlockIndicesPortal.Get(bi);
for (vtkm::IdComponent d = 0; d < numDims; ++d)
{
diyCoords[d] = static_cast<int>(currentCoords[d]);
}
vtkmdiyLocalBlockGids[static_cast<size_t>(bi)] =
RegularDecomposer::coords_to_gid(diyCoords, diyDivisions);
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
"BlocksPerDimension set. Using information provided by caller.");
diyBounds =
vtkm::filter::scalar_topology::internal::ComputeBlockIndices(input,
this->BlocksPerDimension,
this->LocalBlockIndices,
diyDivisions,
vtkmdiyLocalBlockGids);
}
int numDims = diyBounds.min.dimension();
int globalNumberOfBlocks =
std::accumulate(diyDivisions.cbegin(), diyDivisions.cend(), 1, std::multiplies<int>{});
// Record time to compute the local block ids
timingsStream << " " << std::setw(38) << std::left << "Compute Block Ids and Local Links"
@ -840,16 +878,28 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
timer.Start();
// 1.1.3 Setup the block data for DIY and add it to master
// Note: globalPointDimensions is defined outside the loop since it is needed later.
// It may be set multiple times in the loop, but always to the same value.
vtkm::Id3 globalPointDimensions;
for (vtkm::Id bi = 0; bi < input.GetNumberOfPartitions(); bi++)
{
// Get the input block and associated cell set information
auto currBlock = input.GetPartition(static_cast<vtkm::Id>(bi));
vtkm::Id3 pointDimensions, globalPointIndexStart;
currBlock.GetCellSet().CastAndCallForTypes<VTKM_DEFAULT_CELL_SET_LIST_STRUCTURED>(
vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(),
pointDimensions,
globalPointDimensions,
globalPointIndexStart);
// Create the local data block structure and set extents
auto newBlock = new DistributedContourTreeBlockData();
// Copy global block id into the local data block for use in the hierarchical augmentation
newBlock->GlobalBlockId = vtkmdiyLocalBlockGids[bi];
newBlock->LocalBlockNo = bi;
newBlock->BlockOrigin = spatialDecomp.LocalBlockOrigins.ReadPortal().Get(bi);
newBlock->BlockSize = spatialDecomp.LocalBlockSizes.ReadPortal().Get(bi);
newBlock->BlockOrigin = globalPointIndexStart;
newBlock->BlockSize = pointDimensions;
// Save local tree information for fan out; TODO/FIXME: Try to avoid copy
newBlock->ContourTrees.push_back(this->LocalContourTrees[bi]);
@ -868,15 +918,10 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
auto transformedIndex = vtkm::cont::make_ArrayHandleTransform(
permutedSortOrder,
vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler(
this->MultiBlockSpatialDecomposition.LocalBlockOrigins.ReadPortal().Get(
static_cast<vtkm::Id>(bi)),
this->MultiBlockSpatialDecomposition.LocalBlockSizes.ReadPortal().Get(
static_cast<vtkm::Id>(bi)),
this->MultiBlockSpatialDecomposition.GlobalSize));
globalPointIndexStart, pointDimensions, globalPointDimensions));
vtkm::cont::Algorithm::Copy(transformedIndex, localGlobalMeshIndex);
// ... get data values
auto currBlock = input.GetPartition(static_cast<vtkm::Id>(bi));
auto currField =
currBlock.GetField(this->GetActiveFieldName(), this->GetActiveFieldAssociation());
vtkm::cont::ArrayHandle<FieldType> fieldData;
@ -942,16 +987,15 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
RegularDecomposer::BoolVector wrap(3, false);
RegularDecomposer::CoordinateVector ghosts(3, 1);
RegularDecomposer decomposer(static_cast<int>(numDims),
spatialDecomp.GetVTKmDIYBounds(),
static_cast<int>(spatialDecomp.GetGlobalNumberOfBlocks()),
diyBounds,
globalNumberOfBlocks,
shareFace,
wrap,
ghosts,
diyDivisions);
// Define which blocks live on which rank so that vtkmdiy can manage them
vtkmdiy::DynamicAssigner assigner(
comm, static_cast<int>(size), static_cast<int>(spatialDecomp.GetGlobalNumberOfBlocks()));
vtkmdiy::DynamicAssigner assigner(comm, static_cast<int>(size), globalNumberOfBlocks);
for (vtkm::Id bi = 0; bi < input.GetNumberOfPartitions(); bi++)
{
assigner.set_rank(static_cast<int>(rank), vtkmdiyLocalBlockGids[static_cast<size_t>(bi)]);
@ -984,7 +1028,7 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
timer.Start();
// 1.3 Perform fan-in reduction
const vtkm::worklet::contourtree_distributed::ComputeDistributedContourTreeFunctor<FieldType>
computeDistributedContourTreeFunctor(this->MultiBlockSpatialDecomposition.GlobalSize,
computeDistributedContourTreeFunctor(globalPointDimensions,
this->UseBoundaryExtremaOnly,
this->TimingsLogLevel,
this->TreeLogLevel);
@ -1061,17 +1105,22 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
grafter(&(this->LocalMeshes[static_cast<std::size_t>(blockData->LocalBlockNo)]),
blockData->ContourTrees[0],
&(blockData->InteriorForests[0]));
auto currBlock = input.GetPartition(blockData->LocalBlockNo);
vtkm::cont::DataSet currBlock = input.GetPartition(blockData->LocalBlockNo);
auto currField =
currBlock.GetField(this->GetActiveFieldName(), this->GetActiveFieldAssociation());
vtkm::cont::ArrayHandle<FieldType> fieldData;
vtkm::cont::ArrayCopy(currField.GetData(), fieldData);
vtkm::Id3 pointDimensions, globalPointIndexStart;
// globalPointDimensions already defined in parent scope
currBlock.GetCellSet().CastAndCallForTypes<VTKM_DEFAULT_CELL_SET_LIST_STRUCTURED>(
vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(),
pointDimensions,
globalPointDimensions,
globalPointIndexStart);
auto localToGlobalIdRelabeler = vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler(
this->MultiBlockSpatialDecomposition.LocalBlockOrigins.ReadPortal().Get(
blockData->LocalBlockNo),
this->MultiBlockSpatialDecomposition.LocalBlockSizes.ReadPortal().Get(
blockData->LocalBlockNo),
this->MultiBlockSpatialDecomposition.GlobalSize);
globalPointIndexStart, pointDimensions, globalPointDimensions);
grafter.GraftInteriorForests(
0, blockData->HierarchicalTree, fieldData, &localToGlobalIdRelabeler);
@ -1155,9 +1204,30 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
// Add the information to the output data set
blockHierarchcialTree.AddToVTKMDataSet(hierarchicalTreeOutputDataSet[blockData->LocalBlockNo]);
// Save information required to set up DIY
vtkm::cont::ArrayHandle<vtkm::Id> vtkmGlobalBlockIdAH;
vtkmGlobalBlockIdAH.Allocate(1);
auto vtkmGlobalBlockIdWP = vtkmGlobalBlockIdAH.WritePortal();
vtkmGlobalBlockIdWP.Set(0, blockData->GlobalBlockId);
vtkm::cont::Field vtkmGlobalBlockIdField(
"vtkmGlobalBlockId", vtkm::cont::Field::Association::WholeMesh, vtkmGlobalBlockIdAH);
hierarchicalTreeOutputDataSet[blockData->LocalBlockNo].AddField(vtkmGlobalBlockIdField);
vtkm::cont::ArrayHandle<vtkm::Id> vtkmBlocksPerDimensionAH;
vtkmBlocksPerDimensionAH.Allocate(3);
auto vtkmBlocksPerDimensionWP = vtkmBlocksPerDimensionAH.WritePortal();
vtkmBlocksPerDimensionWP.Set(0, this->BlocksPerDimension[0]);
vtkmBlocksPerDimensionWP.Set(1, this->BlocksPerDimension[1]);
vtkmBlocksPerDimensionWP.Set(2, this->BlocksPerDimension[2]);
vtkm::cont::Field vtkmBlocksPerDimensionField("vtkmBlocksPerDimension",
vtkm::cont::Field::Association::WholeMesh,
vtkmBlocksPerDimensionAH);
hierarchicalTreeOutputDataSet[blockData->LocalBlockNo].AddField(vtkmBlocksPerDimensionField);
// Copy cell set from input data set. This is mainly to ensure that the output data set
// has a defined cell set. Without one, serialization for DIY does not work properly.
// Having the extents of the input data set may also help in other use cases.
// For example, the ComputeVolume method gets information from this cell set as will
// the branch decomposition filter.
hierarchicalTreeOutputDataSet[blockData->LocalBlockNo].SetCellSet(
input.GetPartition(blockData->LocalBlockNo).GetCellSet());

52
vtkm/filter/scalar_topology/ContourTreeUniformDistributed.h
View File

@ -56,7 +56,6 @@
#include <vtkm/Types.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/filter/scalar_topology/internal/SpatialDecomposition.h>
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/ContourTree.h>
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/DataSetMesh.h>
#include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/BoundaryTree.h>
@ -113,14 +112,20 @@ public:
/// @param[in] treeLogLevel Set the vtkm::cont:LogLevel to be used to record metadata information
/// about the various trees computed as part of the hierarchical contour tree compute
VTKM_CONT
ContourTreeUniformDistributed(
vtkm::Id3 blocksPerDim, // TODO/FIXME: Possibly pass SpatialDecomposition object instead
vtkm::Id3 globalSize,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockIndices,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockOrigins,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockSizes,
vtkm::cont::LogLevel timingsLogLevel = vtkm::cont::LogLevel::Perf,
vtkm::cont::LogLevel treeLogLevel = vtkm::cont::LogLevel::Info);
VTKM_DEPRECATED(
1.9,
"Use default constructor and set PointSize, GlobalPointIndexStart, and GlobalPointSize in "
"CellSetStructured. Optionally use `SetBlockIndices` accessor (if information is available).")
ContourTreeUniformDistributed(vtkm::Id3 blocksPerDim,
vtkm::Id3 globalSize,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockIndices,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockOrigins,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockSizes,
vtkm::cont::LogLevel timingsLogLevel = vtkm::cont::LogLevel::Perf,
vtkm::cont::LogLevel treeLogLevel = vtkm::cont::LogLevel::Info);
ContourTreeUniformDistributed(vtkm::cont::LogLevel timingsLogLevel = vtkm::cont::LogLevel::Perf,
vtkm::cont::LogLevel treeLogLevel = vtkm::cont::LogLevel::Info);
VTKM_CONT void SetUseBoundaryExtremaOnly(bool useBoundaryExtremaOnly)
{
@ -141,6 +146,13 @@ public:
this->AugmentHierarchicalTree = augmentHierarchicalTree;
}
VTKM_CONT void SetBlockIndices(vtkm::Id3 blocksPerDim,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockIndices)
{
this->BlocksPerDimension = blocksPerDim;
vtkm::cont::ArrayCopy(localBlockIndices, this->LocalBlockIndices);
}
VTKM_CONT bool GetAugmentHierarchicalTree() { return this->AugmentHierarchicalTree; }
VTKM_CONT void SetSaveDotFiles(bool saveDotFiles) { this->SaveDotFiles = saveDotFiles; }
@ -213,8 +225,11 @@ private:
/// Log level to be used for outputting metadata about the trees. Default is vtkm::cont::LogLevel::Info
vtkm::cont::LogLevel TreeLogLevel = vtkm::cont::LogLevel::Info;
/// Information about the spatial decomposition
vtkm::filter::scalar_topology::internal::SpatialDecomposition MultiBlockSpatialDecomposition;
/// Information about block decomposition TODO/FIXME: Remove need for this information
// ... Number of blocks along each dimension
vtkm::Id3 BlocksPerDimension;
// ... Index of the local blocks in x,y,z, i.e., in i,j,k mesh coordinates
vtkm::cont::ArrayHandle<vtkm::Id3> LocalBlockIndices;
/// Intermediate results (one per local data block)...
/// ... local mesh information needed at end of fan out
@ -240,28 +255,23 @@ class VTKM_DEPRECATED(1.8, "Use vtkm::filter::scalar_topology::ContourTreeUnifor
using scalar_topology::ContourTreeUniformDistributed::ContourTreeUniformDistributed;
ContourTreeUniformDistributed(vtkm::Id3 blocksPerDim,
vtkm::Id3 globalSize,
vtkm::Id3,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockIndices,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockOrigins,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockSizes,
const vtkm::cont::ArrayHandle<vtkm::Id3>&,
const vtkm::cont::ArrayHandle<vtkm::Id3>&,
bool useBoundaryExtremaOnly = true,
bool useMarchingCubes = false,
bool augmentHierarchicalTree = false,
bool saveDotFiles = false,
vtkm::cont::LogLevel timingsLogLevel = vtkm::cont::LogLevel::Perf,
vtkm::cont::LogLevel treeLogLevel = vtkm::cont::LogLevel::Info)
: vtkm::filter::scalar_topology::ContourTreeUniformDistributed(blocksPerDim,
globalSize,
localBlockIndices,
localBlockOrigins,
localBlockSizes,
timingsLogLevel,
treeLogLevel)
: vtkm::filter::scalar_topology::ContourTreeUniformDistributed(timingsLogLevel, treeLogLevel)
{
this->SetUseBoundaryExtremaOnly(useBoundaryExtremaOnly);
this->SetUseMarchingCubes(useMarchingCubes);
this->SetAugmentHierarchicalTree(augmentHierarchicalTree);
this->SetSaveDotFiles(saveDotFiles);
this->SetBlockIndices(blocksPerDim, localBlockIndices);
this->SetOutputFieldName("resultData");
}
};

83
vtkm/filter/scalar_topology/DistributedBranchDecompositionFilter.cxx
View File

@ -36,16 +36,11 @@ namespace scalar_topology
// not need to pass it sperately (or check if it can already be derived from
// information stored in PartitionedDataSet)
VTKM_CONT DistributedBranchDecompositionFilter::DistributedBranchDecompositionFilter(
vtkm::Id3 blocksPerDim,
vtkm::Id3 globalSize,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockIndices,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockOrigins,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockSizes)
: MultiBlockSpatialDecomposition(blocksPerDim,
globalSize,
localBlockIndices,
localBlockOrigins,
localBlockSizes)
vtkm::Id3,
vtkm::Id3,
const vtkm::cont::ArrayHandle<vtkm::Id3>&,
const vtkm::cont::ArrayHandle<vtkm::Id3>&,
const vtkm::cont::ArrayHandle<vtkm::Id3>&)
{
}
@ -80,52 +75,59 @@ VTKM_CONT vtkm::cont::PartitionedDataSet DistributedBranchDecompositionFilter::D
BranchDecompositionBlock::Destroy);
timingsStream << " " << std::setw(60) << std::left
<< "Create DIY Master (Branch Decomposition)"
<< "Create DIY Master and Assigner (Branch Decomposition)"
<< ": " << timer.GetElapsedTime() << " seconds" << std::endl;
timer.Start();
// Compute global ids (gids) for our local blocks
using RegularDecomposer = vtkmdiy::RegularDecomposer<vtkmdiy::DiscreteBounds>;
int globalNumberOfBlocks =
static_cast<int>(this->MultiBlockSpatialDecomposition.GetGlobalNumberOfBlocks());
int numDims = static_cast<int>(this->MultiBlockSpatialDecomposition.NumberOfDimensions());
// TODO/FIXME: Is there a better way to set this up?
auto firstDS = input.GetPartition(0);
vtkm::Id3 firstPointDimensions, firstGlobalPointDimensions, firstGlobalPointIndexStart;
firstDS.GetCellSet().CastAndCallForTypes<VTKM_DEFAULT_CELL_SET_LIST_STRUCTURED>(
vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(),
firstPointDimensions,
firstGlobalPointDimensions,
firstGlobalPointIndexStart);
int numDims = firstGlobalPointDimensions[2] > 1 ? 3 : 2;
auto vtkmBlocksPerDimensionRP = input.GetPartition(0)
.GetField("vtkmBlocksPerDimension")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Id>>()
.ReadPortal();
// ... compute division vector for global domain
using RegularDecomposer = vtkmdiy::RegularDecomposer<vtkmdiy::DiscreteBounds>;
RegularDecomposer::DivisionsVector diyDivisions(numDims);
vtkmdiy::DiscreteBounds diyBounds(numDims);
int globalNumberOfBlocks = 1;
for (vtkm::IdComponent d = 0; d < static_cast<vtkm::IdComponent>(numDims); ++d)
{
diyDivisions[d] = static_cast<int>(this->MultiBlockSpatialDecomposition.BlocksPerDimension[d]);
}
// ... compute coordinates of local blocks
std::vector<int> vtkmdiyLocalBlockGids(static_cast<size_t>(input.GetNumberOfPartitions()));
auto localBlockIndicesPortal =
this->MultiBlockSpatialDecomposition.LocalBlockIndices.ReadPortal();
for (vtkm::Id bi = 0; bi < input.GetNumberOfPartitions(); bi++)
{
RegularDecomposer::DivisionsVector diyCoords(static_cast<size_t>(numDims));
auto currentCoords = localBlockIndicesPortal.Get(bi);
for (vtkm::IdComponent d = 0; d < numDims; ++d)
{
diyCoords[d] = static_cast<int>(currentCoords[d]);
}
vtkmdiyLocalBlockGids[static_cast<size_t>(bi)] =
RegularDecomposer::coords_to_gid(diyCoords, diyDivisions);
diyDivisions[d] = static_cast<int>(vtkmBlocksPerDimensionRP.Get(d));
globalNumberOfBlocks *= diyDivisions[d];
diyBounds.min[d] = 0;
diyBounds.max[d] = static_cast<int>(firstGlobalPointDimensions[d]);
}
// Record time to compute the local block ids
timingsStream << " " << std::setw(60) << std::left
<< "Compute Block Ids and Local Links (Branch Decomposition)"
<< "Get DIY Information (Branch Decomposition)"
<< ": " << timer.GetElapsedTime() << " seconds" << std::endl;
timer.Start();
// Initialize branch decomposition computation from data in PartitionedDataSet blocks
vtkmdiy::DynamicAssigner assigner(comm, size, globalNumberOfBlocks);
for (vtkm::Id localBlockIndex = 0; localBlockIndex < input.GetNumberOfPartitions();
++localBlockIndex)
{
int globalBlockId = vtkmdiyLocalBlockGids[localBlockIndex];
const vtkm::cont::DataSet& ds = input.GetPartition(localBlockIndex);
int globalBlockId = static_cast<int>(
vtkm::cont::ArrayGetValue(0,
ds.GetField("vtkmGlobalBlockId")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Id>>()));
BranchDecompositionBlock* newBlock =
new BranchDecompositionBlock(localBlockIndex, globalBlockId, ds);
// NOTE: Use dummy link to make DIY happy. The dummy link is never used, since all
@ -134,6 +136,9 @@ VTKM_CONT vtkm::cont::PartitionedDataSet DistributedBranchDecompositionFilter::D
// NOTE: Since we passed a "Destroy" function to DIY master, it will own the local data
// blocks and delete them when done.
branch_decomposition_master.add(globalBlockId, newBlock, new vtkmdiy::Link());
// Tell assigner that this block lives on this rank so that DIY can manage blocks
assigner.set_rank(rank, globalBlockId);
}
// Log time to copy the data to the HyperSweepBlock data objects
@ -146,19 +151,15 @@ VTKM_CONT vtkm::cont::PartitionedDataSet DistributedBranchDecompositionFilter::D
RegularDecomposer::BoolVector wrap(3, false);
RegularDecomposer::CoordinateVector ghosts(3, 1);
RegularDecomposer decomposer(numDims,
this->MultiBlockSpatialDecomposition.GetVTKmDIYBounds(),
diyBounds,
static_cast<int>(globalNumberOfBlocks),
shareFace,
wrap,
ghosts,
diyDivisions);
// Define which blocks live on which rank so that vtkmdiy can manage them
vtkmdiy::DynamicAssigner assigner(comm, size, globalNumberOfBlocks);
for (vtkm::Id bi = 0; bi < input.GetNumberOfPartitions(); bi++)
{
assigner.set_rank(rank, vtkmdiyLocalBlockGids[static_cast<size_t>(bi)]);
}
timingsStream << " " << std::setw(60) << std::left
@ -196,7 +197,7 @@ VTKM_CONT vtkm::cont::PartitionedDataSet DistributedBranchDecompositionFilter::D
ds.GetField("DependentVolume").GetData().AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Id>>();
// Get global size and compute total volume from it
const auto& globalSize = this->MultiBlockSpatialDecomposition.GlobalSize;
const auto& globalSize = firstGlobalPointDimensions;
vtkm::Id totalVolume = globalSize[0] * globalSize[1] * globalSize[2];
// Compute local best up and down paths by volume

17
vtkm/filter/scalar_topology/DistributedBranchDecompositionFilter.h
View File

@ -45,8 +45,6 @@
#include <vtkm/filter/NewFilterField.h>
#include <vtkm/filter/scalar_topology/vtkm_filter_scalar_topology_export.h>
#include <vtkm/filter/scalar_topology/internal/SpatialDecomposition.h>
namespace vtkm
{
namespace filter
@ -59,20 +57,17 @@ class VTKM_FILTER_SCALAR_TOPOLOGY_EXPORT DistributedBranchDecompositionFilter
: public vtkm::filter::NewFilter
{
public:
VTKM_CONT DistributedBranchDecompositionFilter(
vtkm::Id3 blocksPerDim,
vtkm::Id3 globalSize,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockIndices,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockOrigins,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockSizes);
VTKM_CONT DistributedBranchDecompositionFilter() = default;
VTKM_CONT DistributedBranchDecompositionFilter(vtkm::Id3,
vtkm::Id3,
const vtkm::cont::ArrayHandle<vtkm::Id3>&,
const vtkm::cont::ArrayHandle<vtkm::Id3>&,
const vtkm::cont::ArrayHandle<vtkm::Id3>&);
private:
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet&) override;
VTKM_CONT vtkm::cont::PartitionedDataSet DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& inData) override;
/// Information about the spatial decomposition
vtkm::filter::scalar_topology::internal::SpatialDecomposition MultiBlockSpatialDecomposition;
};
} // namespace scalar_topology

2
vtkm/filter/scalar_topology/internal/CMakeLists.txt
View File

@ -10,8 +10,8 @@
set(headers
BranchDecompositionBlock.h
ComputeBlockIndices.h
ComputeDistributedBranchDecompositionFunctor.h
SpatialDecomposition.h
)
#-----------------------------------------------------------------------------

255
vtkm/filter/scalar_topology/internal/ComputeBlockIndices.cxx Normal file
View File

@ -0,0 +1,255 @@
//============================================================================
// 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) 2018, The Regents of the University of California, through
// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
// from the U.S. Dept. of Energy). All rights reserved.
//
// 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 the University of California, Lawrence Berkeley National
// Laboratory, U.S. Dept. of Energy 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 THE COPYRIGHT HOLDERS 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 THE COPYRIGHT OWNER 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.
//
//=============================================================================
//
// This code is an extension of the algorithm presented in the paper:
// Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
// Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
// Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
// (LDAV), October 2016, Baltimore, Maryland.
//
// The PPP2 algorithm and software were jointly developed by
// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
// Oliver Ruebel (LBNL)
//==============================================================================
#include <algorithm>
#include <vtkm/cont/CellSetStructured.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/EnvironmentTracker.h>
#include <vtkm/filter/scalar_topology/internal/ComputeBlockIndices.h>
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/meshtypes/contourtreemesh/CopyIntoCombinedArrayWorklet.h>
namespace // anonymous namespace for local helper classes
{
struct OriginsBlock
{
std::vector<int> Origins;
OriginsBlock(const std::vector<int>& origins)
: Origins(origins)
{
}
static void Destroy(void* b) { delete static_cast<OriginsBlock*>(b); }
};
struct MergeOriginsFunctor
{
void operator()(OriginsBlock* b,
const vtkmdiy::ReduceProxy& rp, // communication proxy
const vtkmdiy::RegularSwapPartners& // partners of the current block (unused)
) const
{
const auto selfid = rp.gid();
std::vector<int> incoming;
rp.incoming(incoming);
for (const int ingid : incoming)
{
if (ingid != selfid)
{
std::vector<int> incoming_origins;
rp.dequeue(ingid, incoming_origins);
std::vector<int> merged_origins;
std::merge(incoming_origins.begin(),
incoming_origins.end(),
b->Origins.begin(),
b->Origins.end(),
std::inserter(merged_origins, merged_origins.begin()));
auto last = std::unique(merged_origins.begin(), merged_origins.end());
merged_origins.erase(last, merged_origins.end());
std::swap(merged_origins, b->Origins);
}
}
for (int cc = 0; cc < rp.out_link().size(); ++cc)
{
auto target = rp.out_link().target(cc);
if (target.gid != selfid)
{
rp.enqueue(target, b->Origins);
}
}
}
};
}
namespace vtkm
{
namespace filter
{
namespace scalar_topology
{
namespace internal
{
VTKM_CONT vtkmdiy::DiscreteBounds ComputeBlockIndices(const vtkm::cont::PartitionedDataSet& input,
DiscreteBoundsDivisionVector& diyDivisions,
std::vector<int>& diyLocalBlockGids)
{
auto firstDS = input.GetPartition(0);
vtkm::Id3 dummy1, firstGlobalPointDimensions, dummy2;
firstDS.GetCellSet().CastAndCallForTypes<VTKM_DEFAULT_CELL_SET_LIST_STRUCTURED>(
vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(),
dummy1,
firstGlobalPointDimensions,
dummy2);
int numDims = firstGlobalPointDimensions[2] > 1 ? 3 : 2;
diyDivisions.clear();
vtkmdiy::DiscreteBounds diyBounds(numDims);
std::vector<DiscreteBoundsDivisionVector> diyBlockCoords(input.GetNumberOfPartitions());
for (vtkm::IdComponent d = 0; d < numDims; ++d)
{
// Set bounds for this dimension
diyBounds.min[d] = 0;
diyBounds.max[d] = static_cast<int>(firstGlobalPointDimensions[d]);
// Get the list of origins along current coordinate axis
std::vector<int> local_origins;
for (vtkm::Id ds_no = 0; ds_no < input.GetNumberOfPartitions(); ++ds_no)
{
vtkm::Id3 globalPointIndexStart;
input.GetPartition(ds_no)
.GetCellSet()
.CastAndCallForTypes<VTKM_DEFAULT_CELL_SET_LIST_STRUCTURED>(
vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(),
dummy1,
dummy2,
globalPointIndexStart);
local_origins.push_back(static_cast<int>(globalPointIndexStart[d]));
}
// Sort and remove duplicates
OriginsBlock* origins_block = new OriginsBlock(local_origins);
std::sort(origins_block->Origins.begin(), origins_block->Origins.end());
auto last = std::unique(origins_block->Origins.begin(), origins_block->Origins.end());
origins_block->Origins.erase(last, origins_block->Origins.end());
// Create global list of origins across all ranks
auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
auto rank = comm.rank();
auto size = comm.size();
vtkmdiy::Master master(comm, 1, -1, 0, OriginsBlock::Destroy);
master.add(rank, origins_block, new vtkmdiy::Link);
vtkmdiy::ContiguousAssigner assigner(size, size);
vtkmdiy::DiscreteBounds bounds(1);
bounds.min[0] = 0;
bounds.max[0] = size - 1;
vtkmdiy::RegularDecomposer<vtkmdiy::DiscreteBounds> decomposer(1, bounds, size);
vtkmdiy::RegularSwapPartners partners(decomposer, 2, true);
vtkmdiy::reduce(master, assigner, partners, MergeOriginsFunctor{});
// Number of blocks/divisions along axis is number of unique origins along this axis
diyDivisions.push_back(static_cast<int>(origins_block->Origins.size()));
// Block index aling this axis is the index of the origin in that list
for (vtkm::Id ds_no = 0; ds_no < input.GetNumberOfPartitions(); ++ds_no)
{
diyBlockCoords[ds_no].push_back(static_cast<int>(std::find(origins_block->Origins.begin(),
origins_block->Origins.end(),
local_origins[ds_no]) -
origins_block->Origins.begin()));
}
}
// Compute global block IDs
diyLocalBlockGids.clear();
for (vtkm::Id ds_no = 0; ds_no < input.GetNumberOfPartitions(); ++ds_no)
{
diyLocalBlockGids.push_back(vtkmdiy::RegularDecomposer<vtkmdiy::DiscreteBounds>::coords_to_gid(
diyBlockCoords[ds_no], diyDivisions));
}
return diyBounds;
}
VTKM_CONT vtkmdiy::DiscreteBounds ComputeBlockIndices(
const vtkm::cont::PartitionedDataSet& input,
vtkm::Id3 blocksPerDim,
const vtkm::cont::ArrayHandle<vtkm::Id3>& blockIndices,
DiscreteBoundsDivisionVector& diyDivisions,
std::vector<int>& diyLocalBlockGids)
{
auto firstDS = input.GetPartition(0);
vtkm::Id3 dummy1, firstGlobalPointDimensions, dummy2;
firstDS.GetCellSet().CastAndCallForTypes<VTKM_DEFAULT_CELL_SET_LIST_STRUCTURED>(
vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(),
dummy1,
firstGlobalPointDimensions,
dummy2);
int numDims = firstGlobalPointDimensions[2] > 1 ? 3 : 2;
diyDivisions.clear();
vtkmdiy::DiscreteBounds diyBounds(numDims);
for (vtkm::IdComponent d = 0; d < numDims; ++d)
{
// Set bounds for this dimension
diyBounds.min[d] = 0;
diyBounds.max[d] = static_cast<int>(firstGlobalPointDimensions[d]);
diyDivisions.push_back(static_cast<int>(blocksPerDim[d]));
}
// Compute global block IDs
diyLocalBlockGids.clear();
auto blockIndicesPortal = blockIndices.ReadPortal();
for (vtkm::Id ds_no = 0; ds_no < input.GetNumberOfPartitions(); ++ds_no)
{
auto currBlockIndices = blockIndicesPortal.Get(ds_no);
DiscreteBoundsDivisionVector diyBlockCoords(numDims);
for (vtkm::IdComponent d = 0; d < numDims; ++d)
{
diyBlockCoords[d] = static_cast<int>(currBlockIndices[d]);
}
diyLocalBlockGids.push_back(vtkmdiy::RegularDecomposer<vtkmdiy::DiscreteBounds>::coords_to_gid(
diyBlockCoords, diyDivisions));
}
return diyBounds;
}
} // namespace internal
} // namespace scalar_topology
} // namespace filter
} // namespace vtkm

98
vtkm/filter/scalar_topology/internal/ComputeBlockIndices.h Normal file
View File

@ -0,0 +1,98 @@
//============================================================================
// 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) 2018, The Regents of the University of California, through
// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
// from the U.S. Dept. of Energy). All rights reserved.
//
// 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 the University of California, Lawrence Berkeley National
// Laboratory, U.S. Dept. of Energy 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 THE COPYRIGHT HOLDERS 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 THE COPYRIGHT OWNER 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.
//
//=============================================================================
//
// This code is an extension of the algorithm presented in the paper:
// Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
// Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
// Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
// (LDAV), October 2016, Baltimore, Maryland.
//
// The PPP2 algorithm and software were jointly developed by
// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
// Oliver Ruebel (LBNL)
//==============================================================================
#ifndef vtk_m_filter_scalar_topology_internal_ComputeBlockIndices_h_
#define vtk_m_filter_scalar_topology_internal_ComputeBlockIndices_h_
#include <vector>
#include <vtkm/Types.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/PartitionedDataSet.h>
// DIY includes
// clang-format off
VTKM_THIRDPARTY_PRE_INCLUDE
#include <vtkm/thirdparty/diy/Configure.h>
#include <vtkm/thirdparty/diy/diy.h>
VTKM_THIRDPARTY_POST_INCLUDE
// clang-format on
namespace vtkm
{
namespace filter
{
namespace scalar_topology
{
namespace internal
{
using DiscreteBoundsDivisionVector =
vtkmdiy::RegularDecomposer<vtkmdiy::DiscreteBounds>::DivisionsVector;
VTKM_CONT vtkmdiy::DiscreteBounds ComputeBlockIndices(const vtkm::cont::PartitionedDataSet& input,
DiscreteBoundsDivisionVector& diyDivisions,
std::vector<int>& diyLocalBlockGids);
VTKM_CONT vtkmdiy::DiscreteBounds ComputeBlockIndices(
const vtkm::cont::PartitionedDataSet& input,
vtkm::Id3 blocksPerDim,
const vtkm::cont::ArrayHandle<vtkm::Id3>& blockIndices,
DiscreteBoundsDivisionVector& diyDivisions,
std::vector<int>& diyLocalBlockGids);
} // namespace internal
} // namespace scalar_topology
} // namespace filter
} // namespace vtkm
#endif

172
vtkm/filter/scalar_topology/internal/SpatialDecomposition.h
View File

@ -1,172 +0,0 @@
//============================================================================
// 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) 2018, The Regents of the University of California, through
// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
// from the U.S. Dept. of Energy). All rights reserved.
//
// 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 the University of California, Lawrence Berkeley National
// Laboratory, U.S. Dept. of Energy 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 THE COPYRIGHT HOLDERS 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 THE COPYRIGHT OWNER 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.
//
//=============================================================================
//
// This code is an extension of the algorithm presented in the paper:
// Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
// Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
// Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
// (LDAV), October 2016, Baltimore, Maryland.
//
// The PPP2 algorithm and software were jointly developed by
// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
// Oliver Ruebel (LBNL)
//==============================================================================
#ifndef vtk_m_filter_scalar_topology_internal_SpatialDecomposition_h
#define vtk_m_filter_scalar_topology_internal_SpatialDecomposition_h
#include <vtkm/Types.h>
#include <vtkm/cont/BoundsCompute.h>
#include <vtkm/cont/BoundsGlobalCompute.h>
#include <vtkm/cont/EnvironmentTracker.h>
#include <vtkm/cont/PartitionedDataSet.h>
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h>
// clang-format off
VTKM_THIRDPARTY_PRE_INCLUDE
#include <vtkm/thirdparty/diy/diy.h>
VTKM_THIRDPARTY_POST_INCLUDE
// clang-format on
namespace vtkm
{
namespace filter
{
namespace scalar_topology
{
namespace internal
{
// --- Helper class to store the spatial decomposition defined by the PartitionedDataSet input data
class SpatialDecomposition
{
public:
VTKM_CONT
SpatialDecomposition(vtkm::Id3 blocksPerDim,
vtkm::Id3 globalSize,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockIndices,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockOrigins,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockSizes)
: BlocksPerDimension(blocksPerDim)
, GlobalSize(globalSize)
, LocalBlockIndices(localBlockIndices)
, LocalBlockOrigins(localBlockOrigins)
, LocalBlockSizes(localBlockSizes)
{
}
inline vtkmdiy::DiscreteBounds GetVTKmDIYBounds() const
{
if (this->NumberOfDimensions() == 2)
{
vtkmdiy::DiscreteBounds domain(2);
domain.min[0] = domain.min[1] = 0;
domain.max[0] = static_cast<int>(this->GlobalSize[0]);
domain.max[1] = static_cast<int>(this->GlobalSize[1]);
return domain;
}
else
{
vtkmdiy::DiscreteBounds domain(3);
domain.min[0] = domain.min[1] = domain.min[2] = 0;
domain.max[0] = static_cast<int>(this->GlobalSize[0]);
domain.max[1] = static_cast<int>(this->GlobalSize[1]);
domain.max[2] = static_cast<int>(this->GlobalSize[2]);
return domain;
}
}
inline vtkm::Id NumberOfDimensions() const { return GlobalSize[2] > 1 ? 3 : 2; }
inline vtkm::Id GetGlobalNumberOfBlocks() const
{
return BlocksPerDimension[0] * BlocksPerDimension[1] * BlocksPerDimension[2];
}
inline vtkm::Id GetLocalNumberOfBlocks() const { return LocalBlockSizes.GetNumberOfValues(); }
inline static vtkm::Bounds GetGlobalBounds(const vtkm::cont::PartitionedDataSet& input)
{
// Get the spatial bounds of a multi -block data set
vtkm::Bounds bounds = vtkm::cont::BoundsGlobalCompute(input);
return bounds;
}
inline static vtkm::Bounds GetLocalBounds(const vtkm::cont::PartitionedDataSet& input)
{
// Get the spatial bounds of a multi -block data set
vtkm::Bounds bounds = vtkm::cont::BoundsCompute(input);
return bounds;
}
inline static vtkm::Id GetGlobalNumberOfBlocks(const vtkm::cont::PartitionedDataSet& input)
{
auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
vtkm::Id localSize = input.GetNumberOfPartitions();
vtkm::Id globalSize = 0;
#ifdef VTKM_ENABLE_MPI
vtkmdiy::mpi::all_reduce(comm, localSize, globalSize, std::plus<vtkm::Id>{});
#else
globalSize = localSize;
#endif
return globalSize;
}
// Number of blocks along each dimension
vtkm::Id3 BlocksPerDimension;
// Size of the global mesh
vtkm::Id3 GlobalSize;
// Index of the local blocks in x,y,z, i.e., in i,j,k mesh coordinates
vtkm::cont::ArrayHandle<vtkm::Id3> LocalBlockIndices;
// Origin of the local blocks in mesh index space
vtkm::cont::ArrayHandle<vtkm::Id3> LocalBlockOrigins;
// Size of each local block in x, y,z
vtkm::cont::ArrayHandle<vtkm::Id3> LocalBlockSizes;
};
} // internal
} // namespace scalar_topology
} // namespace filter
} // namespace vtkm
#endif

49
vtkm/filter/scalar_topology/testing/TestingContourTreeUniformDistributedFilter.h
View File

@ -203,12 +203,22 @@ inline vtkm::cont::DataSet CreateSubDataSet(const vtkm::cont::DataSet& ds,
{
vtkm::Id2 dimensions{ blockSize[0], blockSize[1] };
vtkm::cont::DataSet dataSet = dsb.Create(dimensions);
vtkm::cont::CellSetStructured<2> cellSet;
cellSet.SetPointDimensions(dimensions);
cellSet.SetGlobalPointDimensions(vtkm::Id2{ globalSize[0], globalSize[1] });
cellSet.SetGlobalPointIndexStart(vtkm::Id2{ blockOrigin[0], blockOrigin[1] });
dataSet.SetCellSet(cellSet);
dataSet.AddField(permutedField);
return dataSet;
}
else
{
vtkm::cont::DataSet dataSet = dsb.Create(blockSize);
vtkm::cont::CellSetStructured<3> cellSet;
cellSet.SetPointDimensions(blockSize);
cellSet.SetGlobalPointDimensions(globalSize);
cellSet.SetGlobalPointIndexStart(blockOrigin);
dataSet.SetCellSet(cellSet);
dataSet.AddField(permutedField);
return dataSet;
}
@ -241,7 +251,8 @@ inline vtkm::cont::PartitionedDataSet RunContourTreeDUniformDistributed(
int numberOfRanks,
bool augmentHierarchicalTree,
bool computeHierarchicalVolumetricBranchDecomposition,
vtkm::Id3& globalSize)
vtkm::Id3& globalSize,
bool passBlockIndices = true)
{
// Get dimensions of data set
vtkm::cont::CastAndCall(
@ -267,15 +278,9 @@ inline vtkm::cont::PartitionedDataSet RunContourTreeDUniformDistributed(
// Created partitioned (split) data set
vtkm::cont::PartitionedDataSet pds;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockIndices;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockOrigins;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockSizes;
localBlockIndices.Allocate(blocksOnThisRank);
localBlockOrigins.Allocate(blocksOnThisRank);
localBlockSizes.Allocate(blocksOnThisRank);
auto localBlockIndicesPortal = localBlockIndices.WritePortal();
auto localBlockOriginsPortal = localBlockOrigins.WritePortal();
auto localBlockSizesPortal = localBlockSizes.WritePortal();
for (vtkm::Id blockNo = 0; blockNo < blocksOnThisRank; ++blockNo)
{
@ -283,20 +288,17 @@ inline vtkm::cont::PartitionedDataSet RunContourTreeDUniformDistributed(
std::tie(blockIndex, blockOrigin, blockSize) =
ComputeBlockExtents(globalSize, blocksPerAxis, startBlockNo + blockNo);
pds.AppendPartition(CreateSubDataSet(ds, blockOrigin, blockSize, fieldName));
localBlockOriginsPortal.Set(blockNo, blockOrigin);
localBlockSizesPortal.Set(blockNo, blockSize);
localBlockIndicesPortal.Set(blockNo, blockIndex);
}
// Run the contour tree analysis
vtkm::filter::scalar_topology::ContourTreeUniformDistributed filter(
blocksPerAxis,
globalSize,
localBlockIndices,
localBlockOrigins,
localBlockSizes,
vtkm::cont::LogLevel::UserVerboseLast,
vtkm::cont::LogLevel::UserVerboseLast);
vtkm::cont::LogLevel::UserVerboseLast, vtkm::cont::LogLevel::UserVerboseLast);
if (passBlockIndices)
{
filter.SetBlockIndices(blocksPerAxis, localBlockIndices);
}
filter.SetUseMarchingCubes(useMarchingCubes);
// Freudenthal: Only use boundary extrema; MC: use all points on boundary
@ -310,8 +312,7 @@ inline vtkm::cont::PartitionedDataSet RunContourTreeDUniformDistributed(
{
using vtkm::filter::scalar_topology::DistributedBranchDecompositionFilter;
DistributedBranchDecompositionFilter bd_filter(
blocksPerAxis, globalSize, localBlockIndices, localBlockOrigins, localBlockSizes);
DistributedBranchDecompositionFilter bd_filter;
result = bd_filter.Execute(result);
}
@ -387,7 +388,8 @@ inline vtkm::cont::PartitionedDataSet RunContourTreeDUniformDistributed(
int rank = 0,
int numberOfRanks = 1,
bool augmentHierarchicalTree = false,
bool computeHierarchicalVolumetricBranchDecomposition = false)
bool computeHierarchicalVolumetricBranchDecomposition = false,
bool passBlockIndices = true)
{
vtkm::Id3 globalSize;
@ -399,7 +401,8 @@ inline vtkm::cont::PartitionedDataSet RunContourTreeDUniformDistributed(
numberOfRanks,
augmentHierarchicalTree,
computeHierarchicalVolumetricBranchDecomposition,
globalSize);
globalSize,
passBlockIndices);
}
inline void TestContourTreeUniformDistributed8x9(int nBlocks, int rank = 0, int size = 1)
@ -575,7 +578,8 @@ inline void TestContourTreeFile(std::string ds_filename,
int rank = 0,
int size = 1,
bool augmentHierarchicalTree = false,
bool computeHierarchicalVolumetricBranchDecomposition = false)
bool computeHierarchicalVolumetricBranchDecomposition = false,
bool passBlockIndices = true)
{
if (rank == 0)
{
@ -611,7 +615,8 @@ inline void TestContourTreeFile(std::string ds_filename,
size,
augmentHierarchicalTree,
computeHierarchicalVolumetricBranchDecomposition,
globalSize);
globalSize,
passBlockIndices);
if (rank == 0)
{

10
vtkm/filter/scalar_topology/testing/UnitTestContourTreeUniformDistributedFilter.cxx
View File

@ -112,6 +112,16 @@ public:
1,
true,
false);
TestContourTreeFile(Testing::DataPath("rectilinear/vanc.vtk"),
"var",
Testing::RegressionImagePath("vanc.augment_hierarchical_tree.ct_txt"),
4,
false,
0,
1,
true,
false,
false);
}
};
}

41
vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h
View File

@ -328,6 +328,47 @@ struct GetPointDimensions
}
};
struct GetLocalAndGlobalPointDimensions
{
void operator()(const vtkm::cont::CellSetStructured<2>& cells,
vtkm::Id3& pointDimensions,
vtkm::Id3& globalPointDimensions,
vtkm::Id3& globalPointIndexStart) const
{
vtkm::Id2 pointDimensions2D = cells.GetPointDimensions();
pointDimensions[0] = pointDimensions2D[0];
pointDimensions[1] = pointDimensions2D[1];
pointDimensions[2] = 1;
vtkm::Id2 globalPointDimensions2D = cells.GetGlobalPointDimensions();
globalPointDimensions[0] = globalPointDimensions2D[0];
globalPointDimensions[1] = globalPointDimensions2D[1];
globalPointDimensions[2] = 1;
vtkm::Id2 pointIndexStart2D = cells.GetGlobalPointIndexStart();
globalPointIndexStart[0] = pointIndexStart2D[0];
globalPointIndexStart[1] = pointIndexStart2D[1];
globalPointIndexStart[2] = 0;
}
void operator()(const vtkm::cont::CellSetStructured<3>& cells,
vtkm::Id3& pointDimensions,
vtkm::Id3& globalPointDimensions,
vtkm::Id3& globalPointIndexStart) const
{
pointDimensions = cells.GetPointDimensions();
globalPointDimensions = cells.GetGlobalPointDimensions();
globalPointIndexStart = cells.GetGlobalPointIndexStart();
}
//@}
/// Raise ErrorBadValue if the input cell set is not a vtkm::cont::CellSetStructured<2> or <3>
template <typename T>
void operator()(const T&, vtkm::Id3&, vtkm::Id3&, vtkm::Id3&) const
{
throw vtkm::cont::ErrorBadValue("Expected 2D or 3D structured cell cet! ");
}
};
} // namespace contourtree_augmented
} // worklet
} // vtkm

49
vtkm/filter/scalar_topology/worklet/contourtree_distributed/MultiBlockContourTreeHelper.h
View File

@ -53,15 +53,12 @@
#ifndef vtk_m_worklet_contourtree_distributed_multiblockcontourtreehelper_h
#define vtk_m_worklet_contourtree_distributed_multiblockcontourtreehelper_h
#include <vtkm/filter/scalar_topology/internal/SpatialDecomposition.h>
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h>
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/meshtypes/ContourTreeMesh.h>
#include <vtkm/Types.h>
#include <vtkm/cont/BoundsCompute.h>
#include <vtkm/cont/BoundsGlobalCompute.h>
//#include <vtkm/cont/AssignerPartitionedDataSet.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/PartitionedDataSet.h>
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/data_set_mesh/IdRelabeler.h>
@ -79,19 +76,21 @@ class MultiBlockContourTreeHelper
public:
VTKM_CONT
MultiBlockContourTreeHelper(vtkm::Id3 blocksPerDim,
vtkm::Id3 globalSize,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockIndices,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockOrigins,
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockSizes)
: MultiBlockSpatialDecomposition(blocksPerDim,
globalSize,
localBlockIndices,
localBlockOrigins,
localBlockSizes)
const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockIndices)
: BlocksPerDimension(blocksPerDim)
, LocalBlockIndices(localBlockIndices)
, LocalContourTrees(static_cast<size_t>(localBlockIndices.GetNumberOfValues()))
, LocalSortOrders(static_cast<size_t>(localBlockIndices.GetNumberOfValues()))
{
}
VTKM_CONT
MultiBlockContourTreeHelper(const vtkm::cont::PartitionedDataSet& input)
: BlocksPerDimension(-1, -1, -1)
, LocalBlockIndices()
, LocalContourTrees(static_cast<size_t>(input.GetNumberOfPartitions()))
, LocalSortOrders(static_cast<size_t>(input.GetNumberOfPartitions()))
{
vtkm::Id localNumBlocks = this->GetLocalNumberOfBlocks();
LocalContourTrees.resize(static_cast<std::size_t>(localNumBlocks));
LocalSortOrders.resize(static_cast<std::size_t>(localNumBlocks));
}
VTKM_CONT
@ -118,25 +117,12 @@ public:
inline vtkm::Id GetLocalNumberOfBlocks() const
{
return this->MultiBlockSpatialDecomposition.GetLocalNumberOfBlocks();
return static_cast<vtkm::Id>(this->LocalContourTrees.size());
}
inline vtkm::Id GetGlobalNumberOfBlocks() const
{
return this->MultiBlockSpatialDecomposition.GetGlobalNumberOfBlocks();
}
inline static vtkm::Id GetGlobalNumberOfBlocks(const vtkm::cont::PartitionedDataSet& input)
{
auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
vtkm::Id localSize = input.GetNumberOfPartitions();
vtkm::Id globalSize = 0;
#ifdef VTKM_ENABLE_MPI
vtkmdiy::mpi::all_reduce(comm, localSize, globalSize, std::plus<vtkm::Id>{});
#else
globalSize = localSize;
#endif
return globalSize;
return this->BlocksPerDimension[0] * this->BlocksPerDimension[1] * this->BlocksPerDimension[2];
}
// Used to compute the local contour tree mesh in after DoExecute. I.e., the function is
@ -198,7 +184,8 @@ public:
}
}
vtkm::filter::scalar_topology::internal::SpatialDecomposition MultiBlockSpatialDecomposition;
vtkm::Id3 BlocksPerDimension;
vtkm::cont::ArrayHandle<vtkm::Id3> LocalBlockIndices;
std::vector<vtkm::worklet::contourtree_augmented::ContourTree> LocalContourTrees;
std::vector<vtkm::worklet::contourtree_augmented::IdArrayType> LocalSortOrders;
}; // end MultiBlockContourTreeHelper

52
vtkm/internal/ConnectivityStructuredInternals.h
View File

@ -37,15 +37,24 @@ public:
VTKM_EXEC_CONT
void SetPointDimensions(vtkm::Id dimensions) { this->PointDimensions = dimensions; }
VTKM_EXEC_CONT
void SetGlobalPointDimensions(vtkm::Id dimensions) { this->GlobalPointDimensions = dimensions; }
VTKM_EXEC_CONT
void SetGlobalPointIndexStart(vtkm::Id start) { this->GlobalPointIndexStart = start; }
VTKM_EXEC_CONT
vtkm::Id GetPointDimensions() const { return this->PointDimensions; }
VTKM_EXEC_CONT
vtkm::Id GetGlobalPointDimensions() const { return this->GlobalPointDimensions; }
VTKM_EXEC_CONT
vtkm::Id GetCellDimensions() const { return this->PointDimensions - 1; }
VTKM_EXEC_CONT
vtkm::Id GetGlobalCellDimensions() const { return this->GlobalPointDimensions - 1; }
VTKM_EXEC_CONT
SchedulingRangeType GetSchedulingRange(vtkm::TopologyElementTagCell) const
{
@ -138,12 +147,15 @@ public:
void PrintSummary(std::ostream& out) const
{
out << " UniformConnectivity<1> ";
out << "this->PointDimensions[" << this->PointDimensions << "] ";
out << "PointDimensions[" << this->PointDimensions << "] ";
out << "GlobalPointDimensions[" << this->GlobalPointDimensions << "] ";
out << "GlobalPointIndexStart[" << this->GlobalPointIndexStart << "] ";
out << "\n";
}
private:
vtkm::Id PointDimensions = 0;
vtkm::Id GlobalPointDimensions = 0;
vtkm::Id GlobalPointIndexStart = 0;
};
@ -157,15 +169,24 @@ public:
VTKM_EXEC_CONT
void SetPointDimensions(vtkm::Id2 dims) { this->PointDimensions = dims; }
VTKM_EXEC_CONT
void SetGlobalPointDimensions(vtkm::Id2 dims) { this->GlobalPointDimensions = dims; }
VTKM_EXEC_CONT
void SetGlobalPointIndexStart(vtkm::Id2 start) { this->GlobalPointIndexStart = start; }
VTKM_EXEC_CONT
const vtkm::Id2& GetPointDimensions() const { return this->PointDimensions; }
VTKM_EXEC_CONT
const vtkm::Id2& GetGlobalPointDimensions() const { return this->GlobalPointDimensions; }
VTKM_EXEC_CONT
vtkm::Id2 GetCellDimensions() const { return this->PointDimensions - vtkm::Id2(1); }
VTKM_EXEC_CONT
vtkm::Id2 GetGlobalCellDimensions() const { return this->GlobalPointDimensions - vtkm::Id2(1); }
VTKM_EXEC_CONT
vtkm::Id GetNumberOfPoints() const { return vtkm::ReduceProduct(this->GetPointDimensions()); }
@ -303,12 +324,18 @@ public:
void PrintSummary(std::ostream& out) const
{
out << " UniformConnectivity<2> ";
out << "pointDim[" << this->PointDimensions[0] << " " << this->PointDimensions[1] << "] ";
out << "PointDimensions[" << this->PointDimensions[0] << " " << this->PointDimensions[1]
<< "] ";
out << "GlobalPointDimensions[" << this->GlobalPointDimensions[0] << " "
<< this->GlobalPointDimensions[1] << "] ";
out << "GlobalPointIndexStart[" << this->GlobalPointIndexStart[0] << " "
<< this->GlobalPointIndexStart[1] << "] ";
out << std::endl;
}
private:
vtkm::Id2 PointDimensions = { 0, 0 };
vtkm::Id2 GlobalPointDimensions = { 0, 0 };
vtkm::Id2 GlobalPointIndexStart = { 0, 0 };
};
@ -327,15 +354,28 @@ public:
this->CellDim01 = (dims[0] - 1) * (dims[1] - 1);
}
VTKM_EXEC_CONT
void SetGlobalPointDimensions(vtkm::Id3 dims)
{
this->GlobalPointDimensions = dims;
this->GlobalCellDimensions = dims - vtkm::Id3(1);
}
VTKM_EXEC_CONT
void SetGlobalPointIndexStart(vtkm::Id3 start) { this->GlobalPointIndexStart = start; }
VTKM_EXEC_CONT
const vtkm::Id3& GetPointDimensions() const { return this->PointDimensions; }
VTKM_EXEC_CONT
const vtkm::Id3& GetGlobalPointDimensions() const { return this->GlobalPointDimensions; }
VTKM_EXEC_CONT
const vtkm::Id3& GetCellDimensions() const { return this->CellDimensions; }
VTKM_EXEC_CONT
const vtkm::Id3& GetGlobalCellDimensions() const { return this->GlobalCellDimensions; }
VTKM_EXEC_CONT
vtkm::Id GetNumberOfPoints() const { return vtkm::ReduceProduct(this->PointDimensions); }
@ -466,8 +506,12 @@ public:
void PrintSummary(std::ostream& out) const
{
out << " UniformConnectivity<3> ";
out << "pointDim[" << this->PointDimensions[0] << " " << this->PointDimensions[1] << " "
out << "PointDimensions[" << this->PointDimensions[0] << " " << this->PointDimensions[1] << " "
<< this->PointDimensions[2] << "] ";
out << "GlobalPointDimensions[" << this->GlobalPointDimensions[0] << " "
<< this->GlobalPointDimensions[1] << " " << this->GlobalPointDimensions[2] << "] ";
out << "GlobalPointIndexStart[" << this->GlobalPointIndexStart[0] << " "
<< this->GlobalPointIndexStart[1] << " " << this->GlobalPointIndexStart[2] << "] ";
out << std::endl;
}
@ -509,6 +553,8 @@ public:
private:
vtkm::Id3 PointDimensions = { 0, 0, 0 };
vtkm::Id3 GlobalPointDimensions = { 0, 0, 0 };
vtkm::Id3 GlobalCellDimensions = { 0, 0, 0 };
vtkm::Id3 GlobalPointIndexStart = { 0, 0, 0 };
vtkm::Id3 CellDimensions = { 0, 0, 0 };
vtkm::Id CellDim01 = 0;

32
vtkm/worklet/testing/UnitTestContourTreeUniformDistributed.cxx
View File

@ -55,7 +55,6 @@
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/scalar_topology/internal/SpatialDecomposition.h>
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/DataSetMesh.h>
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/PrintVectors.h>
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/meshtypes/ContourTreeMesh.h>
@ -154,7 +153,6 @@ void TestHierarchicalHyperSweeper()
"misc/8x9test_HierarchicalAugmentedTree_Block2.dat",
"misc/8x9test_HierarchicalAugmentedTree_Block3.dat" };
vtkm::Id3 globalSize{ 9, 8, 1 };
vtkm::Id3 blocksPerDim{ 2, 2, 1 };
vtkm::Id3 sizes[numBlocks] = { { 5, 4, 1 }, { 5, 5, 1 }, { 5, 4, 1 }, { 5, 5, 1 } };
vtkm::Id3 origins[numBlocks] = { { 0, 0, 0 }, { 0, 3, 0 }, { 4, 0, 0 }, { 4, 3, 0 } };
vtkm::Id3 blockIndices[numBlocks] = { { 0, 0, 0 }, { 0, 1, 0 }, { 1, 0, 0 }, { 1, 1, 0 } };
@ -174,14 +172,6 @@ void TestHierarchicalHyperSweeper()
vtkm::cont::make_ArrayHandle<vtkm::Id>({ 6, 9, 18, 24, 46, 72, 2 })
};
// Create spatial decomposition
vtkm::filter::scalar_topology::internal::SpatialDecomposition spatialDecomp(
blocksPerDim,
globalSize,
vtkm::cont::make_ArrayHandle(blockIndices, numBlocks, vtkm::CopyFlag::Off),
vtkm::cont::make_ArrayHandle(origins, numBlocks, vtkm::CopyFlag::Off),
vtkm::cont::make_ArrayHandle(sizes, numBlocks, vtkm::CopyFlag::Off));
// Load trees
vtkm::worklet::contourtree_distributed::HierarchicalContourTree<vtkm::FloatDefault>
hct[numBlocks];
@ -209,22 +199,21 @@ void TestHierarchicalHyperSweeper()
RegularDecomposer::CoordinateVector ghosts(3, 1);
RegularDecomposer::DivisionsVector diyDivisions{ 2, 2, 1 }; // HARDCODED FOR TEST
int numDims = static_cast<int>(globalSize[2] > 1 ? 3 : 2);
RegularDecomposer decomposer(numDims,
spatialDecomp.GetVTKmDIYBounds(),
static_cast<int>(spatialDecomp.GetGlobalNumberOfBlocks()),
shareFace,
wrap,
ghosts,
diyDivisions);
int numDims = 2;
vtkmdiy::DiscreteBounds diyBounds(2);
diyBounds.min[0] = diyBounds.min[1] = 0;
diyBounds.max[0] = static_cast<int>(globalSize[0]);
diyBounds.max[1] = static_cast<int>(globalSize[1]);
RegularDecomposer decomposer(
numDims, diyBounds, numBlocks, shareFace, wrap, ghosts, diyDivisions);
// ... coordinates of local blocks
auto localBlockIndicesPortal = spatialDecomp.LocalBlockIndices.ReadPortal();
std::vector<int> vtkmdiyLocalBlockGids(numBlocks);
for (vtkm::Id bi = 0; bi < numBlocks; bi++)
{
RegularDecomposer::DivisionsVector diyCoords(static_cast<size_t>(numDims));
auto currentCoords = localBlockIndicesPortal.Get(bi);
auto currentCoords = blockIndices[bi];
for (vtkm::IdComponent d = 0; d < numDims; ++d)
{
diyCoords[d] = static_cast<int>(currentCoords[d]);
@ -234,8 +223,7 @@ void TestHierarchicalHyperSweeper()
}
// Define which blocks live on which rank so that vtkmdiy can manage them
vtkmdiy::DynamicAssigner assigner(
comm, comm.size(), static_cast<int>(spatialDecomp.GetGlobalNumberOfBlocks()));
vtkmdiy::DynamicAssigner assigner(comm, comm.size(), numBlocks);
for (vtkm::Id bi = 0; bi < numBlocks; bi++)
{
assigner.set_rank(static_cast<int>(rank),