support SPR with multiple processes

.github/workflows/cmake.yml

@@ -27,6 +27,13 @@ jobs:
         sudo apt-get install -yq cmake
         cmake --version
 
+    # There is a known bug using MPICH ubuntu-24.04 (ubuntu-latest as of 11/10/2025)
+    # https://bugs.launchpad.net/ubuntu/+source/mpich/+bug/2072338
+    - name: Install mpi
+      run: |
+        sudo apt-get update -yq
+        sudo apt-get install -yq openmpi-bin libopenmpi-dev
+
     ## Kokkos
     - name: Kokkos Checkout repo
       uses: actions/checkout@v4

CMakeLists.txt

@@ -9,6 +9,7 @@ option(MeshFields_USE_Cabana "Build with the Cabana storage backend" OFF)
 
 find_package(Kokkos REQUIRED)
 find_package(Omega_h REQUIRED)
+find_package(MPI REQUIRED)
 
 # Verify Omega_h was built with Kokkos enabled
 if(NOT Omega_h_USE_Kokkos)
@@ -234,6 +235,18 @@ add_test_property(OmegahSprAdapt
   PASS_REGULAR_EXPRESSION
   "afterAdapt nTri: 14439; solution_1 min -4084.78 max 213.322")
 
+test_func(OmegahSprAdapt_p2
+  ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 ${MPIEXEC_PREFLAGS}
+  ./OmegahSprAdapt ${CMAKE_SOURCE_DIR}/meshes/thwaites_basal_effectiveStrain.osh/
+  thwaitesAdapted_p2 #output prefix
+  0.1 # adapt ratio
+  1.0 # min size
+  16.0 # max size
+)
+add_test_property(OmegahSprAdapt_p2
+  PASS_REGULAR_EXPRESSION
+  "afterAdapt nTri: 14468; solution_1 min -4098.9 max 213.322")
+
 #Code Coverage set up -------------------------------------------------------
 
 option(meshfields_ENABLE_COVERAGE_BUILD "Do a coverage build" OFF)

src/MeshField_Integrate.hpp

