diff --git a/src/pcms/coupler/field_communicator.hpp b/src/pcms/coupler/field_communicator.hpp index 3d1014a0..b1874507 100644 --- a/src/pcms/coupler/field_communicator.hpp +++ b/src/pcms/coupler/field_communicator.hpp @@ -35,7 +35,11 @@ class FieldCommunicator serializer_(std::move(serializer)) { PCMS_ALWAYS_ASSERT(&layout_comm.GetLayout() == &field.GetLayout()); - comm_buffer_.resize(layout_comm.GetMsgSize()); + // The exchange plan is per DOF holder; the field payload carries + // num_components values per holder, so the buffer is scaled accordingly. + comm_buffer_.resize( + layout_comm.GetMsgSize() * + static_cast(field.GetLayout().GetNumComponents())); comm_ = layout_comm_.GetChannel().CreateComm(name, layout_comm.GetMPIComm()); layout_comm_.SetOutMessageLayout(comm_); diff --git a/src/pcms/coupler/field_layout_communicator.h b/src/pcms/coupler/field_layout_communicator.h index bc0179f3..68d3572b 100644 --- a/src/pcms/coupler/field_layout_communicator.h +++ b/src/pcms/coupler/field_layout_communicator.h @@ -41,7 +41,19 @@ class FieldLayoutCommunicator template void SetOutMessageLayout(redev::BidirectionalComm& comm) { - comm.SetOutMessageLayout(plan_.dest_ranks, plan_.offsets); + // plan_.offsets is a per-holder CSR layout. Each holder carries + // num_components field values, so the field message's offsets are scaled + // by num_components while the destination ranks are unchanged. + const int num_comp = layout_.GetNumComponents(); + if (num_comp == 1) { + comm.SetOutMessageLayout(plan_.dest_ranks, plan_.offsets); + return; + } + redev::LOs scaled_offsets(plan_.offsets.size()); + for (size_t i = 0; i < plan_.offsets.size(); ++i) { + scaled_offsets[i] = plan_.offsets[i] * num_comp; + } + comm.SetOutMessageLayout(plan_.dest_ranks, scaled_offsets); } void UpdateLayout(); diff --git a/src/pcms/coupler/field_serializer.h b/src/pcms/coupler/field_serializer.h index 13d9ba6f..87246456 100644 --- a/src/pcms/coupler/field_serializer.h +++ b/src/pcms/coupler/field_serializer.h @@ -20,21 +20,22 @@ class FieldSerializer { auto data = field.GetDOFHolderDataHost(); auto owned = layout.GetOwnedHost(); - // owned/permutation/buffer index the flat node-major wire order - // (dof * num_components + component); read the field values through the - // 2D [dof][comp] view. + // 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(data.extent(0)); const LO num_comp = static_cast(data.extent(1)); for (LO i = 0; i < num_dof; ++i) { - for (LO c = 0; c < num_comp; ++c) { - const LO flat = i * num_comp + c; - if (owned[flat]) - buffer[permutation[flat]] = data(i, c); + if (owned[i]) { + for (LO c = 0; c < num_comp; ++c) { + buffer[permutation[i] * num_comp + c] = data(i, c); + } } } } - return static_cast(data.size()); + return static_cast(buffer.size()); } virtual void Deserialize( @@ -42,14 +43,19 @@ class FieldSerializer Rank1View buffer, Rank1View permutation) const { - Kokkos::View sorted("sorted", permutation.size()); + const LO num_dof = layout.GetNumOwnedDofHolder(); + const LO num_comp = layout.GetNumComponents(); + Kokkos::View sorted("sorted", layout.OwnedSize()); auto owned = layout.GetOwnedHost(); - for (LO i = 0; i < static_cast(sorted.size()); ++i) { - if (owned[i]) - sorted[i] = buffer[permutation[i]]; + for (LO i = 0; i < num_dof; ++i) { + if (owned[i]) { + for (LO c = 0; c < num_comp; ++c) { + sorted[i * num_comp + c] = buffer[permutation[i] * num_comp + c]; + } + } } - field.SetDOFHolderDataHost(Rank2View( - sorted.data(), layout.GetNumOwnedDofHolder(), layout.GetNumComponents())); + field.SetDOFHolderDataHost( + Rank2View(sorted.data(), num_dof, num_comp)); } virtual ~FieldSerializer() noexcept = default; diff --git a/src/pcms/coupler/serializer/xgc.h b/src/pcms/coupler/serializer/xgc.h index 1204f706..ce09cf74 100644 --- a/src/pcms/coupler/serializer/xgc.h +++ b/src/pcms/coupler/serializer/xgc.h @@ -35,14 +35,15 @@ class XGCFieldSerializer : public FieldSerializer 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(data.extent(0)); const LO num_comp = static_cast(data.extent(1)); for (LO i = 0; i < num_dof; ++i) { - for (LO c = 0; c < num_comp; ++c) { - const LO flat = i * num_comp + c; - if (owned[flat]) { - buffer[permutation[flat]] = data(i, c); + if (owned[i]) { + for (LO c = 0; c < num_comp; ++c) { + buffer[permutation[i] * num_comp + c] = data(i, c); } } } @@ -72,9 +73,11 @@ class XGCFieldSerializer : public FieldSerializer } } if (rank_participates_) { - for (LO i = 0; i < static_cast(current.size()); ++i) { + for (LO i = 0; i < num_dof; ++i) { if (owned[i]) { - full_data[i] = buffer[permutation[i]]; + for (LO c = 0; c < num_comp; ++c) { + full_data[i * num_comp + c] = buffer[permutation[i] * num_comp + c]; + } } } } diff --git a/src/pcms/discretization/discretization/omega_h.hpp b/src/pcms/discretization/discretization/omega_h.hpp index deeeab90..3fcb7eab 100644 --- a/src/pcms/discretization/discretization/omega_h.hpp +++ b/src/pcms/discretization/discretization/omega_h.hpp @@ -28,6 +28,9 @@ class OmegaHDiscretization : public Discretization Rank1View GetEntityClassificationIds(int entity_dim) const override; + // The Omega_h mesh this discretization is defined on. + Omega_h::Mesh& GetMesh() const noexcept { return mesh_; } + private: static constexpr int max_entity_dim_ = Region; diff --git a/src/pcms/field/evaluator/uniform_grid.h b/src/pcms/field/evaluator/uniform_grid.h index 979d1743..ed061458 100644 --- a/src/pcms/field/evaluator/uniform_grid.h +++ b/src/pcms/field/evaluator/uniform_grid.h @@ -186,7 +186,7 @@ class UniformGridPointEvaluator : public PointEvaluator private: std::shared_ptr> layout_; - const UniformGrid& grid_; + UniformGrid grid_; UniformGridFieldLocalizationHint hint_; Real fill_value_; }; @@ -304,7 +304,7 @@ class UniformGridEvaluatorFactory : public FieldEvaluatorFactory private: std::shared_ptr> layout_; - const UniformGrid& grid_; + UniformGrid grid_; }; using UniformGridEvaluatorFactory2D = UniformGridEvaluatorFactory<2>; diff --git a/src/pcms/field/evaluator/uniform_grid_spline.h b/src/pcms/field/evaluator/uniform_grid_spline.h index 9953d028..3fb693ad 100644 --- a/src/pcms/field/evaluator/uniform_grid_spline.h +++ b/src/pcms/field/evaluator/uniform_grid_spline.h @@ -171,7 +171,7 @@ class UniformGridSplinePointEvaluator2D private: std::shared_ptr> layout_; - const UniformGrid<2>& grid_; + UniformGrid<2> grid_; UniformGridFieldLocalizationHint<2> hint_; Real fill_value_; Kokkos::View x_coords_; @@ -287,7 +287,7 @@ class UniformGridSplineEvaluatorFactory2D : public FieldEvaluatorFactory private: std::shared_ptr> layout_; - const UniformGrid<2>& grid_; + UniformGrid<2> grid_; }; } // namespace pcms diff --git a/src/pcms/field/layout/omega_h_lagrange.cpp b/src/pcms/field/layout/omega_h_lagrange.cpp index 93ebd0bd..a2749a6e 100644 --- a/src/pcms/field/layout/omega_h_lagrange.cpp +++ b/src/pcms/field/layout/omega_h_lagrange.cpp @@ -3,6 +3,7 @@ #include "pcms/utility/mesh_geometry.h" #include "pcms/utility/omega_h_array_utils.h" #include "pcms/utility/profile.h" +#include #include namespace pcms @@ -48,6 +49,35 @@ Omega_h::Write BuildGids(Omega_h::Mesh& mesh, int entity_dim, return gids; } +// Some meshes (e.g. certain XGC meshes) are stored without any geometric +// classification (no class_id/class_dim tags). The layout and its +// discretization read these tags on every mesh dimension, so derive a default +// classification when it is missing. classify_elements + +// finalize_classification populate a geometrically meaningful class_dim on all +// dimensions but do not create class_id (that normally comes from a geometric +// model we don't have here), so fill any missing class_id with 0. +// Classification is only consulted for coupling DOF matching +// (field_exchange_planner); standalone field transfer never reads it, so a +// default id is sufficient to keep such meshes usable. +void EnsureClassification(Omega_h::Mesh& mesh) +{ + if (mesh.has_tag(mesh.dim(), "class_id") && + mesh.has_tag(mesh.dim(), "class_dim")) { + return; + } + if (!mesh.has_tag(mesh.dim(), "class_dim")) { + Omega_h::classify_elements(&mesh); + Omega_h::finalize_classification(&mesh); + } + for (int d = 0; d <= mesh.dim(); ++d) { + if (!mesh.has_tag(d, "class_id")) { + mesh.add_tag( + d, "class_id", 1, + Omega_h::Read(mesh.nents(d), 0, "class_id")); + } + } +} + Kokkos::View BuildOwned(Omega_h::Mesh& mesh, int entity_dim) { @@ -87,6 +117,7 @@ OmegaHLagrangeLayout::OmegaHLagrangeLayout(Omega_h::Mesh& mesh, int order, PCMS_FUNCTION_TIMER; int entity_dim = EntityDimForOrder(order_, mesh_.dim()); + EnsureClassification(mesh_); gids_ = BuildGids(mesh_, entity_dim, global_id_name_); gids_host_ = Omega_h::HostWrite(gids_); coords_ = get_entity_centroids(mesh_, entity_dim); @@ -130,6 +161,7 @@ OmegaHLagrangeLayout::OmegaHLagrangeLayout( PCMS_FUNCTION_TIMER; int entity_dim = EntityDimForOrder(order_, mesh_.dim()); + EnsureClassification(mesh_); gids_ = BuildGids(mesh_, entity_dim, global_id_name_); gids_host_ = Omega_h::HostWrite(gids_); coords_ = get_entity_centroids(mesh_, entity_dim); diff --git a/src/pcms/localization/adj_search.cpp b/src/pcms/localization/adj_search.cpp index e1b937d7..6f88a2ce 100644 --- a/src/pcms/localization/adj_search.cpp +++ b/src/pcms/localization/adj_search.cpp @@ -160,7 +160,10 @@ OMEGA_H_INLINE void add_support_if_unvisited_and_within_cutoff( if (!visited.notVisited(candidate_id)) return; - visited.push_back(candidate_id); + // If the visited buffer is full, stop here so the BFS terminates. Without + // this, unrecorded vertices are treated as unvisited and re-queued forever. + if (!visited.push_back(candidate_id)) + return; Omega_h::Real candidate_coords[max_dim]; for (Omega_h::LO k = 0; k < dim; ++k) { diff --git a/src/pcms/localization/queue_visited.hpp b/src/pcms/localization/queue_visited.hpp index 421281c5..8d768eec 100644 --- a/src/pcms/localization/queue_visited.hpp +++ b/src/pcms/localization/queue_visited.hpp @@ -56,7 +56,7 @@ class Track ~Track() {} OMEGA_H_INLINE - void push_back(const int& item); + bool push_back(const int& item); OMEGA_H_INLINE int size(); @@ -107,15 +107,18 @@ bool Queue::isFull() const } OMEGA_H_INLINE -void Track::push_back(const int& item) +bool Track::push_back(const int& item) { if (count == MAX_SIZE_TRACK) { - printf("track is full %d\n", count); - return; + // 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; tracking_array[last] = item; count++; + return true; } OMEGA_H_INLINE diff --git a/test/test_omega_h_lagrange_field.cpp b/test/test_omega_h_lagrange_field.cpp index fa681b6d..259a0165 100644 --- a/test/test_omega_h_lagrange_field.cpp +++ b/test/test_omega_h_lagrange_field.cpp @@ -210,6 +210,30 @@ TEST_CASE("OmegaHLagrangeField order-1: serialize / deserialize round-trip") pcms::test::CheckSerializeDeserialize(*factory.GetLayout(), field.GetData()); } +TEST_CASE( + "OmegaHLagrangeField order-1: multi-component serialize / deserialize " + "round-trip") +{ + auto lib = Omega_h::Library{}; + auto mesh = MakeBox2D(lib.world()); + const int nc = 3; + // The MeshFields backend is scalar-only; multi-component needs OmegaH. + auto factory = pcms::LagrangeFunctionSpace::FromMesh( + mesh, 1, nc, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto field = factory.CreateField(pcms::FieldMetadata{}); + + const int n = factory.GetLayout()->GetNumOwnedDofHolder(); + std::vector data(static_cast(n) * nc); + for (int i = 0; i < n; ++i) + for (int c = 0; c < nc; ++c) + data[static_cast(i) * nc + c] = i + 0.25 * c; + field.GetData().SetDOFHolderDataHost( + pcms::Rank2View(data.data(), n, nc)); + + pcms::test::CheckSerializeDeserialize(*factory.GetLayout(), field.GetData()); +} + // ---- Order-0 field tests ---------------------------------------------------- TEST_CASE("OmegaHLagrangeField order-0: set/get DOF data round-trip")