Add mfem adapter

CMakeLists.txt

@@ -44,6 +44,7 @@ option(PCMS_ENABLE_CLIENT "enable the coupling client implementation" ON)
 
 option(PCMS_ENABLE_XGC "enable xgc field adapter" ON)
 option(PCMS_ENABLE_OMEGA_H "enable Omega_h field adapter" OFF)
+option(PCMS_ENABLE_MFEM "enable MFEM field adapter" OFF)
 option(PCMS_ENABLE_C "Enable pcms C api" ON)
 option(PCMS_ENABLE_Python "Enable pcms Python api" OFF)
 
@@ -136,6 +137,14 @@ if(PCMS_ENABLE_OMEGA_H)
     message(FATAL_ERROR "Omega_h must be built with MPI enabled.")
   endif()
 endif()
+
+if(PCMS_ENABLE_MFEM)
+  find_package(MFEM REQUIRED)
+  message(STATUS "Found MFEM: ${MFEM_DIR} (found version ${MFEM_VERSION})")
+  if(NOT MFEM_USE_MPI)
+    message(FATAL_ERROR "MFEM must be built with MPI enabled.")
+  endif()
+endif()
 # adios2 adds C and Fortran depending on how it was built
 find_package(ADIOS2 CONFIG 2.10.2 REQUIRED)
 find_package(Kokkos CONFIG 4.5 REQUIRED)

config.cmake.in

@@ -29,6 +29,10 @@ if(@PCMS_ENABLE_OMEGA_H@)
     find_dependency(Omega_h CONFIG HINTS @Omega_h_DIR@)
 endif()
 
+if(@PCMS_ENABLE_MFEM@)
+    find_dependency(MFEM CONFIG HINTS @MFEM_DIR@)
+endif()
+
 include("${CMAKE_CURRENT_LIST_DIR}/pcms_utility-targets.cmake")
 include("${CMAKE_CURRENT_LIST_DIR}/pcms_discretization-targets.cmake")
 include("${CMAKE_CURRENT_LIST_DIR}/pcms_localization-targets.cmake")

src/pcms/configuration.h.in

@@ -2,6 +2,7 @@
 #cmakedefine PCMS_ENABLE_CLIENT
 #cmakedefine PCMS_ENABLE_XGC
 #cmakedefine PCMS_ENABLE_OMEGA_H
+#cmakedefine PCMS_ENABLE_MFEM
 #cmakedefine PCMS_ENABLE_C
 #cmakedefine PCMS_ENABLE_PRINT
 #cmakedefine PCMS_ENABLE_Fortran

src/pcms/coupler/field_communicator.hpp

@@ -51,8 +51,7 @@ class FieldCommunicator
     PCMS_FUNCTION_TIMER;
     PCMS_ALWAYS_ASSERT(layout_comm_.GetChannel().InSendCommunicationPhase());
     auto buffer = make_array_view(comm_buffer_);
-    serializer_->Serialize(field_.GetData(), field_.GetLayout(), buffer,
-                           layout_comm_.GetPermutationArray());
+    serializer_->Serialize(field_, buffer, layout_comm_.GetPermutationArray());
     comm_.Send(buffer.data_handle(), mode);
   }
 
@@ -64,8 +63,7 @@ class FieldCommunicator
     // mode because we make an immediate call to deserialize after a call to
     // receive.
     auto data = comm_.Recv(redev::Mode::Synchronous);
-    serializer_->Deserialize(field_.GetData(), field_.GetLayout(),
-                             make_const_array_view(data),
+    serializer_->Deserialize(field_, make_const_array_view(data),
                              layout_comm_.GetPermutationArray());
   }
 

src/pcms/coupler/field_serializer.h

