Conservative transfer 3D

CMakeLists.txt

@@ -49,6 +49,23 @@ option(PCMS_ENABLE_Python "Enable pcms Python api" OFF)
 
 option(PCMS_ENABLE_PRINT "PCMS print statements enabled" ON)
 
+# Per-thread on-stack ring-buffer capacities for the adjacency-based
+# intersection/adjacency BFS (queue_visited.hpp). Defaults are sized for 3D
+# tetrahedral meshes on the serial/OpenMP backends. Lower them for device
+# (GPU) builds, where each search thread's stack must hold these arrays.
+set(PCMS_INTERSECTION_QUEUE_SIZE 1024 CACHE STRING
+    "Max BFS queue size per target element in the intersection search")
+set(PCMS_INTERSECTION_TRACK_SIZE 2048 CACHE STRING
+    "Max visited-set size per target element in the intersection search")
+set(PCMS_INTERSECTION_ABS_TOL 1e-18 CACHE STRING
+    "Absolute measure tolerance for accepting a clipped intersection")
+set(PCMS_INTERSECTION_REL_TOL 1e-12 CACHE STRING
+    "Relative measure tolerance for accepting a clipped intersection")
+set(PCMS_PROJECTION_KSP_RTOL 1e-14 CACHE STRING
+    "Default relative tolerance for the projection mass-matrix solve")
+set(PCMS_PROJECTION_KSP_ATOL 1e-14 CACHE STRING
+    "Default absolute tolerance for the projection mass-matrix solve")
+
 option(PETSC_LINK_STATIC "Use pkg-config --static results for PETSc"
        ${_pcms_link_petsc_static_default})
 

src/pcms/configuration.h.in

@@ -7,3 +7,11 @@
 #cmakedefine PCMS_ENABLE_Fortran
 #cmakedefine PCMS_ENABLE_MESHFIELDS
 #cmakedefine PCMS_ENABLE_PETSC
+
+#cmakedefine PCMS_INTERSECTION_QUEUE_SIZE @PCMS_INTERSECTION_QUEUE_SIZE@
+#cmakedefine PCMS_INTERSECTION_TRACK_SIZE @PCMS_INTERSECTION_TRACK_SIZE@
+
+#cmakedefine PCMS_INTERSECTION_ABS_TOL @PCMS_INTERSECTION_ABS_TOL@
+#cmakedefine PCMS_INTERSECTION_REL_TOL @PCMS_INTERSECTION_REL_TOL@
+#cmakedefine PCMS_PROJECTION_KSP_RTOL @PCMS_PROJECTION_KSP_RTOL@
+#cmakedefine PCMS_PROJECTION_KSP_ATOL @PCMS_PROJECTION_KSP_ATOL@

src/pcms/field/layout/omega_h_lagrange.cpp

