Make the field API handle owned data by default

src/pcms/capi/client.cpp

@@ -78,7 +78,7 @@ class EmptyFieldFactory : public FieldFactory
           throw pcms_error(
             "EmptyFieldFactory::CreateField: requires SimpleFieldData");
         }
-        if (fd->GetDOFAccess()->GetOwnedValues().size() !=
+        if (fd->GetDOFAccess()->GetValues().size() !=
             detail::ExpectedFlatFieldDataSize(*layout_)) {
           throw pcms_error(
             "EmptyFieldFactory::CreateField: field data size does not match "

src/pcms/coupler/field_communicator.hpp

@@ -54,6 +54,11 @@ class FieldCommunicator
     auto data = comm_.Recv(redev::Mode::Synchronous);
     serializer_.Deserialize(field_, make_const_array_view(data),
                             layout_comm_.GetPermutationArray());
+    // Deserialize writes owned DOF holders only, so ghost copies still hold
+    // their previous values. Anything that evaluates the field reads ghosts
+    // (an owned element can have vertices another rank owns), so bring them up
+    // to date.
+    field_.Synchronize();
   }
 
 private:

src/pcms/coupler/field_exchange_planner.cpp

@@ -29,20 +29,6 @@ struct OutMsg
   redev::LOs offset;
 };
 
-// Returns the mesh entity dimension for DOF local_index based on the
-// entity offsets array.  ent_offsets[d]..ent_offsets[d+1] is the range of
-// DOF indices belonging to mesh entity dimension d.
-static int GetMeshEntityDim(LO local_index, const EntOffsetsArray& ent_offsets)
-{
-  for (int d = 0; d < ent_offsets_len - 1; ++d) {
-    if (local_index >= static_cast<LO>(ent_offsets[d]) &&
-        local_index < static_cast<LO>(ent_offsets[d + 1])) {
-      return d;
-    }
-  }
-  return ent_offsets_len - 1;
-}
-
 static size_t GetMessageBlockIndex(
   LO permutation_entry, Rank1View<const redev::LO, HostMemorySpace> offsets)
 {
@@ -189,6 +175,7 @@ static ReversePartitionMap2 BuildReversePartitionMap(
 {
   PCMS_FUNCTION_TIMER;
   auto owned = layout.GetOwnedHost();
+  auto owned_indices = layout.GetOwnedIndicesHost();
   auto class_dims = layout.GetDOFHolderClassificationDimensionsHost();
   auto class_ids = layout.GetDOFHolderClassificationIdsHost();
   auto coords = layout.GetDOFHolderCoordinates().GetValues();
@@ -209,28 +196,33 @@ static ReversePartitionMap2 BuildReversePartitionMap(
     Kokkos::create_mirror_view_and_copy(HostMemorySpace(), coords_device);
 
   ReversePartitionMap2 reverse_partition;
-  LO n = static_cast<LO>(owned.extent(0));
+  const LO n_owned = static_cast<LO>(coords.extent(0));
   std::array<Real, 3> coord{};
   auto overlap_mask_view = overlap_mask.GetMask(layout);
 
-  for (LO local_index = 0; local_index < n; ++local_index) {
+  for (LO owned_index = 0; owned_index < n_owned; ++owned_index) {
+    // Skip holders not owned by this rank. This is required for xgc, whose
+    // layout marks holders as not owned without having ghosts. The ownership
+    // mask is indexed with ghosts.
+    const LO local_index =
+      owned_indices.size() == 0 ? owned_index : owned_indices(owned_index);
     if (!owned(local_index))
       continue;
 
-    if (!overlap_mask_view[local_index])
+    if (!overlap_mask_view[owned_index])
       continue;
 
     for (int d = 0; d < mesh_dim; ++d)
-      coord[d] = coords_host(local_index, d);
+      coord[d] = coords_host(owned_index, d);
     for (int d = mesh_dim; d < 3; ++d)
       coord[d] = 0.0;
 
-    int mesh_ent_dim = GetMeshEntityDim(local_index, ent_offsets);
-    LO class_dim = class_dims[local_index];
-    LO class_id = class_ids[local_index];
+    int mesh_ent_dim = GetMeshEntityDim(owned_index, ent_offsets);
+    LO class_dim = class_dims[owned_index];
+    LO class_id = class_ids[owned_index];
 
     auto dr = std::visit(GetRank{class_id, class_dim, coord}, partition);
-    reverse_partition[dr].indices.emplace_back(local_index);
+    reverse_partition[dr].indices.emplace_back(owned_index);
 
     for (size_t e = static_cast<size_t>(mesh_ent_dim) + 1; e < ent_offsets_len;
          ++e) {
@@ -296,28 +288,25 @@ void GenericFieldExchangePlanner::FillGidMessage(
                      plan.msg_size + header_size);
 
   auto gids = layout.GetGidsHost();
-  auto owned = layout.GetOwnedHost();
   auto ent_offsets = layout.GetEntOffsets();
   auto offsets = Rank1View<const redev::LO, HostMemorySpace>(
     plan.offsets.data(), plan.offsets.size());
 
   std::vector<EntOffsetsArray> per_rank_offsets(plan.dest_ranks.size());
 
-  for (LO local_index = 0; local_index < static_cast<LO>(gids.size());
-       ++local_index) {
-    LO perm_index = plan.permutation[local_index];
+  for (LO owned_index = 0; owned_index < static_cast<LO>(gids.size());
+       ++owned_index) {
+    LO perm_index = plan.permutation[owned_index];
     // Owned holders outside the overlap region carry the sentinel and have no
     // slot in the message.
     if (perm_index < 0)
       continue;
-    // A holder with a valid permutation slot must be owned.
-    PCMS_ALWAYS_ASSERT(owned[local_index]);
     auto block_index = GetMessageBlockIndex(perm_index, offsets);
     const auto gid_index =
       perm_index + static_cast<LO>((block_index + 1) * ent_offsets_len);
-    gid_message(gid_index) = gids(local_index);
+    gid_message(gid_index) = gids(owned_index);
 
-    int mesh_ent_dim = GetMeshEntityDim(local_index, ent_offsets);
+    int mesh_ent_dim = GetMeshEntityDim(owned_index, ent_offsets);
     for (size_t e = static_cast<size_t>(mesh_ent_dim) + 1; e < ent_offsets_len;
          ++e) {
       per_rank_offsets[block_index][e] += 1;

src/pcms/coupler/field_layout_communicator.cpp

@@ -24,9 +24,9 @@ FieldLayoutCommunicator::FieldLayoutCommunicator(
     name_(name),
     redev_(redev),
     planner_(std::move(planner)),
-    overlap_mask_(overlap_mask ? std::make_unique<OverlapMask>(*overlap_mask)
-                               : std::make_unique<OverlapMask>(
-                                   layout.GetGidsHost().size())),
+    overlap_mask_(overlap_mask
+                    ? std::make_unique<OverlapMask>(*overlap_mask)
+                    : std::make_unique<OverlapMask>(layout.GetNumDofHolder())),
     own_mpi_comm_(own_mpi_comm)
 {
   gid_comm_ = channel.CreateComm<GO>(name_ + "_gids", mpi_comm_);
@@ -111,8 +111,8 @@ void FieldLayoutCommunicator::UpdateLayoutNull()
   PCMS_FUNCTION_TIMER;
   // This rank exchanges nothing, but it still unpacks on receive (backends
   // that keep replicas synchronize them there), so every holder is unmapped.
-  plan_.permutation = MakeUnmappedPermutation(
-    static_cast<size_t>(layout_.GetNumOwnedDofHolder()));
+  plan_.permutation =
+    MakeUnmappedPermutation(static_cast<size_t>(layout_.GetNumDofHolder()));
   if (redev_.GetProcessType() == redev::ProcessType::Client) {
     channel_.BeginSendCommunicationPhase();
     channel_.EndSendCommunicationPhase();

src/pcms/coupler/field_serializer.h

@@ -13,14 +13,15 @@ namespace pcms
 // Packs a field's owned values into, and unpacks them from, an exchange
 // buffer, through the field's DOFAccess. Any field with DOFAccess is
 // exchanged this way. Backends that need extra work on receive (e.g.,
-// synchronizing replicas) do it in DOFAccess::SetOwnedValues.
+// synchronizing replicas) do it in DOFAccess::SetValues.
 template <typename T>
 class FieldSerializer
 {
 public:
   int Serialize(const Field<T>& field, Rank1View<T, HostMemorySpace> buffer,
                 Rank1View<const LO, HostMemorySpace> permutation) const
   {
+    // Only owned (rank-exclusive) DOF holders are serialized.
     const auto data = field.GetDOFHolderDataHost().GetValues();
     if (buffer.size() > 0) {
       const LO num_dof = static_cast<LO>(data.extent(0));
@@ -40,15 +41,14 @@ class FieldSerializer
                    Rank1View<const LO, HostMemorySpace> permutation) const
   {
     const auto& layout = field.GetLayout();
-    const LO num_dof = layout.GetNumOwnedDofHolder();
+    const LO num_dof = layout.GetNumDofHolder();
     const LO num_comp = layout.GetNumComponents();
-    // We initialize with local data so that
-    // data that doesn't get received from deserialized
-    // doesn't get zeroed out. Long term FIXME
+    // We initialize with the current data so that data that doesn't get
+    // received from deserialize doesn't get zeroed out. Long term FIXME
     // is to make Get/SetData operate over only part of the
     // arrays
     const auto current = field.GetDOFHolderDataHost().GetValues();
-    Kokkos::View<T*, HostMemorySpace> sorted("sorted", layout.OwnedSize());
+    Kokkos::View<T*, HostMemorySpace> sorted("sorted", layout.Size());
     for (LO i = 0; i < num_dof; ++i) {
       for (LO c = 0; c < num_comp; ++c) {
         sorted[i * num_comp + c] = current(i, c);

src/pcms/coupler/overlap_mask.h

@@ -15,6 +15,9 @@ struct OverlapMask
   Kokkos::View<bool*, HostMemorySpace> is_overlap_;
   std::function<int8_t(int, int)> in_overlap_func_;
 
+  // Masks are indexed by owned DOF holder, in owned order, matching the
+  // owned-indexed exchange plan.
+
   // Default constructor: all DOF holders are in overlap
   OverlapMask(size_t size) : is_overlap_("overlap_info", size)
   {
@@ -30,15 +33,15 @@ struct OverlapMask
   }
 
   // Construct from a precomputed per-DOF-holder host mask (e.g. an MFEM
-  // subdomain selected by element attribute). Indexed by local DOF holder.
+  // subdomain selected by element attribute), one entry per owned DOF holder.
   OverlapMask(size_t size, Kokkos::View<bool*, HostMemorySpace> is_overlap_host)
     : is_overlap_("overlap_info", size)
   {
     PCMS_ALWAYS_ASSERT(is_overlap_host.extent(0) == size);
     Kokkos::deep_copy(is_overlap_, is_overlap_host);
   }
 
-  // Construct from Omega_h host array
+  // Construct from Omega_h host array, one entry per owned DOF holder
   OverlapMask(size_t size, Omega_h::HostRead<Omega_h::I8> is_overlap_host)
     : is_overlap_("overlap_info", size)
   {
@@ -47,7 +50,7 @@ struct OverlapMask
     }
   }
 
-  // Construct from Omega_h device array
+  // Construct from Omega_h device array, one entry per owned DOF holder
   OverlapMask(size_t size, Omega_h::Read<Omega_h::I8> is_overlap_device)
     : is_overlap_("overlap_info", size)
   {
@@ -57,11 +60,14 @@ struct OverlapMask
     }
   }
 
-  // Get the mask, evaluating the function against the layout if one was given
+  // Get the mask, one entry per owned DOF holder, evaluating the function
+  // against the layout's classification if one was given.
   Kokkos::View<const bool*, HostMemorySpace> GetMask(
     const FieldLayout& layout) const
   {
     if (!in_overlap_func_) {
+      PCMS_ALWAYS_ASSERT(is_overlap_.extent(0) ==
+                         static_cast<size_t>(layout.GetNumDofHolder()));
       return is_overlap_;
     }
     auto class_dims = layout.GetDOFHolderClassificationDimensionsHost();

src/pcms/field/CMakeLists.txt

@@ -53,6 +53,7 @@ if (PCMS_ENABLE_OMEGA_H)
             layout/omega_h_entity.h
             layout/omega_h_lagrange.h
             evaluator/omega_h_lagrange.h
+            data/omega_h.h
             layout/uniform_grid.h
             uniform_grid_binary_field.h
             evaluator/uniform_grid.h