diff --git a/include/GMGPolar/utils.h b/include/GMGPolar/utils.h index af7f123c..2804e51a 100644 --- a/include/GMGPolar/utils.h +++ b/include/GMGPolar/utils.h @@ -43,31 +43,39 @@ void GMGPolar::injection(int current } template -void GMGPolar::extrapolatedProlongation(int current_level, - HostVector result, - HostConstVector x) const +void GMGPolar::extrapolatedProlongation( + int current_level, HostVector result_host, HostConstVector 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 current_grid(levels_[current_level].grid()); - PolarGrid previous_grid(levels_[current_level - 1].grid()); + PolarGrid current_grid = levels_[current_level].grid(); + PolarGrid previous_grid = levels_[current_level - 1].grid(); interpolation_->applyExtrapolatedProlongation(current_grid, previous_grid, result, x); + + Kokkos::deep_copy(result_host, result); } template void GMGPolar::extrapolatedRestriction(int current_level, - HostVector result, - HostConstVector x) const + HostVector result_host, + HostConstVector 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 current_grid(levels_[current_level].grid()); - PolarGrid next_grid(levels_[current_level + 1].grid()); + PolarGrid current_grid = levels_[current_level].grid(); + PolarGrid next_grid = levels_[current_level + 1].grid(); interpolation_->applyExtrapolatedRestriction(current_grid, next_grid, result, x); + + Kokkos::deep_copy(result_host, result); } template diff --git a/include/Interpolation/interpolation.h b/include/Interpolation/interpolation.h index 00810375..7198458f 100644 --- a/include/Interpolation/interpolation.h +++ b/include/Interpolation/interpolation.h @@ -27,18 +27,18 @@ class Interpolation const PolarGrid& fine_grid, HostVector fine_result, HostConstVector coarse_values) const; - void applyExtrapolatedProlongation(const PolarGrid& coarse_grid, - const PolarGrid& fine_grid, HostVector fine_result, - HostConstVector coarse_values) const; + void applyExtrapolatedProlongation(const PolarGrid& coarse_grid, + const PolarGrid& fine_grid, Vector fine_result, + ConstVector coarse_values) const; /* Scaled full weighting (FW) restriction operator. */ void applyRestriction(const PolarGrid& fine_grid, const PolarGrid& coarse_grid, HostVector coarse_result, HostConstVector fine_values) const; - void applyExtrapolatedRestriction(const PolarGrid& fine_grid, - const PolarGrid& coarse_grid, HostVector coarse_result, - HostConstVector fine_values) const; + void applyExtrapolatedRestriction(const PolarGrid& fine_grid, + const PolarGrid& coarse_grid, Vector coarse_result, + ConstVector fine_values) const; /* Bicubic FMG interpolator 1/16 * [-1, 9, 9, -1] */ void applyFMGInterpolation(const PolarGrid& coarse_grid, diff --git a/src/Interpolation/extrapolated_prolongation.cpp b/src/Interpolation/extrapolated_prolongation.cpp index 9b001dee..22c75ee6 100644 --- a/src/Interpolation/extrapolated_prolongation.cpp +++ b/src/Interpolation/extrapolated_prolongation.cpp @@ -48,10 +48,10 @@ using namespace gmgpolar; */ static KOKKOS_INLINE_FUNCTION void fineNodeExtrapolatedProlongation(const int i_r, const int i_theta, - const PolarGrid& coarse_grid, - const PolarGrid& fine_grid, - HostVector& fine_result, - HostConstVector& coarse_values) + const PolarGrid& coarse_grid, + const PolarGrid& fine_grid, + Vector& fine_result, + ConstVector& coarse_values) { const int i_r_coarse = i_r / 2; const int i_theta_coarse = i_theta / 2; @@ -85,10 +85,9 @@ static KOKKOS_INLINE_FUNCTION void fineNodeExtrapolatedProlongation(const int i_ } } -void Interpolation::applyExtrapolatedProlongation(const PolarGrid& coarse_grid, - const PolarGrid& fine_grid, - HostVector fine_result, - HostConstVector coarse_values) const +void Interpolation::applyExtrapolatedProlongation(const PolarGrid& coarse_grid, + const PolarGrid& fine_grid, + Vector fine_result, ConstVector coarse_values) const { assert(std::ssize(coarse_values) == coarse_grid.numberOfNodes()); assert(std::ssize(fine_result) == fine_grid.numberOfNodes()); @@ -99,7 +98,7 @@ void Interpolation::applyExtrapolatedProlongation(const PolarGrid>( // Rank of the index space + Kokkos::MDRangePolicy>( // 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 ), @@ -111,7 +110,7 @@ void Interpolation::applyExtrapolatedProlongation(const PolarGrid>( // Rank of the index space + Kokkos::MDRangePolicy>( // 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 ), diff --git a/src/Interpolation/extrapolated_restriction.cpp b/src/Interpolation/extrapolated_restriction.cpp index 47e59b4f..b3f8b2de 100644 --- a/src/Interpolation/extrapolated_restriction.cpp +++ b/src/Interpolation/extrapolated_restriction.cpp @@ -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& fine_grid, - const PolarGrid& coarse_grid, - HostVector& coarse_result, - HostConstVector& fine_values) + const PolarGrid& fine_grid, + const PolarGrid& coarse_grid, + Vector& coarse_result, + ConstVector& fine_values) { const int i_r = i_r_coarse * 2; const int i_theta = i_theta_coarse * 2; @@ -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& fine_grid, - const PolarGrid& coarse_grid, - HostVector coarse_result, - HostConstVector fine_values) const +void Interpolation::applyExtrapolatedRestriction(const PolarGrid& fine_grid, + const PolarGrid& coarse_grid, + Vector coarse_result, ConstVector fine_values) const { assert(std::ssize(fine_values) == fine_grid.numberOfNodes()); assert(std::ssize(coarse_result) == coarse_grid.numberOfNodes()); @@ -64,7 +63,7 @@ void Interpolation::applyExtrapolatedRestriction(const PolarGrid>( // Rank of the index space + Kokkos::MDRangePolicy>( // 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 ), @@ -77,7 +76,7 @@ void Interpolation::applyExtrapolatedRestriction(const PolarGrid>( // Rank of the index space + Kokkos::MDRangePolicy>( // 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 ), diff --git a/tests/Interpolation/extrapolated_prolongation.cpp b/tests/Interpolation/extrapolated_prolongation.cpp index 47195d9a..137c8ab5 100644 --- a/tests/Interpolation/extrapolated_prolongation.cpp +++ b/tests/Interpolation/extrapolated_prolongation.cpp @@ -47,20 +47,27 @@ TEST(ExtrapolatedProlongationTest, ExtrapolatedProlongationMatchesStencil) std::vector 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 fine_grid(fine_radii, fine_angles); - PolarGrid coarse_grid = coarseningGrid(fine_grid); + PolarGrid fine_grid(fine_radii, fine_angles); + PolarGrid coarse_grid = coarseningGrid(fine_grid); + + PolarGrid h_fine_grid(fine_grid); + PolarGrid h_coarse_grid(coarse_grid); Interpolation I(/*DirBC*/ true); - HostVector coarse_values = generate_random_sample_data(coarse_grid, 1234, 0.0, 1.0); - HostVector fine_result("fine_result", fine_grid.numberOfNodes()); + HostVector h_coarse_values = generate_random_sample_data(h_coarse_grid, 1234, 0.0, 1.0); + Vector 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 << ")"; } } diff --git a/tests/Interpolation/extrapolated_restriction.cpp b/tests/Interpolation/extrapolated_restriction.cpp index 403b4b3b..d6b7273d 100644 --- a/tests/Interpolation/extrapolated_restriction.cpp +++ b/tests/Interpolation/extrapolated_restriction.cpp @@ -46,21 +46,28 @@ TEST(ExtrapolatedRestrictionTest, ExtrapolatedRestrictionMatchesStencil) std::vector 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 fine_grid(fine_radii, fine_angles); - PolarGrid coarse_grid = coarseningGrid(fine_grid); + PolarGrid fine_grid(fine_radii, fine_angles); + PolarGrid coarse_grid = coarseningGrid(fine_grid); + + PolarGrid h_fine_grid(fine_grid); + PolarGrid h_coarse_grid(coarse_grid); Interpolation I(/*DirBC*/ true); - HostVector fine_values = generate_random_sample_data(fine_grid, 9012, 0.0, 1.0); - HostVector coarse_result("coarse_result", coarse_grid.numberOfNodes()); + HostVector h_fine_values = generate_random_sample_data(h_fine_grid, 9012, 0.0, 1.0); + Vector 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 << ")"; } }