@@ -14,22 +14,15 @@ template <typename T>
 class FieldSerializer
 {
 public:
-  virtual int Serialize(const FieldData<T>& field, const FieldLayout& layout,
+  virtual int Serialize(const Field<T>& field,
                         Rank1View<T, HostMemorySpace> buffer,
                         Rank1View<const LO, HostMemorySpace> permutation) const
   {
     auto data = field.GetDOFHolderDataHost();
-    auto owned = layout.GetOwnedHost();
-    // The exchange plan is per DOF holder: owned[i] and permutation[i] are
-    // indexed by holder. All num_components components of a holder share its
-    // location, so they occupy one contiguous block permutation[i]*num_comp in
-    // the wire buffer.
     if (buffer.size() > 0) {
       const LO num_dof = static_cast<LO>(data.extent(0));
       const LO num_comp = static_cast<LO>(data.extent(1));
       for (LO i = 0; i < num_dof; ++i) {
-        // A negative permutation entry marks a holder outside the exchange
-        // (non-owned, or owned but outside the overlap region); it has no slot.
         if (permutation[i] >= 0) {
           for (LO c = 0; c < num_comp; ++c) {
             buffer[permutation[i] * num_comp + c] = data(i, c);
@@ -41,17 +34,25 @@ class FieldSerializer
   }
 
   virtual void Deserialize(
-    FieldData<T>& field, const FieldLayout& layout,
-    Rank1View<const T, HostMemorySpace> buffer,
+    Field<T>& field, Rank1View<const T, HostMemorySpace> buffer,
     Rank1View<const LO, HostMemorySpace> permutation) const
   {
+    const auto& layout = field.GetLayout();
     const LO num_dof = layout.GetNumOwnedDofHolder();
     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
+    // is to make Get/SetData operate over only part of the
+    // arrays
+    auto current = field.GetDOFHolderDataHost();
     Kokkos::View<T*, HostMemorySpace> sorted("sorted", layout.OwnedSize());
     for (LO i = 0; i < num_dof; ++i) {
-      // A negative permutation entry marks a holder outside the exchange (owned
-      // but outside the overlap region); no data was received for it, so its
-      // zero-initialized `sorted` slot is left as-is.
+      for (LO c = 0; c < num_comp; ++c) {
+        sorted[i * num_comp + c] = current(i, c);
+      }
+    }
+    for (LO i = 0; i < num_dof; ++i) {
       if (permutation[i] >= 0) {
         for (LO c = 0; c < num_comp; ++c) {
           sorted[i * num_comp + c] = buffer[permutation[i] * num_comp + c];

src/pcms/coupler/overlap_mask.h

@@ -3,6 +3,7 @@
 
 #include "pcms/field/field_layout.h"
 #include "pcms/utility/arrays.h"
+#include "pcms/utility/assert.h"
 #include <Omega_h_array.hpp>
 #include <functional>
 
@@ -27,6 +28,15 @@ struct OverlapMask
     Kokkos::deep_copy(is_overlap_, true);
   }
 
+  // Construct from a precomputed per-DOF-holder host mask (e.g. an MFEM
+  // subdomain selected by element attribute). Indexed by local 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
   OverlapMask(size_t size, Omega_h::HostRead<Omega_h::I8> is_overlap_host)
     : is_overlap_("overlap_info", size)

src/pcms/coupler/serializer/xgc.h

@@ -20,23 +20,20 @@ class XGCFieldSerializer : public FieldSerializer<T>
   {
   }
 
-  int Serialize(const FieldData<T>& field, const FieldLayout& layout,
-                Rank1View<T, HostMemorySpace> buffer,
+  int Serialize(const Field<T>& field, Rank1View<T, HostMemorySpace> buffer,
                 Rank1View<const LO, HostMemorySpace> permutation) const override
   {
     if (!rank_participates_) {
       return 0;
     }
 
-    auto const* xgc_field = dynamic_cast<const XGCFieldData<T>*>(&field);
+    auto const* xgc_field =
+      dynamic_cast<const XGCFieldData<T>*>(&field.GetData());
     if (!xgc_field) {
       throw pcms_error("XGCFieldSerializer::Serialize: incompatible FieldData");
     }
 
     auto data = xgc_field->GetDOFHolderDataHost();
-    auto owned = layout.GetOwnedHost();
-    // Per-holder plan: owned[i]/permutation[i] index holders; a holder's
-    // num_components values form one contiguous block in the wire buffer.
     if (buffer.size() > 0) {
       const LO num_dof = static_cast<LO>(data.extent(0));
       const LO num_comp = static_cast<LO>(data.extent(1));
@@ -54,17 +51,17 @@ class XGCFieldSerializer : public FieldSerializer<T>
   }
 
   void Deserialize(
-    FieldData<T>& field, const FieldLayout& layout,
-    Rank1View<const T, HostMemorySpace> buffer,
+    Field<T>& field, Rank1View<const T, HostMemorySpace> buffer,
     Rank1View<const LO, HostMemorySpace> permutation) const override
   {
-    auto* xgc_field = dynamic_cast<XGCFieldData<T>*>(&field);
+    auto* xgc_field = dynamic_cast<XGCFieldData<T>*>(&field.GetData());
     if (!xgc_field) {
       throw pcms_error(
         "XGCFieldSerializer::Deserialize: incompatible FieldData");
     }
 
     auto current = xgc_field->GetDOFHolderDataHost();
+    const auto& layout = field.GetLayout();
     const LO num_dof = static_cast<LO>(current.extent(0));
     const LO num_comp = static_cast<LO>(current.extent(1));
     std::vector<T> full_data(current.size());

src/pcms/field/CMakeLists.txt

@@ -60,6 +60,18 @@ if (PCMS_ENABLE_OMEGA_H)
     )
 endif ()
 
+if (PCMS_ENABLE_MFEM)
+    list(APPEND PCMS_FIELD_SOURCES
+            layout/mfem.cpp
+            data/mfem.cpp
+    )
+    list(APPEND PCMS_FIELD_HEADERS
+            layout/mfem.hpp
+            data/mfem.hpp
+            function_space/mfem.hpp
+    )
+endif ()
+
 if (PCMS_ENABLE_MESHFIELDS)
     list(APPEND PCMS_FIELD_SOURCES
             layout/mesh_fields.cpp
@@ -103,6 +115,10 @@ if(PCMS_ENABLE_OMEGA_H)
     target_link_libraries(pcms_field PRIVATE Kokkos::kokkoskernels)
 endif()
 
+if(PCMS_ENABLE_MFEM)
+    target_link_libraries(pcms_field PUBLIC mfem)
+endif()
+
 if (PCMS_ENABLE_MESHFIELDS)
     target_link_libraries(pcms_field PUBLIC meshfields::meshfields)
 endif ()

src/pcms/field/data/mfem.cpp

@@ -0,0 +1,90 @@
+#include "pcms/field/data/mfem.hpp"
+
+#include "pcms/utility/arrays.h"
+#include "pcms/utility/assert.h"
+#include "pcms/utility/profile.h"
+
+namespace pcms
+{
+
+MFEMVertexFieldData::MFEMVertexFieldData(mfem::ParFiniteElementSpace& pfes,
+                                         mfem::ParGridFunction& gf,
+                                         FieldMetadata metadata)
+  : pfes_(pfes),
+    gf_(gf),
+    metadata_(metadata),
+    host_data_("mfem_field_data_host", pfes.GetNDofs(), 1),
+    device_data_("mfem_field_data_device", pfes.GetNDofs(), 1)
+{
+  PCMS_ALWAYS_ASSERT(pfes_.GetVDim() == 1);
+}
+
+const FieldMetadata& MFEMVertexFieldData::GetMetadata() const
+{
+  return metadata_;
+}
+
+Rank2View<const Real, HostMemorySpace>
+MFEMVertexFieldData::GetDOFHolderDataHost() const
+{
+  PCMS_FUNCTION_TIMER;
+
+  // Pull owner values into a true-DOF vector via the parallel restriction so
+  // every shared vertex reports its owner's value.
+  mfem::Vector true_values(pfes_.GetTrueVSize());
+  gf_.GetTrueDofs(true_values);
+
+  mfem::Array<int> vdofs;
+  const int nv = pfes_.GetNDofs();
+  for (int v = 0; v < nv; ++v) {
+    pfes_.GetVertexDofs(v, vdofs);
+    PCMS_ALWAYS_ASSERT(vdofs.Size() == 1);
+    const int lt = pfes_.GetLocalTDofNumber(vdofs[0]);
+    host_data_(v, 0) = (lt >= 0) ? static_cast<Real>(true_values[lt]) : Real{0};
+  }
+
+  return MakeConstRank2View(host_data_);
+}
+
+void MFEMVertexFieldData::SetDOFHolderDataHost(
+  Rank2View<const Real, HostMemorySpace> data)
+{
+  PCMS_FUNCTION_TIMER;
+  const int nv = pfes_.GetNDofs();
+  PCMS_ALWAYS_ASSERT(static_cast<int>(data.extent(0)) == nv);
+  PCMS_ALWAYS_ASSERT(data.extent(1) == 1);
+
+  // Scatter received owner values into a true-DOF vector, then distribute to
+  // all local DOFs (including shared, non-owned vertices) via the parallel
+  // prolongation.
+  mfem::Vector true_values(pfes_.GetTrueVSize());
+  mfem::Array<int> vdofs;
+  for (int v = 0; v < nv; ++v) {
+    pfes_.GetVertexDofs(v, vdofs);
+    PCMS_ALWAYS_ASSERT(vdofs.Size() == 1);
+    const int lt = pfes_.GetLocalTDofNumber(vdofs[0]);
+    if (lt >= 0) {
+      true_values[lt] = static_cast<double>(data(v, 0));
+    }
+  }
+
+  gf_.SetFromTrueDofs(true_values);
+}
+
+Rank2View<const Real, DeviceMemorySpace> MFEMVertexFieldData::GetDOFHolderData()
+  const
+{
+  GetDOFHolderDataHost();
+  DeepCopyMismatchLayouts(device_data_, host_data_);
+  return MakeConstRank2View(device_data_);
+}
+
+void MFEMVertexFieldData::SetDOFHolderData(
+  Rank2View<const Real, DeviceMemorySpace> data)
+{
+  CopyDeviceRank2ViewToDeviceView(device_data_, data);
+  DeepCopyMismatchLayouts(host_data_, device_data_);
+  SetDOFHolderDataHost(MakeConstRank2View(host_data_));
+}
+
+} // namespace pcms

src/pcms/field/data/mfem.hpp

@@ -0,0 +1,39 @@
+#ifndef PCMS_FIELD_DATA_MFEM_HPP
+#define PCMS_FIELD_DATA_MFEM_HPP
+
+#include "pcms/field/field_data.h"
+#include "pcms/field/field_metadata.h"
+#include "pcms/utility/arrays.h"
+
+#include <mfem.hpp>
+
+namespace pcms
+{
+
+// FieldData backend for an MFEM order-1 H1 (vertex) scalar field.
+class MFEMVertexFieldData : public FieldData<Real>
+{
+public:
+  MFEMVertexFieldData(mfem::ParFiniteElementSpace& pfes,
+                      mfem::ParGridFunction& gf, FieldMetadata metadata = {});
+
+  const FieldMetadata& GetMetadata() const override;
+
+  Rank2View<const Real, HostMemorySpace> GetDOFHolderDataHost() const override;
+  void SetDOFHolderDataHost(
+    Rank2View<const Real, HostMemorySpace> data) override;
+
+  Rank2View<const Real, DeviceMemorySpace> GetDOFHolderData() const override;
+  void SetDOFHolderData(Rank2View<const Real, DeviceMemorySpace> data) override;
+
+private:
+  mfem::ParFiniteElementSpace& pfes_;
+  mfem::ParGridFunction& gf_;
+  FieldMetadata metadata_;
+  mutable Kokkos::View<Real**, HostMemorySpace> host_data_;
+  mutable Kokkos::View<Real**, DeviceMemorySpace> device_data_;
+};
+
+} // namespace pcms
+
+#endif // PCMS_FIELD_DATA_MFEM_HPP