diff --git a/CMakeLists.txt b/CMakeLists.txt index 9d407a53..65417731 100644 --- a/CMakeLists.txt +++ b/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) @@ -119,6 +120,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) diff --git a/src/pcms/capi/client.cpp b/src/pcms/capi/client.cpp index 3d320094..4faada7b 100644 --- a/src/pcms/capi/client.cpp +++ b/src/pcms/capi/client.cpp @@ -1,5 +1,6 @@ #include "client.h" #include "pcms.h" +#include "pcms/field/function_space.h" #include "pcms/field/function_space/xgc.h" #include "pcms/field/data/xgc.h" #include "pcms/field/layout/xgc.h" @@ -23,11 +24,11 @@ namespace detail template struct XGCFieldRegistration { - XGCFunctionSpace function_space; + XGCFieldFactory function_space; MPI_Comm plane_comm; Rank1View data; - XGCFieldRegistration(XGCFunctionSpace fs, MPI_Comm comm, + XGCFieldRegistration(XGCFieldFactory fs, MPI_Comm comm, Rank1View d) : function_space(std::move(fs)), plane_comm(comm), data(d) { @@ -251,7 +252,7 @@ void pcms_create_xgc_field_adapter_t( { PCMS_ALWAYS_ASSERT((size > 0) ? (data != nullptr) : true); auto function_space = - pcms::XGCFunctionSpace(reverse_classification, in_overlap, size); + pcms::XGCFieldFactory(reverse_classification, in_overlap, size); pcms::Rank1View data_view( reinterpret_cast(data), size); field_adapter.emplace>( diff --git a/src/pcms/coupler/field_exchange_planner.cpp b/src/pcms/coupler/field_exchange_planner.cpp index 83ed0eb5..78cbedc0 100644 --- a/src/pcms/coupler/field_exchange_planner.cpp +++ b/src/pcms/coupler/field_exchange_planner.cpp @@ -12,6 +12,10 @@ namespace pcms namespace { +// Sentinel in a receive permutation: this local DOF's GID was not present in +// the received message. The deserializer skips it (negative => not received). +constexpr redev::LO kUnreceivedDof = -1; + struct PartitionMapping { std::vector indices; @@ -134,8 +138,16 @@ static redev::LOs ConstructPermutation( const auto start = ent_offsets[e]; const auto end = ent_offsets[e + 1]; - for (size_t i = start; i < end; ++i) - permutation.push_back(gid_to_buffer_index[e][local_gids[i]]); + for (size_t i = start; i < end; ++i) { + // A local DOF whose GID was not in the received message (e.g. it lies + // outside the sender's overlap mask) gets the sentinel kUnreceivedDof. + // The deserializer must skip these and preserve the field's existing + // value, rather than reading buffer[0]. Use find() rather than + // operator[] so missing keys are not silently inserted as 0. + const auto it = gid_to_buffer_index[e].find(local_gids[i]); + permutation.push_back( + it != gid_to_buffer_index[e].end() ? it->second : kUnreceivedDof); + } } REDEV_ALWAYS_ASSERT(permutation.size() == local_gids.size()); @@ -275,7 +287,8 @@ ExchangePlan GenericFieldExchangePlanner::BuildReceivePlan( void GenericFieldExchangePlanner::FillGidMessage( const FieldLayout& layout, const ExchangePlan& plan, - Rank1View gid_message) const + Rank1View gid_message, + const OverlapMask* overlap_mask) const { PCMS_FUNCTION_TIMER; PCMS_ALWAYS_ASSERT(static_cast(gid_message.size()) == plan.msg_size); @@ -286,11 +299,19 @@ void GenericFieldExchangePlanner::FillGidMessage( auto offsets = Rank1View( plan.offsets.data(), plan.offsets.size()); + // Participation must match BuildReversePartitionMap (owned AND in overlap); + // otherwise owned-but-non-overlap DOFs inflate the per-rank counts written + // here and corrupt the message. A default all-true mask is used when none is + // provided so behavior is unchanged for callers without an overlap mask. + OverlapMask default_mask(static_cast(gids.size())); + auto overlap = + (overlap_mask ? *overlap_mask : default_mask).GetMask(layout); + std::vector per_rank_offsets(plan.dest_ranks.size()); for (LO local_index = 0; local_index < static_cast(gids.size()); ++local_index) { - if (!owned[local_index]) + if (!owned[local_index] || !overlap[local_index]) continue; LO perm_index = plan.permutation[local_index]; diff --git a/src/pcms/coupler/field_exchange_planner.h b/src/pcms/coupler/field_exchange_planner.h index 4ea2fa54..d6278c93 100644 --- a/src/pcms/coupler/field_exchange_planner.h +++ b/src/pcms/coupler/field_exchange_planner.h @@ -32,7 +32,8 @@ class FieldExchangePlanner virtual void FillGidMessage( const FieldLayout& layout, const ExchangePlan& plan, - Rank1View gid_message) const = 0; + Rank1View gid_message, + const OverlapMask* overlap_mask = nullptr) const = 0; virtual ~FieldExchangePlanner() noexcept = default; }; @@ -51,7 +52,8 @@ class GenericFieldExchangePlanner : public FieldExchangePlanner void FillGidMessage( const FieldLayout& layout, const ExchangePlan& plan, - Rank1View gid_message) const override; + Rank1View gid_message, + const OverlapMask* overlap_mask = nullptr) const override; }; } // namespace pcms diff --git a/src/pcms/coupler/field_layout_communicator.cpp b/src/pcms/coupler/field_layout_communicator.cpp index 713ba98c..3a0d50c7 100644 --- a/src/pcms/coupler/field_layout_communicator.cpp +++ b/src/pcms/coupler/field_layout_communicator.cpp @@ -78,7 +78,8 @@ void FieldLayoutCommunicator::UpdateLayout() std::vector gid_message(plan_.msg_size); planner_->FillGidMessage( layout_, plan_, - Rank1View(gid_message.data(), gid_message.size())); + Rank1View(gid_message.data(), gid_message.size()), + overlap_mask_.get()); channel_.BeginSendCommunicationPhase(); gid_comm_.Send(gid_message.data()); diff --git a/src/pcms/coupler/field_serializer.h b/src/pcms/coupler/field_serializer.h index fba5af57..2ec7082f 100644 --- a/src/pcms/coupler/field_serializer.h +++ b/src/pcms/coupler/field_serializer.h @@ -34,10 +34,19 @@ class FieldSerializer Rank1View buffer, Rank1View permutation) const { + // Seed from the field's current values so DOFs that were not received are + // preserved rather than overwritten. This matters for masked coupling: a + // permutation entry < 0 (kUnreceivedDof) means the sender did not include + // that DOF's GID, so its existing value must be left untouched instead of + // reading buffer[0]. + auto current = field.GetDOFHolderDataHost(); Kokkos::View sorted("sorted", permutation.size()); + for (LO i = 0; i < static_cast(sorted.size()); ++i) { + sorted[i] = current[i]; + } auto owned = layout.GetOwnedHost(); for (LO i = 0; i < static_cast(sorted.size()); ++i) { - if (owned[i]) + if (owned[i] && permutation[i] >= 0) sorted[i] = buffer[permutation[i]]; } field.SetDOFHolderDataHost(make_const_array_view(sorted)); diff --git a/src/pcms/coupler/overlap_mask.h b/src/pcms/coupler/overlap_mask.h index db795b45..fef2d753 100644 --- a/src/pcms/coupler/overlap_mask.h +++ b/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 #include @@ -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 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 is_overlap_host) : is_overlap_("overlap_info", size) diff --git a/src/pcms/coupler/serializer/xgc.h b/src/pcms/coupler/serializer/xgc.h index 8a6f6cfe..d31b9b46 100644 --- a/src/pcms/coupler/serializer/xgc.h +++ b/src/pcms/coupler/serializer/xgc.h @@ -64,7 +64,9 @@ class XGCFieldSerializer : public FieldSerializer } if (rank_participates_) { for (LO i = 0; i < static_cast(current.size()); ++i) { - if (owned[i]) { + // permutation[i] < 0 (kUnreceivedDof) => this DOF's GID was not in the + // received message; preserve its current value (and avoid buffer[-1]). + if (owned[i] && permutation[i] >= 0) { full_data[i] = buffer[permutation[i]]; } } diff --git a/src/pcms/field/CMakeLists.txt b/src/pcms/field/CMakeLists.txt index 49f9b8ff..3c4f0286 100644 --- a/src/pcms/field/CMakeLists.txt +++ b/src/pcms/field/CMakeLists.txt @@ -9,6 +9,7 @@ set(PCMS_FIELD_HEADERS layout/empty.h out_of_bounds_policy.h field_data.h + field_factory.h function_space.h point_evaluator.h data/simple.h @@ -54,6 +55,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.h + data/mfem.h + function_space/mfem.h + ) +endif () + if (PCMS_ENABLE_MESHFIELDS) list(APPEND PCMS_FIELD_SOURCES layout/mesh_fields.cpp @@ -97,6 +110,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 () diff --git a/src/pcms/field/data/mfem.cpp b/src/pcms/field/data/mfem.cpp new file mode 100644 index 00000000..f1f667e3 --- /dev/null +++ b/src/pcms/field/data/mfem.cpp @@ -0,0 +1,89 @@ +#include "pcms/field/data/mfem.h" + +#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()), + device_data_("mfem_field_data_device", pfes.GetNDofs()) +{ + PCMS_ALWAYS_ASSERT(pfes_.GetVDim() == 1); +} + +const FieldMetadata& MFEMVertexFieldData::GetMetadata() const +{ + return metadata_; +} + +Rank1View +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 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) = (lt >= 0) ? static_cast(true_values[lt]) : Real{0}; + } + + return make_const_array_view(host_data_); +} + +void MFEMVertexFieldData::SetDOFHolderDataHost( + Rank1View data) +{ + PCMS_FUNCTION_TIMER; + PCMS_ALWAYS_ASSERT(static_cast(data.size()) == pfes_.GetNDofs()); + + // 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 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]); + if (lt >= 0) { + true_values[lt] = static_cast(data[v]); + } + } + + gf_.SetFromTrueDofs(true_values); +} + +Rank1View +MFEMVertexFieldData::GetDOFHolderData() const +{ + GetDOFHolderDataHost(); + Kokkos::deep_copy(device_data_, host_data_); + return make_const_array_view(device_data_); +} + +void MFEMVertexFieldData::SetDOFHolderData( + Rank1View data) +{ + PCMS_ALWAYS_ASSERT(static_cast(data.size()) == pfes_.GetNDofs()); + CopyDeviceRank1ViewToHostView(host_data_, data); + SetDOFHolderDataHost(make_const_array_view(host_data_)); +} + +} // namespace pcms diff --git a/src/pcms/field/data/mfem.h b/src/pcms/field/data/mfem.h new file mode 100644 index 00000000..d5ca1c52 --- /dev/null +++ b/src/pcms/field/data/mfem.h @@ -0,0 +1,47 @@ +#ifndef PCMS_FIELD_DATA_MFEM_H +#define PCMS_FIELD_DATA_MFEM_H + +#include "pcms/field/field_data.h" +#include "pcms/field/field_metadata.h" +#include "pcms/utility/arrays.h" + +#include + +namespace pcms +{ + +// FieldData backend for an MFEM order-1 H1 (vertex) scalar field. The +// coefficient store is the live mfem::ParGridFunction: Get/Set operate on it +// directly so the coupler and the MFEM solver share state. +// +// DOF-holder ordering is the local vertex ordering, matching MFEMLayout. To +// stay consistent across process boundaries, Get/Set round-trip through the +// FE space true DOFs: Get reads owner values via the parallel restriction and +// Set distributes received owner values via the parallel prolongation, filling +// shared (non-owned) vertices on this rank. +class MFEMVertexFieldData : public FieldData +{ +public: + MFEMVertexFieldData(mfem::ParFiniteElementSpace& pfes, + mfem::ParGridFunction& gf, FieldMetadata metadata = {}); + + const FieldMetadata& GetMetadata() const override; + + Rank1View GetDOFHolderDataHost() const override; + void SetDOFHolderDataHost( + Rank1View data) override; + + Rank1View GetDOFHolderData() const override; + void SetDOFHolderData(Rank1View data) override; + +private: + mfem::ParFiniteElementSpace& pfes_; + mfem::ParGridFunction& gf_; + FieldMetadata metadata_; + mutable Kokkos::View host_data_; + mutable Kokkos::View device_data_; +}; + +} // namespace pcms + +#endif // PCMS_FIELD_DATA_MFEM_H diff --git a/src/pcms/field/data/xgc.h b/src/pcms/field/data/xgc.h index 364b6042..2f9c8871 100644 --- a/src/pcms/field/data/xgc.h +++ b/src/pcms/field/data/xgc.h @@ -25,7 +25,7 @@ class XGCFieldData : public FieldData layout_->GetFullDataSize()); } - // Self-allocating constructor: XGCFunctionSpace::CreateFieldImpl uses this + // Self-allocating constructor: XGCFieldFactory::CreateFieldImpl uses this // to produce a field with internally-managed storage. XGCFieldData(std::shared_ptr layout, FieldMetadata metadata) diff --git a/src/pcms/field/field.h b/src/pcms/field/field.h index 98f0301f..54fb67cf 100644 --- a/src/pcms/field/field.h +++ b/src/pcms/field/field.h @@ -2,7 +2,6 @@ #define PCMS_COUPLING_FIELD_H #include "field_data.h" -#include "field_evaluator_factory.h" #include "field_layout.h" #include "pcms/utility/arrays.h" #include "pcms/utility/memory_spaces.h" @@ -13,14 +12,9 @@ namespace pcms { -class FunctionSpace; +class FieldFactory; -// Field is a composed per-field object: it owns coefficient data and holds -// a shared reference to the evaluator factory so the function space stays alive -// as long as the field does. -// -// Fields are typically created via LagrangeFunctionSpace::CreateField(). -// They are move-only (unique_ptr member). +// Field is FieldLayout (topology / coupling identity) plus owned FieldData template class Field { @@ -59,22 +53,18 @@ class Field class CtorKey { CtorKey() = default; - friend class FunctionSpace; + friend class FieldFactory; }; Field(CtorKey, std::shared_ptr layout, - std::shared_ptr> evaluator_factory, std::unique_ptr> data) - : layout_(std::move(layout)), - evaluator_factory_(std::move(evaluator_factory)), - data_(std::move(data)) + : layout_(std::move(layout)), data_(std::move(data)) { } - friend class FunctionSpace; + friend class FieldFactory; std::shared_ptr layout_; - std::shared_ptr> evaluator_factory_; std::unique_ptr> data_; }; diff --git a/src/pcms/field/field_factory.h b/src/pcms/field/field_factory.h new file mode 100644 index 00000000..8ea5bc8d --- /dev/null +++ b/src/pcms/field/field_factory.h @@ -0,0 +1,95 @@ +#ifndef PCMS_FIELD_FACTORY_H +#define PCMS_FIELD_FACTORY_H + +#include "field.h" +#include "field_data.h" +#include "field_layout.h" +#include "field_metadata.h" +#include "pcms/utility/common.h" +#include "pcms/utility/types.h" +#include + +namespace pcms +{ + +namespace detail +{ + +inline size_t ExpectedFlatFieldDataSize(const FieldLayout& layout) +{ + return static_cast(layout.GetNumOwnedDofHolder()) * + static_cast(layout.GetNumComponents()); +} + +} // namespace detail + +// Compile-time gate: true only for the five supported field value types. +template +inline constexpr bool is_supported_field_type_v = + std::is_same_v || std::is_same_v || + std::is_same_v || std::is_same_v || + std::is_same_v; + +// FieldFactory is an abstract type that can produce +// Fields over a known layout. This is the entire surface the coupler needs +// (GetLayout + CreateField). A FieldFactory makes discrete fields (values on +// DOF holders); it makes no claim about evaluating them as functions. +// +// Communication-only adapters (MFEM, XGC) are FieldFactories. Evaluatable +// spaces are FunctionSpaces, which derive from FieldFactory +class FieldFactory +{ +public: + virtual std::shared_ptr GetLayout() const noexcept = 0; + + virtual ~FieldFactory() noexcept = default; + + // Create a new field with freshly allocated data for this factory. + // Compile-time error for unsupported T; runtime error for T unsupported by + // the concrete backend. + template + [[nodiscard]] Field CreateField(FieldMetadata metadata = {}) const; + + // Expert API: wrap externally constructed field data into a Field for this + // factory. The concrete factory validates backend-specific field-data type + // and storage size compatibility. + template + [[nodiscard]] Field CreateField(std::unique_ptr> data) const; + +protected: + template + static Field WrapField(std::shared_ptr layout, + std::unique_ptr> data) + { + return Field(typename Field::CtorKey{}, std::move(layout), + std::move(data)); + } + + virtual FieldVariant CreateFieldImpl(Type value_type, + FieldMetadata metadata) const = 0; + + virtual FieldVariant CreateFieldImpl(FieldDataVariant data) const = 0; +}; + +template +Field FieldFactory::CreateField(FieldMetadata metadata) const +{ + static_assert(is_supported_field_type_v, + "T is not a supported field type"); + return std::get>(CreateFieldImpl(TypeEnumFromType(), metadata)); +} + +template +Field FieldFactory::CreateField(std::unique_ptr> data) const +{ + static_assert(is_supported_field_type_v, + "T is not a supported field type"); + if (!data) { + throw pcms_error("FieldFactory::CreateField: data must not be null"); + } + return std::get>(CreateFieldImpl(FieldDataVariant{std::move(data)})); +} + +} // namespace pcms + +#endif // PCMS_FIELD_FACTORY_H diff --git a/src/pcms/field/function_space.h b/src/pcms/field/function_space.h index a5e6f867..3badf34e 100644 --- a/src/pcms/field/function_space.h +++ b/src/pcms/field/function_space.h @@ -5,7 +5,7 @@ #include "evaluation_request.h" #include "field.h" #include "field_data.h" -#include "field_evaluator_factory.h" +#include "field_factory.h" #include "field_layout.h" #include "field_metadata.h" #include "out_of_bounds_policy.h" @@ -20,33 +20,8 @@ namespace pcms { -namespace detail -{ - -inline size_t ExpectedFlatFieldDataSize(const FieldLayout& layout) -{ - return static_cast(layout.GetNumOwnedDofHolder()) * - static_cast(layout.GetNumComponents()); -} - -} // namespace detail - -// Compile-time gate: true only for the five supported field value types. -template -inline constexpr bool is_supported_field_type_v = - std::is_same_v || std::is_same_v || - std::is_same_v || std::is_same_v || - std::is_same_v; - -// FunctionSpace is an abstract interface representing an evaluatable field -// space: layout, evaluation rules, and coordinate interpretation. -// -// Concrete implementations (e.g. LagrangeFunctionSpace) provide backends for -// specific discretizations or mesh types. -// -// FunctionSpace is used as the parameter type for operation objects such as -// Interpolator, so that operations are not coupled to a specific backend. -class FunctionSpace +// A FunctionSpace allows you to construct fields and evaluators for those fields +class FunctionSpace : public FieldFactory { public: virtual std::shared_ptr GetDiscretization() @@ -55,24 +30,8 @@ class FunctionSpace return GetLayout()->GetDiscretization(); } - virtual std::shared_ptr GetLayout() const noexcept = 0; - virtual CoordinateSystem GetCoordinateSystem() const noexcept = 0; - virtual ~FunctionSpace() noexcept = default; - - // Create a new field with freshly allocated data for this function space. - // Compile-time error for unsupported T; runtime error for T unsupported by - // the concrete backend. - template - [[nodiscard]] Field CreateField(FieldMetadata metadata = {}) const; - - // Expert API: wrap externally constructed field data into a Field for this - // function space. The concrete function space validates backend-specific - // field-data type and storage size compatibility. - template - [[nodiscard]] Field CreateField(std::unique_ptr> data) const; - // Create a point evaluator for the given evaluation request. // Compile-time error for unsupported T; runtime error for T or capability // unsupported by the concrete backend. @@ -87,44 +46,10 @@ class FunctionSpace const EvaluationRequest& request) const; protected: - template - static Field WrapField(std::shared_ptr layout, - std::unique_ptr> data, - std::shared_ptr> - evaluator_factory = nullptr) - { - return Field(typename Field::CtorKey{}, std::move(layout), - std::move(evaluator_factory), std::move(data)); - } - - virtual FieldVariant CreateFieldImpl(Type value_type, - FieldMetadata metadata) const = 0; - - virtual FieldVariant CreateFieldImpl(FieldDataVariant data) const = 0; - virtual PointEvaluatorVariant CreatePointEvaluatorImpl( Type value_type, const EvaluationRequest& request) const = 0; }; -template -Field FunctionSpace::CreateField(FieldMetadata metadata) const -{ - static_assert(is_supported_field_type_v, - "T is not a supported field type"); - return std::get>(CreateFieldImpl(TypeEnumFromType(), metadata)); -} - -template -Field FunctionSpace::CreateField(std::unique_ptr> data) const -{ - static_assert(is_supported_field_type_v, - "T is not a supported field type"); - if (!data) { - throw pcms_error("FunctionSpace::CreateField: data must not be null"); - } - return std::get>(CreateFieldImpl(FieldDataVariant{std::move(data)})); -} - template std::unique_ptr> FunctionSpace::CreatePointEvaluator( const EvaluationRequest& request) const diff --git a/src/pcms/field/function_space/lagrange.cpp b/src/pcms/field/function_space/lagrange.cpp index 8328e266..99c8f286 100644 --- a/src/pcms/field/function_space/lagrange.cpp +++ b/src/pcms/field/function_space/lagrange.cpp @@ -251,8 +251,7 @@ FieldVariant LagrangeFunctionSpace::CreateFieldImpl( [this](auto&& fd) -> FieldVariant { using FD = std::decay_t; using T = typename FD::element_type::value_type; - return WrapField(layout_, std::forward(fd), - evaluator_factory_); + return WrapField(layout_, std::forward(fd)); }, std::move(field_data)); } @@ -269,8 +268,7 @@ FieldVariant LagrangeFunctionSpace::CreateFieldImpl(FieldDataVariant data) const "LagrangeFunctionSpace: int8_t is not a supported field type"); } else { ValidateLagrangeWrappedFieldData(*layout_, *fd); - return WrapField(layout_, std::forward(fd), - evaluator_factory_); + return WrapField(layout_, std::forward(fd)); } }, std::move(data)); diff --git a/src/pcms/field/function_space/mfem.h b/src/pcms/field/function_space/mfem.h new file mode 100644 index 00000000..d1beafb3 --- /dev/null +++ b/src/pcms/field/function_space/mfem.h @@ -0,0 +1,79 @@ +#ifndef PCMS_FUNCTION_SPACE_MFEM_H +#define PCMS_FUNCTION_SPACE_MFEM_H + +#include "pcms/field/coordinate_system.h" +#include "pcms/field/data/mfem.h" +#include "pcms/field/field_factory.h" +#include "pcms/field/layout/mfem.h" +#include "pcms/utility/assert.h" +#include "pcms/utility/common.h" + +#include + +#include +#include + +namespace pcms +{ + +// Field factory for an MFEM order-1 H1 (vertex) scalar +// field. The produced field's data is bound to the live mfem::ParGridFunction, +// so coupled get/set operate on the solver's state. +// Lifetime: the factory and any Field it produces must not outlive +// pmesh / pfes / gf. +class MFEMFieldFactory : public FieldFactory +{ +public: + MFEMFieldFactory(mfem::ParMesh& pmesh, mfem::ParFiniteElementSpace& pfes, + mfem::ParGridFunction& gf, + CoordinateSystem coordinate_system) + : layout_( + std::make_shared(pmesh, pfes, coordinate_system)), + pfes_(pfes), + gf_(gf) + { + } + + [[nodiscard]] std::shared_ptr GetLayout() + const noexcept override + { + return layout_; + } + +protected: + [[nodiscard]] FieldVariant CreateFieldImpl( + Type value_type, FieldMetadata metadata) const override + { + if (value_type != Type::Real) { + throw pcms_error("MFEM adapter only supports Real (double) fields"); + } + return WrapField( + layout_, std::make_unique(pfes_, gf_, metadata)); + } + + [[nodiscard]] FieldVariant CreateFieldImpl( + FieldDataVariant data) const override + { + auto* real = std::get_if>>(&data); + if (real == nullptr || *real == nullptr || + dynamic_cast(real->get()) == nullptr) { + throw pcms_error( + "MFEMFieldFactory::CreateField: requires MFEMVertexFieldData"); + } + if ((*real)->GetDOFHolderDataHost().size() != + detail::ExpectedFlatFieldDataSize(*layout_)) { + throw pcms_error("MFEMFieldFactory::CreateField: field data size does " + "not match layout"); + } + return WrapField(layout_, std::move(*real)); + } + +private: + std::shared_ptr layout_; + mfem::ParFiniteElementSpace& pfes_; + mfem::ParGridFunction& gf_; +}; + +} // namespace pcms + +#endif // PCMS_FUNCTION_SPACE_MFEM_H diff --git a/src/pcms/field/function_space/polynomial_reconstruction.cpp b/src/pcms/field/function_space/polynomial_reconstruction.cpp index 060b9c24..d26beec1 100644 --- a/src/pcms/field/function_space/polynomial_reconstruction.cpp +++ b/src/pcms/field/function_space/polynomial_reconstruction.cpp @@ -93,8 +93,7 @@ FieldVariant PolynomialReconstructionFunctionSpace::CreateFieldImpl( "supported"); } else { return WrapField( - layout_, std::make_unique>(layout_, metadata), - evaluator_factory_); + layout_, std::make_unique>(layout_, metadata)); } }); } @@ -120,7 +119,7 @@ FieldVariant PolynomialReconstructionFunctionSpace::CreateFieldImpl( "PolynomialReconstructionFunctionSpace::CreateField: field data size " "does not match layout"); } - return WrapField(layout_, std::move(fd), evaluator_factory_); + return WrapField(layout_, std::move(fd)); } PointEvaluatorVariant diff --git a/src/pcms/field/function_space/spline.cpp b/src/pcms/field/function_space/spline.cpp index f4530fc8..167048b8 100644 --- a/src/pcms/field/function_space/spline.cpp +++ b/src/pcms/field/function_space/spline.cpp @@ -47,8 +47,7 @@ FieldVariant SplineFunctionSpace::CreateFieldImpl(Type value_type, throw pcms_error("SplineFunctionSpace: only double (Real) is supported"); } else { return WrapField( - layout_, std::make_unique>(layout_, metadata), - evaluator_factory_); + layout_, std::make_unique>(layout_, metadata)); } }); } @@ -70,7 +69,7 @@ FieldVariant SplineFunctionSpace::CreateFieldImpl(FieldDataVariant data) const "SplineFunctionSpace::CreateField: field data size does not match " "layout"); } - return WrapField(layout_, std::move(fd), evaluator_factory_); + return WrapField(layout_, std::move(fd)); } PointEvaluatorVariant SplineFunctionSpace::CreatePointEvaluatorImpl( diff --git a/src/pcms/field/function_space/xgc.h b/src/pcms/field/function_space/xgc.h index 164abc7a..c1a3564a 100644 --- a/src/pcms/field/function_space/xgc.h +++ b/src/pcms/field/function_space/xgc.h @@ -1,10 +1,12 @@ -#ifndef PCMS_XGC_FUNCTION_SPACE_H -#define PCMS_XGC_FUNCTION_SPACE_H +#ifndef PCMS_XGC_FIELD_FACTORY_H +#define PCMS_XGC_FIELD_FACTORY_H +#include "pcms/field/coordinate_system.h" #include "pcms/field/data/xgc.h" #include "pcms/field/field.h" -#include "pcms/field/function_space.h" +#include "pcms/field/field_factory.h" #include "pcms/field/field_metadata.h" +#include "pcms/utility/assert.h" #include "pcms/utility/common.h" #include #include @@ -12,12 +14,13 @@ namespace pcms { -class XGCFunctionSpace : public FunctionSpace +// Field factory for XGC fields. +class XGCFieldFactory : public FieldFactory { public: - XGCFunctionSpace(const ReverseClassificationVertex& reverse_classification, - std::function in_overlap, - LO num_plane_nodes) + XGCFieldFactory(const ReverseClassificationVertex& reverse_classification, + std::function in_overlap, + LO num_plane_nodes) : layout_(std::make_shared( reverse_classification, std::move(in_overlap), num_plane_nodes)) { @@ -29,11 +32,17 @@ class XGCFunctionSpace : public FunctionSpace return layout_; } - [[nodiscard]] CoordinateSystem GetCoordinateSystem() const noexcept override + [[nodiscard]] CoordinateSystem GetCoordinateSystem() const noexcept { return CoordinateSystem::XGC; } + [[nodiscard]] std::shared_ptr GetXGCLayout() + const noexcept + { + return layout_; + } + protected: [[nodiscard]] FieldVariant CreateFieldImpl( Type value_type, FieldMetadata metadata) const override @@ -55,12 +64,12 @@ class XGCFunctionSpace : public FunctionSpace PCMS_ALWAYS_ASSERT(fd != nullptr); if (dynamic_cast*>(fd.get()) == nullptr) { throw pcms_error( - "XGCFunctionSpace::CreateField: requires XGCFieldData"); + "XGCFieldFactory::CreateField: requires XGCFieldData"); } if (fd->GetDOFHolderDataHost().size() != static_cast(layout_->GetFullDataSize())) { throw pcms_error( - "XGCFunctionSpace::CreateField: field data size does not match " + "XGCFieldFactory::CreateField: field data size does not match " "layout"); } return WrapField(layout_, std::move(fd)); @@ -68,23 +77,10 @@ class XGCFunctionSpace : public FunctionSpace std::move(data)); } - [[nodiscard]] PointEvaluatorVariant CreatePointEvaluatorImpl( - Type /*value_type*/, const EvaluationRequest& /*request*/) const override - { - throw pcms_error("XGCFunctionSpace does not support point evaluation"); - } - -public: - [[nodiscard]] std::shared_ptr GetXGCLayout() - const noexcept - { - return layout_; - } - private: std::shared_ptr layout_; }; } // namespace pcms -#endif // PCMS_XGC_FUNCTION_SPACE_H +#endif // PCMS_XGC_FIELD_FACTORY_H diff --git a/src/pcms/field/layout/mfem.cpp b/src/pcms/field/layout/mfem.cpp new file mode 100644 index 00000000..245bc801 --- /dev/null +++ b/src/pcms/field/layout/mfem.cpp @@ -0,0 +1,189 @@ +#include "pcms/field/layout/mfem.h" + +#include "pcms/discretization/discretization/point_cloud.hpp" +#include "pcms/utility/arrays.h" +#include "pcms/utility/assert.h" +#include "pcms/utility/profile.h" + +#include + +namespace pcms +{ + +namespace +{ + +std::shared_ptr MakeVertexDiscretization( + int dim, Kokkos::View coords_host) +{ + return std::make_shared(dim, coords_host); +} + +} // namespace + +MFEMLayout::MFEMLayout(mfem::ParMesh& pmesh, + mfem::ParFiniteElementSpace& pfes, + CoordinateSystem coordinate_system) + : pmesh_(pmesh), + pfes_(pfes), + dim_(pmesh.SpaceDimension()), + coordinate_system_(coordinate_system) +{ + PCMS_FUNCTION_TIMER; + + AssertVertexScalarSpace(); + + const int nv = pmesh_.GetNV(); + + // Host coordinates, one row per vertex. GetVertex returns a pointer to the + // contiguous spatial coordinates of the vertex, independent of MFEM's + // internal storage ordering. + Kokkos::View coords_host("mfem_coords_host", nv, + dim_); + for (int v = 0; v < nv; ++v) { + const double* x = pmesh_.GetVertex(v); + for (int d = 0; d < dim_; ++d) { + coords_host(v, d) = static_cast(x[d]); + } + } + coords_ = Kokkos::View("mfem_coords", nv, dim_); + Kokkos::deep_copy(coords_, coords_host); + + // Global vertex ids. + mfem::Array vertex_gids; + pmesh_.GetGlobalVertexIndices(vertex_gids); + gids_host_ = Kokkos::View("mfem_gids", nv); + + // Ownership: a vertex DOF holder is owned by this rank iff it maps to a + // local true DOF. Shared vertices are owned by exactly one rank, so every + // global DOF holder has a single owner across the communicator. + owned_host_ = Kokkos::View("mfem_owned", nv); + + // Placeholder classification: vertex entity dim (0) and local-index ids. + // RCB partitioning ignores these; they exist to satisfy the FieldLayout and + // OverlapMask contracts. + classification_dims_host_ = + Kokkos::View("mfem_class_dims", nv); + classification_ids_host_ = + Kokkos::View("mfem_class_ids", nv); + + mfem::Array vdofs; + for (int v = 0; v < nv; ++v) { + pfes_.GetVertexDofs(v, vdofs); + PCMS_ALWAYS_ASSERT(vdofs.Size() == 1); + const int dof = vdofs[0]; + + gids_host_(v) = static_cast(vertex_gids[v]); + owned_host_(v) = (pfes_.GetLocalTDofNumber(dof) >= 0); + classification_dims_host_(v) = 0; + classification_ids_host_(v) = v; + } + + discretization_ = MakeVertexDiscretization(dim_, coords_host); +} + +void MFEMLayout::AssertVertexScalarSpace() const +{ + // Order-1 H1, single scalar component: exactly one DOF per vertex. + PCMS_ALWAYS_ASSERT(pfes_.GetVDim() == 1); + PCMS_ALWAYS_ASSERT(pfes_.GetNDofs() == pmesh_.GetNV()); +} + +std::shared_ptr MFEMLayout::GetDiscretization() + const noexcept +{ + return discretization_; +} + +int MFEMLayout::GetNumComponents() const +{ + return 1; +} + +LO MFEMLayout::GetNumOwnedDofHolder() const +{ + // Number of local DOF holders (all vertices on this rank). The owned mask + // distinguishes the single-owner subset used for communication. + return static_cast(pmesh_.GetNV()); +} + +GO MFEMLayout::GetNumGlobalDofHolder() const +{ + // Each globally unique vertex corresponds to one true DOF. + return static_cast(pfes_.GlobalTrueVSize()); +} + +Rank1View MFEMLayout::GetOwnedHost() const +{ + return make_const_array_view(owned_host_); +} + +GlobalIDView MFEMLayout::GetGidsHost() const +{ + return GlobalIDView(gids_host_.data(), gids_host_.size()); +} + +CoordinateView MFEMLayout::GetDOFHolderCoordinates() const +{ + return CoordinateView{coordinate_system_, + MakeConstRank2View(coords_)}; +} + +bool MFEMLayout::IsDistributed() const +{ + return true; +} + +EntOffsetsArray MFEMLayout::GetEntOffsets() const +{ + EntOffsetsArray offsets{}; + offsets.fill(0); + const auto n = static_cast(GetNumOwnedDofHolder()); + // Vertex DOF holders occupy entity dimension 0. + for (int i = 1; i < ent_offsets_len; ++i) { + offsets[i] = n; + } + return offsets; +} + +int MFEMLayout::GetDimension() const +{ + return dim_; +} + +Rank1View +MFEMLayout::GetDOFHolderClassificationDimensionsHost() const +{ + return make_const_array_view(classification_dims_host_); +} + +Rank1View +MFEMLayout::GetDOFHolderClassificationIdsHost() const +{ + return make_const_array_view(classification_ids_host_); +} + +Kokkos::View MFEMLayout::OverlapMaskFromAttribute( + mfem::ParMesh& pmesh, int attribute) +{ + PCMS_FUNCTION_TIMER; + + const int nv = pmesh.GetNV(); + Kokkos::View overlap("mfem_overlap_mask", nv); + Kokkos::deep_copy(overlap, false); + + mfem::Array verts; + for (int e = 0; e < pmesh.GetNE(); ++e) { + if (pmesh.GetAttribute(e) != attribute) { + continue; + } + pmesh.GetElementVertices(e, verts); + for (int j = 0; j < verts.Size(); ++j) { + overlap(verts[j]) = true; + } + } + + return overlap; +} + +} // namespace pcms diff --git a/src/pcms/field/layout/mfem.h b/src/pcms/field/layout/mfem.h new file mode 100644 index 00000000..f780e64e --- /dev/null +++ b/src/pcms/field/layout/mfem.h @@ -0,0 +1,80 @@ +#ifndef PCMS_FIELD_LAYOUT_MFEM_H +#define PCMS_FIELD_LAYOUT_MFEM_H + +#include "pcms/field/coordinate_system.h" +#include "pcms/field/field_layout.h" + +#include + +#include + +namespace pcms +{ + +// Field layout for an MFEM order-1 H1 (vertex) scalar field. +// +// Each DOF holder is a mesh vertex; DOF-holder ordering follows the MFEM local +// vertex ordering. This is the minimal layout matching the original MFEM +// adapter scope: vertex-only, order-1, single component. +// +// Partitioning relies on DOF-holder coordinates (redev::RCBPtn). The +// classification arrays required by FieldLayout are filled with placeholders +// (entity dim 0, id = local index) because MFEM has no Omega_h-style geometric +// classification; an RCB partition ignores them. +class MFEMLayout : public FieldLayout +{ +public: + MFEMLayout(mfem::ParMesh& pmesh, mfem::ParFiniteElementSpace& pfes, + CoordinateSystem coordinate_system); + + std::shared_ptr GetDiscretization() + const noexcept override; + + int GetNumComponents() const override; + LO GetNumOwnedDofHolder() const override; + GO GetNumGlobalDofHolder() const override; + + Rank1View GetOwnedHost() const override; + GlobalIDView GetGidsHost() const override; + CoordinateView GetDOFHolderCoordinates() const override; + + [[nodiscard]] bool IsDistributed() const override; + EntOffsetsArray GetEntOffsets() const override; + int GetDimension() const override; + + Rank1View + GetDOFHolderClassificationDimensionsHost() const override; + + Rank1View GetDOFHolderClassificationIdsHost() + const override; + + mfem::ParMesh& GetMesh() const noexcept { return pmesh_; } + mfem::ParFiniteElementSpace& GetFESpace() const noexcept { return pfes_; } + + // Build a per-vertex overlap mask from MFEM element attributes: a vertex is + // in the overlap if it is incident to at least one element with the given + // attribute. Mirrors the create_mask strategy used in the mfem-pcms-example. + // The returned host array is indexed by local vertex (DOF holder) and can be + // passed to pcms::OverlapMask. + static Kokkos::View OverlapMaskFromAttribute( + mfem::ParMesh& pmesh, int attribute); + +private: + void AssertVertexScalarSpace() const; + + mfem::ParMesh& pmesh_; + mfem::ParFiniteElementSpace& pfes_; + int dim_; + CoordinateSystem coordinate_system_; + + Kokkos::View coords_; + Kokkos::View owned_host_; + Kokkos::View gids_host_; + Kokkos::View classification_dims_host_; + Kokkos::View classification_ids_host_; + std::shared_ptr discretization_; +}; + +} // namespace pcms + +#endif // PCMS_FIELD_LAYOUT_MFEM_H diff --git a/src/pcms/transfer/field_compatibility.hpp b/src/pcms/transfer/field_compatibility.hpp new file mode 100644 index 00000000..c60178a4 --- /dev/null +++ b/src/pcms/transfer/field_compatibility.hpp @@ -0,0 +1,41 @@ +#ifndef PCMS_TRANSFER_FIELD_COMPATIBILITY_H +#define PCMS_TRANSFER_FIELD_COMPATIBILITY_H + +#include "pcms/discretization/discretization.h" +#include "pcms/field/field.h" +#include "pcms/field/field_layout.h" +#include "pcms/utility/common.h" +#include + +namespace pcms +{ +namespace detail +{ + +// Verify that a field handed to a TransferOperator::Apply belongs to the same +// space the operator was constructed from. +// Throws pcms_error on mismatch. +template +inline void CheckTransferFieldLayout(const Field& field, + const FieldLayout& expected, + const char* role) +{ + const FieldLayout& actual = field.GetLayout(); + if (&actual == &expected) { + return; + } + auto actual_disc = actual.GetDiscretization(); + auto expected_disc = expected.GetDiscretization(); + if (actual_disc && expected_disc && + actual_disc->SameEntities(*expected_disc)) { + return; + } + throw pcms_error(std::string("TransferOperator: ") + role + + " field does not match the function space the operator was " + "constructed from"); +} + +} // namespace detail +} // namespace pcms + +#endif // PCMS_TRANSFER_FIELD_COMPATIBILITY_H diff --git a/src/pcms/transfer/interpolator.h b/src/pcms/transfer/interpolator.h index d60b31bb..1fc24e8d 100644 --- a/src/pcms/transfer/interpolator.h +++ b/src/pcms/transfer/interpolator.h @@ -6,6 +6,7 @@ #include "pcms/field/function_space.h" #include "pcms/field/out_of_bounds_policy.h" #include "pcms/field/point_evaluator.h" +#include "pcms/transfer/field_compatibility.hpp" #include "pcms/utility/arrays.h" #include "pcms/utility/memory_spaces.h" #include "pcms/utility/profile.h" @@ -41,6 +42,8 @@ class Interpolator : public TransferOperator .GetCoordinates() .extent(0))), n_comp_(target_space.GetLayout()->GetNumComponents()), + source_layout_(source_space.GetLayout()), + target_layout_(target_space.GetLayout()), evaluator_(source_space.CreatePointEvaluator( EvaluationRequest::FromFunctionSpace(target_space, policy))) { @@ -51,6 +54,8 @@ class Interpolator : public TransferOperator void Apply(const Field& source, Field& target) const override { PCMS_FUNCTION_TIMER; + detail::CheckTransferFieldLayout(source, *source_layout_, "source"); + detail::CheckTransferFieldLayout(target, *target_layout_, "target"); const LO num_points = num_points_; const int n_comp = n_comp_; Kokkos::View output("interp_output", num_points, @@ -72,6 +77,8 @@ class Interpolator : public TransferOperator private: LO num_points_; int n_comp_; + std::shared_ptr source_layout_; + std::shared_ptr target_layout_; std::unique_ptr> evaluator_; }; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 8412bac2..5ea45764 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -401,6 +401,10 @@ if(Catch2_FOUND) test_point_evaluator.cpp) endif() + if(PCMS_ENABLE_MFEM) + list(APPEND PCMS_UNIT_TEST_SOURCES test_mfem_adapter.cpp) + endif() + if(PCMS_ENABLE_MESHFIELDS) list(APPEND PCMS_UNIT_TEST_SOURCES test_load_vector.cpp) @@ -417,6 +421,10 @@ if(Catch2_FOUND) target_link_libraries(unit_tests PUBLIC meshfields::meshfields) endif() + if(PCMS_ENABLE_MFEM) + target_link_libraries(unit_tests PUBLIC mfem) + endif() + target_link_libraries(unit_tests PUBLIC Catch2::Catch2 pcms::core pcms::transfer Kokkos::kokkoskernels) target_include_directories(unit_tests PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) @@ -452,6 +460,34 @@ else() message(WARNING "Catch2 not found. Disabling Unit Tests") endif() +if(PCMS_ENABLE_MFEM) + add_executable(test_mfem_coupling test_mfem_coupling.cpp) + target_link_libraries(test_mfem_coupling PUBLIC pcms::core mfem) + if(HOST_NPROC GREATER_EQUAL 2) + dual_mpi_test( + TESTNAME + test_mfem_coupling + TIMEOUT + 60 + NAME1 + server + EXE1 + $ + PROCS1 + 1 + ARGS1 + -1 + NAME2 + client + EXE2 + $ + PROCS2 + 1 + ARGS2 + 0) + endif() +endif() + if(PCMS_ENABLE_C) find_package(Kokkos REQUIRED) if(PCMS_ENABLE_XGC) diff --git a/test/test_mfem_adapter.cpp b/test/test_mfem_adapter.cpp new file mode 100644 index 00000000..87a6e712 --- /dev/null +++ b/test/test_mfem_adapter.cpp @@ -0,0 +1,139 @@ +#include + +#include +#include +#include +#include + +#include + +#include + +namespace +{ + +double LinearField(const mfem::Vector& x) +{ + return 1.0 + 2.0 * x[0] + 3.0 * x[1]; +} + +} // namespace + +TEST_CASE("MFEM vertex-scalar field adapter") +{ + // Common setup (re-run for each SECTION): 2x2 Cartesian quad mesh + // (9 vertices), order-1 H1 scalar space. + auto serial = mfem::Mesh::MakeCartesian2D(2, 2, mfem::Element::QUADRILATERAL); + mfem::ParMesh pmesh(MPI_COMM_WORLD, serial); + mfem::H1_FECollection fec(1, pmesh.Dimension()); + mfem::ParFiniteElementSpace pfes(&pmesh, &fec); + mfem::ParGridFunction gf(&pfes); + gf = 0.0; + + SECTION("layout reports vertex-scalar properties") + { + pcms::MFEMLayout layout(pmesh, pfes, pcms::CoordinateSystem::Cartesian); + + REQUIRE(layout.GetNumComponents() == 1); + REQUIRE(layout.GetDimension() == 2); + REQUIRE(layout.GetNumOwnedDofHolder() == pmesh.GetNV()); + REQUIRE(layout.IsDistributed()); + + auto gids = layout.GetGidsHost(); + REQUIRE(static_cast(gids.size()) == pmesh.GetNV()); + + auto coords = layout.GetDOFHolderCoordinates().GetCoordinates(); + REQUIRE(static_cast(coords.extent(0)) == pmesh.GetNV()); + REQUIRE(static_cast(coords.extent(1)) == 2); + + // Owned mask covers exactly the global unique vertices across ranks. + auto owned = layout.GetOwnedHost(); + int local_owned = 0; + for (size_t i = 0; i < owned.size(); ++i) { + if (owned[i]) + ++local_owned; + } + int total_owned = 0; + MPI_Allreduce(&local_owned, &total_owned, 1, MPI_INT, MPI_SUM, + MPI_COMM_WORLD); + REQUIRE(total_owned == static_cast(layout.GetNumGlobalDofHolder())); + } + + SECTION("field data round-trips through the grid function") + { + mfem::FunctionCoefficient coeff(LinearField); + gf.ProjectCoefficient(coeff); + + pcms::MFEMLayout layout(pmesh, pfes, pcms::CoordinateSystem::Cartesian); + pcms::MFEMVertexFieldData data(pfes, gf); + + auto host = data.GetDOFHolderDataHost(); + REQUIRE(static_cast(host.size()) == pmesh.GetNV()); + + std::vector captured(host.size()); + for (size_t v = 0; v < host.size(); ++v) { + captured[v] = host[v]; + } + + // Write the captured values back and confirm the grid function is restored. + gf = 0.0; + Kokkos::View in("in", captured.size()); + for (size_t v = 0; v < captured.size(); ++v) { + in(v) = captured[v]; + } + data.SetDOFHolderDataHost(pcms::make_const_array_view(in)); + + auto host2 = data.GetDOFHolderDataHost(); + auto owned = layout.GetOwnedHost(); + for (size_t v = 0; v < host2.size(); ++v) { + if (owned[v]) { + REQUIRE(host2[v] == captured[v]); + } + } + } + + SECTION("serializer identity round-trip on a single rank") + { + int nproc = 0; + MPI_Comm_size(MPI_COMM_WORLD, &nproc); + if (nproc != 1) { + SUCCEED("serializer identity round-trip is only checked on one rank"); + return; + } + + mfem::FunctionCoefficient coeff(LinearField); + gf.ProjectCoefficient(coeff); + + auto fs = pcms::MFEMFieldFactory( + pmesh, pfes, gf, pcms::CoordinateSystem::Cartesian); + auto field = fs.CreateField(); + const auto& layout = field.GetLayout(); + + const auto n = static_cast(layout.GetNumOwnedDofHolder()); + + // Identity permutation: buffer index equals DOF-holder index. + Kokkos::View perm("perm", n); + for (size_t i = 0; i < n; ++i) { + perm(i) = static_cast(i); + } + Kokkos::View buffer("buffer", n); + + pcms::FieldSerializer serializer; + serializer.Serialize(field.GetData(), layout, + pcms::make_array_view(buffer), + pcms::make_const_array_view(perm)); + + // Deserialize into a fresh field bound to a second grid function. + mfem::ParGridFunction gf2(&pfes); + gf2 = 0.0; + pcms::MFEMVertexFieldData data2(pfes, gf2); + serializer.Deserialize(data2, layout, pcms::make_const_array_view(buffer), + pcms::make_const_array_view(perm)); + + auto restored = data2.GetDOFHolderDataHost(); + auto original = field.GetDOFHolderDataHost(); + for (size_t v = 0; v < restored.size(); ++v) { + REQUIRE(restored[v] == original[v]); + } + } +} diff --git a/test/test_mfem_coupling.cpp b/test/test_mfem_coupling.cpp new file mode 100644 index 00000000..9b2508b5 --- /dev/null +++ b/test/test_mfem_coupling.cpp @@ -0,0 +1,212 @@ +// Two-application coupling test for the MFEM field adapter. +// +// A client application owns an MFEM order-1 vertex scalar field and sends it, +// restricted to a masked overlap domain, to a rendezvous "coupler" server. +// The overlap domain is selected from MFEM element attributes following the +// create_mask strategy in the mfem-pcms-example: a vertex participates if it is +// incident to an element with the target attribute. Routing uses a single-rank +// RCB (coordinate) partition. +// +// Usage: test_mfem_coupling + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +using pcms::Real; + +namespace +{ + +constexpr int TargetAttribute = 2; + +// Seeded into receiver DOFs outside the overlap before receiving. Distinct from +// any sent value (gid + 1 >= 1), so a correct masked receive must leave it +// untouched. +constexpr Real OutsideOverlapSentinel = -7.0; + +// Build a 4x4 quad mesh and tag the left half (centroid x < 0.5) with the +// target attribute; the rest keep attribute 1. The same construction runs on +// both apps so vertex global ids and coordinates match. +mfem::Mesh MakeAttributedMesh() +{ + auto mesh = mfem::Mesh::MakeCartesian2D(4, 4, mfem::Element::QUADRILATERAL); + for (int e = 0; e < mesh.GetNE(); ++e) { + mfem::Vector center; + mesh.GetElementCenter(e, center); + mesh.SetAttribute(e, center[0] < 0.5 ? TargetAttribute : 1); + } + mesh.SetAttributes(); + return mesh; +} + +redev::Partition MakeRCBPartition(int dim) +{ + redev::LOs ranks(1); + std::iota(ranks.begin(), ranks.end(), 0); + redev::Reals cuts = {0}; + return redev::Partition{redev::RCBPtn{dim, ranks, cuts}}; +} + +// Expected sent value at a vertex: its global id + 1 (strictly positive so it +// is distinguishable from the receiver's initial zero state). +Real ExpectedValue(pcms::GO gid) +{ + return static_cast(gid + 1); +} + +int RunClient(MPI_Comm comm) +{ + auto serial = MakeAttributedMesh(); + mfem::ParMesh pmesh(comm, serial); + mfem::H1_FECollection fec(1, pmesh.Dimension()); + mfem::ParFiniteElementSpace pfes(&pmesh, &fec); + mfem::ParGridFunction gf(&pfes); + gf = 0.0; + + pcms::Coupler cpl("mfem_overlap_coupler", comm, false, redev::Partition{}); + auto* app = cpl.AddApplication("mfem_app"); + + auto fs = pcms::MFEMFieldFactory(pmesh, pfes, gf, + pcms::CoordinateSystem::Cartesian); + auto layout = fs.GetLayout(); + + auto overlap_view = + pcms::MFEMLayout::OverlapMaskFromAttribute(pmesh, TargetAttribute); + app->SetLayoutOverlapMask( + "field", std::make_unique( + static_cast(layout->GetNumOwnedDofHolder()), + overlap_view)); + app->AddLayout("field", layout); + + auto handle = app->AddField("field", fs.CreateField()); + + // Seed the field so each vertex holds (gid + 1). + auto gids = layout->GetGidsHost(); + const auto n = static_cast(layout->GetNumOwnedDofHolder()); + Kokkos::View values("client_values", n); + for (size_t v = 0; v < n; ++v) { + values(v) = ExpectedValue(gids[v]); + } + handle.GetField().SetDOFHolderDataHost(pcms::make_const_array_view(values)); + + app->SendPhase([&]() { handle.Send(); }); + return 0; +} + +int RunServer(MPI_Comm comm) +{ + auto serial = MakeAttributedMesh(); + mfem::ParMesh pmesh(comm, serial); + mfem::H1_FECollection fec(1, pmesh.Dimension()); + mfem::ParFiniteElementSpace pfes(&pmesh, &fec); + mfem::ParGridFunction gf(&pfes); + gf = 0.0; + + pcms::Coupler cpl("mfem_overlap_coupler", comm, true, + MakeRCBPartition(pmesh.SpaceDimension())); + auto* app = cpl.AddApplication("mfem_app"); + + auto fs = pcms::MFEMFieldFactory(pmesh, pfes, gf, + pcms::CoordinateSystem::Cartesian); + auto layout = fs.GetLayout(); + app->AddLayout("field", layout); + auto handle = app->AddField("field", fs.CreateField()); + + // Seed the whole receiver field with a sentinel. The masked send only carries + // the overlap DOFs, so a correct receive must overwrite only those and leave + // every DOF outside the overlap holding the sentinel. + { + const auto n = static_cast(layout->GetNumOwnedDofHolder()); + Kokkos::View seed("server_seed", n); + for (size_t v = 0; v < n; ++v) { + seed(v) = OutsideOverlapSentinel; + } + handle.GetField().SetDOFHolderDataHost(pcms::make_const_array_view(seed)); + } + + app->ReceivePhase([&]() { handle.Receive(); }); + + // Verify: every overlap vertex received the expected value. The overlap set + // is recomputed from the identical mesh's attributes. + auto overlap = + pcms::MFEMLayout::OverlapMaskFromAttribute(pmesh, TargetAttribute); + auto gids = layout->GetGidsHost(); + auto received = handle.GetField().GetDOFHolderDataHost(); + + int overlap_count = 0; + int mismatches = 0; + int corrupted_outside = 0; + for (size_t v = 0; v < received.size(); ++v) { + if (overlap(v)) { + ++overlap_count; + const Real expected = ExpectedValue(gids[v]); + if (received[v] != expected) { + ++mismatches; + std::cerr << "Mismatch at vertex " << v << ": expected " << expected + << " got " << received[v] << "\n"; + } + } else if (received[v] != OutsideOverlapSentinel) { + // A DOF outside the overlap was overwritten on receive. + ++corrupted_outside; + } + } + + const int nv = pmesh.GetNV(); + std::cout << "MFEM overlap coupling: " << overlap_count << " / " << nv + << " overlap vertices received; " << corrupted_outside + << " DOFs outside the overlap corrupted\n"; + + if (overlap_count == 0 || overlap_count == nv) { + std::cerr << "Overlap mask is trivial (count=" << overlap_count + << ", nv=" << nv << "); test is not meaningful\n"; + return 1; + } + if (mismatches != 0) { + std::cerr << "MFEM overlap coupling FAILED with " << mismatches + << " mismatches\n"; + return 1; + } + if (corrupted_outside != 0) { + std::cerr << "MFEM overlap coupling FAILED: " << corrupted_outside + << " DOFs outside the overlap were overwritten on receive " + "(expected to be preserved)\n"; + return 1; + } + std::cout << "MFEM overlap coupling PASSED\n"; + return 0; +} + +} // namespace + +int main(int argc, char** argv) +{ + MPI_Init(&argc, &argv); + int rc = 0; + { + Kokkos::ScopeGuard kokkos(argc, argv); + if (argc != 2) { + std::cerr << "Usage: " << argv[0] << " \n"; + MPI_Finalize(); + return EXIT_FAILURE; + } + const int role = std::atoi(argv[1]); + try { + rc = (role == -1) ? RunServer(MPI_COMM_WORLD) : RunClient(MPI_COMM_WORLD); + } catch (const std::exception& e) { + std::cerr << "Exception: " << e.what() << "\n"; + rc = 1; + } + } + MPI_Finalize(); + return rc; +} diff --git a/test/test_proxy_coupling_xgc_server.cpp b/test/test_proxy_coupling_xgc_server.cpp index 4a76b748..a46a3873 100644 --- a/test/test_proxy_coupling_xgc_server.cpp +++ b/test/test_proxy_coupling_xgc_server.cpp @@ -51,7 +51,7 @@ void xgc_coupler(MPI_Comm comm, Omega_h::Mesh& mesh, std::string_view cpn_file) // FIXME: The current C/Fortran proxy API couples layout registration to // field registration, so each XGC plane is registered as a separate layout // communicator even though the layouts are geometrically identical. - auto function_space = pcms::XGCFunctionSpace( + auto function_space = pcms::XGCFieldFactory( rc, ts::IsModelEntInOverlap{}, static_cast(mesh.nverts())); auto field = function_space.CreateField( std::make_unique>( diff --git a/test/test_proxy_coupling_xgc_server_overlap.cpp b/test/test_proxy_coupling_xgc_server_overlap.cpp index 0161ee4e..08176372 100644 --- a/test/test_proxy_coupling_xgc_server_overlap.cpp +++ b/test/test_proxy_coupling_xgc_server_overlap.cpp @@ -54,7 +54,7 @@ void xgc_coupler_with_overlap(MPI_Comm comm, Omega_h::Mesh& mesh, application->SetLayoutOverlapMask(ss.str(), std::move(overlap_mask)); - auto function_space = pcms::XGCFunctionSpace( + auto function_space = pcms::XGCFieldFactory( rc, ts::IsModelEntInOverlap{}, static_cast(mesh.nverts())); application->AddLayout(ss.str(), function_space.GetLayout()); diff --git a/test/test_xgc_field_data.cpp b/test/test_xgc_field_data.cpp index de3f8411..9d448e49 100644 --- a/test/test_xgc_field_data.cpp +++ b/test/test_xgc_field_data.cpp @@ -104,11 +104,11 @@ TEST_CASE("XGC FieldData serializer preserves inactive entries") } } -TEST_CASE("XGCFunctionSpace creates fields and rejects evaluator access") +TEST_CASE("XGCFieldFactory creates fields") { static constexpr int data_size = 16; auto rc = create_dummy_rc(data_size); - pcms::XGCFunctionSpace function_space(rc, in_overlap, data_size); + pcms::XGCFieldFactory function_space(rc, in_overlap, data_size); std::vector data(data_size); std::iota(data.begin(), data.end(), 0.0); @@ -119,14 +119,4 @@ TEST_CASE("XGCFunctionSpace creates fields and rejects evaluator access") REQUIRE(&field.GetLayout() == function_space.GetLayout().get()); REQUIRE(function_space.GetCoordinateSystem() == pcms::CoordinateSystem::XGC); - // XGC does not support point evaluation; CreatePointEvaluator throws. - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - pcms::Rank2View - empty_coords{nullptr, 0, 2}; - REQUIRE_THROWS_AS(function_space.CreatePointEvaluator( - pcms::EvaluationRequest::FromCoordinates( - pcms::CoordinateView{ - pcms::CoordinateSystem::XGC, empty_coords})), - pcms::pcms_error); }