From 99900576c558b65a105a884cd0863ddc03325411 Mon Sep 17 00:00:00 2001 From: "Gunther H. Weber" Date: Wed, 30 Nov 2022 19:57:43 -0800 Subject: [PATCH 1/3] Add ComputeGlobalPointSize --- .../ContourTreeApp.cxx | 75 +++++++++++++++++++ .../ContourTreeAppDataIO.h | 15 ++-- 2 files changed, 85 insertions(+), 5 deletions(-) diff --git a/examples/contour_tree_distributed/ContourTreeApp.cxx b/examples/contour_tree_distributed/ContourTreeApp.cxx index 293ce2dbb..9c454855e 100644 --- a/examples/contour_tree_distributed/ContourTreeApp.cxx +++ b/examples/contour_tree_distributed/ContourTreeApp.cxx @@ -154,6 +154,80 @@ private: std::vector mCLOptions; }; +void ComputeGlobalPointSize(vtkm::cont::PartitionedDataSet& pds) +{ + // Compute GlobalPointDimensions as maximim of GlobalPointIndexStart + PointDimensions + // for each dimension across all blocks + + // Compute GlobalPointDimensions for all data sets on this MPI rank + std::vector globalPointDimensionsThisRank; + using ds_const_iterator = vtkm::cont::PartitionedDataSet::const_iterator; + for (ds_const_iterator ds_it = pds.cbegin(); ds_it != pds.cend(); ++ds_it) + { + ds_it->GetCellSet().CastAndCallForTypes( + [&globalPointDimensionsThisRank](const auto& css) { + globalPointDimensionsThisRank.resize(css.Dimension, -1); + for (vtkm::IdComponent d = 0; d < css.Dimension; ++d) + { + globalPointDimensionsThisRank[d] = + std::max(globalPointDimensionsThisRank[d], + css.GetGlobalPointIndexStart()[d] + css.GetPointDimensions()[d]); + } + }); + } + + // Perform global reduction to find GlobalPointDimensions across all ranks + std::vector globalPointDimensions; + auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator(); + vtkmdiy::mpi::all_reduce( + comm, globalPointDimensionsThisRank, globalPointDimensions, vtkmdiy::mpi::maximum{}); + + // Set this information in all cell sets + using ds_iterator = vtkm::cont::PartitionedDataSet::iterator; + for (ds_iterator ds_it = pds.begin(); ds_it != pds.end(); ++ds_it) + { +#if 0 + // This does not work, i.e., it does not really change the cell set for the DataSet + ds_it->GetCellSet().CastAndCallForTypes( + [&globalPointDimensions](auto& css) + { + typename std::remove_reference_t::SchedulingRangeType gpd; + for (vtkm::IdComponent d = 0; d < css.Dimension; ++d) + { + gpd[d] = globalPointDimensions[d]; + } + css.SetGlobalPointDimensions(gpd); + }); +#else + if (globalPointDimensions.size() == 2) + { + vtkm::cont::CellSetStructured<2> cs; + ds_it->GetCellSet().AsCellSet>(cs); + vtkm::Id2 gpd{ globalPointDimensions[0], globalPointDimensions[1] }; + cs.SetGlobalPointDimensions(gpd); + ds_it->SetCellSet(cs); + } + else if (globalPointDimensions.size() == 3) + { + + vtkm::cont::CellSetStructured<3> cs; + ds_it->GetCellSet().AsCellSet>(cs); + vtkm::Id3 gpd{ globalPointDimensions[0], globalPointDimensions[1], globalPointDimensions[2] }; + cs.SetGlobalPointDimensions(gpd); + ds_it->SetCellSet(cs); + } + else + { + throw vtkm::cont::ErrorBadValue("StructuredCellSet dimension must be 2 or 3"); + } +#endif + } + + for (ds_const_iterator ds_it = pds.cbegin(); ds_it != pds.cend(); ++ds_it) + { + ds_it->PrintSummary(std::cout); + } +} // Compute and render an isosurface for a uniform grid example int main(int argc, char* argv[]) @@ -627,6 +701,7 @@ int main(int argc, char* argv[]) filter.SetActiveField("values"); // Execute the contour tree analysis + ComputeGlobalPointSize(useDataSet); auto result = filter.Execute(useDataSet); currTime = totalTime.GetElapsedTime(); diff --git a/examples/contour_tree_distributed/ContourTreeAppDataIO.h b/examples/contour_tree_distributed/ContourTreeAppDataIO.h index acbf32b1c..bef3ebec6 100644 --- a/examples/contour_tree_distributed/ContourTreeAppDataIO.h +++ b/examples/contour_tree_distributed/ContourTreeAppDataIO.h @@ -346,7 +346,8 @@ bool read3DHDF5File(const int& mpi_rank, ds = dsb.Create(v_dims, v_origin, v_spacing); vtkm::cont::CellSetStructured<3> cs; cs.SetPointDimensions(v_dims); - cs.SetGlobalPointDimensions(globalSize); + // NOTE: Comment out for test purposes + //cs.SetGlobalPointDimensions(globalSize); cs.SetGlobalPointIndexStart(vtkm::Id3{ v_origin[0], v_origin[1], v_origin[2] }); ds.SetCellSet(cs); @@ -587,7 +588,8 @@ bool readPreSplitFiles(const int& rank, 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] }); + // NOTE: Comment out for test purposes + //cs.SetGlobalPointDimensions(vtkm::Id2{ globalSize[0], globalSize[1] }); cs.SetGlobalPointIndexStart(vtkm::Id2{ offset[0], offset[1] }); ds.SetCellSet(cs); } @@ -604,7 +606,8 @@ bool readPreSplitFiles(const int& rank, ds = dsb.Create(v_dims, v_origin, v_spacing); vtkm::cont::CellSetStructured<3> cs; cs.SetPointDimensions(v_dims); - cs.SetGlobalPointDimensions(globalSize); + // NOTE: Comment out for test purposes + //cs.SetGlobalPointDimensions(globalSize); cs.SetGlobalPointIndexStart(vtkm::Id3{ offset[0], offset[1], offset[2] }); ds.SetCellSet(cs); } @@ -815,7 +818,8 @@ bool readSingleBlockFile(const int& rank, ds = dsb.Create(vdims, origin, spacing); vtkm::cont::CellSetStructured<2> cs; cs.SetPointDimensions(vdims); - cs.SetGlobalPointDimensions(vtkm::Id2{ globalSize[0], globalSize[1] }); + // NOTE: Comment out for test purposes + //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)); @@ -832,7 +836,8 @@ bool readSingleBlockFile(const int& rank, ds = dsb.Create(vdims, origin, spacing); vtkm::cont::CellSetStructured<3> cs; cs.SetPointDimensions(vdims); - cs.SetGlobalPointDimensions(globalSize); + // NOTE: Comment out for test purposes + //cs.SetGlobalPointDimensions(globalSize); cs.SetGlobalPointIndexStart(vtkm::Id3(0, 0, blockStart / blockSliceSize)); ds.SetCellSet(cs); localBlockIndicesPortal.Set(localBlockIndex, vtkm::Id3(0, 0, blockIndex)); From 8b87ae192fff74b154b8e680d0ca512f4049c168 Mon Sep 17 00:00:00 2001 From: "Gunther H. Weber" Date: Thu, 1 Dec 2022 12:04:59 -0800 Subject: [PATCH 2/3] Updates based on review --- .../ContourTreeApp.cxx | 42 ++++--------------- 1 file changed, 9 insertions(+), 33 deletions(-) diff --git a/examples/contour_tree_distributed/ContourTreeApp.cxx b/examples/contour_tree_distributed/ContourTreeApp.cxx index 9c454855e..b45bce9e7 100644 --- a/examples/contour_tree_distributed/ContourTreeApp.cxx +++ b/examples/contour_tree_distributed/ContourTreeApp.cxx @@ -154,9 +154,9 @@ private: std::vector mCLOptions; }; -void ComputeGlobalPointSize(vtkm::cont::PartitionedDataSet& pds) +void ComputeGlobalPointDimensions(vtkm::cont::PartitionedDataSet& pds) { - // Compute GlobalPointDimensions as maximim of GlobalPointIndexStart + PointDimensions + // Compute GlobalPointDimensions as maximum of GlobalPointIndexStart + PointDimensions // for each dimension across all blocks // Compute GlobalPointDimensions for all data sets on this MPI rank @@ -186,47 +186,23 @@ void ComputeGlobalPointSize(vtkm::cont::PartitionedDataSet& pds) using ds_iterator = vtkm::cont::PartitionedDataSet::iterator; for (ds_iterator ds_it = pds.begin(); ds_it != pds.end(); ++ds_it) { -#if 0 // This does not work, i.e., it does not really change the cell set for the DataSet ds_it->GetCellSet().CastAndCallForTypes( - [&globalPointDimensions](auto& css) - { + [&globalPointDimensions, &ds_it](auto& css) { typename std::remove_reference_t::SchedulingRangeType gpd; for (vtkm::IdComponent d = 0; d < css.Dimension; ++d) { gpd[d] = globalPointDimensions[d]; } css.SetGlobalPointDimensions(gpd); + // Why is the following necessary? Shouldn't it be sufficient to update the + // CellSet through the reference? + ds_it->SetCellSet(css); }); -#else - if (globalPointDimensions.size() == 2) - { - vtkm::cont::CellSetStructured<2> cs; - ds_it->GetCellSet().AsCellSet>(cs); - vtkm::Id2 gpd{ globalPointDimensions[0], globalPointDimensions[1] }; - cs.SetGlobalPointDimensions(gpd); - ds_it->SetCellSet(cs); - } - else if (globalPointDimensions.size() == 3) - { - - vtkm::cont::CellSetStructured<3> cs; - ds_it->GetCellSet().AsCellSet>(cs); - vtkm::Id3 gpd{ globalPointDimensions[0], globalPointDimensions[1], globalPointDimensions[2] }; - cs.SetGlobalPointDimensions(gpd); - ds_it->SetCellSet(cs); - } - else - { - throw vtkm::cont::ErrorBadValue("StructuredCellSet dimension must be 2 or 3"); - } -#endif } - for (ds_const_iterator ds_it = pds.cbegin(); ds_it != pds.cend(); ++ds_it) - { - ds_it->PrintSummary(std::cout); - } + // Debug + // pds.PrintSummary(std::cout); } // Compute and render an isosurface for a uniform grid example @@ -701,7 +677,7 @@ int main(int argc, char* argv[]) filter.SetActiveField("values"); // Execute the contour tree analysis - ComputeGlobalPointSize(useDataSet); + ComputeGlobalPointDimensions(useDataSet); auto result = filter.Execute(useDataSet); currTime = totalTime.GetElapsedTime(); From fae6ef8c93fc75314c763443012e230ac5703d85 Mon Sep 17 00:00:00 2001 From: "Gunther H. Weber" Date: Fri, 27 Jan 2023 13:46:54 -0800 Subject: [PATCH 3/3] Add code to shift logical extents --- .../ContourTreeApp.cxx | 57 ++++++++++++++++++- 1 file changed, 56 insertions(+), 1 deletion(-) diff --git a/examples/contour_tree_distributed/ContourTreeApp.cxx b/examples/contour_tree_distributed/ContourTreeApp.cxx index b45bce9e7..f28156ff9 100644 --- a/examples/contour_tree_distributed/ContourTreeApp.cxx +++ b/examples/contour_tree_distributed/ContourTreeApp.cxx @@ -202,7 +202,61 @@ void ComputeGlobalPointDimensions(vtkm::cont::PartitionedDataSet& pds) } // Debug - // pds.PrintSummary(std::cout); + pds.PrintSummary(std::cout); +} +void ShiftLogicalOriginToZero(vtkm::cont::PartitionedDataSet& pds) +{ + // Shift the logical origin (minimum of LocalPointIndexStart) to zero + // along each dimension + + // Compute minimum global point index start for all data sets on this MPI rank + std::vector minimumGlobalPointIndexStartThisRank; + using ds_const_iterator = vtkm::cont::PartitionedDataSet::const_iterator; + for (ds_const_iterator ds_it = pds.cbegin(); ds_it != pds.cend(); ++ds_it) + { + ds_it->GetCellSet().CastAndCallForTypes( + [&minimumGlobalPointIndexStartThisRank](const auto& css) { + minimumGlobalPointIndexStartThisRank.resize(css.Dimension, + std::numeric_limits::max()); + for (vtkm::IdComponent d = 0; d < css.Dimension; ++d) + { + minimumGlobalPointIndexStartThisRank[d] = + std::min(minimumGlobalPointIndexStartThisRank[d], css.GetGlobalPointIndexStart()[d]); + } + }); + } + + // Perform global reduction to find GlobalPointDimensions across all ranks + std::vector minimumGlobalPointIndexStart; + auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator(); + vtkmdiy::mpi::all_reduce(comm, + minimumGlobalPointIndexStartThisRank, + minimumGlobalPointIndexStart, + vtkmdiy::mpi::minimum{}); + + // Shift all cell sets so that minimum global point index start + // along each dimension is zero + using ds_iterator = vtkm::cont::PartitionedDataSet::iterator; + for (ds_iterator ds_it = pds.begin(); ds_it != pds.end(); ++ds_it) + { + // This does not work, i.e., it does not really change the cell set for the DataSet + ds_it->GetCellSet().CastAndCallForTypes( + [&minimumGlobalPointIndexStart, &ds_it](auto& css) { + auto pointIndexStart = css.GetGlobalPointIndexStart(); + typename std::remove_reference_t::SchedulingRangeType shiftedPointIndexStart; + for (vtkm::IdComponent d = 0; d < css.Dimension; ++d) + { + shiftedPointIndexStart[d] = pointIndexStart[d] - minimumGlobalPointIndexStart[d]; + } + css.SetGlobalPointIndexStart(shiftedPointIndexStart); + // Why is the following necessary? Shouldn't it be sufficient to update the + // CellSet through the reference? + ds_it->SetCellSet(css); + }); + } + + // Debug + //pds.PrintSummary(std::cout); } // Compute and render an isosurface for a uniform grid example @@ -677,6 +731,7 @@ int main(int argc, char* argv[]) filter.SetActiveField("values"); // Execute the contour tree analysis + ShiftLogicalOriginToZero(useDataSet); ComputeGlobalPointDimensions(useDataSet); auto result = filter.Execute(useDataSet);