diff --git a/examples/contour_tree_augmented/CMakeLists.txt b/examples/contour_tree_augmented/CMakeLists.txt index 8e94647f0..35d6be675 100644 --- a/examples/contour_tree_augmented/CMakeLists.txt +++ b/examples/contour_tree_augmented/CMakeLists.txt @@ -59,7 +59,7 @@ find_package(VTKm REQUIRED QUIET) # Serial #################################### add_executable(ContourTree_Augmented ContourTreeApp.cxx) -target_link_libraries(ContourTree_Augmented vtkm_filter) +target_link_libraries(ContourTree_Augmented vtkm_filter vtkm_io) vtkm_add_target_information(ContourTree_Augmented DROP_UNUSED_SYMBOLS MODIFY_CUDA_FLAGS DEVICE_SOURCES ContourTreeApp.cxx) @@ -78,7 +78,7 @@ endif() #################################### if (VTKm_ENABLE_MPI) add_executable(ContourTree_Augmented_MPI ContourTreeApp.cxx) - target_link_libraries(ContourTree_Augmented_MPI vtkm_filter) + target_link_libraries(ContourTree_Augmented_MPI vtkm_filter vtkm_io) vtkm_add_target_information(ContourTree_Augmented_MPI MODIFY_CUDA_FLAGS DEVICE_SOURCES ContourTreeApp.cxx) diff --git a/examples/contour_tree_augmented/ContourTreeApp.cxx b/examples/contour_tree_augmented/ContourTreeApp.cxx index 0859a4595..7efffd5a6 100644 --- a/examples/contour_tree_augmented/ContourTreeApp.cxx +++ b/examples/contour_tree_augmented/ContourTreeApp.cxx @@ -64,9 +64,13 @@ #include #include #include +#include +#include #include #include #include +#include + #include #include #include @@ -228,18 +232,21 @@ int main(int argc, char* argv[]) int numThreads = tbb::task_scheduler_init::default_num_threads(); if (parser.hasOption("--numThreads")) { - // Print warning about mismatch between the --numThreads and -d/--device opton - VTKM_LOG_IF_S(vtkm::cont::LogLevel::Warn, - device != vtkm::cont::DeviceAdapterTagTBB{}, - "WARNING: Mismatch between --numThreads and -d/--device option." - "numThreads option requires the use of TBB as device. " - "Ignoring the numThread option."); + bool deviceIsTBB = (device.GetName() == "TBB"); // Set the number of threads to be used for TBB - if (device == vtkm::cont::DeviceAdapterTagTBB{}) + if (deviceIsTBB) { numThreads = std::stoi(parser.getOption("--numThreads")); tbb::task_scheduler_init schedulerInit(numThreads); } + // Print warning about mismatch between the --numThreads and -d/--device option + else + { + VTKM_LOG_S(vtkm::cont::LogLevel::Warn, + "WARNING: Mismatch between --numThreads and -d/--device option." + "numThreads option requires the use of TBB as device. " + "Ignoring the numThread option."); + } } #endif @@ -432,84 +439,99 @@ int main(int argc, char* argv[]) /////////////////////////////////////////////// // Read the input data /////////////////////////////////////////////// - std::ifstream inFile(filename); - if (inFile.bad()) - return 0; - - // Read the dimensions of the mesh, i.e,. number of elementes in x, y, and z - std::vector dims; - std::string line; - getline(inFile, line); - std::istringstream linestream(line); - std::size_t dimVertices; - while (linestream >> dimVertices) - { - dims.push_back(dimVertices); - } - - // Compute the number of vertices, i.e., xdim * ydim * zdim - unsigned short nDims = static_cast(dims.size()); - std::size_t numVertices = static_cast( - std::accumulate(dims.begin(), dims.end(), std::size_t(1), std::multiplies())); - - // Print the mesh metadata - if (rank == 0) - { - VTKM_LOG_S(vtkm::cont::LogLevel::Info, - std::endl - << " ---------------- Input Mesh Properties --------------" - << std::endl - << " Number of dimensions: " - << nDims - << std::endl - << " Number of mesh vertices: " - << numVertices - << std::endl); - } - - // Check for fatal input errors - // Check the the number of dimensiosn is either 2D or 3D - bool invalidNumDimensions = (nDims < 2 || nDims > 3); - // Check if marching cubes is enabled for non 3D data - bool invalidMCOption = (useMarchingCubes && nDims != 3); - // Log any errors if found on rank 0 - VTKM_LOG_IF_S(vtkm::cont::LogLevel::Error, - invalidNumDimensions && (rank == 0), - "The input mesh is " << nDims << "D. " - "The input data must be either 2D or 3D."); - VTKM_LOG_IF_S( - vtkm::cont::LogLevel::Error, - invalidMCOption && (rank == 0), - "The input mesh is " << nDims << "D. " - << "Contour tree using marching cubes is only supported for 3D data."); - // If we found any errors in the setttings than finalize MPI and exit the execution - if (invalidNumDimensions || invalidMCOption) + vtkm::Float64 dataReadTime = 0; + vtkm::Float64 buildDatasetTime = 0; + std::vector::size_type nDims = 0; + vtkm::cont::DataSet inDataSet; + std::vector values; + std::vector dims; + if (filename.compare(filename.length() - 3, 3, "bov") == 0) { + std::cout << "Reading BOV file" << std::endl; + vtkm::io::BOVDataSetReader reader(filename); + inDataSet = reader.ReadDataSet(); + nDims = 3; + currTime = totalTime.GetElapsedTime(); + dataReadTime = currTime - prevTime; + prevTime = currTime; #ifdef WITH_MPI + // Copy the data into the values array so we can construct a multiblock dataset + // TODO All we should need to do to implement BOV support is to copy the values + // in the values vector and copy the dimensions in the dims vector + vtkm::Id nRows, nCols, nSlices; + vtkm::filter::GetRowsColsSlices temp; + temp(inDataSet.GetCellSet(), nRows, nCols, nSlices); + dims[0] = nRows; + dims[1] = nCols; + dims[2] = nSlices; + auto tempField = inDataSet.GetField("values").GetData(); + values.resize(static_cast(tempField.GetNumberOfValues())); + auto tempFieldHandle = tempField.AsVirtual().ReadPortal(); + for (vtkm::Id i = 0; i < tempField.GetNumberOfValues(); i++) + { + values[static_cast(i)] = static_cast(tempFieldHandle.Get(i)); + } + VTKM_LOG_S(vtkm::cont::LogLevel::Error, + "BOV reader not yet support in MPI mode by this example"); MPI_Finalize(); -#endif return EXIT_SUCCESS; +#endif } - - // Read data - std::vector values(numVertices); - for (std::size_t vertex = 0; vertex < numVertices; ++vertex) + else // Read ASCII data input { - inFile >> values[vertex]; - } + std::cout << "Reading ASCII file" << std::endl; + std::ifstream inFile(filename); + if (inFile.bad()) + return 0; - // finish reading the data - inFile.close(); + // Read the dimensions of the mesh, i.e,. number of elementes in x, y, and z + std::string line; + getline(inFile, line); + std::istringstream linestream(line); + vtkm::Id dimVertices; + while (linestream >> dimVertices) + { + dims.push_back(dimVertices); + } - currTime = totalTime.GetElapsedTime(); - vtkm::Float64 dataReadTime = currTime - prevTime; - prevTime = currTime; + // Compute the number of vertices, i.e., xdim * ydim * zdim + nDims = static_cast(dims.size()); + std::size_t numVertices = static_cast( + std::accumulate(dims.begin(), dims.end(), std::size_t(1), std::multiplies())); - vtkm::cont::DataSetBuilderUniform dsb; -#ifndef WITH_MPI // construct regular, single-block VTK-M input dataset - vtkm::cont::DataSet inDataSet; // Single block dataset - { + // Check for fatal input errors + // Check the the number of dimensiosn is either 2D or 3D + bool invalidNumDimensions = (nDims < 2 || nDims > 3); + // Log any errors if found on rank 0 + VTKM_LOG_IF_S(vtkm::cont::LogLevel::Error, + invalidNumDimensions && (rank == 0), + "The input mesh is " << nDims << "D. The input data must be either 2D or 3D."); + // If we found any errors in the setttings than finalize MPI and exit the execution + if (invalidNumDimensions) + { +#ifdef WITH_MPI + MPI_Finalize(); +#endif + return EXIT_SUCCESS; + } + + // Read data + values.resize(numVertices); + for (std::size_t vertex = 0; vertex < numVertices; ++vertex) + { + inFile >> values[vertex]; + } + + // finish reading the data + inFile.close(); + + currTime = totalTime.GetElapsedTime(); + dataReadTime = currTime - prevTime; + prevTime = currTime; + +#ifndef WITH_MPI // We only need the inDataSet if are not using MPI otherwise we'll constructe a multi-block dataset // build the input dataset + vtkm::cont::DataSetBuilderUniform dsb; // 2D data if (nDims == 2) { @@ -528,9 +550,41 @@ int main(int argc, char* argv[]) inDataSet = dsb.Create(vdims); } inDataSet.AddPointField("values", values); +#endif + } // END ASCII Read + + // Print the mesh metadata + if (rank == 0) + { + VTKM_LOG_S(vtkm::cont::LogLevel::Info, + std::endl + << " ---------------- Input Mesh Properties --------------" + << std::endl + << " Number of dimensions: " + << nDims); } -#else // Create a multi-block dataset for multi-block DIY-paralle processing - vtkm::cont::PartitionedDataSet inDataSet; // Partitioned variant of the input dataset + + // Check if marching cubes is enabled for non 3D data + bool invalidMCOption = (useMarchingCubes && nDims != 3); + VTKM_LOG_IF_S( + vtkm::cont::LogLevel::Error, + invalidMCOption && (rank == 0), + "The input mesh is " << nDims << "D. " + << "Contour tree using marching cubes is only supported for 3D data."); + + // If we found any errors in the setttings than finalize MPI and exit the execution + if (invalidMCOption) + { +#ifdef WITH_MPI + MPI_Finalize(); +#endif + return EXIT_SUCCESS; + } + +#ifndef WITH_MPI // construct regular, single-block VTK-M input dataset + vtkm::cont::DataSet useDataSet = inDataSet; // Single block dataset +#else // Create a multi-block dataset for multi-block DIY-paralle processing + vtkm::cont::PartitionedDataSet useDataSet; // Partitioned variant of the input dataset vtkm::Id3 blocksPerDim = nDims == 3 ? vtkm::Id3(1, 1, numBlocks) : vtkm::Id3(1, numBlocks, 1); // Decompose the data into vtkm::Id3 globalSize = nDims == 3 ? vtkm::Id3(static_cast(dims[0]), @@ -584,6 +638,7 @@ int main(int argc, char* argv[]) } vtkm::Id currBlockSize = (vtkm::Id)((blockEnd - blockStart) / blockSliceSize); + vtkm::cont::DataSetBuilderUniform dsb; vtkm::cont::DataSet ds; // 2D data @@ -625,13 +680,13 @@ int main(int argc, char* argv[]) (values.begin() + blockEnd)); ds.AddPointField("values", subValues); - inDataSet.AppendPartition(ds); + useDataSet.AppendPartition(ds); } } -#endif // WITH_MPI construct input dataset +#endif // WITH_MPI construct input dataset currTime = totalTime.GetElapsedTime(); - vtkm::Float64 buildDatasetTime = currTime - prevTime; + buildDatasetTime = currTime - prevTime; prevTime = currTime; // Convert the mesh of values into contour tree, pairs of vertex ids @@ -645,7 +700,7 @@ int main(int argc, char* argv[]) // Execute the contour tree analysis. NOTE: If MPI is used the result will be // a vtkm::cont::PartitionedDataSet instead of a vtkm::cont::DataSet - auto result = filter.Execute(inDataSet); + auto result = filter.Execute(useDataSet); currTime = totalTime.GetElapsedTime(); vtkm::Float64 computeContourTreeTime = currTime - prevTime; @@ -710,7 +765,7 @@ int main(int argc, char* argv[]) bool dataFieldIsSorted = true; #else vtkm::cont::ArrayHandle dataField; - inDataSet.GetField(0).GetData().CopyTo(dataField); + useDataSet.GetField(0).GetData().CopyTo(dataField); bool dataFieldIsSorted = false; #endif @@ -916,6 +971,12 @@ int main(int argc, char* argv[]) << ": " << ct.Hyperarcs.GetNumberOfValues() << std::endl); + // Print hyperstructure statistics + VTKM_LOG_S(vtkm::cont::LogLevel::Info, + std::endl + << ct.PrintHyperStructureStatistics(false) + << std::endl); + // Flush ouput streams just to make sure everything has been logged (in particular when using MPI) std::cout << std::flush; std::cerr << std::flush; diff --git a/vtkm/worklet/contourtree_augmented/ContourTree.h b/vtkm/worklet/contourtree_augmented/ContourTree.h index 647244b04..ba130e972 100644 --- a/vtkm/worklet/contourtree_augmented/ContourTree.h +++ b/vtkm/worklet/contourtree_augmented/ContourTree.h @@ -57,6 +57,8 @@ #include #include #include +#include +#include // local includes #include @@ -161,15 +163,14 @@ public: inline void Init(vtkm::Id dataSize); // debug routine - inline void DebugPrint(const char* message, const char* fileName, long lineNum); + inline void DebugPrint(const char* message, const char* fileName, long lineNum) const; // print contents inline void PrintContent() const; // print routines - //void PrintDotHyperStructure(); - inline void PrintDotSuperStructure(); - //void PrintDotRegularStructure(); + inline void PrintDotSuperStructure() const; + inline std::string PrintHyperStructureStatistics(bool print = true) const; }; // class ContourTree @@ -217,7 +218,7 @@ inline void ContourTree::PrintContent() const PrintIndices("Augmentarcs", this->Augmentarcs); } -void ContourTree::DebugPrint(const char* message, const char* fileName, long lineNum) +void ContourTree::DebugPrint(const char* message, const char* fileName, long lineNum) const { // DebugPrint() #ifdef DEBUG_PRINT std::cout << "---------------------------" << std::endl; @@ -237,13 +238,7 @@ void ContourTree::DebugPrint(const char* message, const char* fileName, long lin #endif } // DebugPrint() -// print routines -// void ContourTree::PrintDotHyperStructure() -// { // PrintDotHyperStructure() -// } // PrintDotHyperStructure() - - -void ContourTree::PrintDotSuperStructure() +void ContourTree::PrintDotSuperStructure() const { // PrintDotSuperStructure() // print the header information printf("digraph G\n\t{\n"); @@ -330,9 +325,80 @@ void ContourTree::PrintDotSuperStructure() printf("\t}\n"); } // PrintDotSuperStructure() -// void ContourTree::PrintDotRegularStructure() -// { // PrintDotRegularStructure() -// } // PrintDotRegularStructure() +std::string ContourTree::PrintHyperStructureStatistics(bool print) const +{ // PrintHyperStructureStatistics() + // arrays for collecting statistics + std::vector minPath; + std::vector maxPath; + std::vector supernodeCount; + std::vector hypernodeCount; + auto whenTransferredPortal = this->WhenTransferred.ReadPortal(); + auto hypernodesPortal = this->Hypernodes.ReadPortal(); + + // set an initial iteration number to negative to get it started + long whichIteration = -1; + + // loop through the hypernodes + for (vtkm::Id hypernode = 0; hypernode < this->Hypernodes.GetNumberOfValues(); hypernode++) + { // per hypernode + // retrieve corresponding supernode ID + vtkm::Id supernodeID = hypernodesPortal.Get(hypernode); + // and the iteration of transfer + vtkm::Id iterationNo = MaskedIndex(whenTransferredPortal.Get(supernodeID)); + + // if it doesn't match, we've hit a boundary + if (whichIteration != iterationNo) + { // new iteration + // initialise the next iteration + // this one is larger than the maximum possible to force minimum + minPath.push_back(static_cast(this->Supernodes.GetNumberOfValues() + 1)); + maxPath.push_back(0); + supernodeCount.push_back(0); + hypernodeCount.push_back(0); + // and increment the iteration ID + whichIteration++; + } // new iteration + + // now compute the new path length - default to off the end + vtkm::Id pathLength = static_cast(this->Supernodes.GetNumberOfValues() - supernodeID); + // for all except the last, take the next one + if (hypernode != this->Hypernodes.GetNumberOfValues() - 1) + { + pathLength = hypernodesPortal.Get(hypernode + 1) - supernodeID; + } + // update the statistics + if (pathLength < minPath[static_cast(whichIteration)]) + { + minPath[static_cast(whichIteration)] = pathLength; + } + if (pathLength > maxPath[static_cast(whichIteration)]) + { + maxPath[static_cast(whichIteration)] = pathLength; + } + supernodeCount[static_cast(whichIteration)] += pathLength; + hypernodeCount[static_cast(whichIteration)]++; + } // per hypernode + + // now print out the statistics + std::stringstream resultString; + for (std::size_t iteration = 0; iteration < minPath.size(); iteration++) + { // per iteration + double averagePath = static_cast(supernodeCount[iteration]) / + static_cast(hypernodeCount[iteration]); + resultString << "Iteration: " << iteration << " Hyper: " << hypernodeCount[iteration] + << " Super: " << supernodeCount[iteration] << " Min: " << minPath[iteration] + << " Avg: " << averagePath << " Max: " << maxPath[iteration] << std::endl; + } // per iteration + resultString << "Total Hypernodes: " << this->Hypernodes.GetNumberOfValues() + << " Supernodes: " << this->Supernodes.GetNumberOfValues() << std::endl; + if (print) + { + std::cout << resultString.str() << std::endl; + } + + return resultString.str(); +} // PrintHyperStructureStatistics() + } // namespace contourtree_augmented } // worklet