Refactoring coarsening out of base solver

- Fix bug where coarsening across processor boundaries might segfault - Fix bug where dorfler marking would segfault for empty estimate array (occurs for empty partitions or sometimes coarsening)
awslabs · Oct 25, 2023 · 304f330 · 304f330
1 parent 368c643
commit 304f330
Show file tree

Hide file tree

Showing 3 changed files with 64 additions and 33 deletions.
diff --git a/palace/drivers/basesolver.cpp b/palace/drivers/basesolver.cpp
@@ -216,26 +216,8 @@ void BaseSolver::SolveEstimateMarkRefine(
     {
       // Perform a Dörfler style marking looking for the largest number of derefinement
       // opportunities to represent a fraction of the derefinable error.
-      const auto &derefinement_table = mesh.back()->pncmesh->GetDerefinementTable();
-      Vector coarse_error(derefinement_table.Size());
-      mfem::Array<int> row;
-      for (int i = 0; i < derefinement_table.Size(); i++)
-      {
-        // Sum the error for all sub elements that can be combined.
-        derefinement_table.GetRow(i, row);
-        coarse_error[i] =
-            std::sqrt(std::accumulate(row.begin(), row.end(), 0.0,
-                                      [&indicators](double s, int i) {
-                                        return s += std::pow(indicators.Local()[i], 2.0);
-                                      }));
-      }
-
-      // Given the coarse errors, we use the Dörfler marking strategy to identify the
-      // smallest set of original elements that make up (1 - θ) of the total error. The
-      // complement of this set is then the largest number of elements that make up θ of the
-      // total error.
-      const double threshold =
-          utils::ComputeDorflerThreshold(comm, 1 - param.update_fraction, coarse_error);
+      const double threshold = utils::ComputeDorflerCoarseningThreshold(
+          *mesh.back(), param.update_fraction, indicators.Local());
 
       const auto initial_elem_count = mesh.back()->GetGlobalNE();
       constexpr int aggregate_operation = 3;  // sum of squares

diff --git a/palace/utils/dorfler.cpp b/palace/utils/dorfler.cpp
@@ -6,6 +6,7 @@
 #include <algorithm>
 #include <limits>
 #include <numeric>
+#include <mfem.hpp>
 #include "utils/communication.hpp"
 
 namespace palace::utils
@@ -17,7 +18,6 @@ double ComputeDorflerThreshold(MPI_Comm comm, double fraction, const Vector &e)
   e.HostRead();  // Pull the data out to the host
   std::vector<double> estimates(e.begin(), e.end());
   std::sort(estimates.begin(), estimates.end());
-  MFEM_ASSERT(estimates.size() > 0, "Estimates must be non-empty");
 
   // Accumulate the squares of the estimates.
   std::vector<double> sum(estimates.size());
@@ -32,19 +32,27 @@ double ComputeDorflerThreshold(MPI_Comm comm, double fraction, const Vector &e)
   }
 
   // The pivot is the first point which leaves (1-θ) of the total sum after it.
-  auto pivot = std::lower_bound(sum.begin(), sum.end(), (1 - fraction) * sum.back());
+  const double local_total = sum.size() > 0 ? sum.back() : 0.0;
+  auto pivot = std::lower_bound(sum.begin(), sum.end(), (1 - fraction) * local_total);
   auto index = std::distance(sum.begin(), pivot);
-  double error_threshold = estimates[index];
+  double error_threshold = estimates.size() > 0 ? estimates[index] : 0.0;
 
-  // Compute the number of elements, and error, marked by threshold value e.
-  auto marked = [&estimates, &sum](double e) -> std::pair<std::size_t, double>
+  // Compute the number of elements, and amount of error, marked by threshold value e.
+  auto marked = [&estimates, &sum, &local_total](double e) -> std::pair<std::size_t, double>
   {
-    const auto lb = std::lower_bound(estimates.begin(), estimates.end(), e);
-    const auto elems_marked = std::distance(lb, estimates.end());
-    const double error_unmarked =
-        lb != estimates.begin() ? sum[sum.size() - elems_marked - 1] : 0;
-    const double error_marked = sum.back() - error_unmarked;
-    return {elems_marked, error_marked};
+    if (local_total > 0)
+    {
+      const auto lb = std::lower_bound(estimates.begin(), estimates.end(), e);
+      const auto elems_marked = std::distance(lb, estimates.end());
+      const double error_unmarked =
+          lb != estimates.begin() ? sum[sum.size() - elems_marked - 1] : 0;
+      const double error_marked = local_total - error_unmarked;
+      return {elems_marked, error_marked};
+    }
+    else
+    {
+      return {0, 0.0};
+    }
   };
 
   // Each processor will compute a different threshold: if a given processor has lots of low
@@ -70,14 +78,14 @@ double ComputeDorflerThreshold(MPI_Comm comm, double fraction, const Vector &e)
     double min_marked;
     double max_marked;
   } error;
