Skip to content

Put applyExtrapolated Restriction and Prolongation on GPU #286

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
EmilyBourne wants to merge 7 commits into
base: main
Choose a base branch
from
Open

Put applyExtrapolated Restriction and Prolongation on GPU #286

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
7f1f500
Put applyExtrapolatedRestriction and applyExtrapolatedProlongation on…
EmilyBourne May 28, 2026
dc04e5a
Clang format
EmilyBourne May 29, 2026
d4ebff8
Const doesn't need to be copied
EmilyBourne May 29, 2026
42e054b
Don't copy const
EmilyBourne May 29, 2026
5aea55d
Correct MemorySpace
EmilyBourne May 29, 2026
4031e7d
Fix memory spaces
EmilyBourne May 29, 2026
8e8dc27
Clang format
EmilyBourne May 29, 2026
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
26 changes: 17 additions & 9 deletions include/GMGPolar/utils.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,31 +43,39 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::injection(int current
}

template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolatedProlongation(int current_level,
HostVector<double> result,
HostConstVector<double> x) const
void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolatedProlongation(
int current_level, HostVector<double> result_host, HostConstVector<double> x_host) const
{
auto result = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), result_host);
auto x = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), x_host);

assert(current_level < number_of_levels_ && 1 <= current_level);
if (!interpolation_)
throw std::runtime_error("Interpolation not initialized.");

PolarGrid<Kokkos::HostSpace> current_grid(levels_[current_level].grid());
PolarGrid<Kokkos::HostSpace> previous_grid(levels_[current_level - 1].grid());
PolarGrid<DefaultMemorySpace> current_grid = levels_[current_level].grid();
PolarGrid<DefaultMemorySpace> previous_grid = levels_[current_level - 1].grid();
interpolation_->applyExtrapolatedProlongation(current_grid, previous_grid, result, x);

Kokkos::deep_copy(result_host, result);
}

template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolatedRestriction(int current_level,
HostVector<double> result,
HostConstVector<double> x) const
HostVector<double> result_host,
HostConstVector<double> x_host) const
{
auto result = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), result_host);
auto x = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), x_host);
assert(current_level < number_of_levels_ - 1 && 0 <= current_level);
if (!interpolation_)
throw std::runtime_error("Interpolation not initialized.");

PolarGrid<Kokkos::HostSpace> current_grid(levels_[current_level].grid());
PolarGrid<Kokkos::HostSpace> next_grid(levels_[current_level + 1].grid());
PolarGrid<DefaultMemorySpace> current_grid = levels_[current_level].grid();
PolarGrid<DefaultMemorySpace> next_grid = levels_[current_level + 1].grid();
interpolation_->applyExtrapolatedRestriction(current_grid, next_grid, result, x);

Kokkos::deep_copy(result_host, result);
}

template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
Expand Down
12 changes: 6 additions & 6 deletions include/Interpolation/interpolation.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,18 +27,18 @@ class Interpolation
const PolarGrid<Kokkos::HostSpace>& fine_grid, HostVector<double> fine_result,
HostConstVector<double> coarse_values) const;

void applyExtrapolatedProlongation(const PolarGrid<Kokkos::HostSpace>& coarse_grid,
const PolarGrid<Kokkos::HostSpace>& fine_grid, HostVector<double> fine_result,
HostConstVector<double> coarse_values) const;
void applyExtrapolatedProlongation(const PolarGrid<DefaultMemorySpace>& coarse_grid,
const PolarGrid<DefaultMemorySpace>& fine_grid, Vector<double> fine_result,
ConstVector<double> coarse_values) const;

/* Scaled full weighting (FW) restriction operator. */
void applyRestriction(const PolarGrid<Kokkos::HostSpace>& fine_grid,
const PolarGrid<Kokkos::HostSpace>& coarse_grid, HostVector<double> coarse_result,
HostConstVector<double> fine_values) const;

void applyExtrapolatedRestriction(const PolarGrid<Kokkos::HostSpace>& fine_grid,
const PolarGrid<Kokkos::HostSpace>& coarse_grid, HostVector<double> coarse_result,
HostConstVector<double> fine_values) const;
void applyExtrapolatedRestriction(const PolarGrid<DefaultMemorySpace>& fine_grid,
const PolarGrid<DefaultMemorySpace>& coarse_grid, Vector<double> coarse_result,
ConstVector<double> fine_values) const;

/* Bicubic FMG interpolator 1/16 * [-1, 9, 9, -1] */
void applyFMGInterpolation(const PolarGrid<Kokkos::HostSpace>& coarse_grid,
Expand Down
19 changes: 9 additions & 10 deletions src/Interpolation/extrapolated_prolongation.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,10 +48,10 @@ using namespace gmgpolar;
*/

static KOKKOS_INLINE_FUNCTION void fineNodeExtrapolatedProlongation(const int i_r, const int i_theta,
const PolarGrid<Kokkos::HostSpace>& coarse_grid,
const PolarGrid<Kokkos::HostSpace>& fine_grid,
HostVector<double>& fine_result,
HostConstVector<double>& coarse_values)
const PolarGrid<DefaultMemorySpace>& coarse_grid,
const PolarGrid<DefaultMemorySpace>& fine_grid,
Vector<double>& fine_result,
ConstVector<double>& coarse_values)
{
const int i_r_coarse = i_r / 2;
const int i_theta_coarse = i_theta / 2;
Expand Down Expand Up @@ -85,10 +85,9 @@ static KOKKOS_INLINE_FUNCTION void fineNodeExtrapolatedProlongation(const int i_
}
}