@@ -196,9 +196,12 @@ auto getJacobianDeterminants(FieldElement &fes,
  * - `post`: Called after the integration process ends.
  * - `atPoints`: Called for each integration point to perform user-defined
  *   computations.
+ * - `parallelReduce`: Called after `atPoints` to reduce process-local
+ *   integrations into a global mesh integration.
  *
  * To use this class, derive from it and implement the `atPoints` method.
- * Optionally, override `pre` and `post` for additional setup or cleanup.
+ * Optionally, override `parallelReduce`, `pre` and `post` for distributed
+ * memory support, additional setup, or cleanup.
  */
 class Integrator {
 public:
@@ -227,9 +230,18 @@ class Integrator {
    */
   virtual void atPoints(Kokkos::View<Real **> p, Kokkos::View<Real *> w,
                         Kokkos::View<Real *> dV) = 0;
+
+  /** \brief User callback: distributed memory accumulation.
+   *
+   * \details If needed, this function supports the use of
+   * inter-process communication (e.g., MPI) to reduce
+   * process-local integrations into a global
+   * mesh integration.
+   */
+  virtual void parallelReduce(){};
+
   /** \brief Run the Integrator over the local field elements.
-   * \param fes FieldElement
-   * FIXME make the sensible
+   *  \param fes FieldElement
    * */
   template <typename FieldElement> void process(FieldElement &fes) {
     constexpr auto topo = decltype(FieldElement::elm2dof)::getTopology();
@@ -239,9 +251,8 @@ class Integrator {
     auto weights = getIntegrationPointWeights(fes, ip);
     auto dV = getJacobianDeterminants(fes, localCoords, ip.size());
     atPoints(localCoords, weights, dV);
+    parallelReduce();
     post();
-    // TODO support distributed meshes by running a parallel reduction with user
-    // functor
   }
 
 protected:

src/MeshField_SPR_ErrorEstimator.hpp

@@ -19,9 +19,12 @@ static double getNaN() { return std::numeric_limits<double>::quiet_NaN(); }
 /* common base for Scalar Integrator. */
 class SInt : public MeshField::Integrator {
 public:
-  SInt(int order) : MeshField::Integrator(order), r(0) {}
+  SInt(int order, Omega_h::CommPtr comm_in)
+      : MeshField::Integrator(order), r(0), comm(comm_in) {}
   void reset() { r = 0; }
   MeshField::Real r;
+  Omega_h::CommPtr comm;
+  void parallelReduce() { r = comm->allreduce(r, OMEGA_H_SUM); }
 };
 
 template <typename ShapeField> class Estimation {
@@ -83,8 +86,8 @@ template <typename EstimationT, typename OmegahMeshField>
 class SelfProduct : public SInt {
 public:
   SelfProduct(EstimationT &estimation_in, OmegahMeshField &omf_in)
-      : SInt(estimation_in.integration_order), estimation(estimation_in),
-        omf(omf_in) {}
+      : SInt(estimation_in.integration_order, estimation_in.mesh.comm()),
+        estimation(estimation_in), omf(omf_in) {}
   void atPoints(Kokkos::View<MeshField::Real **> p,
                 Kokkos::View<MeshField::Real *> w,
                 Kokkos::View<MeshField::Real *> dV) {
@@ -155,8 +158,9 @@ template <typename EstimationT, typename OmegahMeshField>
 class Error : public SInt {
 public:
   Error(EstimationT &estimation_in, OmegahMeshField &omf_in)
-      : SInt(estimation_in.integration_order), estimation(estimation_in),
-        omf(omf_in), errorNorm("errorNorm", estimation_in.mesh.nelems()) {}
+      : SInt(estimation_in.integration_order, estimation_in.mesh.comm()),
+        estimation(estimation_in), omf(omf_in),
+        errorNorm("errorNorm", estimation_in.mesh.nelems()) {}
   void atPoints(Kokkos::View<MeshField::Real **> p,
                 Kokkos::View<MeshField::Real *> w,
                 Kokkos::View<MeshField::Real *> dV) {
@@ -270,7 +274,8 @@ template <typename EstimationT, typename OmegahMeshField, typename FieldElement>
 std::tuple<Kokkos::View<MeshField::Real *>, MeshField::Real>
 getSprSizeField(EstimationT &e, OmegahMeshField &omf, FieldElement &coordFe) {
   Error errorIntegrator(e, omf);
-  errorIntegrator.process(coordFe);
+  errorIntegrator.process(
+      coordFe); // FIXME - needs to sync results across processes
   computeSizeFactor(e, omf, coordFe, errorIntegrator);
   getElementSizeField(e, errorIntegrator);
   return {averageToVertex(e.mesh, e.element_size), errorIntegrator.r};

test/testSprThwaitesAdapt.cpp

@@ -157,6 +157,7 @@ int main(int argc, char **argv) {
   auto world = lib.world();
   Omega_h::Mesh mesh(&lib);
   Omega_h::binary::read(argv[1], world, &mesh);
+
   const auto prefix = std::string(argv[2]);
   // adaptRatio = 0.1 is used in scorec/core:cws/sprThwaites test/spr_test.cc
   const MeshField::Real adaptRatio = std::stof(argv[3]);
@@ -166,13 +167,21 @@ int main(int argc, char **argv) {
       argv[5]); // a value of 8 or 16 roughly maintains the boundary shape
                 // in the downstream parts of the domain
                 // where there is significant coarsening
-  std::cout << "input mesh: " << argv[1] << " outputMeshPrefix: " << prefix
-            << " adaptRatio: " << adaptRatio << " max_size: " << max_size
-            << " min_size: " << min_size << "\n";
+  if (!mesh.comm()->rank()) {
+    std::cout << "input mesh: " << argv[1] << " outputMeshPrefix: " << prefix
+              << " adaptRatio: " << adaptRatio << " max_size: " << max_size
+              << " min_size: " << min_size << "\n";
+  }
   const auto outname = prefix + "_adaptRatio_" + std::string(argv[3]) +
                        "_maxSz_" + std::to_string(max_size) + "_minSz_" +
                        std::to_string(min_size);
 
+  // equally distribute elements among processes; required when loading a serial
+  // mesh
+  mesh.balance();
+  // Omega_h::project_by_fit used by spr requires ghosts
+  mesh.set_parting(Omega_h_Parting::OMEGA_H_GHOSTED);
+
   auto effectiveStrain = getEffectiveStrainRate(mesh);
   auto recoveredStrain = recoverLinearStrain(mesh, effectiveStrain);
   mesh.add_tag<Real>(VERT, "recoveredStrain", 1, recoveredStrain, false,
@@ -203,10 +212,10 @@ int main(int argc, char **argv) {
                      Omega_h::ArrayType::VectorND);
 
   { // write vtk
-    const std::string vtkFileName = "beforeAdapt" + outname + ".vtk";
+    const std::string vtkFileName = "beforeAdapt_" + outname + ".vtk";
     Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2);
     const std::string vtkFileName_edges =
-        "beforeAdapt" + outname + "_edges.vtk";
+        "beforeAdapt_" + outname + "_edges.vtk";
     Omega_h::vtk::write_parallel(vtkFileName_edges, &mesh, 1);
   }
 
@@ -223,12 +232,14 @@ int main(int argc, char **argv) {
   printTriCount(&mesh, "afterAdapt");
 
   const auto sol = mesh.get_array<Omega_h::Real>(VERT, "solution_1");
-  const auto sol_min = Omega_h::get_min(sol);
-  const auto sol_max = Omega_h::get_max(sol);
-  std::cout << "solution_1 min " << sol_min << " max " << sol_max << '\n';
+  const auto sol_min = Omega_h::get_min(mesh.comm(), sol);
+  const auto sol_max = Omega_h::get_max(mesh.comm(), sol);
+  if (!mesh.comm()->rank()) {
+    std::cout << "solution_1 min " << sol_min << " max " << sol_max << '\n';
+  }
 
   { // write vtk and osh for adapted mesh
-    const std::string outfilename = "afterAdapt" + outname;
+    const std::string outfilename = "afterAdapt_" + outname;
     Omega_h::vtk::write_parallel(outfilename + ".vtk", &mesh, 2);
     Omega_h::binary::write(outfilename + ".osh", &mesh);
     std::cout << "wrote adapted mesh: " << outfilename + ".osh"