From 41d86536b66a5a489d4cabdcb05557bef00168c0 Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Wed, 2 Sep 2026 21:17:09 -0400 Subject: [PATCH] support multi component uniform grid interpolator --- src/pcms/field/data/simple.h | 25 ++-- src/pcms/field/evaluator/uniform_grid.h | 40 ++++--- src/pcms/transfer/linear_interpolant.hpp | 35 ++++-- test/test_uniform_grid_field.cpp | 139 +++++++++++++++++++++++ 4 files changed, 192 insertions(+), 47 deletions(-) diff --git a/src/pcms/field/data/simple.h b/src/pcms/field/data/simple.h index 6deda15a..93332245 100644 --- a/src/pcms/field/data/simple.h +++ b/src/pcms/field/data/simple.h @@ -14,9 +14,8 @@ namespace pcms { -// SimpleFieldData is a generic concrete FieldData backed by a flat -// Kokkos::View. It works for any backend whose DOF data -// is a flat coefficient array (OmegaH, UniformGrid, PointCloud, etc.). +// SimpleFieldData is a generic concrete FieldData for backends whose DOF +// can be stored in a simple 2D array format. // // Ownership of the layout is shared — the layout is typically held by the // factory that created this field data object. @@ -29,9 +28,11 @@ class SimpleFieldData : public FieldData : layout_(std::move(layout)), metadata_(metadata), host_data_("simple_field_data", - static_cast(layout_->OwnedSize())), + static_cast(layout_->GetNumOwnedDofHolder()), + static_cast(layout_->GetNumComponents())), device_data_("simple_field_data_device", - static_cast(layout_->OwnedSize())) + static_cast(layout_->GetNumOwnedDofHolder()), + static_cast(layout_->GetNumComponents())) { } @@ -39,10 +40,8 @@ class SimpleFieldData : 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,9 +53,7 @@ class SimpleFieldData : public FieldData Rank2View GetDOFHolderData() const override { - return Rank2View( - device_data_.data(), layout_->GetNumOwnedDofHolder(), - layout_->GetNumComponents()); + return MakeConstRank2View(device_data_); } void SetDOFHolderData(Rank2View values) override @@ -69,8 +66,8 @@ class SimpleFieldData : public FieldData private: std::shared_ptr layout_; FieldMetadata metadata_; - mutable Kokkos::View host_data_; - Kokkos::View device_data_; + mutable Kokkos::View host_data_; + Kokkos::View device_data_; }; } // namespace pcms diff --git a/src/pcms/field/evaluator/uniform_grid.h b/src/pcms/field/evaluator/uniform_grid.h index ed061458..320b4978 100644 --- a/src/pcms/field/evaluator/uniform_grid.h +++ b/src/pcms/field/evaluator/uniform_grid.h @@ -127,7 +127,9 @@ class UniformGridPointEvaluator : public PointEvaluator "set_dimensions_view", Kokkos::RangePolicy(0, Dim), KOKKOS_LAMBDA(const unsigned d) { - dimensions_view(d) = cell_divisions[d] + 1; + // cell_dim_indices uses GetDimensionedIndex(), which returns indices in + // the opposite order of the coordinates, so we reverse the order here. + dimensions_view(d) = cell_divisions[Dim - 1 - d] + 1; }); RealMatView parametric_coords("parametric_coords", num_points, Dim); @@ -140,7 +142,10 @@ class UniformGridPointEvaluator : public PointEvaluator Real coord = coordinates(i, d); Real cell_min = cell_bbox.center[d] - cell_bbox.half_width[d]; Real cell_max = cell_bbox.center[d] + cell_bbox.half_width[d]; - parametric_coords(i, d) = (coord - cell_min) / (cell_max - cell_min); + // cell_dim_indices is filled from GetDimensionedIndex(), which + // returns indices in the opposite order of the coordinates + parametric_coords(i, Dim - 1 - d) = + (coord - cell_min) / (cell_max - cell_min); } }); @@ -153,18 +158,16 @@ class UniformGridPointEvaluator : public PointEvaluator cell_indices_interp(i, d) = cell_dim_indices(i, d); }); - // n_comp == 1 for now (multi-component path would need per-component calls) - PCMS_ALWAYS_ASSERT( - n_comp == 1 && - "UniformGridPointEvaluator: multi-component order-1 not yet supported"); - - RealVecView values_interp("values_interp", - static_cast(dof_data.size())); + RealMatView values_interp( + "values_interp", static_cast(layout_->GetNumOwnedDofHolder()), + static_cast(n_comp)); Kokkos::parallel_for( "copy_dof_data_to_values_interp", - Kokkos::RangePolicy( - 0, static_cast(dof_data.size())), - KOKKOS_LAMBDA(const LO i) { values_interp(i) = dof_data(i, 0); }); + Kokkos::MDRangePolicy>( + {0, 0}, {static_cast(layout_->GetNumOwnedDofHolder()), n_comp}), + KOKKOS_LAMBDA(const LO i, const int c) { + values_interp(i, c) = dof_data(i, c); + }); auto interpolator = RegularGridInterpolator( parametric_coords, values_interp, cell_indices_interp, dimensions_view); @@ -174,12 +177,12 @@ class UniformGridPointEvaluator : public PointEvaluator Real fv = fill_value_; Kokkos::parallel_for( "evaluate_order1_device", - Kokkos::RangePolicy(0, num_points), - KOKKOS_LAMBDA(const LO i) { + Kokkos::MDRangePolicy>({0, 0}, {num_points, n_comp}), + KOKKOS_LAMBDA(const LO i, const int c) { if (is_out_of_bounds(i) && mode == OutOfBoundsMode::FILL) { - values(i, 0) = fv; + values(i, c) = fv; } else { - values(i, 0) = interpolated(i); + values(i, c) = interpolated(i, c); } }); } @@ -235,11 +238,6 @@ class UniformGridEvaluatorFactory : public FieldEvaluatorFactory "UniformGridEvaluatorFactory: nearest-boundary evaluation is not " "supported"); } - if (layout_->GetOrder() == 1 && layout_->GetNumComponents() != 1) { - throw pcms_error( - "UniformGridEvaluatorFactory: order-1 multi-component evaluation is " - "not implemented"); - } auto coordinates = coords.GetValues(); LO num_points = static_cast(coordinates.extent(0)); diff --git a/src/pcms/transfer/linear_interpolant.hpp b/src/pcms/transfer/linear_interpolant.hpp index c058d48f..bf0afe35 100644 --- a/src/pcms/transfer/linear_interpolant.hpp +++ b/src/pcms/transfer/linear_interpolant.hpp @@ -68,7 +68,8 @@ void basis_function(const RealVecView& parametric_coord, KOKKOS_INLINE_FUNCTION double linear_interpolant(const IntVecView& dimensions, const double* linear_basis_each_dir, - const IntVecView& indices, const RealVecView& values) + const IntVecView& indices, const RealMatView& values, + const int component) { int dim = dimensions.extent(0); double sum = 0; @@ -91,50 +92,60 @@ double linear_interpolant(const IntVecView& dimensions, } int idx = calculateIndex(dimensions, ids); - double corner_values = values(idx); + double corner_values = values(idx, component); sum += temp * corner_values; } return sum; } +// Multilinear interpolation over a regular grid. Supports one or more +// components per grid point; the scalar case is num_components == 1. +// `values` is [num_grid_points][num_components] and the result is +// [num_query_points][num_components]. class RegularGridInterpolator { private: const RealMatView parametric_coords; - const RealVecView values; + const RealMatView values; // 2D: [num_grid_points][num_components] const IntMatView indices; const IntVecView dimensions; + const int num_components; public: RegularGridInterpolator(const RealMatView& parametric_coords_, - const RealVecView& values_, + const RealMatView& values_, const IntMatView& indices_, const IntVecView& dimensions_) : parametric_coords(parametric_coords_), values(values_), indices(indices_), - dimensions(dimensions_){}; + dimensions(dimensions_), + num_components(values_.extent(1)) + { + } - RealVecView linear_interpolation() + RealMatView linear_interpolation() { - int dim = dimensions.extent(0); int N = parametric_coords.extent(0); - RealVecView interpolated_values("approximated values", N); + RealMatView interpolated_values("approximated values", N, num_components); auto parametric_coords_ = parametric_coords; auto dimensions_ = dimensions; auto values_ = values; auto indices_ = indices; + int n_comp = num_components; + Kokkos::parallel_for( - "linear interpolation function", N, KOKKOS_LAMBDA(int j) { + "linear interpolation", N, KOKKOS_LAMBDA(int j) { double linear_basis_each_dir[MAX_DIM] = {0.0}; auto parametric_coord_each_point = Kokkos::subview(parametric_coords_, j, Kokkos::ALL()); auto index_each_point = Kokkos::subview(indices_, j, Kokkos::ALL()); basis_function(parametric_coord_each_point, linear_basis_each_dir); - auto approx_value = linear_interpolant( - dimensions_, linear_basis_each_dir, index_each_point, values_); - interpolated_values(j) = approx_value; + for (int c = 0; c < n_comp; ++c) { + interpolated_values(j, c) = linear_interpolant( + dimensions_, linear_basis_each_dir, index_each_point, values_, c); + } }); return interpolated_values; diff --git a/test/test_uniform_grid_field.cpp b/test/test_uniform_grid_field.cpp index 699ec26a..eebef1b5 100644 --- a/test/test_uniform_grid_field.cpp +++ b/test/test_uniform_grid_field.cpp @@ -659,3 +659,142 @@ TEST_CASE("UniformGrid workflow") VerifyMaskFieldValues(grid, mask_field); } + +TEST_CASE("UniformGrid Order-1 multi-component field evaluation") +{ + // Create a 2x2 cell uniform grid + pcms::UniformGrid<2> grid; + grid.bot_left = {0.0, 0.0}; + grid.edge_length = {10.0, 10.0}; + grid.divisions = {2, 2}; + + // Components: x, y, x+y (separable) and x*y (non-separable). Off-center + // query points plus the non-separable component catch axis-mixing bugs that + // symmetric cell-center samples would hide. + const int num_components = 4; + + auto layout = std::make_shared>( + grid, num_components, pcms::CoordinateSystem::Cartesian); + auto field_space = pcms::LagrangeFunctionSpace::FromUniformGrid( + grid, num_components, pcms::CoordinateSystem::Cartesian); + auto field = field_space->CreateFunction(); + pcms::UniformGridEvaluatorFactory<2> eval_factory(layout); + + // Vertex layout: + // v6---v7---v8 + // | 2 | 3 | + // v3---v4---v5 + // | 0 | 1 | + // v0---v1---v2 + std::vector data; + for (int j = 0; j <= grid.divisions[1]; ++j) { + for (int i = 0; i <= grid.divisions[0]; ++i) { + pcms::Real x = + grid.bot_left[0] + i * (grid.edge_length[0] / grid.divisions[0]); + pcms::Real y = + grid.bot_left[1] + j * (grid.edge_length[1] / grid.divisions[1]); + data.push_back(x); // Component 0: x + data.push_back(y); // Component 1: y + data.push_back(x + y); // Component 2: x + y + data.push_back(x * y); // Component 3: x * y (non-separable) + } + } + + // Reshape data to [num_vertices][num_components] + pcms::LO num_vertices = (grid.divisions[0] + 1) * (grid.divisions[1] + 1); + field.SetDOFHolderDataHost( + pcms::Rank2View( + data.data(), num_vertices, num_components)); + + // Mix of cell centers and off-center points. + const int num_points = 4; + std::vector eval_coords = { + 2.5, 2.5, // cell 0 center + 1.0, 3.0, // cell 0 off-center + 6.0, 1.5, // cell 1 off-center + 7.5, 7.5 // cell 3 center + }; + auto device_coords = pcms::test::CreateDeviceCoordinateView( + eval_coords, pcms::CoordinateSystem::Cartesian); + auto evaluator = eval_factory.CreatePointEvaluator( + pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); + + Kokkos::View results_device( + "results_device", num_points, num_components); + evaluator->Evaluate(field, pcms::MakeRank2View(results_device)); + auto results_host = Kokkos::create_mirror_view_and_copy( + pcms::HostMemorySpace(), results_device); + + // Expected [x, y, x+y, x*y] at each query point. + const pcms::Real expected[num_points][num_components] = { + {2.5, 2.5, 5.0, 6.25}, // (2.5, 2.5) + {1.0, 3.0, 4.0, 3.0}, // (1.0, 3.0) + {6.0, 1.5, 7.5, 9.0}, // (6.0, 1.5) + {7.5, 7.5, 15.0, 56.25}, // (7.5, 7.5) + }; + + for (int p = 0; p < num_points; ++p) { + for (int c = 0; c < num_components; ++c) { + INFO("point " << p << " component " << c); + REQUIRE(std::abs(results_host(p, c) - expected[p][c]) < 1e-10); + } + } +} + +TEST_CASE("UniformGrid Order-0 multi-component field (regression)") +{ + // Verify Order-0 (piecewise constant) still works with multiple components + pcms::UniformGrid<2> grid; + grid.bot_left = {0.0, 0.0}; + grid.edge_length = {10.0, 10.0}; + grid.divisions = {2, 2}; + + const int num_components = 2; + + auto layout = std::make_shared>( + grid, num_components, pcms::CoordinateSystem::Cartesian, 0); // Order 0 + auto field_space = pcms::LagrangeFunctionSpace::FromUniformGrid( + grid, num_components, pcms::CoordinateSystem::Cartesian, 0); + auto field = field_space->CreateFunction(); + pcms::UniformGridEvaluatorFactory<2> eval_factory(layout); + + // Set cell values (4 cells, 2 components each) + std::vector data = { + 1.0, 10.0, // Cell 0 + 2.0, 20.0, // Cell 1 + 3.0, 30.0, // Cell 2 + 4.0, 40.0 // Cell 3 + }; + + field.SetDOFHolderDataHost( + pcms::Rank2View(data.data(), 4, + num_components)); + + // Evaluate at cell centers + std::vector eval_coords = { + 2.5, 2.5, // Cell 0 center + 7.5, 2.5, // Cell 1 center + 2.5, 7.5, // Cell 2 center + 7.5, 7.5 // Cell 3 center + }; + auto device_coords = pcms::test::CreateDeviceCoordinateView( + eval_coords, pcms::CoordinateSystem::Cartesian); + auto evaluator = eval_factory.CreatePointEvaluator( + pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); + + Kokkos::View results_device( + "results_device", 4, num_components); + evaluator->Evaluate(field, pcms::MakeRank2View(results_device)); + auto results_host = Kokkos::create_mirror_view_and_copy( + pcms::HostMemorySpace(), results_device); + + // For Order-0, values should match cell values exactly + REQUIRE(results_host(0, 0) == 1.0); + REQUIRE(results_host(0, 1) == 10.0); + REQUIRE(results_host(1, 0) == 2.0); + REQUIRE(results_host(1, 1) == 20.0); + REQUIRE(results_host(2, 0) == 3.0); + REQUIRE(results_host(2, 1) == 30.0); + REQUIRE(results_host(3, 0) == 4.0); + REQUIRE(results_host(3, 1) == 40.0); +}