-  error.total = sum.back();
+  error.total = local_total;
   std::tie(elements.max_marked, error.max_marked) = marked(min_threshold);
   std::tie(elements.min_marked, error.min_marked) = marked(max_threshold);
   Mpi::GlobalSum(3, &elements.total, comm);
   Mpi::GlobalSum(3, &error.total, comm);
   const double max_indicator = [&]()
   {
-    double max_indicator = estimates.back();
+    double max_indicator = estimates.size() > 0 ? estimates.back() : 0.0;
     Mpi::GlobalMax(1, &max_indicator, comm);
     return max_indicator;
   }();
@@ -163,4 +171,32 @@ double ComputeDorflerThreshold(MPI_Comm comm, double fraction, const Vector &e)
   return error_threshold;
 }
 
+double ComputeDorflerCoarseningThreshold(const mfem::ParMesh &mesh, double fraction,
+                                         const Vector &e)
+{
+  MFEM_VERIFY(mesh.Nonconforming(), "Can only perform coarsening on a Nonconforming mesh!");
+  const auto &derefinement_table = mesh.pncmesh->GetDerefinementTable();
+  mfem::Array<double> elem_error(e.Size());
+  for (int i = 0; i < e.Size(); i++)
+  {
+    elem_error[i] = e[i];
+  }
+  mesh.pncmesh->SynchronizeDerefinementData(elem_error, derefinement_table);
+  Vector coarse_error(derefinement_table.Size());
+  mfem::Array<int> row;
+  for (int i = 0; i < derefinement_table.Size(); i++)
+  {
+    derefinement_table.GetRow(i, row);
+    coarse_error[i] = std::sqrt(std::accumulate(
+        row.begin(), row.end(), 0.0,
+        [&elem_error](double s, int i) { return s += std::pow(elem_error[i], 2.0); }));
+  }
+
+  // Given the coarse errors, we use the Dörfler marking strategy to identify the
+  // smallest set of original elements that make up (1 - θ) of the total error. The
+  // complement of this set is then the largest number of elements that make up θ of the
+  // total error.
+  return ComputeDorflerThreshold(mesh.GetComm(), 1 - fraction, coarse_error);
+}
+
 }  // namespace palace::utils
diff --git a/palace/utils/dorfler.hpp b/palace/utils/dorfler.hpp
@@ -7,6 +7,12 @@
 #include "linalg/vector.hpp"
 #include "utils/communication.hpp"
 
+namespace mfem
+{
+
+class ParMesh;
+
+}  // namespace mfem
 namespace palace::utils
 {
 
@@ -18,6 +24,13 @@ namespace palace::utils
 //            Numer. Anal. (1996).
 double ComputeDorflerThreshold(MPI_Comm comm, double fraction, const Vector &e);
 
+// Given a nonconforming mesh, target fraction and error estimates, compute a threshold
+// value that will mark the largest number of elements that make up the specified fraction
+// of the total coarsening opportunities. This is analogous to ComputeDorflerThreshold, but
+// operates only the list of available derefinement opportunities within the mesh.
+double ComputeDorflerCoarseningThreshold(const mfem::ParMesh &mesh, double fraction,
+                                         const Vector &e);
+
 }  // namespace palace::utils
 
 #endif  // PALACE_UTILS_DORFLER_HPP