void Interpolation::applyExtrapolatedProlongation(const PolarGrid<Kokkos::HostSpace>& coarse_grid,
const PolarGrid<Kokkos::HostSpace>& fine_grid,
HostVector<double> fine_result,
HostConstVector<double> coarse_values) const
void Interpolation::applyExtrapolatedProlongation(const PolarGrid<DefaultMemorySpace>& coarse_grid,
const PolarGrid<DefaultMemorySpace>& fine_grid,
Vector<double> fine_result, ConstVector<double> coarse_values) const
{
assert(std::ssize(coarse_values) == coarse_grid.numberOfNodes());
assert(std::ssize(fine_result) == fine_grid.numberOfNodes());
Expand All @@ -99,7 +98,7 @@ void Interpolation::applyExtrapolatedProlongation(const PolarGrid<Kokkos::HostSp
// The For loop matches circular access pattern */
Kokkos::parallel_for(
"Interpolation: Extrapolated Prolongation (Circular)",
Kokkos::MDRangePolicy<Kokkos::DefaultHostExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
{0, 0}, // Starting point of the index space
{fine_grid.numberSmootherCircles(), fine_grid.ntheta()} // Ending point of the index space
),
Expand All @@ -111,7 +110,7 @@ void Interpolation::applyExtrapolatedProlongation(const PolarGrid<Kokkos::HostSp
/* For loop matches radial access pattern */
Kokkos::parallel_for(
"Interpolation: Extrapolated Prolongation (Radial)",
Kokkos::MDRangePolicy<Kokkos::DefaultHostExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
{0, fine_grid.numberSmootherCircles()}, // Starting point of the index space
{fine_grid.ntheta(), fine_grid.nr()} // Ending point of the index space
),
Expand Down
19 changes: 9 additions & 10 deletions src/Interpolation/extrapolated_restriction.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,10 +23,10 @@ using namespace gmgpolar;
*/

static KOKKOS_INLINE_FUNCTION void coarseNodeExtrapolatedRestriction(const int i_r_coarse, const int i_theta_coarse,
const PolarGrid<Kokkos::HostSpace>& fine_grid,
const PolarGrid<Kokkos::HostSpace>& coarse_grid,
HostVector<double>& coarse_result,
HostConstVector<double>& fine_values)
const PolarGrid<DefaultMemorySpace>& fine_grid,
const PolarGrid<DefaultMemorySpace>& coarse_grid,
Vector<double>& coarse_result,
ConstVector<double>& fine_values)
{
const int i_r = i_r_coarse * 2;
const int i_theta = i_theta_coarse * 2;
Expand All @@ -50,10 +50,9 @@ static KOKKOS_INLINE_FUNCTION void coarseNodeExtrapolatedRestriction(const int i
coarse_result[coarse_grid.index(i_r_coarse, i_theta_coarse)] = value;
}

void Interpolation::applyExtrapolatedRestriction(const PolarGrid<Kokkos::HostSpace>& fine_grid,
const PolarGrid<Kokkos::HostSpace>& coarse_grid,
HostVector<double> coarse_result,
HostConstVector<double> fine_values) const
void Interpolation::applyExtrapolatedRestriction(const PolarGrid<DefaultMemorySpace>& fine_grid,
const PolarGrid<DefaultMemorySpace>& coarse_grid,
Vector<double> coarse_result, ConstVector<double> fine_values) const
{
assert(std::ssize(fine_values) == fine_grid.numberOfNodes());
assert(std::ssize(coarse_result) == coarse_grid.numberOfNodes());
Expand All @@ -64,7 +63,7 @@ void Interpolation::applyExtrapolatedRestriction(const PolarGrid<Kokkos::HostSpa
// The For loop matches circular access pattern */
Kokkos::parallel_for(
"Interpolation: Extrapolated Restriction (Circular)",
Kokkos::MDRangePolicy<Kokkos::DefaultHostExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
{0, 0}, // Starting point of the index space
{coarse_grid.numberSmootherCircles(), coarse_grid.ntheta()} // Ending point of the index space
),
Expand All @@ -77,7 +76,7 @@ void Interpolation::applyExtrapolatedRestriction(const PolarGrid<Kokkos::HostSpa
/* For loop matches radial access pattern */
Kokkos::parallel_for(
"Interpolation: Extrapolated Restriction (Radial)",
Kokkos::MDRangePolicy<Kokkos::DefaultHostExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
{0, coarse_grid.numberSmootherCircles()}, // Starting point of the index space
{coarse_grid.ntheta(), coarse_grid.nr()} // Ending point of the index space
),
Expand Down
19 changes: 13 additions & 6 deletions tests/Interpolation/extrapolated_prolongation.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,20 +47,27 @@ TEST(ExtrapolatedProlongationTest, ExtrapolatedProlongationMatchesStencil)
std::vector<double> fine_angles = {
0, M_PI / 16, M_PI / 8, M_PI / 2, M_PI, M_PI + M_PI / 16, M_PI + M_PI / 8, M_PI + M_PI / 2, 2 * M_PI};

PolarGrid<Kokkos::HostSpace> fine_grid(fine_radii, fine_angles);
PolarGrid<Kokkos::HostSpace> coarse_grid = coarseningGrid(fine_grid);
PolarGrid<DefaultMemorySpace> fine_grid(fine_radii, fine_angles);
PolarGrid<DefaultMemorySpace> coarse_grid = coarseningGrid(fine_grid);

PolarGrid<Kokkos::HostSpace> h_fine_grid(fine_grid);
PolarGrid<Kokkos::HostSpace> h_coarse_grid(coarse_grid);

Interpolation I(/*DirBC*/ true);

HostVector<double> coarse_values = generate_random_sample_data(coarse_grid, 1234, 0.0, 1.0);
HostVector<double> fine_result("fine_result", fine_grid.numberOfNodes());
HostVector<double> h_coarse_values = generate_random_sample_data(h_coarse_grid, 1234, 0.0, 1.0);
Vector<double> fine_result("fine_result", fine_grid.numberOfNodes());

auto coarse_values = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_coarse_values);

I.applyExtrapolatedProlongation(coarse_grid, fine_grid, fine_result, coarse_values);

auto h_fine_result = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), fine_result);

for (int i_r = 0; i_r < fine_grid.nr(); ++i_r) {
for (int i_theta = 0; i_theta < fine_grid.ntheta(); ++i_theta) {
double expected = expected_extrapolated_value(coarse_grid, fine_grid, coarse_values, i_r, i_theta);
double got = fine_result[fine_grid.index(i_r, i_theta)];
double expected = expected_extrapolated_value(h_coarse_grid, h_fine_grid, h_coarse_values, i_r, i_theta);
double got = h_fine_result[h_fine_grid.index(i_r, i_theta)];
ASSERT_NEAR(expected, got, 1e-10) << "Mismatch at (" << i_r << ", " << i_theta << ")";
}
}
Expand Down
21 changes: 14 additions & 7 deletions tests/Interpolation/extrapolated_restriction.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,21 +46,28 @@ TEST(ExtrapolatedRestrictionTest, ExtrapolatedRestrictionMatchesStencil)
std::vector<double> fine_angles = {
0, M_PI / 16, M_PI / 8, M_PI / 2, M_PI, M_PI + M_PI / 16, M_PI + M_PI / 8, M_PI + M_PI / 2, 2 * M_PI};

PolarGrid<Kokkos::HostSpace> fine_grid(fine_radii, fine_angles);
PolarGrid<Kokkos::HostSpace> coarse_grid = coarseningGrid(fine_grid);
PolarGrid<DefaultMemorySpace> fine_grid(fine_radii, fine_angles);
PolarGrid<DefaultMemorySpace> coarse_grid = coarseningGrid(fine_grid);

PolarGrid<Kokkos::HostSpace> h_fine_grid(fine_grid);
PolarGrid<Kokkos::HostSpace> h_coarse_grid(coarse_grid);

Interpolation I(/*DirBC*/ true);

HostVector<double> fine_values = generate_random_sample_data(fine_grid, 9012, 0.0, 1.0);
HostVector<double> coarse_result("coarse_result", coarse_grid.numberOfNodes());
HostVector<double> h_fine_values = generate_random_sample_data(h_fine_grid, 9012, 0.0, 1.0);
Vector<double> coarse_result("coarse_result", coarse_grid.numberOfNodes());

auto fine_values = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_fine_values);

I.applyExtrapolatedRestriction(fine_grid, coarse_grid, coarse_result, fine_values);

auto h_coarse_result = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), coarse_result);

for (int i_r_coarse = 0; i_r_coarse < coarse_grid.nr(); ++i_r_coarse) {
for (int i_theta_coarse = 0; i_theta_coarse < coarse_grid.ntheta(); ++i_theta_coarse) {
double expected = expected_extrapolated_restriction_value(fine_grid, coarse_grid, fine_values, i_r_coarse,
i_theta_coarse);
double got = coarse_result[coarse_grid.index(i_r_coarse, i_theta_coarse)];
double expected = expected_extrapolated_restriction_value(h_fine_grid, h_coarse_grid, h_fine_values,
i_r_coarse, i_theta_coarse);
double got = h_coarse_result[h_coarse_grid.index(i_r_coarse, i_theta_coarse)];
ASSERT_NEAR(expected, got, 1e-10) << "Mismatch at (" << i_r_coarse << ", " << i_theta_coarse << ")";
}
}
Expand Down
Loading