diff --git a/examples/contour_tree_augmented/ContourTreeApp.cxx b/examples/contour_tree_augmented/ContourTreeApp.cxx index 96aca16a8..2ac66b8a1 100644 --- a/examples/contour_tree_augmented/ContourTreeApp.cxx +++ b/examples/contour_tree_augmented/ContourTreeApp.cxx @@ -548,14 +548,8 @@ int main(int argc, char* argv[]) static_cast(dims[1]), static_cast(0)); vtkm::cont::ArrayHandle localBlockIndices; - vtkm::cont::ArrayHandle localBlockOrigins; - vtkm::cont::ArrayHandle 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 origin(0, blockIndex * blockSize); vtkm::Vec 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(dims[0]), 0)); + localBlockIndicesPortal.Set(localBlockIndex, vtkm::Id3(0, blockIndex, 0)); } // 3D data else { vtkm::Id3 vdims; - vdims[0] = static_cast(dims[1]); - vdims[1] = static_cast(dims[0]); + vdims[0] = static_cast(dims[0]); + vdims[1] = static_cast(dims[1]); vdims[2] = static_cast(currBlockSize); - vtkm::Vec origin(0, 0, (blockIndex * blockSize)); + vtkm::Vec origin(0, 0, blockIndex * blockSize); vtkm::Vec 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(dims[0]), static_cast(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 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"); diff --git a/examples/contour_tree_distributed/ContourTreeApp.cxx b/examples/contour_tree_distributed/ContourTreeApp.cxx index 2a0fabd47..a2cb3c4d7 100644 --- a/examples/contour_tree_distributed/ContourTreeApp.cxx +++ b/examples/contour_tree_distributed/ContourTreeApp.cxx @@ -401,14 +401,8 @@ int main(int argc, char* argv[]) vtkm::Id3 globalSize; vtkm::Id3 blocksPerDim; vtkm::cont::ArrayHandle localBlockIndices; - vtkm::cont::ArrayHandle localBlockOrigins; - vtkm::cont::ArrayHandle 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(offset[1]) }; const vtkm::Vec 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(offset[2]) }; vtkm::Vec 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(blockIndex[0]), static_cast(blockIndex[1]), static_cast(nDims == 3 ? blockIndex[2] : 0) }); - localBlockOriginsPortal.Set(blockNo, - vtkm::Id3{ static_cast(offset[0]), - static_cast(offset[1]), - static_cast(nDims == 3 ? offset[2] : 0) }); - localBlockSizesPortal.Set(blockNo, - vtkm::Id3{ static_cast(dims[0]), - static_cast(dims[1]), - static_cast(nDims == 3 ? dims[2] : 0) }); if (blockNo == 0) { @@ -753,7 +749,6 @@ int main(int argc, char* argv[]) : vtkm::Id3(static_cast(dims[0]), static_cast(dims[1]), static_cast(1)); - std::cout << blocksPerDim << " " << globalSize << std::endl; { vtkm::Id lastDimSize = (nDims == 2) ? static_cast(dims[1]) : static_cast(dims[2]); @@ -801,12 +796,12 @@ int main(int argc, char* argv[]) vtkm::Vec origin(0, blockIndex * blockSize); vtkm::Vec 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(dims[0]), currBlockSize, 0)); } // 3D data else @@ -818,14 +813,12 @@ int main(int argc, char* argv[]) vtkm::Vec origin(0, 0, (blockIndex * blockSize)); vtkm::Vec 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(dims[0]), - static_cast(dims[1]), - currBlockSize)); } std::vector 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) diff --git a/vtkm/cont/CellSetStructured.h b/vtkm/cont/CellSetStructured.h index 26475b46e..6399c4ee4 100644 --- a/vtkm/cont/CellSetStructured.h +++ b/vtkm/cont/CellSetStructured.h @@ -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); } }; diff --git a/vtkm/cont/PartitionedDataSet.cxx b/vtkm/cont/PartitionedDataSet.cxx index 9dbcb6e0d..761bde302 100644 --- a/vtkm/cont/PartitionedDataSet.cxx +++ b/vtkm/cont/PartitionedDataSet.cxx @@ -15,6 +15,13 @@ #include #include #include +#include + +// clang-format off +VTKM_THIRDPARTY_PRE_INCLUDE +#include +VTKM_THIRDPARTY_POST_INCLUDE +// clang-format on namespace vtkm { @@ -73,6 +80,19 @@ vtkm::Id PartitionedDataSet::GetNumberOfPartitions() const return static_cast(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{}); + return globalSize; +#else + return GetNumberOfPartitions(); +#endif +} + VTKM_CONT const vtkm::cont::DataSet& PartitionedDataSet::GetPartition(vtkm::Id blockId) const { diff --git a/vtkm/cont/PartitionedDataSet.h b/vtkm/cont/PartitionedDataSet.h index 7152939ac..e2f4f465a 100644 --- a/vtkm/cont/PartitionedDataSet.h +++ b/vtkm/cont/PartitionedDataSet.h @@ -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& 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& 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& partitions); diff --git a/vtkm/filter/scalar_topology/CMakeLists.txt b/vtkm/filter/scalar_topology/CMakeLists.txt index 44cee269d..01a938d26 100644 --- a/vtkm/filter/scalar_topology/CMakeLists.txt +++ b/vtkm/filter/scalar_topology/CMakeLists.txt @@ -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 diff --git a/vtkm/filter/scalar_topology/ContourTreeUniformAugmented.cxx b/vtkm/filter/scalar_topology/ContourTreeUniformAugmented.cxx index 720bfa9d2..6bf2b7268 100644 --- a/vtkm/filter/scalar_topology/ContourTreeUniformAugmented.cxx +++ b/vtkm/filter/scalar_topology/ContourTreeUniformAugmented.cxx @@ -51,6 +51,7 @@ //============================================================================== #include +#include #include #include @@ -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& localBlockIndices, - const vtkm::cont::ArrayHandle& localBlockOrigins, - const vtkm::cont::ArrayHandle& localBlockSizes) + const vtkm::cont::ArrayHandle& localBlockIndices) { if (this->MultiBlockTreeHelper) { @@ -96,7 +94,7 @@ void ContourTreeAugmented::SetSpatialDecomposition( } this->MultiBlockTreeHelper = std::make_unique( - 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(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(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(bi)); auto currField = currBlock.GetField(this->GetActiveFieldName(), this->GetActiveFieldAssociation()); + + vtkm::Id3 pointDimensions, globalPointDimensions, globalPointIndexStart; + currBlock.GetCellSet().CastAndCallForTypes( + vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(), + pointDimensions, + globalPointDimensions, + globalPointIndexStart); + //const vtkm::cont::ArrayHandle &fieldData = currField.GetData().Cast >(); vtkm::cont::ArrayHandle fieldData; vtkm::cont::ArrayCopy(currField.GetData(), fieldData); auto currContourTreeMesh = vtkm::worklet::contourtree_distributed::MultiBlockContourTreeHelper:: - ComputeLocalContourTreeMesh( - this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.LocalBlockOrigins.ReadPortal() - .Get(static_cast(bi)), - this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.LocalBlockSizes.ReadPortal().Get( - static_cast(bi)), - this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.GlobalSize, - fieldData, - MultiBlockTreeHelper->LocalContourTrees[bi], - MultiBlockTreeHelper->LocalSortOrders[bi], - compRegularStruct); + ComputeLocalContourTreeMesh(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(); @@ -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(bi)); - localDataBlocks[bi]->BlockSize = localBlocksSizesPortal.Get(static_cast(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; - 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(spatialDecomp.NumberOfDimensions()); - ++d) + RegularDecomposer::DivisionsVector diyDivisions; + std::vector vtkmdiyLocalBlockGids; + vtkmdiy::DiscreteBounds diyBounds(0); + if (this->MultiBlockTreeHelper->BlocksPerDimension[0] == -1) { - diyDivisions[d] = static_cast(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 vtkmdiyLocalBlockGids(static_cast(input.GetNumberOfPartitions())); - for (vtkm::Id bi = 0; bi < input.GetNumberOfPartitions(); bi++) + else { - RegularDecomposer::DivisionsVector diyCoords(static_cast(numDims)); - auto currentCoords = localBlockIndicesPortal.Get(bi); - for (vtkm::IdComponent d = 0; d < numDims; ++d) - { - diyCoords[d] = static_cast(currentCoords[d]); - } - vtkmdiyLocalBlockGids[static_cast(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{}); // Add my local blocks to the vtkmdiy master. for (std::size_t bi = 0; bi < static_cast(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(numDims), - spatialDecomp.GetVTKmDIYBounds(), - static_cast(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(size), static_cast(spatialDecomp.GetGlobalNumberOfBlocks())); + vtkmdiy::DynamicAssigner assigner(comm, static_cast(size), globalNumberOfBlocks); for (vtkm::Id bi = 0; bi < input.GetNumberOfPartitions(); bi++) { assigner.set_rank(static_cast(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::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 diff --git a/vtkm/filter/scalar_topology/ContourTreeUniformAugmented.h b/vtkm/filter/scalar_topology/ContourTreeUniformAugmented.h index 0b5dc342c..466b75c90 100644 --- a/vtkm/filter/scalar_topology/ContourTreeUniformAugmented.h +++ b/vtkm/filter/scalar_topology/ContourTreeUniformAugmented.h @@ -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& 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& localBlockIndices, - const vtkm::cont::ArrayHandle& localBlockOrigins, - const vtkm::cont::ArrayHandle& localBlockSizes); + const vtkm::cont::ArrayHandle&, + const vtkm::cont::ArrayHandle&) + { + SetBlockIndices(blocksPerDim, localBlockIndices); + } //@{ /// Get the contour tree computed by the filter diff --git a/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.cxx b/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.cxx index 3ccbdbc2d..d237dd8b0 100644 --- a/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.cxx +++ b/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.cxx @@ -51,6 +51,7 @@ //============================================================================== // vtkm includes +#include #include #include @@ -61,7 +62,7 @@ #include // distributed contour tree includes -#include +#include #include #include #include @@ -165,10 +166,10 @@ void SaveHierarchicalTreeDot( //----------------------------------------------------------------------------- ContourTreeUniformDistributed::ContourTreeUniformDistributed( vtkm::Id3 blocksPerDim, - vtkm::Id3 globalSize, + vtkm::Id3, // globalSize, -> Now in CellSetStructured const vtkm::cont::ArrayHandle& localBlockIndices, - const vtkm::cont::ArrayHandle& localBlockOrigins, - const vtkm::cont::ArrayHandle& localBlockSizes, + const vtkm::cont::ArrayHandle&, // localBlockOrigins, -> Use from CellSetStructured + const vtkm::cont::ArrayHandle&, // 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(localBlockSizes.GetNumberOfValues())) - , LocalContourTrees(static_cast(localBlockSizes.GetNumberOfValues())) - , LocalBoundaryTrees(static_cast(localBlockSizes.GetNumberOfValues())) - , LocalInteriorForests(static_cast(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 void ContourTreeUniformDistributed::ComputeLocalTreeImpl( const vtkm::Id blockIndex, - const vtkm::cont::DataSet&, // input, + const vtkm::cont::DataSet& ds, // input, const vtkm::cont::ArrayHandle& 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::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( @@ -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(input.GetGlobalNumberOfPartitions())); + this->LocalContourTrees.resize(static_cast(input.GetGlobalNumberOfPartitions())); + this->LocalBoundaryTrees.resize(static_cast(input.GetGlobalNumberOfPartitions())); + this->LocalInteriorForests.resize(static_cast(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::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 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; - 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(numDims); ++d) + RegularDecomposer::DivisionsVector diyDivisions; + std::vector vtkmdiyLocalBlockGids; + vtkmdiy::DiscreteBounds diyBounds(0); + if (this->BlocksPerDimension[0] == -1) { - diyDivisions[d] = static_cast(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 vtkmdiyLocalBlockGids(static_cast(input.GetNumberOfPartitions())); - for (vtkm::Id bi = 0; bi < input.GetNumberOfPartitions(); bi++) + else { - RegularDecomposer::DivisionsVector diyCoords(static_cast(numDims)); - auto currentCoords = localBlockIndicesPortal.Get(bi); - for (vtkm::IdComponent d = 0; d < numDims; ++d) - { - diyCoords[d] = static_cast(currentCoords[d]); - } - vtkmdiyLocalBlockGids[static_cast(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{}); // 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(bi)); + vtkm::Id3 pointDimensions, globalPointIndexStart; + currBlock.GetCellSet().CastAndCallForTypes( + 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(bi)), - this->MultiBlockSpatialDecomposition.LocalBlockSizes.ReadPortal().Get( - static_cast(bi)), - this->MultiBlockSpatialDecomposition.GlobalSize)); + globalPointIndexStart, pointDimensions, globalPointDimensions)); vtkm::cont::Algorithm::Copy(transformedIndex, localGlobalMeshIndex); // ... get data values - auto currBlock = input.GetPartition(static_cast(bi)); auto currField = currBlock.GetField(this->GetActiveFieldName(), this->GetActiveFieldAssociation()); vtkm::cont::ArrayHandle fieldData; @@ -942,16 +987,15 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( RegularDecomposer::BoolVector wrap(3, false); RegularDecomposer::CoordinateVector ghosts(3, 1); RegularDecomposer decomposer(static_cast(numDims), - spatialDecomp.GetVTKmDIYBounds(), - static_cast(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(size), static_cast(spatialDecomp.GetGlobalNumberOfBlocks())); + vtkmdiy::DynamicAssigner assigner(comm, static_cast(size), globalNumberOfBlocks); for (vtkm::Id bi = 0; bi < input.GetNumberOfPartitions(); bi++) { assigner.set_rank(static_cast(rank), vtkmdiyLocalBlockGids[static_cast(bi)]); @@ -984,7 +1028,7 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( timer.Start(); // 1.3 Perform fan-in reduction const vtkm::worklet::contourtree_distributed::ComputeDistributedContourTreeFunctor - 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(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 fieldData; vtkm::cont::ArrayCopy(currField.GetData(), fieldData); + + vtkm::Id3 pointDimensions, globalPointIndexStart; + // globalPointDimensions already defined in parent scope + currBlock.GetCellSet().CastAndCallForTypes( + 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 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 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()); diff --git a/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.h b/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.h index d173bb60c..5f915d250 100644 --- a/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.h +++ b/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.h @@ -56,7 +56,6 @@ #include #include -#include #include #include #include @@ -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& localBlockIndices, - const vtkm::cont::ArrayHandle& localBlockOrigins, - const vtkm::cont::ArrayHandle& 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& localBlockIndices, + const vtkm::cont::ArrayHandle& localBlockOrigins, + const vtkm::cont::ArrayHandle& 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& 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 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& localBlockIndices, - const vtkm::cont::ArrayHandle& localBlockOrigins, - const vtkm::cont::ArrayHandle& localBlockSizes, + const vtkm::cont::ArrayHandle&, + const vtkm::cont::ArrayHandle&, 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"); } }; diff --git a/vtkm/filter/scalar_topology/DistributedBranchDecompositionFilter.cxx b/vtkm/filter/scalar_topology/DistributedBranchDecompositionFilter.cxx index 22fc3b6a9..f2adc22cc 100644 --- a/vtkm/filter/scalar_topology/DistributedBranchDecompositionFilter.cxx +++ b/vtkm/filter/scalar_topology/DistributedBranchDecompositionFilter.cxx @@ -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& localBlockIndices, - const vtkm::cont::ArrayHandle& localBlockOrigins, - const vtkm::cont::ArrayHandle& localBlockSizes) - : MultiBlockSpatialDecomposition(blocksPerDim, - globalSize, - localBlockIndices, - localBlockOrigins, - localBlockSizes) + vtkm::Id3, + vtkm::Id3, + const vtkm::cont::ArrayHandle&, + const vtkm::cont::ArrayHandle&, + const vtkm::cont::ArrayHandle&) { } @@ -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; - int globalNumberOfBlocks = - static_cast(this->MultiBlockSpatialDecomposition.GetGlobalNumberOfBlocks()); - int numDims = static_cast(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::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(), + firstPointDimensions, + firstGlobalPointDimensions, + firstGlobalPointIndexStart); + int numDims = firstGlobalPointDimensions[2] > 1 ? 3 : 2; + auto vtkmBlocksPerDimensionRP = input.GetPartition(0) + .GetField("vtkmBlocksPerDimension") + .GetData() + .AsArrayHandle>() + .ReadPortal(); // ... compute division vector for global domain + using RegularDecomposer = vtkmdiy::RegularDecomposer; RegularDecomposer::DivisionsVector diyDivisions(numDims); + vtkmdiy::DiscreteBounds diyBounds(numDims); + int globalNumberOfBlocks = 1; + for (vtkm::IdComponent d = 0; d < static_cast(numDims); ++d) { - diyDivisions[d] = static_cast(this->MultiBlockSpatialDecomposition.BlocksPerDimension[d]); - } - - // ... compute coordinates of local blocks - std::vector vtkmdiyLocalBlockGids(static_cast(input.GetNumberOfPartitions())); - - auto localBlockIndicesPortal = - this->MultiBlockSpatialDecomposition.LocalBlockIndices.ReadPortal(); - for (vtkm::Id bi = 0; bi < input.GetNumberOfPartitions(); bi++) - { - RegularDecomposer::DivisionsVector diyCoords(static_cast(numDims)); - auto currentCoords = localBlockIndicesPortal.Get(bi); - for (vtkm::IdComponent d = 0; d < numDims; ++d) - { - diyCoords[d] = static_cast(currentCoords[d]); - } - vtkmdiyLocalBlockGids[static_cast(bi)] = - RegularDecomposer::coords_to_gid(diyCoords, diyDivisions); + diyDivisions[d] = static_cast(vtkmBlocksPerDimensionRP.Get(d)); + globalNumberOfBlocks *= diyDivisions[d]; + diyBounds.min[d] = 0; + diyBounds.max[d] = static_cast(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( + vtkm::cont::ArrayGetValue(0, + ds.GetField("vtkmGlobalBlockId") + .GetData() + .AsArrayHandle>())); + 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(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(bi)]); } timingsStream << " " << std::setw(60) << std::left @@ -196,7 +197,7 @@ VTKM_CONT vtkm::cont::PartitionedDataSet DistributedBranchDecompositionFilter::D ds.GetField("DependentVolume").GetData().AsArrayHandle>(); // 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 diff --git a/vtkm/filter/scalar_topology/DistributedBranchDecompositionFilter.h b/vtkm/filter/scalar_topology/DistributedBranchDecompositionFilter.h index fc6458694..244bbfa71 100644 --- a/vtkm/filter/scalar_topology/DistributedBranchDecompositionFilter.h +++ b/vtkm/filter/scalar_topology/DistributedBranchDecompositionFilter.h @@ -45,8 +45,6 @@ #include #include -#include - 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& localBlockIndices, - const vtkm::cont::ArrayHandle& localBlockOrigins, - const vtkm::cont::ArrayHandle& localBlockSizes); + VTKM_CONT DistributedBranchDecompositionFilter() = default; + VTKM_CONT DistributedBranchDecompositionFilter(vtkm::Id3, + vtkm::Id3, + const vtkm::cont::ArrayHandle&, + const vtkm::cont::ArrayHandle&, + const vtkm::cont::ArrayHandle&); 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 diff --git a/vtkm/filter/scalar_topology/internal/CMakeLists.txt b/vtkm/filter/scalar_topology/internal/CMakeLists.txt index 6585e20c7..896c847f5 100644 --- a/vtkm/filter/scalar_topology/internal/CMakeLists.txt +++ b/vtkm/filter/scalar_topology/internal/CMakeLists.txt @@ -10,8 +10,8 @@ set(headers BranchDecompositionBlock.h + ComputeBlockIndices.h ComputeDistributedBranchDecompositionFunctor.h - SpatialDecomposition.h ) #----------------------------------------------------------------------------- diff --git a/vtkm/filter/scalar_topology/internal/ComputeBlockIndices.cxx b/vtkm/filter/scalar_topology/internal/ComputeBlockIndices.cxx new file mode 100644 index 000000000..ff881c797 --- /dev/null +++ b/vtkm/filter/scalar_topology/internal/ComputeBlockIndices.cxx @@ -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 +#include +#include +#include +#include +#include + +namespace // anonymous namespace for local helper classes +{ + +struct OriginsBlock +{ + std::vector Origins; + OriginsBlock(const std::vector& origins) + : Origins(origins) + { + } + static void Destroy(void* b) { delete static_cast(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 incoming; + rp.incoming(incoming); + for (const int ingid : incoming) + { + if (ingid != selfid) + { + std::vector incoming_origins; + rp.dequeue(ingid, incoming_origins); + + std::vector 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& diyLocalBlockGids) +{ + auto firstDS = input.GetPartition(0); + vtkm::Id3 dummy1, firstGlobalPointDimensions, dummy2; + firstDS.GetCellSet().CastAndCallForTypes( + vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(), + dummy1, + firstGlobalPointDimensions, + dummy2); + int numDims = firstGlobalPointDimensions[2] > 1 ? 3 : 2; + + diyDivisions.clear(); + vtkmdiy::DiscreteBounds diyBounds(numDims); + std::vector 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(firstGlobalPointDimensions[d]); + + // Get the list of origins along current coordinate axis + std::vector local_origins; + + for (vtkm::Id ds_no = 0; ds_no < input.GetNumberOfPartitions(); ++ds_no) + { + vtkm::Id3 globalPointIndexStart; + input.GetPartition(ds_no) + .GetCellSet() + .CastAndCallForTypes( + vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(), + dummy1, + dummy2, + globalPointIndexStart); + + local_origins.push_back(static_cast(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 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(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(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::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& blockIndices, + DiscreteBoundsDivisionVector& diyDivisions, + std::vector& diyLocalBlockGids) +{ + auto firstDS = input.GetPartition(0); + vtkm::Id3 dummy1, firstGlobalPointDimensions, dummy2; + firstDS.GetCellSet().CastAndCallForTypes( + 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(firstGlobalPointDimensions[d]); + diyDivisions.push_back(static_cast(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(currBlockIndices[d]); + } + + diyLocalBlockGids.push_back(vtkmdiy::RegularDecomposer::coords_to_gid( + diyBlockCoords, diyDivisions)); + } + + return diyBounds; +} + +} // namespace internal +} // namespace scalar_topology +} // namespace filter +} // namespace vtkm diff --git a/vtkm/filter/scalar_topology/internal/ComputeBlockIndices.h b/vtkm/filter/scalar_topology/internal/ComputeBlockIndices.h new file mode 100644 index 000000000..1b492e9b5 --- /dev/null +++ b/vtkm/filter/scalar_topology/internal/ComputeBlockIndices.h @@ -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 + +#include +#include +#include + +// DIY includes +// clang-format off +VTKM_THIRDPARTY_PRE_INCLUDE +#include +#include +VTKM_THIRDPARTY_POST_INCLUDE +// clang-format on + +namespace vtkm +{ +namespace filter +{ +namespace scalar_topology +{ +namespace internal +{ + +using DiscreteBoundsDivisionVector = + vtkmdiy::RegularDecomposer::DivisionsVector; + +VTKM_CONT vtkmdiy::DiscreteBounds ComputeBlockIndices(const vtkm::cont::PartitionedDataSet& input, + DiscreteBoundsDivisionVector& diyDivisions, + std::vector& diyLocalBlockGids); + +VTKM_CONT vtkmdiy::DiscreteBounds ComputeBlockIndices( + const vtkm::cont::PartitionedDataSet& input, + vtkm::Id3 blocksPerDim, + const vtkm::cont::ArrayHandle& blockIndices, + DiscreteBoundsDivisionVector& diyDivisions, + std::vector& diyLocalBlockGids); + +} // namespace internal +} // namespace scalar_topology +} // namespace filter +} // namespace vtkm + +#endif diff --git a/vtkm/filter/scalar_topology/internal/SpatialDecomposition.h b/vtkm/filter/scalar_topology/internal/SpatialDecomposition.h deleted file mode 100644 index 5ab7ab689..000000000 --- a/vtkm/filter/scalar_topology/internal/SpatialDecomposition.h +++ /dev/null @@ -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 -#include -#include -#include -#include - -#include - -// clang-format off -VTKM_THIRDPARTY_PRE_INCLUDE -#include -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& localBlockIndices, - const vtkm::cont::ArrayHandle& localBlockOrigins, - const vtkm::cont::ArrayHandle& 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(this->GlobalSize[0]); - domain.max[1] = static_cast(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(this->GlobalSize[0]); - domain.max[1] = static_cast(this->GlobalSize[1]); - domain.max[2] = static_cast(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{}); -#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 LocalBlockIndices; - // Origin of the local blocks in mesh index space - vtkm::cont::ArrayHandle LocalBlockOrigins; - // Size of each local block in x, y,z - vtkm::cont::ArrayHandle LocalBlockSizes; -}; - - -} // internal -} // namespace scalar_topology -} // namespace filter -} // namespace vtkm - -#endif diff --git a/vtkm/filter/scalar_topology/testing/TestingContourTreeUniformDistributedFilter.h b/vtkm/filter/scalar_topology/testing/TestingContourTreeUniformDistributedFilter.h index 2747866c2..431ef376a 100644 --- a/vtkm/filter/scalar_topology/testing/TestingContourTreeUniformDistributedFilter.h +++ b/vtkm/filter/scalar_topology/testing/TestingContourTreeUniformDistributedFilter.h @@ -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 localBlockIndices; - vtkm::cont::ArrayHandle localBlockOrigins; - vtkm::cont::ArrayHandle 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) { diff --git a/vtkm/filter/scalar_topology/testing/UnitTestContourTreeUniformDistributedFilter.cxx b/vtkm/filter/scalar_topology/testing/UnitTestContourTreeUniformDistributedFilter.cxx index 35e878c32..e37194004 100644 --- a/vtkm/filter/scalar_topology/testing/UnitTestContourTreeUniformDistributedFilter.cxx +++ b/vtkm/filter/scalar_topology/testing/UnitTestContourTreeUniformDistributedFilter.cxx @@ -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); } }; } diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h b/vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h index f9752a05f..1927b3831 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h @@ -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 + 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 diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/MultiBlockContourTreeHelper.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/MultiBlockContourTreeHelper.h index ba9095521..5aad7f3ae 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/MultiBlockContourTreeHelper.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/MultiBlockContourTreeHelper.h @@ -53,15 +53,12 @@ #ifndef vtk_m_worklet_contourtree_distributed_multiblockcontourtreehelper_h #define vtk_m_worklet_contourtree_distributed_multiblockcontourtreehelper_h -#include - #include #include #include #include #include -//#include #include #include #include @@ -79,19 +76,21 @@ class MultiBlockContourTreeHelper public: VTKM_CONT MultiBlockContourTreeHelper(vtkm::Id3 blocksPerDim, - vtkm::Id3 globalSize, - const vtkm::cont::ArrayHandle& localBlockIndices, - const vtkm::cont::ArrayHandle& localBlockOrigins, - const vtkm::cont::ArrayHandle& localBlockSizes) - : MultiBlockSpatialDecomposition(blocksPerDim, - globalSize, - localBlockIndices, - localBlockOrigins, - localBlockSizes) + const vtkm::cont::ArrayHandle& localBlockIndices) + : BlocksPerDimension(blocksPerDim) + , LocalBlockIndices(localBlockIndices) + , LocalContourTrees(static_cast(localBlockIndices.GetNumberOfValues())) + , LocalSortOrders(static_cast(localBlockIndices.GetNumberOfValues())) + { + } + + VTKM_CONT + MultiBlockContourTreeHelper(const vtkm::cont::PartitionedDataSet& input) + : BlocksPerDimension(-1, -1, -1) + , LocalBlockIndices() + , LocalContourTrees(static_cast(input.GetNumberOfPartitions())) + , LocalSortOrders(static_cast(input.GetNumberOfPartitions())) { - vtkm::Id localNumBlocks = this->GetLocalNumberOfBlocks(); - LocalContourTrees.resize(static_cast(localNumBlocks)); - LocalSortOrders.resize(static_cast(localNumBlocks)); } VTKM_CONT @@ -118,25 +117,12 @@ public: inline vtkm::Id GetLocalNumberOfBlocks() const { - return this->MultiBlockSpatialDecomposition.GetLocalNumberOfBlocks(); + return static_cast(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{}); -#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 LocalBlockIndices; std::vector LocalContourTrees; std::vector LocalSortOrders; }; // end MultiBlockContourTreeHelper diff --git a/vtkm/internal/ConnectivityStructuredInternals.h b/vtkm/internal/ConnectivityStructuredInternals.h index 5240850ae..7e2e0717a 100644 --- a/vtkm/internal/ConnectivityStructuredInternals.h +++ b/vtkm/internal/ConnectivityStructuredInternals.h @@ -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; diff --git a/vtkm/worklet/testing/UnitTestContourTreeUniformDistributed.cxx b/vtkm/worklet/testing/UnitTestContourTreeUniformDistributed.cxx index 41e773088..e5c68258f 100644 --- a/vtkm/worklet/testing/UnitTestContourTreeUniformDistributed.cxx +++ b/vtkm/worklet/testing/UnitTestContourTreeUniformDistributed.cxx @@ -55,7 +55,6 @@ #include #include -#include #include #include #include @@ -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({ 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 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(globalSize[2] > 1 ? 3 : 2); - RegularDecomposer decomposer(numDims, - spatialDecomp.GetVTKmDIYBounds(), - static_cast(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(globalSize[0]); + diyBounds.max[1] = static_cast(globalSize[1]); + + RegularDecomposer decomposer( + numDims, diyBounds, numBlocks, shareFace, wrap, ghosts, diyDivisions); // ... coordinates of local blocks - auto localBlockIndicesPortal = spatialDecomp.LocalBlockIndices.ReadPortal(); std::vector vtkmdiyLocalBlockGids(numBlocks); for (vtkm::Id bi = 0; bi < numBlocks; bi++) { RegularDecomposer::DivisionsVector diyCoords(static_cast(numDims)); - auto currentCoords = localBlockIndicesPortal.Get(bi); + auto currentCoords = blockIndices[bi]; for (vtkm::IdComponent d = 0; d < numDims; ++d) { diyCoords[d] = static_cast(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(spatialDecomp.GetGlobalNumberOfBlocks())); + vtkmdiy::DynamicAssigner assigner(comm, comm.size(), numBlocks); for (vtkm::Id bi = 0; bi < numBlocks; bi++) { assigner.set_rank(static_cast(rank),