From d0e32ad5a88f32861ba851c72ad6dde09637bf4f Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Thu, 27 Aug 2026 13:23:05 -0400 Subject: [PATCH 1/2] support meshfields multi-component backend --- src/pcms/field/evaluator/mesh_fields.h | 25 +++----- .../field/evaluator/mesh_fields_backend.h | 61 +++++++++++++++---- src/pcms/field/function_space/lagrange.cpp | 5 -- 3 files changed, 58 insertions(+), 33 deletions(-) diff --git a/src/pcms/field/evaluator/mesh_fields.h b/src/pcms/field/evaluator/mesh_fields.h index 277664ce..b3924a89 100644 --- a/src/pcms/field/evaluator/mesh_fields.h +++ b/src/pcms/field/evaluator/mesh_fields.h @@ -42,8 +42,6 @@ class MeshFieldsPointEvaluator : public PointEvaluator hint_.coordinates_.extent(0) + hint_.num_missing_); PCMS_ALWAYS_ASSERT(values.extent(1) == static_cast(layout_->GetNumComponents())); - // ensure that only scalar fields are supported - PCMS_ALWAYS_ASSERT(layout_->GetNumComponents() == 1); auto const* mesh_field_data = dynamic_cast*>(&field.GetData()); if (!mesh_field_data) { @@ -51,27 +49,26 @@ class MeshFieldsPointEvaluator : public PointEvaluator "MeshFieldsPointEvaluator::Evaluate: incompatible FieldData type"); } - // Use device views directly from hint (no copy needed) auto eval_results = mesh_field_data->GetMeshFieldBackend()->evaluate( hint_.coordinates_d_, hint_.offsets_d_); - // Scatter results directly on device (no host copy) + int ncomp = layout_->GetNumComponents(); + int num_eval = static_cast(eval_results.extent(0)); Kokkos::parallel_for( "CopyEvalResultsToValues", - Kokkos::RangePolicy( - 0, eval_results.extent(0)), - KOKKOS_CLASS_LAMBDA(LO i) { - values(hint_.indices_d_(i), 0) = eval_results(i, 0); + Kokkos::MDRangePolicy>({0, 0}, {num_eval, ncomp}), + KOKKOS_CLASS_LAMBDA(LO i, int c) { + values(hint_.indices_d_(i), c) = eval_results(i, c); }); if (hint_.num_missing_ > 0 && hint_.mode_ == OutOfBoundsMode::FILL) { T fill_val = static_cast(fill_value_); + int num_missing = static_cast(hint_.num_missing_); Kokkos::parallel_for( "FillMissingValues", - Kokkos::RangePolicy( - 0, hint_.num_missing_), - KOKKOS_CLASS_LAMBDA(LO i) { - values(hint_.missing_indices_d_(i), 0) = fill_val; + Kokkos::MDRangePolicy>({0, 0}, {num_missing, ncomp}), + KOKKOS_CLASS_LAMBDA(LO i, int c) { + values(hint_.missing_indices_d_(i), c) = fill_val; }); } } @@ -98,10 +95,6 @@ class MeshFieldsEvaluatorFactory : public FieldEvaluatorFactory if (mesh_.dim() == 3) { throw pcms_error("MeshFieldsEvaluatorFactory does not support 3D meshes"); } - if (layout_->GetNumComponents() != 1) { - throw pcms_error( - "MeshFieldsEvaluatorFactory only supports single-component fields"); - } } const FieldLayout& GetLayout() const override { return *layout_; } diff --git a/src/pcms/field/evaluator/mesh_fields_backend.h b/src/pcms/field/evaluator/mesh_fields_backend.h index 04f81c2a..e59025e8 100644 --- a/src/pcms/field/evaluator/mesh_fields_backend.h +++ b/src/pcms/field/evaluator/mesh_fields_backend.h @@ -27,8 +27,8 @@ class MeshFieldBackend { public: virtual ~MeshFieldBackend() = default; - virtual Kokkos::View evaluate(Kokkos::View localCoords, - Kokkos::View offsets) const = 0; + virtual Kokkos::View evaluate(Kokkos::View localCoords, + Kokkos::View offsets) const = 0; virtual void SetData(Rank1View data, size_t num_nodes, size_t num_components, int dim) = 0; virtual void GetData(Rank1View data, size_t num_nodes, @@ -38,21 +38,23 @@ class MeshFieldBackend // --------------------------------------------------------------------------- // Concrete backend implementation // --------------------------------------------------------------------------- -template +template class MeshFieldBackendImpl : public MeshFieldBackend { public: MeshFieldBackendImpl(Omega_h::Mesh& mesh) : mesh_(mesh), mesh_field_(mesh), - shape_field_(mesh_field_.template CreateLagrangeField()) + shape_field_( + mesh_field_.template CreateLagrangeField()) { } - Kokkos::View evaluate(Kokkos::View localCoords, - Kokkos::View offsets) const override + Kokkos::View evaluate(Kokkos::View localCoords, + Kokkos::View offsets) const override { - auto self = const_cast*>(this); + auto self = + const_cast*>(this); return self->mesh_field_.triangleLocalPointEval(localCoords, offsets, shape_field_.field); } @@ -92,13 +94,47 @@ class MeshFieldBackendImpl : public MeshFieldBackend private: Omega_h::Mesh& mesh_; MeshField::OmegahMeshField mesh_field_; - using FWC = decltype(mesh_field_.template CreateLagrangeField()); + using FWC = + decltype(mesh_field_ + .template CreateLagrangeField()); FWC shape_field_; // FWC = FieldWithController; keeps ctrlr alive }; // --------------------------------------------------------------------------- // Factory function: create a MeshFieldBackend from a layout // --------------------------------------------------------------------------- + +// Helper to dispatch on num_components at runtime for a fixed (Dim, Order). +template +std::shared_ptr> MakeBackendForComponents( + Omega_h::Mesh& mesh, int num_components) +{ + switch (num_components) { + case 1: + return std::make_shared>(mesh); + case 2: + return std::make_shared>(mesh); + case 3: + return std::make_shared>(mesh); + case 4: + return std::make_shared>(mesh); + case 5: + return std::make_shared>(mesh); + case 6: + return std::make_shared>(mesh); + case 7: + return std::make_shared>(mesh); + case 8: + return std::make_shared>(mesh); + case 9: + return std::make_shared>(mesh); + default: + throw pcms_error("MeshFieldBackend: num_components " + + std::to_string(num_components) + + " exceeds maximum supported (9)."); + } +} + template std::shared_ptr> MakeMeshFieldBackend( const MeshFieldsAdapterLayout& layout) @@ -115,18 +151,19 @@ std::shared_ptr> MakeMeshFieldBackend( throw pcms_error("MeshFieldBackend does not support 3D meshes"); } auto nodes_per_dim = layout.GetNodesPerDim(); + int num_components = layout.GetNumComponents(); if (nodes_per_dim[0] == 1 && nodes_per_dim[1] == 0 && nodes_per_dim[2] == 0 && nodes_per_dim[3] == 0) { switch (mesh.dim()) { - case 1: return std::make_shared>(mesh); - case 2: return std::make_shared>(mesh); + case 1: return MakeBackendForComponents(mesh, num_components); + case 2: return MakeBackendForComponents(mesh, num_components); default: break; } } else if (nodes_per_dim[0] == 1 && nodes_per_dim[1] == 1 && nodes_per_dim[2] == 0 && nodes_per_dim[3] == 0) { switch (mesh.dim()) { - case 2: return std::make_shared>(mesh); - case 3: return std::make_shared>(mesh); + case 2: return MakeBackendForComponents(mesh, num_components); + case 3: return MakeBackendForComponents(mesh, num_components); default: break; } } diff --git a/src/pcms/field/function_space/lagrange.cpp b/src/pcms/field/function_space/lagrange.cpp index 8d15a6ec..bd6498ec 100644 --- a/src/pcms/field/function_space/lagrange.cpp +++ b/src/pcms/field/function_space/lagrange.cpp @@ -75,11 +75,6 @@ std::shared_ptr LagrangeFunctionSpace::FromMesh( } if (backend == Backend::MeshFields) { #ifdef PCMS_ENABLE_MESHFIELDS - if (num_components != 1) { - throw pcms_error( - "LagrangeFunctionSpace::FromMesh: MeshFields backend only supports " - "single-component fields"); - } std::array nodes_per_dim{}; switch (order) { case 1: nodes_per_dim = {1, 0, 0, 0}; break; From 28ae1572283faee95121f132f9d15a8b0ac12d3b Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Thu, 27 Aug 2026 13:23:38 -0400 Subject: [PATCH 2/2] debug meshfields layout mismatch --- src/pcms/field/data/mesh_fields.h | 65 ++++++++++++------ test/test_point_evaluator.cpp | 108 ++++++++++++++++++++++++++++-- 2 files changed, 147 insertions(+), 26 deletions(-) diff --git a/src/pcms/field/data/mesh_fields.h b/src/pcms/field/data/mesh_fields.h index 4c27a49d..1760e9aa 100644 --- a/src/pcms/field/data/mesh_fields.h +++ b/src/pcms/field/data/mesh_fields.h @@ -14,6 +14,32 @@ namespace pcms { +namespace +{ +// Functor that flattens component-major 2D data (LayoutLeft) into dof-major +// order. +template +struct FlattenFunctor +{ + Rank2View data_; + Kokkos::View flat_; + LO row_off_; + size_t nc_; + FlattenFunctor(Rank2View d, + Kokkos::View f, LO ro, size_t nc) + : data_(d), flat_(f), row_off_(ro), nc_(nc) + { + } + KOKKOS_INLINE_FUNCTION void operator()(LO local) const + { + LO global_dof = row_off_ + local; + for (size_t c = 0; c < nc_; ++c) { + flat_(local * nc_ + c) = data_(global_dof, c); + } + } +}; +} // namespace + template class MeshFieldsFieldData : public FieldData { @@ -24,9 +50,11 @@ class MeshFieldsFieldData : public FieldData metadata_(metadata), mesh_field_(MakeMeshFieldBackend(*layout_)), host_data_("meshfields_field_data", - static_cast(layout_->OwnedSize())), + static_cast(layout_->GetNumOwnedDofHolder()), + static_cast(layout_->GetNumComponents())), device_data_("meshfields_field_data_device", - static_cast(layout_->OwnedSize())) + static_cast(layout_->GetNumOwnedDofHolder()), + static_cast(layout_->GetNumComponents())) { if (!mesh_field_) { throw pcms_error( @@ -38,10 +66,8 @@ class MeshFieldsFieldData : public FieldData Rank2View GetDOFHolderDataHost() const override { - Kokkos::deep_copy(host_data_, device_data_); - return Rank2View(host_data_.data(), - layout_->GetNumOwnedDofHolder(), - layout_->GetNumComponents()); + DeepCopyMismatchLayouts(host_data_, device_data_); + return MakeConstRank2View(host_data_); } void SetDOFHolderDataHost(Rank2View values) override @@ -54,12 +80,7 @@ class MeshFieldsFieldData : public FieldData Rank2View GetDOFHolderData() const override { - // The Rank2View will wrap the dof-major data with layout left when device - // memory is enabled. This may cause issues in multi component cases. See - // issue #342 - return Rank2View( - device_data_.data(), layout_->GetNumOwnedDofHolder(), - layout_->GetNumComponents()); + return MakeConstRank2View(device_data_); } void SetDOFHolderData(Rank2View values) override @@ -81,18 +102,22 @@ class MeshFieldsFieldData : public FieldData auto nodes_per_dim = layout_->GetNodesPerDim(); auto num_components = layout_->GetNumComponents(); auto& mesh = layout_->GetMesh(); - // data is [dof_holder][component], contiguous node-major, so each mesh - // dimension owns a contiguous block of rows; SetData consumes a flat - // node-major span over that block. + // device_data_ is rank-2 LayoutLeft (component-major), but SetData + // expects a flat dof-major span. size_t row_offset = 0; for (int i = 0; i <= mesh.dim(); ++i) { if (nodes_per_dim[i]) { size_t num_rows = static_cast(mesh.nents(i)) * static_cast(nodes_per_dim[i]); size_t len = num_rows * static_cast(num_components); - Rank1View subspan{ - data.data_handle() + row_offset * static_cast(num_components), - len}; + Kokkos::View flat("sync_flat", len); + Kokkos::parallel_for( + "SyncBackendReorder", + Kokkos::RangePolicy( + 0, static_cast(num_rows)), + FlattenFunctor(data, flat, static_cast(row_offset), + num_components)); + Rank1View subspan(flat.data(), len); mesh_field_->SetData(subspan, nodes_per_dim[i], num_components, i); row_offset += num_rows; } @@ -102,8 +127,8 @@ class MeshFieldsFieldData : public FieldData std::shared_ptr layout_; FieldMetadata metadata_; std::shared_ptr> mesh_field_; - mutable Kokkos::View host_data_; - Kokkos::View device_data_; + mutable Kokkos::View host_data_; + Kokkos::View device_data_; }; } // namespace pcms diff --git a/test/test_point_evaluator.cpp b/test/test_point_evaluator.cpp index e9e235a9..36807bdf 100644 --- a/test/test_point_evaluator.cpp +++ b/test/test_point_evaluator.cpp @@ -361,15 +361,111 @@ TEST_CASE( *evaluator, field_b, pts, OMEGA_H_LAMBDA(Real, Real) { return Real(42); }); } -TEST_CASE("LagrangeFunctionSpace: MeshFields rejects multi-component fields") +TEST_CASE("PointEvaluator: MeshFields order-1 multi-component (2) evaluation") { auto lib = Omega_h::Library{}; auto mesh = - Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1, 1, 0, 10, 10, 0, false); + Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1, 1, 0, 50, 50, 0, false); + + auto factory = pcms::LagrangeFunctionSpace::FromMesh( + mesh, 1, 2, CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::MeshFields); - REQUIRE_THROWS_AS(pcms::LagrangeFunctionSpace::FromMesh( - mesh, 1, 2, CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::MeshFields), - pcms::pcms_error); + auto field = factory->CreateFunction(); + auto layout = factory->GetLayout(); + REQUIRE(layout->GetNumComponents() == 2); + int num_dof = layout->GetNumOwnedDofHolder(); + + // Component 0: f(x,y)=x+y, Component 1: f(x,y)=2x-y + auto dof_coords_mdspan = layout->GetDOFHolderCoordinates().GetValues(); + Kokkos::View dof_dev("dof_dev", num_dof, 2); + pcms::ConvertMismatchLayoutView2D(dof_dev, dof_coords_mdspan); + auto coords_host = + Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), dof_dev); + + std::vector host_data(static_cast(num_dof * 2)); + for (int i = 0; i < num_dof; ++i) { + Real x = coords_host(i, 0); + Real y = coords_host(i, 1); + host_data[static_cast(i) * 2 + 0] = x + y; + host_data[static_cast(i) * 2 + 1] = 2.0 * x - y; + } + field.SetDOFHolderDataHost(pcms::Rank2View( + host_data.data(), num_dof, 2)); + + auto pts = pcms::test::StandardEvalCoords2D(); + int n = static_cast(pts.size()) / 2; + auto device_coords = + pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); + auto evaluator = factory->CreatePointEvaluator( + pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); + + Kokkos::View out("out", n, 2); + evaluator->Evaluate(field, pcms::MakeRank2View(out)); + auto out_host = + Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out); + + for (int i = 0; i < n; ++i) { + Real x = pts[2 * static_cast(i)]; + Real y = pts[2 * static_cast(i) + 1]; + INFO("Point " << i << " (" << x << ", " << y << ")"); + REQUIRE(out_host(i, 0) == Catch::Approx(x + y).margin(1e-8)); + REQUIRE(out_host(i, 1) == Catch::Approx(2.0 * x - y).margin(1e-8)); + } +} + +TEST_CASE("PointEvaluator: MeshFields order-1 multi-component (3) evaluation") +{ + auto lib = Omega_h::Library{}; + auto mesh = + Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1, 1, 0, 50, 50, 0, false); + + auto factory = pcms::LagrangeFunctionSpace::FromMesh( + mesh, 1, 3, CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::MeshFields); + + auto field = factory->CreateFunction(); + auto layout = factory->GetLayout(); + REQUIRE(layout->GetNumComponents() == 3); + int num_dof = layout->GetNumOwnedDofHolder(); + + // Component 0: x, Component 1: y, Component 2: x*y + auto dof_coords_mdspan = layout->GetDOFHolderCoordinates().GetValues(); + Kokkos::View dof_dev("dof_dev", num_dof, 2); + pcms::ConvertMismatchLayoutView2D(dof_dev, dof_coords_mdspan); + auto coords_host = + Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), dof_dev); + + std::vector host_data(static_cast(num_dof * 3)); + for (int i = 0; i < num_dof; ++i) { + Real x = coords_host(i, 0); + Real y = coords_host(i, 1); + host_data[static_cast(i) * 3 + 0] = x; + host_data[static_cast(i) * 3 + 1] = y; + host_data[static_cast(i) * 3 + 2] = x * y; + } + field.SetDOFHolderDataHost(pcms::Rank2View( + host_data.data(), num_dof, 3)); + + auto pts = pcms::test::StandardEvalCoords2D(); + int n = static_cast(pts.size()) / 2; + auto device_coords = + pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); + auto evaluator = factory->CreatePointEvaluator( + pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); + + Kokkos::View out("out", n, 3); + evaluator->Evaluate(field, pcms::MakeRank2View(out)); + auto out_host = + Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out); + + for (int i = 0; i < n; ++i) { + Real x = pts[2 * static_cast(i)]; + Real y = pts[2 * static_cast(i) + 1]; + INFO("Point " << i << " (" << x << ", " << y << ")"); + REQUIRE(out_host(i, 0) == Catch::Approx(x).margin(1e-8)); + REQUIRE(out_host(i, 1) == Catch::Approx(y).margin(1e-8)); + REQUIRE(out_host(i, 2) == Catch::Approx(x * y).margin(1e-8)); + } } #endif // PCMS_ENABLE_MESHFIELDS