diff --git a/examples/contour_tree_distributed/CMakeLists.txt b/examples/contour_tree_distributed/CMakeLists.txt index f5ceeac9d..0f91d470f 100644 --- a/examples/contour_tree_distributed/CMakeLists.txt +++ b/examples/contour_tree_distributed/CMakeLists.txt @@ -59,7 +59,7 @@ find_package(VTKm REQUIRED QUIET) # MPI #################################### if (VTKm_ENABLE_MPI) - add_executable(ContourTree_Distributed ContourTreeApp.cxx) + add_executable(ContourTree_Distributed ContourTreeApp.cxx TreeCompiler.cxx) target_link_libraries(ContourTree_Distributed vtkm_filter vtkm_io MPI::MPI_CXX) vtkm_add_target_information(ContourTree_Distributed MODIFY_CUDA_FLAGS @@ -70,6 +70,7 @@ if (VTKm_ENABLE_MPI) mark_as_advanced(VTKM_EXAMPLE_CONTOURTREE_ENABLE_DEBUG_PRINT) if (VTKM_EXAMPLE_CONTOURTREE_ENABLE_DEBUG_PRINT) target_compile_definitions(ContourTree_Distributed PRIVATE "DEBUG_PRINT") + target_compile_definitions(ContourTree_Distributed PRIVATE "DEBUG_PRINT_CTUD") endif() if (TARGET vtkm::tbb) diff --git a/examples/contour_tree_distributed/ContourTreeApp.cxx b/examples/contour_tree_distributed/ContourTreeApp.cxx index d8e50015e..b14a4f43a 100644 --- a/examples/contour_tree_distributed/ContourTreeApp.cxx +++ b/examples/contour_tree_distributed/ContourTreeApp.cxx @@ -98,8 +98,9 @@ VTKM_THIRDPARTY_POST_INCLUDE #include #include +#include "TreeCompiler.h" + #define PRESPLIT_FILE -using ValueType = vtkm::Float32; namespace ctaug_ns = vtkm::worklet::contourtree_augmented; @@ -154,7 +155,42 @@ private: std::vector mCLOptions; }; +#if 0 +struct PrintArrayContentsFunctor +{ + template + VTKM_CONT void operator()(const vtkm::cont::ArrayHandle& array) const + { + this->PrintArrayPortal(array.GetPortalConstControl()); + } +private: + template + VTKM_CONT void PrintArrayPortal(const PortalType& portal) const + { + for (vtkm::Id index = 0; index < portal.GetNumberOfValues(); index++) + { + // All ArrayPortal objects have ValueType for the type of each value. + using ValueType = typename PortalType::ValueType; + + ValueType value = portal.Get(index); + vtkm::IdComponent numComponents = + vtkm::VecTraits ::GetNumberOfComponents(value); + for (vtkm::IdComponent componentIndex = 0; componentIndex < numComponents; componentIndex ++) + { + std::cout << " " << vtkm::VecTraits::GetComponent(value, componentIndex); + } + //std::cout << std::endl; + } + } +}; + +template +void PrintArrayContents(const VariantArrayType& array) +{ + array.CastAndCall(PrintArrayContentsFunctor()); +} +#endif // Compute and render an isosurface for a uniform grid example int main(int argc, char* argv[]) @@ -471,6 +507,7 @@ int main(int argc, char* argv[]) } // Read data + using ValueType = vtkm::FloatDefault; std::vector values(numVertices); if (filename.compare(filename.length() - 5, 5, ".bdem") == 0) { @@ -534,6 +571,19 @@ int main(int argc, char* argv[]) static_cast(dims[1]), static_cast(nDims == 3 ? dims[2] : 0) }); +#ifdef DEBUG_PRINT_CTUD + std::cout << "blockIndex: " + << vtkm::Id3{ static_cast(std::ceil(static_cast(offset[0]) / + static_cast(dims[0]))), + static_cast(std::ceil(static_cast(offset[1]) / + static_cast(dims[1]))), + static_cast(nDims == 3 + ? std::ceil(static_cast(offset[2]) / + static_cast(dims[2])) + : 0) } + << std::endl; +#endif + if (blockNo == 0) { // FIXME: Hack: Approximate blocks per dim @@ -546,6 +596,9 @@ int main(int argc, char* argv[]) static_cast(dims[2])) : 1) }; +#ifdef DEBUG_PRINT_CTUD + std::cout << "blocksPerDim: " << blocksPerDim << std::endl; +#endif } } @@ -786,10 +839,40 @@ int main(int argc, char* argv[]) useMarchingCubes); filter.SetActiveField("values"); - // Execute the contour tree analysis. NOTE: If MPI is used the result will be - // a vtkm::cont::PartitionedDataSet instead of a vtkm::cont::DataSet + // Execute the contour tree analysis auto result = filter.Execute(useDataSet); +#if 0 + std::cout << "Result dataset has " << result.GetNumberOfPartitions() << " partitions" << std::endl; + + for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no) + { + auto ds = result.GetPartition(ds_no); + for (vtkm::Id f_no = 0; f_no < ds.GetNumberOfFields(); ++f_no) + { + auto field = ds.GetField(f_no); + std::cout << field.GetName() << ": "; + PrintArrayContents(field.GetData()); + std::cout << std::endl; + } + } +#endif + + for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no) + { + TreeCompiler treeCompiler; + treeCompiler.AddHierarchicalTree(result.GetPartition(ds_no)); + char fname[256]; + std::snprintf(fname, + sizeof(fname), + "TreeComplilerOutput_Rank%d_Block%d.dat", + rank, + static_cast(ds_no)); + FILE* out_file = std::fopen(fname, "wb"); + treeCompiler.WriteBinary(out_file); + std::fclose(out_file); + } + currTime = totalTime.GetElapsedTime(); vtkm::Float64 computeContourTreeTime = currTime - prevTime; prevTime = currTime; diff --git a/examples/contour_tree_distributed/TreeCompiler.cxx b/examples/contour_tree_distributed/TreeCompiler.cxx new file mode 100644 index 000000000..5ce7941c6 --- /dev/null +++ b/examples/contour_tree_distributed/TreeCompiler.cxx @@ -0,0 +1,239 @@ +#include "TreeCompiler.h" +#include +#include + +#define PRINT_WIDTH 12 + +// stream output +std::ostream& operator<<(std::ostream& outStream, SupernodeOnSuperarc& node) +{ // stream output + outStream << node.lowEnd << " " << node.highEnd << " " << node.dataValue << " " << node.globalID + << std::endl; + return outStream; +} // stream output + +// stream input +std::istream& operator>>(std::istream& inStream, SupernodeOnSuperarc& node) +{ // stream input + inStream >> node.lowEnd >> node.highEnd >> node.dataValue >> node.globalID; + return inStream; +} // stream input + +// routine to add a known hierarchical tree to it +// note that this DOES NOT finalise - we don't want too many sorts +void TreeCompiler::AddHierarchicalTree(const vtkm::cont::DataSet& addedTree) +{ // TreeCompiler::AddHierarchicalTree() + // Copy relevant tree content to STL arrays + vtkm::cont::VariantArrayHandle dataValues_array = addedTree.GetField("DataValues").GetData(); + std::vector dataValues(dataValues_array.GetNumberOfValues()); + auto dataValues_handle = vtkm::cont::make_ArrayHandle(dataValues, vtkm::CopyFlag::Off); + vtkm::cont::ArrayCopy(dataValues_array.ResetTypes(vtkm::List{}), + dataValues_handle); + dataValues_handle.SyncControlArray(); + + auto regularNodeGlobalIds_array = addedTree.GetField("RegularNodeGlobalIds").GetData(); + std::vector regularNodeGlobalIds(regularNodeGlobalIds_array.GetNumberOfValues()); + auto regularNodeGlobalIds_handle = + vtkm::cont::make_ArrayHandle(regularNodeGlobalIds, vtkm::CopyFlag::Off); + vtkm::cont::ArrayCopy(regularNodeGlobalIds_array.ResetTypes(vtkm::List{}), + regularNodeGlobalIds_handle); + regularNodeGlobalIds_handle + .SyncControlArray(); //Forces values to get updated if copy happened on GPU + + auto superarcs_array = addedTree.GetField("Superarcs").GetData(); + std::vector added_tree_superarcs(superarcs_array.GetNumberOfValues()); + auto superarcs_handle = vtkm::cont::make_ArrayHandle(added_tree_superarcs, vtkm::CopyFlag::Off); + vtkm::cont::ArrayCopy(superarcs_array.ResetTypes(vtkm::List{}), superarcs_handle); + superarcs_handle.SyncControlArray(); //Forces values to get updated if copy happened on GPU + + auto supernodes_array = addedTree.GetField("Supernodes").GetData(); + std::vector added_tree_supernodes(supernodes_array.GetNumberOfValues()); + auto supernodes_handle = vtkm::cont::make_ArrayHandle(added_tree_supernodes, vtkm::CopyFlag::Off); + vtkm::cont::ArrayCopy(supernodes_array.ResetTypes(vtkm::List{}), supernodes_handle); + supernodes_handle.SyncControlArray(); //Forces values to get updated if copy happened on GPU + + auto superparents_array = addedTree.GetField("Superparents").GetData(); + std::vector superparents(superparents_array.GetNumberOfValues()); + auto superparents_handle = vtkm::cont::make_ArrayHandle(superparents, vtkm::CopyFlag::Off); + vtkm::cont::ArrayCopy(superparents_array.ResetTypes(vtkm::List{}), superparents_handle); + superparents_handle.SyncControlArray(); //Forces values to get updated if copy happened on GPU + + // loop through all of the supernodes in the hierarchical tree + for (indexType supernode = 0; supernode < added_tree_supernodes.size(); supernode++) + { // per supernode + // retrieve the regular ID for the supernode + indexType regularId = added_tree_supernodes[supernode]; + indexType globalId = regularNodeGlobalIds[regularId]; + dataType dataVal = dataValues[regularId]; + + // retrieve the supernode at the far end + indexType superTo = added_tree_superarcs[supernode]; + + // now test - if it is NO_SUCH_ELEMENT, there are two possibilities + if (vtkm::worklet::contourtree_augmented::NoSuchElement(superTo)) + { // no Superto + + // retrieve the superparent + indexType superparent = superparents[regularId]; + + + // the root node will have itself as its superparent + if (superparent == supernode) + continue; + else + { // not own superparent - an attachment point + // retrieve the superparent's from & to + indexType regularFrom = added_tree_supernodes[superparent]; + indexType globalFrom = regularNodeGlobalIds[regularFrom]; + indexType superParentTo = added_tree_superarcs[superparent]; + indexType regularTo = + added_tree_supernodes[vtkm::worklet::contourtree_augmented::MaskedIndex(superParentTo)]; + indexType globalTo = regularNodeGlobalIds[regularTo]; + + // test the superTo to see whether we ascend or descend + // note that we will never have NO_SUCH_ELEMENT here + if (vtkm::worklet::contourtree_augmented::IsAscending(superParentTo)) + { // ascending + this->supernodes.push_back(SupernodeOnSuperarc(globalId, dataVal, globalFrom, globalTo)); + } // ascending + else + { // descending + this->supernodes.push_back(SupernodeOnSuperarc(globalId, dataVal, globalTo, globalFrom)); + } // descending + } // not own superparent - an attachment point + } // no Superto + else + { // there is a valid superarc + // retrieve the "to" and convert to global + indexType maskedTo = vtkm::worklet::contourtree_augmented::MaskedIndex(superTo); + indexType regularTo = added_tree_supernodes[maskedTo]; + indexType globalTo = regularNodeGlobalIds[regularTo]; + dataType dataTo = dataValues[regularTo]; + + // test the superTo to see whether we ascend or descend + // note that we will never have NO_SUCH_ELEMENT here + // we add both ends + if (vtkm::worklet::contourtree_augmented::IsAscending(superTo)) + { // ascending + this->supernodes.push_back(SupernodeOnSuperarc(globalId, dataVal, globalId, globalTo)); + this->supernodes.push_back(SupernodeOnSuperarc(globalTo, dataTo, globalId, globalTo)); + } // ascending + else + { // descending + this->supernodes.push_back(SupernodeOnSuperarc(globalId, dataVal, globalTo, globalId)); + this->supernodes.push_back(SupernodeOnSuperarc(globalTo, dataTo, globalTo, globalId)); + } // descending + } // there is a valid superarc + } // per supernode + +} // TreeCompiler::AddHierarchicalTree() + +// routine to compute the actual superarcs +void TreeCompiler::ComputeSuperarcs() +{ // TreeCompiler::ComputeSuperarcs() + // first we sort the vector + std::sort(supernodes.begin(), supernodes.end()); + + // we could do a unique test here, but it's easier just to suppress it inside the loop + + // now we loop through it: note the -1 + // this is because we know a priori that the last one is the last supernode on a superarc + // and would fail the test inside the loop. By putting it in the loop test, we avoid having + // to have an explicit if statement inside the loop + for (indexType supernode = 0; supernode < supernodes.size() - 1; supernode++) + { // loop through supernodes + // this is actually painfully simple: if the (lowEnd, highEnd) don't match the next one, + // then we're at the end of the group and do nothing. Otherwise, we link to the next one + if ((supernodes[supernode].lowEnd != supernodes[supernode + 1].lowEnd) || + (supernodes[supernode].highEnd != supernodes[supernode + 1].highEnd)) + continue; + + // if the supernode matches, then we have a repeat, and can suppress + if (supernodes[supernode].globalID == supernodes[supernode + 1].globalID) + continue; + + // otherwise, add a superarc to the list + superarcs.push_back(Edge(supernodes[supernode].globalID, supernodes[supernode + 1].globalID)); + } // loop through supernodes + + // now sort them + std::sort(superarcs.begin(), superarcs.end()); +} // TreeCompiler::ComputeSuperarcs() + +// routine to print the superarcs +void TreeCompiler::PrintSuperarcs() +{ // TreeCompiler::PrintSuperarcs() + std::cout << "============" << std::endl; + std::cout << "Contour Tree" << std::endl; + + for (indexType superarc = 0; superarc < superarcs.size(); superarc++) + { // per superarc + if (superarcs[superarc].low < superarcs[superarc].high) + { // order by ID not value + std::cout << std::setw(PRINT_WIDTH) << superarcs[superarc].low << " "; + std::cout << std::setw(PRINT_WIDTH) << superarcs[superarc].high << std::endl; + } // order by ID not value + else + { // order by ID not value + std::cout << std::setw(PRINT_WIDTH) << superarcs[superarc].high << " "; + std::cout << std::setw(PRINT_WIDTH) << superarcs[superarc].low << std::endl; + } // order by ID not value + + } // per superarc + +} // TreeCompiler::PrintSuperarcs() + +// routine to write out binary file +void TreeCompiler::WriteBinary(FILE* outFile) +{ // WriteBinary() + // do a bulk write of the entire contents + // no error checking, no type checking, no nothing + fwrite(&(supernodes[0]), sizeof(SupernodeOnSuperarc), supernodes.size(), outFile); +} // WriteBinary() + +// routine to read in binary file and append +void TreeCompiler::ReadBinary(FILE* inFile) +{ // ReadBinary() + // use fseek to jump to the end + fseek(inFile, 0, SEEK_END); + + // use fTell to retrieve the size of the file + long long nBytes = ftell(inFile); + // now rewind + rewind(inFile); + + // compute how many elements are to be read + long long nSupernodes = nBytes / sizeof(SupernodeOnSuperarc); + + // retrieve the current size + long long currentSize = supernodes.size(); + + // resize to add the right number + supernodes.resize(currentSize + nSupernodes); + + // now read directly into the right chunk + fread(&(supernodes[currentSize]), sizeof(SupernodeOnSuperarc), nSupernodes, inFile); + +} // ReadBinary() + +// stream output - just dumps the supernodeonsuperarcs +std::ostream& operator<<(std::ostream& outStream, TreeCompiler& tree) +{ // stream output + for (indexType supernode = 0; supernode < tree.supernodes.size(); supernode++) + outStream << tree.supernodes[supernode]; + return outStream; +} // stream output + +// stream input - reads in the supernodeonsuperarcs & appends them +std::istream& operator>>(std::istream& inStream, TreeCompiler& tree) +{ // stream input + while (!inStream.eof()) + { + SupernodeOnSuperarc tempNode; + inStream >> tempNode; + tree.supernodes.push_back(tempNode); + } + // we will overshoot, so subtract one + tree.supernodes.resize(tree.supernodes.size() - 1); + return inStream; +} // stream input diff --git a/examples/contour_tree_distributed/TreeCompiler.h b/examples/contour_tree_distributed/TreeCompiler.h new file mode 100644 index 000000000..19b279a8e --- /dev/null +++ b/examples/contour_tree_distributed/TreeCompiler.h @@ -0,0 +1,138 @@ +#ifndef _TREECOMPILER_H_ +#define _TREECOMPILER_H_ + +#include +#include +#include +#include + +// FIXME/HACK: Define here for compatibility with PPP TreeCompiler +typedef double dataType; +typedef unsigned long long indexType; + +// small class for storing the contour arcs +class Edge +{ // Edge +public: + indexType low, high; + + // constructor - defaults to -1 + Edge(vtkm::Id Low = -1, vtkm::Id High = -1) + : low(Low) + , high(High) + { + } +}; // Edge + +// comparison operator < +inline bool operator<(const Edge LHS, const Edge RHS) +{ // operator < + if (LHS.low < RHS.low) + return true; + if (LHS.low > RHS.low) + return false; + if (LHS.high < RHS.high) + return true; + if (LHS.high > RHS.high) + return false; + return false; +} // operator < + +// a helper class which stores a single supernode inserted onto a superarc +class SupernodeOnSuperarc +{ // class SupernodeOnSuperarc +public: + // the global ID of the supernode + indexType globalID; + // the data value stored at the supernode + dataType dataValue; + + // the low and high ends of the superarc it is on (may be itself) + indexType lowEnd, highEnd; + + // constructor + SupernodeOnSuperarc(indexType GlobalID = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT, + dataType DataValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT, + indexType LowEnd = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT, + indexType HighEnd = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT) + : globalID(GlobalID) + , dataValue(DataValue) + , lowEnd(LowEnd) + , highEnd(HighEnd) + { // constructor + } // constructor +}; // class SupernodeOnSuperarc + +// overloaded comparison operator +// primary sort is by superarc (low, high), +// then secondary sort on datavalue +// tertiary on globalID to implement simulated simplicity +inline bool operator<(const SupernodeOnSuperarc& left, const SupernodeOnSuperarc& right) +{ // < operator + // simple lexicographic sort + if (left.lowEnd < right.lowEnd) + return true; + if (left.lowEnd > right.lowEnd) + return false; + if (left.highEnd < right.highEnd) + return true; + if (left.highEnd > right.highEnd) + return false; + if (left.dataValue < right.dataValue) + return true; + if (left.dataValue > right.dataValue) + return false; + if (left.globalID < right.globalID) + return true; + if (left.globalID > right.globalID) + return false; + + // fall-through (shouldn't happen, but) + // if they're the same, it's false + return false; +} // < operator + +// stream output +std::ostream& operator<<(std::ostream& outStream, SupernodeOnSuperarc& node); + +// stream input +std::istream& operator>>(std::istream& inStream, SupernodeOnSuperarc& node); + +// the class that compiles the contour tree +class TreeCompiler +{ // class TreeCompiler +public: + // we want a vector of supernodes on superarcs + std::vector supernodes; + + // and a vector of Edges (the output) + std::vector superarcs; + + // default constructor sets it to empty + TreeCompiler() + { // constructor + // clear out the supernode array + supernodes.resize(0); + // and the superarc array + superarcs.resize(0); + } // constructor + + // routine to add a known hierarchical tree to it + // note that this DOES NOT finalise - we don't want too many sorts + void AddHierarchicalTree(const vtkm::cont::DataSet& addedTree); + + // routine to compute the actual superarcs + void ComputeSuperarcs(); + + // routine to print the superarcs + void PrintSuperarcs(); + + // routine to write out binary file + void WriteBinary(FILE* outFile); + + // routine to read in binary file & append to contents + void ReadBinary(FILE* inFile); + +}; // class TreeCompiler + +#endif diff --git a/vtkm/filter/ContourTreeUniformDistributed.hxx b/vtkm/filter/ContourTreeUniformDistributed.hxx index 95ccfd0d7..0e59436c9 100644 --- a/vtkm/filter/ContourTreeUniformDistributed.hxx +++ b/vtkm/filter/ContourTreeUniformDistributed.hxx @@ -283,6 +283,7 @@ void ContourTreeUniformDistributed::ComputeLocalTreeImpl( boundaryTreeMaker.Construct(&localToGlobalIdRelabeler); +#ifdef DEBUG_PRINT_CTUD // Print BRACT for debug vtkm::Id rank = vtkm::cont::EnvironmentTracker::GetCommunicator().rank(); char buffer[256]; @@ -300,7 +301,7 @@ void ContourTreeUniformDistributed::ComputeLocalTreeImpl( this->MultiBlockSpatialDecomposition.LocalBlockSizes.ReadPortal().Get(blockIndex), this->MultiBlockSpatialDecomposition.GlobalSize) << std::endl; - +#endif // TODO: Fix calling conventions so dot print works here // At this point, I'm reasonably certain that the contour tree has been computed regardless of data push/pull // So although it might be logical to print things out earlier, I'll do it here @@ -378,7 +379,7 @@ inline VTKM_CONT void ContourTreeUniformDistributed::PreExecute( template VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( const vtkm::cont::PartitionedDataSet& input, - vtkm::cont::PartitionedDataSet&, // temporarily commented out: output, + vtkm::cont::PartitionedDataSet& result, const vtkm::filter::FieldMetadata&, // dummy parameter for field meta data const vtkm::cont::ArrayHandle&, // dummy parameter to get the type vtkm::filter::PolicyBase) // dummy parameter for policy @@ -440,6 +441,7 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( fieldData, localGlobalMeshIndex); +#ifdef DEBUG_PRINT_CTUD // ... save for debugging in text and .gv/.dot format char buffer[256]; std::snprintf(buffer, @@ -460,6 +462,7 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( std::string(" Initial Step 5 BRACT Mesh"), localDataBlocks[bi]->ContourTreeMeshes.back(), worklet::contourtree_distributed::SHOW_CONTOUR_TREE_MESH_ALL); +#endif } // Setup vtkmdiy to do global binary reduction of neighbouring blocks. See also RecuctionOperation struct for example @@ -533,9 +536,12 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( comm.barrier(); // Be safe! + std::vector hierarchicalTreeOutputDataSet( + localDataBlocks.size()); // DataSets for creating output data master.foreach ([&](vtkm::worklet::contourtree_distributed::DistributedContourTreeBlockData< FieldType>* b, const vtkmdiy::Master::ProxyWithLink&) { +#ifdef DEBUG_PRINT_CTUD VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Fan In Complete" << std::endl << "# of CTs: " << b->ContourTrees.size() << std::endl @@ -564,6 +570,7 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( for (const auto& info : b->InteriorForests) info.PrintContent(os); os << std::endl; +#endif // Fan out auto nRounds = b->ContourTrees.size() - 1; @@ -619,6 +626,24 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( this->MultiBlockSpatialDecomposition.GlobalSize); grafter.GraftInteriorForests(0, hierarchicalTree, fieldData, &localToGlobalIdRelabeler); + // Create data set from output + vtkm::cont::Field dataValuesField( + "DataValues", vtkm::cont::Field::Association::WHOLE_MESH, hierarchicalTree.DataValues); + hierarchicalTreeOutputDataSet[b->BlockIndex].AddField(dataValuesField); + vtkm::cont::Field regularNodeGlobalIdsField("RegularNodeGlobalIds", + vtkm::cont::Field::Association::WHOLE_MESH, + hierarchicalTree.RegularNodeGlobalIds); + hierarchicalTreeOutputDataSet[b->BlockIndex].AddField(regularNodeGlobalIdsField); + vtkm::cont::Field superarcsField( + "Superarcs", vtkm::cont::Field::Association::WHOLE_MESH, hierarchicalTree.Superarcs); + hierarchicalTreeOutputDataSet[b->BlockIndex].AddField(superarcsField); + vtkm::cont::Field supernodesField( + "Supernodes", vtkm::cont::Field::Association::WHOLE_MESH, hierarchicalTree.Supernodes); + hierarchicalTreeOutputDataSet[b->BlockIndex].AddField(supernodesField); + vtkm::cont::Field superparentsField( + "Superparents", vtkm::cont::Field::Association::WHOLE_MESH, hierarchicalTree.Superparents); + hierarchicalTreeOutputDataSet[b->BlockIndex].AddField(superparentsField); + // TODO: GET THIS COMPILING // save the corresponding .gv file // std::string hierarchicalTreeFileName = std::string("Rank_") + std::to_string(static_cast(rank)) + std::string("_Block_") + std::to_string(static_cast(b->BlockIndex)) + "_Round_" + std::to_string(nRounds) + "_Hierarchical_Tree.gv"; @@ -629,56 +654,7 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( // ); }); - // FIXME: The following is still old code and needs to be changed!!! - -#if 0 - // Now run the contour tree algorithm on the last block to compute the final tree - vtkm::Id currNumIterations; - vtkm::worklet::contourtree_augmented::ContourTree currContourTree; - vtkm::worklet::contourtree_augmented::IdArrayType currSortOrder; - vtkm::worklet::ContourTreeAugmented worklet; - vtkm::cont::ArrayHandle currField; - // Construct the mesh boundary exectuion object needed for boundary augmentation - vtkm::Id3 minIdx(0, 0, 0); - vtkm::Id3 maxIdx = this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.GlobalSize; - maxIdx[0] = maxIdx[0] - 1; - maxIdx[1] = maxIdx[1] - 1; - maxIdx[2] = maxIdx[2] > 0 ? (maxIdx[2] - 1) : 0; - // TODO Check setting. In parallel we should need to augment the tree in order to be able to do the merging - const unsigned int computeRegularStructure = 1; - auto meshBoundaryExecObj = localDataBlocks[0]->ContourTreeMeshes.back().GetMeshBoundaryExecutionObject( - this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.GlobalSize[0], - this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.GlobalSize[1], - minIdx, - maxIdx); - // Run the worklet to compute the final contour tree - worklet.Run( - localDataBlocks[0]->ContourTreeMeshes.back().SortedValues, // Unused param. Provide something to keep API happy - localDataBlocks[0]->ContourTreeMeshes.back(), - this->ContourTreeData, - this->MeshSortOrder, - currNumIterations, - computeRegularStructure, - meshBoundaryExecObj); - - // Set the final mesh sort order we need to use - this->MeshSortOrder = localDataBlocks[0]->ContourTreeMeshes.back().GlobalMeshIndex; - // Remeber the number of iterations for the output - this->NumIterations = currNumIterations; - - // Return the sorted values of the contour tree as the result - // TODO the result we return for the parallel and serial case are different right now. This should be made consistent. However, only in the parallel case are we useing the result output - vtkm::cont::DataSet temp; - vtkm::cont::Field rfield(this->GetOutputFieldName(), - vtkm::cont::Field::Association::WHOLE_MESH, - localDataBlocks[0]->ContourTreeMeshes.back().SortedValues); - temp.AddField(rfield); - output = vtkm::cont::PartitionedDataSet(temp); - /* Not necessary? - this->ContourTreeData = MultiBlockTreeHelper->LocalContourTrees[0]; - this->MeshSortOrder = MultiBlockTreeHelper->LocalSortOrders[0]; - */ -#endif + result = vtkm::cont::PartitionedDataSet(hierarchicalTreeOutputDataSet); } diff --git a/vtkm/worklet/contourtree_distributed/ComputeDistributedContourTreeFunctor.h b/vtkm/worklet/contourtree_distributed/ComputeDistributedContourTreeFunctor.h index 1e634cf08..8b86716db 100644 --- a/vtkm/worklet/contourtree_distributed/ComputeDistributedContourTreeFunctor.h +++ b/vtkm/worklet/contourtree_distributed/ComputeDistributedContourTreeFunctor.h @@ -100,10 +100,11 @@ public: ) const { const auto selfid = rp.gid(); - //DEBUGOUT//std::cout << "Entering ComputeDistributedContourTreeFunctor::operator() for GID " << selfid << std::endl; +#ifdef DEBUG_PRINT_CTUD // Get rank (for debug output only) const vtkm::Id rank = vtkm::cont::EnvironmentTracker::GetCommunicator().rank(); +#endif // Here we do the deque first before the send due to the way the iteration is handled in DIY, i.e., in each iteration // A block needs to first collect the data from its neighours and then send the combined block to its neighbours @@ -124,6 +125,9 @@ public: rp.dequeue(ingid, otherBlockSize); vtkm::worklet::contourtree_augmented::ContourTreeMesh otherContourTreeMesh; rp.dequeue(ingid, otherContourTreeMesh); + +#ifdef DEBUG_PRINT_CTUD + char buffer[256]; #if 0 // FIXME: Delete after debugging std::cout << "Our block has extents: " << block->BlockOrigin << " " << block->BlockSize << std::endl; @@ -133,20 +137,19 @@ public: block->ContourTreeMeshes.back().DebugPrint("OUR CTM", __FILE__, __LINE__); otherContourTreeMesh.DebugPrint("OTHER CTM", __FILE__, __LINE__); #endif - #if 0 - char buffer[256]; std::snprintf(buffer, sizeof(buffer), "BeforeCombine_MyMesh_ContourTreeMesh_Rank%lld_Round%d.txt", rank, rp.round()); block->ContourTreeMeshes.back().Save(buffer); std::snprintf(buffer, sizeof(buffer), "BeforeCombine_Other_ContourTreeMesh_Rank%lld_Round%d.txt", rank, rp.round()); otherContourTreeMesh.Save(buffer); +#endif #endif // Merge the two contour tree meshes vtkm::cont::TryExecute( MergeContourTreeMeshFunctor{}, otherContourTreeMesh, block->ContourTreeMeshes.back()); - char buffer[256]; +#ifdef DEBUG_PRINT_CTUD std::snprintf(buffer, sizeof(buffer), "AfterCombine_ContourTreeMesh_Rank%d_Block%d_Round%d.txt", @@ -154,12 +157,12 @@ public: static_cast(block->BlockIndex), rp.round()); block->ContourTreeMeshes.back().Save(buffer); - #if 0 // FIXME: Delete after debugging std::cout << "================== AFTER COMBINING =================" << std::endl; std::cout << "OUR CTM" << std::endl; block->ContourTreeMeshes.back().DebugPrint("OUR CTM", __FILE__, __LINE__); +#endif #endif // Compute the origin and size of the new block @@ -196,13 +199,15 @@ public: block->ContourTrees.back(), currSortOrder, currNumIterations, - 1, // Fully augmented FIXME: Verify + 1, // Fully augmented meshBoundaryExecObj); +#ifdef DEBUG_PRINT_CTUD #if 0 // FIXME: Delete after debugging std::cout << "=================== BEGIN: COMBINED CONTOUR TREE =========================" << std::endl; block->ContourTrees.back().PrintContent(); std::cout << "=================== END: COMBINED CONTOUR TREE =========================" << std::endl; +#endif #endif // Update block extents @@ -216,22 +221,19 @@ public: // last round) then compute contour tree mesh to send and save it. if (rp.round() != 0 && rp.out_link().size() != 0) { +#ifdef DEBUG_PRINT_CTUD + char buffer[256]; #if 0 - // If we have the fully augmented contour tree - block->ContourTreeMeshes.emplace_back( - block->ContourTrees.back().Arcs, block->ContourTreeMeshes.back() - ); - - char buffer[256]; std::snprintf(buffer, sizeof(buffer), "CombinedMeshes_GID%d_Round%d.txt", selfid, rp.round()); block->ContourTreeMeshes.back().Save(buffer); +#endif #endif vtkm::Id3 maxIdx{ block->BlockOrigin[0] + block->BlockSize[0] - 1, block->BlockOrigin[1] + block->BlockSize[1] - 1, block->BlockOrigin[2] + block->BlockSize[2] - 1 }; - char buffer[256]; +#ifdef DEBUG_PRINT_CTUD std::snprintf(buffer, sizeof(buffer), "BRACTInputCTM_Rank%d_Block%d_Round%d.ctm_txt", @@ -256,28 +258,29 @@ public: // TODO: GET THIS COMPILING // save the corresponding .gv file for the contour tree mesh - // std::string contourTreeMeshFileName = std::string("Rank_") + std::to_string(static_cast(rank)) + std::string("_Block_") + std::to_string(static_cast(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_0_Combined_Mesh.gv"); - // std::ofstream contourTreeMeshFile(contourTreeMeshFileName); - // contourTreeMeshFile << vtkm::worklet::contourtree_distributed::ContourTreeMeshDotGraphPrint - // ( std::string("Block ") + std::to_string(static_cast(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 0 Combined Mesh"), - // block->ContourTreeMeshes.back(), worklet::contourtree_distributed::SHOW_CONTOUR_TREE_MESH_ALL); - // - // // and the ones for the contour tree regular and superstructures - // std::string regularStructureFileName = std::string("Rank_") + std::to_string(static_cast(rank)) + std::string("_Block_") + std::to_string(static_cast(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string("_Step_1_Contour_Tree_Regular_Structure.gv"); - // std::ofstream regularStructureFile(regularStructureFileName); - // regularStructureFile << worklet::contourtree_distributed::ContourTreeDotGraphPrint - // ( std::string("Block ") + std::to_string(static_cast(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 1 Contour Tree Regular Structure"), - // block->Meshes.back(), - // block->ContourTrees.back(), - // worklet::contourtree_distributed::SHOW_REGULAR_STRUCTURE|worklet::contourtree_distributed::SHOW_ALL_IDS); - // - // std::string superStructureFileName = std::string("Rank_") + std::to_string(static_cast(rank)) + std::string("_Block_") + std::to_string(static_cast(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string("_Step_2_Contour_Tree_Super_Structure.gv"); - // std::ofstream superStructureFile(superStructureFileName); - // superStructureFile << worklet::contourtree_distributed::ContourTreeDotGraphPrint - // ( std::string("Block ") + std::to_string(static_cast(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 2 Contour Tree Super Structure"), - // block->Meshes.back(), - // block->ContourTrees.back(), - // worklet::contourtree_distributed::SHOW_SUPER_STRUCTURE|worklet::contourtree_distributed::SHOW_HYPER_STRUCTURE|worklet::contourtree_distributed::SHOW_ALL_IDS|worklet::contourtree_distributed::SHOW_ALL_SUPERIDS|worklet::contourtree_distributed::SHOW_ALL_HYPERIDS); +// std::string contourTreeMeshFileName = std::string("Rank_") + std::to_string(static_cast(rank)) + std::string("_Block_") + std::to_string(static_cast(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_0_Combined_Mesh.gv"); +// std::ofstream contourTreeMeshFile(contourTreeMeshFileName); +// contourTreeMeshFile << vtkm::worklet::contourtree_distributed::ContourTreeMeshDotGraphPrint +// ( std::string("Block ") + std::to_string(static_cast(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 0 Combined Mesh"), +// block->ContourTreeMeshes.back(), worklet::contourtree_distributed::SHOW_CONTOUR_TREE_MESH_ALL); +// +// // and the ones for the contour tree regular and superstructures +// std::string regularStructureFileName = std::string("Rank_") + std::to_string(static_cast(rank)) + std::string("_Block_") + std::to_string(static_cast(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string("_Step_1_Contour_Tree_Regular_Structure.gv"); +// std::ofstream regularStructureFile(regularStructureFileName); +// regularStructureFile << worklet::contourtree_distributed::ContourTreeDotGraphPrint +// ( std::string("Block ") + std::to_string(static_cast(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 1 Contour Tree Regular Structure"), +// block->Meshes.back(), +// block->ContourTrees.back(), +// worklet::contourtree_distributed::SHOW_REGULAR_STRUCTURE|worklet::contourtree_distributed::SHOW_ALL_IDS); +// +// std::string superStructureFileName = std::string("Rank_") + std::to_string(static_cast(rank)) + std::string("_Block_") + std::to_string(static_cast(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string("_Step_2_Contour_Tree_Super_Structure.gv"); +// std::ofstream superStructureFile(superStructureFileName); +// superStructureFile << worklet::contourtree_distributed::ContourTreeDotGraphPrint +// ( std::string("Block ") + std::to_string(static_cast(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 2 Contour Tree Super Structure"), +// block->Meshes.back(), +// block->ContourTrees.back(), +// worklet::contourtree_distributed::SHOW_SUPER_STRUCTURE|worklet::contourtree_distributed::SHOW_HYPER_STRUCTURE|worklet::contourtree_distributed::SHOW_ALL_IDS|worklet::contourtree_distributed::SHOW_ALL_SUPERIDS|worklet::contourtree_distributed::SHOW_ALL_HYPERIDS); +#endif // Compute BRACT vtkm::worklet::contourtree_distributed::BoundaryTree boundaryTree; @@ -286,7 +289,6 @@ public: this->GlobalSize, block->BlockOrigin, maxIdx); // Make the BRACT and InteriorForest (i.e., residue) block->InteriorForests.emplace_back(); - std::cout << "Calling BRACT Maker" << std::endl; auto boundaryTreeMaker = vtkm::worklet::contourtree_distributed::BoundaryTreeMaker< vtkm::worklet::contourtree_augmented::ContourTreeMesh, vtkm::worklet::contourtree_augmented::MeshBoundaryContourTreeMeshExec>( @@ -298,6 +300,7 @@ public: // Construct the BRACT and InteriorForest. Since we are working on a ContourTreeMesh we do // not need to provide and IdRelabeler here in order to compute the InteriorForest boundaryTreeMaker.Construct(); +#ifdef DEBUG_PRINT_CTUD os << "================= OUTPUT ===================" << std::endl; os << "+++++++++++++++++ BRACT +++++++++++++++++++" << std::endl; os << boundaryTree.DebugPrint("validate", __FILE__, __LINE__); @@ -319,9 +322,12 @@ public: rp.round()); std::ofstream dotStream(buffer); dotStream << debug_dot; +#endif // Construct contour tree mesh from BRACT block->ContourTreeMeshes.emplace_back( boundaryTree.VertexIndex, boundaryTree.Superarcs, block->ContourTreeMeshes.back()); + +#ifdef DEBUG_PRINT_CTUD block->ContourTreeMeshes.back().PrintContent(os); //char buffer[256]; std::snprintf(buffer, @@ -331,34 +337,35 @@ public: static_cast(block->BlockIndex), block->ContourTreeMeshes.size() - 1); block->ContourTreeMeshes.back().Save(buffer); - } - // TODO: GET THIS COMPILING - // save the Boundary Tree as a dot file - // std::string boundaryTreeFileName = std::string("Rank_") + std::to_string(static_cast(rank)) + std::string("_Block_") + std::to_string(static_cast(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_3_Boundary_Tree.gv"); - // std::ofstream boundaryTreeFile(boundaryTreeFileName); - // boundaryTreeFile << vtkm::worklet::contourtree_distributed::BoundaryTreeDotGraphPrint - // ( std::string("Block ") + std::to_string(static_cast(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 3 Boundary Tree"), - // block->Meshes.back()], - // block->BoundaryTrees.back()); - // - // // and save the Interior Forest as another dot file - // std::string interiorForestFileName = std::string("Rank_") + std::to_string(static_cast(rank)) + std::string("_Block_") + std::to_string(static_cast(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_4_Interior_Forest.gv"); - // std::ofstream interiorForestFile(interiorForestFileName); - // interiorForestFileName << InteriorForestDotGraphPrintFile - // ( std::string("Block ") + std::to_string(static_cast(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 4 Interior Forest"), - // block->InteriorForests.back(), - // block->ContourTrees.back(), - // block->BoundaryTrees.back(), - // block->Meshes.back()); - // - // // save the corresponding .gv file - // std::string boundaryTreeMeshFileName = std::string("Rank_") + std::to_string(static_cast(rank)) + std::string("_Block_") + std::to_string(static_cast(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_5_Boundary_Tree_Mesh.gv"); - // std::ofstream boundaryTreeMeshFile(boundaryTreeMeshFileName); - // boundaryTreeMeshFile << vtkm::worklet::contourtree_distributed::ContourTreeMeshDotGraphPrint - // ( std::string("Block ") + std::to_string(static_cast(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 5 Boundary Tree Mesh"), - // block->ContourTreeMeshes.back(), - // worklet::contourtree_distributed::SHOW_CONTOUR_TREE_MESH_ALL); + // TODO: GET THIS COMPILING + // save the Boundary Tree as a dot file +// std::string boundaryTreeFileName = std::string("Rank_") + std::to_string(static_cast(rank)) + std::string("_Block_") + std::to_string(static_cast(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_3_Boundary_Tree.gv"); +// std::ofstream boundaryTreeFile(boundaryTreeFileName); +// boundaryTreeFile << vtkm::worklet::contourtree_distributed::BoundaryTreeDotGraphPrint +// ( std::string("Block ") + std::to_string(static_cast(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 3 Boundary Tree"), +// block->Meshes.back()], +// block->BoundaryTrees.back()); +// +// // and save the Interior Forest as another dot file +// std::string interiorForestFileName = std::string("Rank_") + std::to_string(static_cast(rank)) + std::string("_Block_") + std::to_string(static_cast(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_4_Interior_Forest.gv"); +// std::ofstream interiorForestFile(interiorForestFileName); +// interiorForestFileName << InteriorForestDotGraphPrintFile +// ( std::string("Block ") + std::to_string(static_cast(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 4 Interior Forest"), +// block->InteriorForests.back(), +// block->ContourTrees.back(), +// block->BoundaryTrees.back(), +// block->Meshes.back()); +// +// // save the corresponding .gv file +// std::string boundaryTreeMeshFileName = std::string("Rank_") + std::to_string(static_cast(rank)) + std::string("_Block_") + std::to_string(static_cast(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_5_Boundary_Tree_Mesh.gv"); +// std::ofstream boundaryTreeMeshFile(boundaryTreeMeshFileName); +// boundaryTreeMeshFile << vtkm::worklet::contourtree_distributed::ContourTreeMeshDotGraphPrint +// ( std::string("Block ") + std::to_string(static_cast(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 5 Boundary Tree Mesh"), +// block->ContourTreeMeshes.back(), +// worklet::contourtree_distributed::SHOW_CONTOUR_TREE_MESH_ALL); +#endif + } // Send our current block (which is either our original block or the one we just combined from the ones we received) to our next neighbour. // Once a rank has send his block (either in its orignal or merged form) it is done with the reduce