@@ -171,6 +171,7 @@ OmegaHLagrangeLayout::OmegaHLagrangeLayout(
   owned_ = BuildOwned(mesh_, entity_dim, owned_mask);
   owned_host_ =
     Kokkos::View<bool*, HostMemorySpace>("owned_host", owned_.size());
+  Kokkos::deep_copy(owned_host_, owned_);
 
   class_ids_ = Omega_h::Read<Omega_h::ClassId>(
     mesh_.get_array<Omega_h::ClassId>(entity_dim, "class_id"));

src/pcms/localization/queue_visited.hpp

@@ -7,16 +7,15 @@
 #include <Omega_h_library.hpp>
 #include <Omega_h_mesh.hpp>
 #include <Omega_h_reduce.hpp>
-#define MAX_SIZE_QUEUE 500
-#define MAX_SIZE_TRACK 800
+#include "pcms/configuration.h"
 
 namespace pcms
 {
 
 class Queue
 {
 private:
-  Omega_h::LO queue_array[MAX_SIZE_QUEUE];
+  Omega_h::LO queue_array[PCMS_INTERSECTION_QUEUE_SIZE];
   int first = 0, last = -1, count = 0;
 
 public:
@@ -45,7 +44,7 @@ class Queue
 class Track
 {
 private:
-  Omega_h::LO tracking_array[MAX_SIZE_TRACK];
+  Omega_h::LO tracking_array[PCMS_INTERSECTION_TRACK_SIZE];
   int first = 0, last = -1, count = 0;
 
 public:
@@ -68,11 +67,11 @@ class Track
 OMEGA_H_INLINE
 void Queue::push_back(const int& item)
 {
-  if (count == MAX_SIZE_QUEUE) {
+  if (count == PCMS_INTERSECTION_QUEUE_SIZE) {
     printf("queue is full %d\n", count);
     return;
   }
-  last = (last + 1) % MAX_SIZE_QUEUE;
+  last = (last + 1) % PCMS_INTERSECTION_QUEUE_SIZE;
   queue_array[last] = item;
   count++;
 }
@@ -84,7 +83,7 @@ void Queue::pop_front()
     printf("queue is empty\n");
     return;
   }
-  first = (first + 1) % MAX_SIZE_QUEUE;
+  first = (first + 1) % PCMS_INTERSECTION_QUEUE_SIZE;
   count--;
 }
 
@@ -103,19 +102,19 @@ bool Queue::isEmpty() const
 OMEGA_H_INLINE
 bool Queue::isFull() const
 {
-  return count == MAX_SIZE_QUEUE;
+  return count == PCMS_INTERSECTION_QUEUE_SIZE;
 }
 
 OMEGA_H_INLINE
 bool Track::push_back(const int& item)
 {
-  if (count == MAX_SIZE_TRACK) {
+  if (count == PCMS_INTERSECTION_TRACK_SIZE) {
     // Visited buffer is full. Report failure so callers can stop expanding the
     // search; otherwise new vertices can never be marked visited and the BFS
     // re-queues them forever (infinite loop).
     return false;
   }
-  last = (last + 1) % MAX_SIZE_TRACK;
+  last = (last + 1) % PCMS_INTERSECTION_TRACK_SIZE;
   tracking_array[last] = item;
   count++;
   return true;
@@ -126,7 +125,7 @@ bool Track::notVisited(const int& item)
 {
   int id;
   for (int i = 0; i < count; ++i) {
-    id = (first + i) % MAX_SIZE_TRACK;
+    id = (first + i) % PCMS_INTERSECTION_TRACK_SIZE;
     if (tracking_array[id] == item) {
       return false;
     }

src/pcms/transfer/conservative_projection_solver.cpp

@@ -1,5 +1,6 @@
 #include <petscksp.h>
 
+#include "pcms/configuration.h"
 #include "pcms/transfer/conservative_projection_solver.hpp"
 #include "pcms/utility/arrays.h"
 #include "pcms/transfer/petsc_utils.hpp"
@@ -65,6 +66,10 @@ GalerkinProjectionSolver::GalerkinProjectionSolver(
     ierr = PCSetType(pc, PCJACOBI);
     CHKERRABORT(PETSC_COMM_SELF, ierr);
   }
+  ierr =
+    KSPSetTolerances(ksp_, PCMS_PROJECTION_KSP_RTOL, PCMS_PROJECTION_KSP_ATOL,
+                     PETSC_DEFAULT, PETSC_DEFAULT);
+  CHKERRABORT(PETSC_COMM_SELF, ierr);
   ierr = KSPSetFromOptions(ksp_);
   CHKERRABORT(PETSC_COMM_SELF, ierr);
   ierr = KSPSetUp(ksp_);

src/pcms/transfer/mass_matrix_integrator.hpp

@@ -18,16 +18,19 @@ template <typename FieldElement>
 class MassMatrixIntegrator : public MeshField::Integrator
 {
 public:
+  // Linear simplex: numNodes = spatial dim + 1 (3 for triangles, 4 for tets).
+  static constexpr int numNodes = FieldElement::MeshEntDim + 1;
+
   MassMatrixIntegrator(Omega_h::Mesh& mesh_in, FieldElement& fe_in,
                        int order = 2)
     : mesh(mesh_in),
       fe(fe_in),
-      subMatrixSize(3 * 3), // FIXME remove hard coded size
-      elmMassMatrix("elmMassMatrix", mesh_in.nelems() * 3 * 3),
+      subMatrixSize(numNodes * numNodes),
+      elmMassMatrix("elmMassMatrix", mesh_in.nelems() * numNodes * numNodes),
       Integrator(order)
   {
     Kokkos::deep_copy(elmMassMatrix, 0);
-    assert(mesh.dim() == 2); // TODO support 1d,2d,3d
+    assert(mesh.dim() == 2 || mesh.dim() == 3);
     assert(mesh.family() == OMEGA_H_SIMPLEX);
   }
   void atPoints(Kokkos::View<MeshField::Real**> p,
@@ -58,12 +61,18 @@ class MassMatrixIntegrator : public MeshField::Integrator
           }
           const auto N = shapeFn.getValues(localCoord);
           const auto wPt = w(pt);
-          const auto dVPt = dV(pt);
+          // Use the unsigned volume element: MeshField returns a signed
+          // Jacobian determinant, which is negative for tetrahedra whose vertex
+          // ordering has negative orientation. A mass matrix integrates against
+          // the positive volume measure, so take the magnitude (a no-op in 2D
+          // where the differential area is already positive).
+          const auto dVPt = Kokkos::fabs(dV(pt));
           //       printf("Shape Functions: %f, %f, %f \n", N[0], N[1], N[2]);
           //       printf("wPt, dVPt: %f, %f \n", wPt, dVPt);
           for (auto i = 0; i < N.size(); i++) {
             for (auto j = 0; j < N.size(); j++) {
-              massMatrix(elm * subMat + i * 3 + j) += N[i] * N[j] * wPt * dVPt;
+              massMatrix(elm * subMat + i * numNodes + j) +=
+                N[i] * N[j] * wPt * dVPt;
             }
           }
         }

src/pcms/transfer/mesh_intersection.cpp

@@ -4,51 +4,66 @@
 
 namespace pcms
 {
+namespace
+{
+// Construct the source-mesh containing-element search appropriate to the
+// spatial dimension: a uniform 20^Dim background grid over the source mesh.
+template <int Dim>
+auto MakeGridPointSearch(Omega_h::Mesh& source_mesh)
+{
+  if constexpr (Dim == 3) {
+    return pcms::GridPointSearch3D(source_mesh, 20, 20, 20);
+  } else {
+    return pcms::GridPointSearch2D(source_mesh, 20, 20);
+  }
+}
+} // namespace
+
+template <int Dim>
 void FindIntersections::adjBasedIntersectSearch(
   const Omega_h::LOs& tgt2src_offsets,
   Omega_h::Write<Omega_h::LO>& nIntersections,
   Omega_h::Write<Omega_h::LO>& tgt2src_indices, bool is_count_only)
 {
-
+  // Element entity dimension equals the spatial dimension (FACE for 2D, REGION
+  // for 3D); measures are triangle areas (2D) or tet volumes (3D).
   const auto& tgt_coords = target_mesh_.coords();
   const auto& src_coords = source_mesh_.coords();
-  const auto& tgt_faces2nodes =
-    target_mesh_.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b;
-  const auto& src_faces2nodes =
-    source_mesh_.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b;
-  const auto& src_elem_areas = measure_elements_real(&source_mesh_);
-  const auto& tgt_elem_areas = measure_elements_real(&target_mesh_);
+  const auto& tgt_elems2nodes = target_mesh_.ask_down(Dim, Omega_h::VERT).ab2b;
+  const auto& src_elems2nodes = source_mesh_.ask_down(Dim, Omega_h::VERT).ab2b;
+  const auto& src_elem_measures = measure_elements_real(&source_mesh_);
+  const auto& tgt_elem_measures = measure_elements_real(&target_mesh_);
   const auto& t2t =
     source_mesh_.ask_dual(); // gives connected element neighbors
   const auto& t2tt = t2t.a2ab;
   const auto& tt2t = t2t.ab2b;
 
-  const auto flat_centroids =
-    pcms::get_entity_centroids(target_mesh_, Omega_h::FACE);
+  const auto flat_centroids = pcms::get_entity_centroids(target_mesh_, Dim);
   // Convert layout_right 1D Omega_h array to 2D Kokkos view with correct layout
-  auto centroids = ConvertCoordsTo2D(flat_centroids, target_mesh_.nfaces(), 2);
+  auto centroids =
+    ConvertCoordsTo2D(flat_centroids, target_mesh_.nelems(), Dim);
 
-  pcms::GridPointSearch2D search_cell(source_mesh_, 20, 20);
+  auto search_cell = MakeGridPointSearch<Dim>(source_mesh_);
   auto results = search_cell(centroids);
   auto owning_cell_ids = search_cell.GetOwningElementIds(results);
 
-  auto nfaces_target = target_mesh_.nfaces();
+  auto nelems_target = target_mesh_.nelems();
   Omega_h::parallel_for(
-    nfaces_target,
+    nelems_target,
     OMEGA_H_LAMBDA(const Omega_h::LO id) {
       Queue queue;
       Track visited;
 
       auto current_cell_id = owning_cell_ids(id);
-      auto current_tgt_elm_area = tgt_elem_areas[id];
+      auto current_tgt_elm_measure = tgt_elem_measures[id];
 
       OMEGA_H_CHECK_PRINTF(current_cell_id >= 0,
                            "ERROR: source cell id not found for given target "
-                           "centroid %d (%f, %f)\n",
-                           id, centroids(id, 0), centroids(id, 1));
+                           "centroid %d\n",
+                           id);
 
       auto tgt_elm_vert_coords =
-        get_vert_coords_of_elem(tgt_coords, tgt_faces2nodes, id);
+        get_vert_coords_of_elem<Dim>(tgt_coords, tgt_elems2nodes, id);
 
       Omega_h::LO start_counter;
       if (!is_count_only) {
@@ -76,22 +91,45 @@ void FindIntersections::adjBasedIntersectSearch(
           auto neighborElmId = tt2t[i];
 
           if (visited.notVisited(neighborElmId)) {
-            visited.push_back(neighborElmId);
-            auto elm_vert_coords = get_vert_coords_of_elem(
-              src_coords, src_faces2nodes, neighborElmId);
-            r3d::Polytope<2> intersection;
+            // If the visited buffer is full, skip this neighbor so the BFS
+            // terminates. Without this, an unrecorded neighbor stays "not
+            // visited" and is re-queued forever (infinite loop). Mirrors the
+            // guard in adj_search.cpp. Raise PCMS_INTERSECTION_TRACK_SIZE if
+            // this fires.
+            if (!visited.push_back(neighborElmId)) {
+              printf("ERROR: visited buffer full "
+                     "(PCMS_INTERSECTION_TRACK_SIZE=%d) for "
+                     "target %d; some intersections may be missed\n",
+                     PCMS_INTERSECTION_TRACK_SIZE, id);
+              continue;
+            }
+            auto elm_vert_coords = get_vert_coords_of_elem<Dim>(
+              src_coords, src_elems2nodes, neighborElmId);
+            r3d::Polytope<Dim> intersection;
             r3d::intersect_simplices(intersection, tgt_elm_vert_coords,
                                      elm_vert_coords);
-            auto intersected_area = r3d::measure(intersection);
-            auto current_src_elm_area = src_elem_areas[neighborElmId];
+            // Take the magnitude: r3d::measure is signed by the orientation of
+            // the target simplex used to initialize the polytope, which can be
+            // negative for tetrahedra. This mirrors the fabs applied in the
+            // sub-simplex decomposition and mass assembly; without it a
+            // negatively-oriented target element would reject all of its real
+            // overlaps and break conservation.
+            auto intersected_measure = Kokkos::fabs(r3d::measure(intersection));
+            auto current_src_elm_measure = src_elem_measures[neighborElmId];
             auto scale =
-              Kokkos::fmax(current_tgt_elm_area, current_src_elm_area);
-            auto eps = Kokkos::fmax(abs_tol, rel_tol * scale);
-            if (intersection.nverts >= 3 && intersected_area >= eps) {
+              Kokkos::fmax(current_tgt_elm_measure, current_src_elm_measure);
+            auto eps = Kokkos::fmax(PCMS_INTERSECTION_ABS_TOL,
+                                    PCMS_INTERSECTION_REL_TOL * scale);
+            // A valid intersection is a non-degenerate simplex-simplex overlap:
+            // at least Dim+1 vertices (a polygon in 2D, a polyhedron in 3D).
+            if (intersection.nverts >= Dim + 1 && intersected_measure >= eps) {
               count++;
 
               OMEGA_H_CHECK_PRINTF(
-                count < 500, "WARNING: count exceeds 500 for target %d", id);
+                count < PCMS_INTERSECTION_QUEUE_SIZE,
+                "intersection count for target %d reached the cap %d; raise "
+                "PCMS_INTERSECTION_QUEUE_SIZE",
+                id, PCMS_INTERSECTION_QUEUE_SIZE);
 
               queue.push_back(neighborElmId);
 
@@ -113,20 +151,32 @@ void FindIntersections::adjBasedIntersectSearch(
     }, // end of lambda
     "count the number of intersections for each target element");
 }
-IntersectionResults intersectTargets(Omega_h::Mesh& source_mesh,
-                                     Omega_h::Mesh& target_mesh)
+
+// Explicit instantiations for the supported spatial dimensions.
+template void FindIntersections::adjBasedIntersectSearch<2>(
+  const Omega_h::LOs&, Omega_h::Write<Omega_h::LO>&,
+  Omega_h::Write<Omega_h::LO>&, bool);
+template void FindIntersections::adjBasedIntersectSearch<3>(
+  const Omega_h::LOs&, Omega_h::Write<Omega_h::LO>&,
+  Omega_h::Write<Omega_h::LO>&, bool);
+
+namespace
+{
+template <int Dim>
+IntersectionResults intersectTargetsImpl(Omega_h::Mesh& source_mesh,
+                                         Omega_h::Mesh& target_mesh)
 {
   FindIntersections intersect(source_mesh, target_mesh);
 
-  auto nfaces_target = target_mesh.nfaces();
+  auto nelems_target = target_mesh.nelems();
 
   Omega_h::Write<Omega_h::LO> nIntersections(
-    nfaces_target, 0, "number of intersections in each target vertex");
+    nelems_target, 0, "number of intersections in each target element");
 
   Omega_h::Write<Omega_h::LO> tgt2src_indices;
 
-  intersect.adjBasedIntersectSearch(Omega_h::LOs(), nIntersections,
-                                    tgt2src_indices, true);
+  intersect.adjBasedIntersectSearch<Dim>(Omega_h::LOs(), nIntersections,
+                                         tgt2src_indices, true);
 
   Kokkos::fence();
   auto tgt2src_offsets = Omega_h::offset_scan(Omega_h::Read(nIntersections),
@@ -139,9 +189,20 @@ IntersectionResults intersectTargets(Omega_h::Mesh& source_mesh,
     ntotal_intersections, 0,
     "indices of the source elements that intersect the given target element");
 
-  intersect.adjBasedIntersectSearch(tgt2src_offsets, nIntersections,
-                                    tgt2src_indices, false);
+  intersect.adjBasedIntersectSearch<Dim>(tgt2src_offsets, nIntersections,
+                                         tgt2src_indices, false);
   return {.tgt2src_offsets = tgt2src_offsets,
           .tgt2src_indices = Omega_h::read(tgt2src_indices)};
 }
+} // namespace
+
+IntersectionResults intersectTargets(Omega_h::Mesh& source_mesh,
+                                     Omega_h::Mesh& target_mesh)
+{
+  OMEGA_H_CHECK(source_mesh.dim() == target_mesh.dim());
+  if (source_mesh.dim() == 3) {
+    return intersectTargetsImpl<3>(source_mesh, target_mesh);
+  }
+  return intersectTargetsImpl<2>(source_mesh, target_mesh);
+}
 } // namespace pcms