Skip to content

Support multi component uniform grid interpolator #376

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Sichao25 wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

Support multi component uniform grid interpolator #376

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
25 changes: 11 additions & 14 deletions src/pcms/field/data/simple.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,9 +14,8 @@
namespace pcms
{

// SimpleFieldData<T> is a generic concrete FieldData<T> backed by a flat
// Kokkos::View<T*, HostMemorySpace>. It works for any backend whose DOF data
// is a flat coefficient array (OmegaH, UniformGrid, PointCloud, etc.).
// SimpleFieldData<T> is a generic concrete FieldData<T> 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.
Expand All @@ -29,20 +28,20 @@ class SimpleFieldData : public FieldData<T>
: layout_(std::move(layout)),
metadata_(metadata),
host_data_("simple_field_data",
static_cast<size_t>(layout_->OwnedSize())),
static_cast<size_t>(layout_->GetNumOwnedDofHolder()),
static_cast<size_t>(layout_->GetNumComponents())),
device_data_("simple_field_data_device",
static_cast<size_t>(layout_->OwnedSize()))
static_cast<size_t>(layout_->GetNumOwnedDofHolder()),
static_cast<size_t>(layout_->GetNumComponents()))
{
}

const FieldMetadata& GetMetadata() const override { return metadata_; }

Rank2View<const T, HostMemorySpace> GetDOFHolderDataHost() const override
{
Kokkos::deep_copy(host_data_, device_data_);
return Rank2View<const T, HostMemorySpace>(host_data_.data(),
layout_->GetNumOwnedDofHolder(),
layout_->GetNumComponents());
DeepCopyMismatchLayouts(host_data_, device_data_);
return MakeConstRank2View(host_data_);
}

void SetDOFHolderDataHost(Rank2View<const T, HostMemorySpace> values) override
Expand All @@ -54,9 +53,7 @@ class SimpleFieldData : public FieldData<T>

Rank2View<const T, DeviceMemorySpace> GetDOFHolderData() const override
{
return Rank2View<const T, DeviceMemorySpace>(
device_data_.data(), layout_->GetNumOwnedDofHolder(),
layout_->GetNumComponents());
return MakeConstRank2View(device_data_);
}

void SetDOFHolderData(Rank2View<const T, DeviceMemorySpace> values) override
Expand All @@ -69,8 +66,8 @@ class SimpleFieldData : public FieldData<T>
private:
std::shared_ptr<const FieldLayout> layout_;
FieldMetadata metadata_;
mutable Kokkos::View<T*, HostMemorySpace> host_data_;
Kokkos::View<T*, DeviceMemorySpace> device_data_;
mutable Kokkos::View<T**, HostMemorySpace> host_data_;
Kokkos::View<T**, DeviceMemorySpace> device_data_;
};

} // namespace pcms
Expand Down
40 changes: 19 additions & 21 deletions src/pcms/field/evaluator/uniform_grid.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,9 @@ class UniformGridPointEvaluator : public PointEvaluator<Real, LayoutPolicy>
"set_dimensions_view",
Kokkos::RangePolicy<DeviceMemorySpace::execution_space>(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);
Expand All @@ -140,7 +142,10 @@ class UniformGridPointEvaluator : public PointEvaluator<Real, LayoutPolicy>
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);
}
});

Expand All @@ -153,18 +158,16 @@ class UniformGridPointEvaluator : public PointEvaluator<Real, LayoutPolicy>
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<size_t>(dof_data.size()));
RealMatView values_interp(
"values_interp", static_cast<size_t>(layout_->GetNumOwnedDofHolder()),
static_cast<size_t>(n_comp));
Kokkos::parallel_for(
"copy_dof_data_to_values_interp",
Kokkos::RangePolicy<DeviceMemorySpace::execution_space>(
0, static_cast<LO>(dof_data.size())),
KOKKOS_LAMBDA(const LO i) { values_interp(i) = dof_data(i, 0); });
Kokkos::MDRangePolicy<Kokkos::Rank<2>>(
{0, 0}, {static_cast<LO>(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);
Expand All @@ -174,12 +177,12 @@ class UniformGridPointEvaluator : public PointEvaluator<Real, LayoutPolicy>
Real fv = fill_value_;
Kokkos::parallel_for(
"evaluate_order1_device",
Kokkos::RangePolicy<DeviceMemorySpace::execution_space>(0, num_points),
KOKKOS_LAMBDA(const LO i) {
Kokkos::MDRangePolicy<Kokkos::Rank<2>>({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);
}
});
}
Expand Down Expand Up @@ -235,11 +238,6 @@ class UniformGridEvaluatorFactory : public FieldEvaluatorFactory<Real>
"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<LO>(coordinates.extent(0));
Expand Down
35 changes: 23 additions & 12 deletions src/pcms/transfer/linear_interpolant.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand Down
139 changes: 139 additions & 0 deletions test/test_uniform_grid_field.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<pcms::UniformGridFieldLayout<2>>(
grid, num_components, pcms::CoordinateSystem::Cartesian);
auto field_space = pcms::LagrangeFunctionSpace::FromUniformGrid(
grid, num_components, pcms::CoordinateSystem::Cartesian);
auto field = field_space->CreateFunction<pcms::Real>();
pcms::UniformGridEvaluatorFactory<2> eval_factory(layout);

// Vertex layout:
// v6---v7---v8
// | 2 | 3 |
// v3---v4---v5
// | 0 | 1 |
// v0---v1---v2
std::vector<pcms::Real> 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<const pcms::Real, pcms::HostMemorySpace>(
data.data(), num_vertices, num_components));

// Mix of cell centers and off-center points.
const int num_points = 4;
std::vector<pcms::Real> 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<pcms::Real**, pcms::DeviceMemorySpace> 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<pcms::UniformGridFieldLayout<2>>(
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::Real>();
pcms::UniformGridEvaluatorFactory<2> eval_factory(layout);

// Set cell values (4 cells, 2 components each)
std::vector<pcms::Real> 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<const pcms::Real, pcms::HostMemorySpace>(data.data(), 4,
num_components));

// Evaluate at cell centers
std::vector<pcms::Real> 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<pcms::Real**, pcms::DeviceMemorySpace> 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);
}
Loading