handling multiple domains

Alpine-DAV · Aug 22, 2024 · edffaea · edffaea
1 parent 778db6b
commit edffaea
Showing 1 changed file with 86 additions and 18 deletions.
diff --git a/src/libs/vtkh/filters/ContourTree.cpp b/src/libs/vtkh/filters/ContourTree.cpp
@@ -19,6 +19,69 @@
 
 namespace caugmented_ns = vtkm::worklet::contourtree_augmented;
 
+#ifdef DEBUG
+void SaveToBDEM(const vtkm::cont::DataSet& ds,
+                const std::string& fieldname,
+                const std::string& filename)
+{
+  std::cout << "Saving block to " << filename << std::endl;
+  vtkm::Id3 blockSize;
+  vtkm::Id3 globalSize;
+  vtkm::Id3 pointIndexStart;
+  vtkm::cont::CastAndCall(ds.GetCellSet(),
+                          vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(),
+                          blockSize,
+                          globalSize,
+                          pointIndexStart);
+  std::ofstream os(filename, std::ios::binary);
+  os << "#GLOBAL_EXTENTS " << globalSize[1] << " " << globalSize[0];
+  if (globalSize[2] > 1)
+  {
+    os << " " << globalSize[2];
+  }
+  os << std::endl;
+  os << "#OFFSET " << pointIndexStart[1] << " " << pointIndexStart[0];
+  if (globalSize[2] > 1)
+  {
+    os << " " << pointIndexStart[2];
+  }
+  os << std::endl;
+  os << blockSize[1] << " " << blockSize[0];
+  if (globalSize[2] > 1)
+  {
+    os << " " << blockSize[2];
+  }
+  os << std::endl;
+  try
+  {
+    auto dataArray =
+      ds.GetField(fieldname).GetData().AsArrayHandle<vtkm::cont::ArrayHandleBasic<vtkm::Float32>>();
+    os.write(reinterpret_cast<const char*>(dataArray.GetReadPointer()),
+             dataArray.GetNumberOfValues() * sizeof(vtkm::Float32));
+  }
+  catch (vtkm::cont::ErrorBadType& e)
+  {
+    std::cerr << "Error: Cannot get as ArrayHandleBasic<vtkm::Float32>: " << e.GetMessage()
+              << std::endl;
+  }
+}
+
+void SaveToBDEM(const vtkm::cont::PartitionedDataSet& pds,
+                const std::string& fieldname,
+                const std::string& basefilename,
+                int rank,
+                int blocksPerRank)
+{
+  for (vtkm::Id i = 0; i < pds.GetNumberOfPartitions(); ++i)
+  {
+    std::string filename = basefilename + '_' + std::to_string(rank*blocksPerRank+i) + ".bdem";
+    std::cout << "Saving partition " << i << " to " << filename << std::endl;
+    SaveToBDEM(pds.GetPartition(i), fieldname, filename);
+  }
+}
+#endif
+
+#ifdef VTKH_PARALLEL
 static void ShiftLogicalOriginToZero(vtkm::cont::PartitionedDataSet& pds)
 {
   // Shift the logical origin (minimum of LocalPointIndexStart) to zero
@@ -120,7 +183,6 @@ static void ComputeGlobalPointSize(vtkm::cont::PartitionedDataSet& pds)
   }
 }
 
-#ifdef VTKH_PARALLEL
 #include <mpi.h>
 
 // This is from VTK-m diy mpi_cast.hpp. Need the make_DIY_MPI_Comm
@@ -442,8 +504,18 @@ void ContourTree::distributed_analysis(const vtkm::cont::PartitionedDataSet& res
 
 void ContourTree::DoExecute()
 { 
+{
+#ifdef VTKH_PARALLEL
+  int mpi_rank = 0;
+
+  MPI_Comm mpi_comm = MPI_Comm_f2c(vtkh::GetMPICommHandle());
+  vtkm::cont::EnvironmentTracker::SetCommunicator(vtkmdiy::mpi::communicator(vtkmdiy::mpi::make_DIY_MPI_Comm(mpi_comm)));
+
+  MPI_Comm_rank(mpi_comm, &mpi_rank);
+#endif
+}
+
   vtkh::DataSet *old_input = this->m_input;
-  const int before_num_domains = this->m_input->GetNumberOfDomains();
 
   // make sure we have a node-centered field
   bool valid_field = false;
@@ -485,7 +557,6 @@ void ContourTree::DoExecute()
   this->m_output = new DataSet();
 
   const int num_domains = this->m_input->GetNumberOfDomains();
-  assert(num_domains == 1);
 
 #ifndef VTKH_PARALLEL
   vtkm::cont::DataSet inDataSet;
@@ -506,11 +577,15 @@ void ContourTree::DoExecute()
 
   vtkm::cont::PartitionedDataSet inDataSet;
 
-  vtkm::Id domain_id;
-  vtkm::cont::DataSet dom;
+  for (int i = 0; i < num_domains; ++i)
+  {
+     vtkm::Id domain_id;
+     vtkm::cont::DataSet dom;
 
-  this->m_input->GetDomain(0, dom, domain_id);
-  inDataSet.AppendPartition(dom);
+     this->m_input->GetDomain(i, dom, domain_id);
+     inDataSet.AppendPartition(dom);
+     this->m_output->AddDomain(dom, domain_id);
+  }
 
 #ifdef DEBUG
   if( mpi_size != 1 )
@@ -533,12 +608,6 @@ void ContourTree::DoExecute()
 
   //Convert the mesh of values into contour tree, pairs of vertex ids
   vtkm::filter::scalar_topology::ContourTreeAugmented filter(useMarchingCubes, computeRegularStructure);
-
-#ifdef DEBUG
-  std::cout << "--- BEGIN_SUMMARY inDataSet" << std::endl;
-  inDataSet.PrintSummary( std::cout );
-  std::cout << "--- END_SUMMARY inDataSet" << std::endl;
-#endif
 #else // VTKH_PARALLEL
   // TODO SET VTKM to do vtkm::cont::LogLovel::Info
   // vtkm::filter::scalar_topology::ContourTreeUniformDistributed filter(vtkm::cont::LogLevel::Info, vtkm::cont::LogLevel::Info);
@@ -556,15 +625,14 @@ void ContourTree::DoExecute()
     filterField = &filter;
   }
 
-#ifdef DEBUG
-  std::cout << "--- BEGIN_SUMMARY inDataSet:mpi_rank=" << mpi_rank << std::endl;
-  inDataSet.PrintSummary( std::cout );
-  std::cout << "--- END_SUMMARY inDataSet:mpi_rank=" << mpi_rank << std::endl;
-#endif
 #endif // VTKH_PARALLEL
 
   filterField->SetActiveField(m_field_name);
 
+#ifdef DEBUG
+  SaveToBDEM(inDataSet, m_field_name, "debug-ascent", mpi_rank, num_domains);
+#endif
+
   auto result = filterField->Execute(inDataSet);
 
   m_iso_values.resize(m_levels);