diff --git a/include/DirectSolver/DirectSolverGive/directSolverGive.h b/include/DirectSolver/DirectSolverGive/directSolverGive.h index a21063ce..458fc977 100644 --- a/include/DirectSolver/DirectSolverGive/directSolverGive.h +++ b/include/DirectSolver/DirectSolverGive/directSolverGive.h @@ -13,7 +13,7 @@ class DirectSolverGive : public DirectSolver bool DirBC_Interior); // Note: The rhs (right-hand side) vector gets overwritten during the solution process. - void solveInPlace(HostVector solution) override; + void solveInPlace(Vector solution) override; private: #ifdef GMGPOLAR_USE_MUMPS diff --git a/include/DirectSolver/DirectSolverGive/directSolverGive.inl b/include/DirectSolver/DirectSolverGive/directSolverGive.inl index 05f09d13..b9d856b3 100644 --- a/include/DirectSolver/DirectSolverGive/directSolverGive.inl +++ b/include/DirectSolver/DirectSolverGive/directSolverGive.inl @@ -14,9 +14,8 @@ DirectSolverGive::DirectSolverGive(const PolarGrid -void DirectSolverGive::solveInPlace(HostVector h_solution) +void DirectSolverGive::solveInPlace(Vector solution) { - auto solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_solution); // Adjusts the right-hand side vector to account for symmetry corrections. // This transforms the system matrixA * solution = rhs into the equivalent system: // symmetric_DBc(matrixA) * solution = rhs - applySymmetryShift(rhs). @@ -25,6 +24,4 @@ void DirectSolverGive::solveInPlace(HostVector h_solutio applySymmetryShift(solution); // Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver. system_solver_.solveInPlace(solution); - - Kokkos::deep_copy(h_solution, solution); } diff --git a/include/DirectSolver/DirectSolverTake/directSolverTake.h b/include/DirectSolver/DirectSolverTake/directSolverTake.h index 60042499..ffc2301a 100644 --- a/include/DirectSolver/DirectSolverTake/directSolverTake.h +++ b/include/DirectSolver/DirectSolverTake/directSolverTake.h @@ -13,7 +13,7 @@ class DirectSolverTake : public DirectSolver bool DirBC_Interior); // Note: The rhs (right-hand side) vector gets overwritten during the solution process. - void solveInPlace(HostVector solution) override; + void solveInPlace(Vector solution) override; private: #ifdef GMGPOLAR_USE_MUMPS diff --git a/include/DirectSolver/DirectSolverTake/directSolverTake.inl b/include/DirectSolver/DirectSolverTake/directSolverTake.inl index 9bdea7e9..144ad98f 100644 --- a/include/DirectSolver/DirectSolverTake/directSolverTake.inl +++ b/include/DirectSolver/DirectSolverTake/directSolverTake.inl @@ -14,9 +14,8 @@ DirectSolverTake::DirectSolverTake(const PolarGrid -void DirectSolverTake::solveInPlace(HostVector h_solution) +void DirectSolverTake::solveInPlace(Vector solution) { - auto solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_solution); // Adjusts the right-hand side vector to account for symmetry corrections. // This transforms the system matrixA * solution = rhs into the equivalent system: // symmetric_DBc(matrixA) * solution = rhs - applySymmetryShift(rhs). @@ -25,6 +24,4 @@ void DirectSolverTake::solveInPlace(HostVector h_solutio applySymmetryShift(solution); // Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver. system_solver_.solveInPlace(solution); - - Kokkos::deep_copy(h_solution, solution); } diff --git a/include/DirectSolver/directSolver.h b/include/DirectSolver/directSolver.h index 881dd6bf..1dc9e392 100644 --- a/include/DirectSolver/directSolver.h +++ b/include/DirectSolver/directSolver.h @@ -37,7 +37,7 @@ class DirectSolver virtual ~DirectSolver() = default; // Note: The rhs (right-hand side) vector gets overwritten during the solution process. - virtual void solveInPlace(HostVector solution) = 0; + virtual void solveInPlace(Vector solution) = 0; protected: const PolarGrid& grid_; diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h index c5c216cc..9d104085 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h @@ -68,7 +68,7 @@ class ExtrapolatedSmootherGive : public ExtrapolatedSmoother // Parallel implementation using OpenMP: // Scedule every 2nd/4th line in parallel to avoid race conditions arising from the A-Give distribution. // Sceduling every 3rd line in parallel would also be possible, but is less natural for the 2 coloring. - void extrapolatedSmoothing(HostVector x, HostConstVector rhs, HostVector temp) override; + void extrapolatedSmoothing(Vector x, ConstVector rhs, Vector temp) override; private: /* ------------------- */ diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.inl index a670b3e8..c355f861 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.inl @@ -42,14 +42,9 @@ ExtrapolatedSmootherGive::ExtrapolatedSmootherGive(const PolarGr // - The system is then solved in-place in temp, and the results // are copied back to x. template -void ExtrapolatedSmootherGive::extrapolatedSmoothing(HostVector h_x, - HostConstVector h_rhs, - HostVector h_temp) +void ExtrapolatedSmootherGive::extrapolatedSmoothing(Vector x, ConstVector rhs, + Vector temp) { - auto x = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_x); - auto rhs = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_rhs); - auto temp = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_temp); - assert(x.size() == rhs.size()); assert(temp.size() == rhs.size()); @@ -91,7 +86,4 @@ void ExtrapolatedSmootherGive::extrapolatedSmoothing(HostVector< applyAscOrthoWhiteRadialSection(x, rhs, temp); solveWhiteRadialSection(x, temp); - - Kokkos::deep_copy(h_x, x); - Kokkos::deep_copy(h_temp, temp); } diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.h b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.h index d8552bbe..a7964365 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.h +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.h @@ -66,7 +66,7 @@ class ExtrapolatedSmootherTake : public ExtrapolatedSmoother // Performs one full coupled extrapolated smoothing sweep: // BC -> WC -> BR -> WR // using temp as RHS workspace. - void extrapolatedSmoothing(HostVector x, HostConstVector rhs, HostVector temp) override; + void extrapolatedSmoothing(Vector x, ConstVector rhs, Vector temp) override; private: /* ------------------- */ diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.inl index bc3e9687..961195d2 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.inl @@ -43,14 +43,9 @@ ExtrapolatedSmootherTake::ExtrapolatedSmootherTake(const PolarGr // are copied back to x. template -void ExtrapolatedSmootherTake::extrapolatedSmoothing(HostVector h_x, - HostConstVector h_rhs, - HostVector h_temp) +void ExtrapolatedSmootherTake::extrapolatedSmoothing(Vector x, ConstVector rhs, + Vector temp) { - auto x = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_x); - auto rhs = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_rhs); - auto temp = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_temp); - assert(x.size() == rhs.size()); assert(temp.size() == rhs.size()); @@ -65,7 +60,4 @@ void ExtrapolatedSmootherTake::extrapolatedSmoothing(HostVector< applyAscOrthoWhiteRadialSection(x, rhs, temp); solveWhiteRadialSection(x, temp); - - Kokkos::deep_copy(h_x, x); - Kokkos::deep_copy(h_temp, temp); } diff --git a/include/ExtrapolatedSmoother/extrapolatedSmoother.h b/include/ExtrapolatedSmoother/extrapolatedSmoother.h index 38ca430f..9f1fd652 100644 --- a/include/ExtrapolatedSmoother/extrapolatedSmoother.h +++ b/include/ExtrapolatedSmoother/extrapolatedSmoother.h @@ -36,7 +36,7 @@ class ExtrapolatedSmoother } virtual ~ExtrapolatedSmoother() = default; - virtual void extrapolatedSmoothing(HostVector x, HostConstVector rhs, HostVector temp) = 0; + virtual void extrapolatedSmoothing(Vector x, ConstVector rhs, Vector temp) = 0; protected: const PolarGrid& grid_; diff --git a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h index a647a905..78f05d34 100644 --- a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h +++ b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h @@ -2,9 +2,9 @@ template void GMGPolar::extrapolated_multigrid_F_Cycle(int level_depth, - HostVector solution, - HostConstVector rhs, - HostVector residual) + HostVector h_solution, + HostConstVector h_rhs, + HostVector h_residual) { assert(level_depth == 0); assert(number_of_levels_ >= 2); @@ -25,6 +25,12 @@ void GMGPolar::extrapolated_multigri /* ------------ */ auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now(); + auto residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_residual); + auto solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_solution); + auto rhs = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_rhs); // const + auto next_level_residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.residual()); + auto next_level_solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.solution()); + for (int i = 0; i < pre_smoothing_steps_; i++) { if (level.level_depth() == 0 && !full_grid_smoothing_) { level.extrapolatedSmoothing(solution, rhs, residual); @@ -38,21 +44,23 @@ void GMGPolar::extrapolated_multigri t_avg_MGC_preSmoothing_ += std::chrono::duration(end_MGC_preSmoothing - start_MGC_preSmoothing).count(); /* -------------------------------------------- */ - /* Compute the restricted extrapolated residual */ + /* Compute the restricted extrapolated h_residual */ /* -------------------------------------------- */ auto start_MGC_residual = std::chrono::high_resolution_clock::now(); // P_ex^T (f_l - A_l*u_l) level.computeResidual(residual, rhs, solution); - extrapolatedRestriction(level.level_depth(), next_level.residual(), residual); + extrapolatedRestriction(level.level_depth(), next_level_residual, residual); // f_{l-1} - A_{l-1}* Inject(u_l) - injection(level.level_depth(), next_level.solution(), solution); - next_level.computeResidual(next_level.error_correction(), next_level.rhs(), next_level.solution()); + injection(level.level_depth(), next_level_solution, solution); + next_level.computeResidual(next_level.error_correction(), next_level.rhs(), next_level_solution); // res_ex = 4/3 * P_ex^T (f_l - A_l*u_l) - 1/3 * (f_{l-1} - A_{l-1}* Inject(u_l)) - linear_combination(next_level.residual(), 4.0 / 3.0, HostConstVector(next_level.error_correction()), + linear_combination(next_level_residual, 4.0 / 3.0, ConstVector(next_level.error_correction()), -1.0 / 3.0); + Kokkos::deep_copy(next_level.residual(), next_level_residual); + Kokkos::deep_copy(next_level.solution(), next_level_solution); auto end_MGC_residual = std::chrono::high_resolution_clock::now(); t_avg_MGC_residual_ += std::chrono::duration(end_MGC_residual - start_MGC_residual).count(); @@ -66,12 +74,13 @@ void GMGPolar::extrapolated_multigri assign(next_level.error_correction(), 0.0); /* Solve for the error by recursively calling the multigrid cycle. */ - multigrid_F_Cycle(next_level.level_depth(), next_level.error_correction(), next_level.residual(), + auto h_next_level_error_correction = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, next_level.error_correction()); + multigrid_F_Cycle(next_level.level_depth(), h_next_level_error_correction, next_level.residual(), next_level.solution()); /* Don't do a second recursion on the coarsest level since the DirectSolver is exact. */ if (next_level.level_depth() != number_of_levels_ - 1) { - multigrid_V_Cycle(next_level.level_depth(), next_level.error_correction(), next_level.residual(), + multigrid_V_Cycle(next_level.level_depth(), h_next_level_error_correction, next_level.residual(), next_level.solution()); } @@ -80,12 +89,13 @@ void GMGPolar::extrapolated_multigri /* -------------------------- */ // Use 'residual' instead of 'level.error_correction()' as a temporary buffer. // Note: 'level.error_correction()' has size 0 at level depth = 0. - extrapolatedProlongation(next_level.level_depth(), residual, next_level.error_correction()); + Kokkos::deep_copy(next_level.error_correction(), h_next_level_error_correction); + extrapolatedProlongation(next_level.level_depth(), residual, next_level.error_correction()); /* ----------------------------------- */ /* Compute the corrected approximation */ /* ----------------------------------- */ - add(solution, HostConstVector(residual)); + add(solution, ConstVector(residual)); /* ------------- */ /* Postsmoothing */ @@ -100,6 +110,8 @@ void GMGPolar::extrapolated_multigri level.smoothing(solution, rhs, residual); } } + Kokkos::deep_copy(h_residual, residual); + Kokkos::deep_copy(h_solution, solution); auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now(); t_avg_MGC_postSmoothing_ += std::chrono::duration(end_MGC_postSmoothing - start_MGC_postSmoothing).count(); diff --git a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h index 87c916aa..a97212f2 100644 --- a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h +++ b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h @@ -2,9 +2,9 @@ template void GMGPolar::extrapolated_multigrid_V_Cycle(int level_depth, - HostVector solution, - HostConstVector rhs, - HostVector residual) + HostVector h_solution, + HostConstVector h_rhs, + HostVector h_residual) { assert(level_depth == 0); assert(number_of_levels_ >= 2); @@ -25,6 +25,12 @@ void GMGPolar::extrapolated_multigri /* ------------ */ auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now(); + auto residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_residual); + auto solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_solution); + auto rhs = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_rhs); // const + auto next_level_residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.residual()); + auto next_level_solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.solution()); + for (int i = 0; i < pre_smoothing_steps_; i++) { if (level.level_depth() == 0 && !full_grid_smoothing_) { level.extrapolatedSmoothing(solution, rhs, residual); @@ -44,15 +50,17 @@ void GMGPolar::extrapolated_multigri // P_ex^T (f_l - A_l*u_l) level.computeResidual(residual, rhs, solution); - extrapolatedRestriction(level.level_depth(), next_level.residual(), residual); + extrapolatedRestriction(level.level_depth(), next_level_residual, residual); // f_{l-1} - A_{l-1}* Inject(u_l) - injection(level.level_depth(), next_level.solution(), solution); - next_level.computeResidual(next_level.error_correction(), next_level.rhs(), next_level.solution()); + injection(level.level_depth(), next_level_solution, solution); + next_level.computeResidual(next_level.error_correction(), next_level.rhs(), next_level_solution); // res_ex = 4/3 * P_ex^T (f_l - A_l*u_l) - 1/3 * (f_{l-1} - A_{l-1}* Inject(u_l)) - linear_combination(next_level.residual(), 4.0 / 3.0, HostConstVector(next_level.error_correction()), + linear_combination(next_level_residual, 4.0 / 3.0, ConstVector(next_level.error_correction()), -1.0 / 3.0); + Kokkos::deep_copy(next_level.residual(), next_level_residual); + Kokkos::deep_copy(next_level.solution(), next_level_solution); auto end_MGC_residual = std::chrono::high_resolution_clock::now(); t_avg_MGC_residual_ += std::chrono::duration(end_MGC_residual - start_MGC_residual).count(); @@ -66,7 +74,8 @@ void GMGPolar::extrapolated_multigri assign(next_level.error_correction(), 0.0); /* Solve for the error by recursively calling the multigrid cycle. */ - multigrid_V_Cycle(next_level.level_depth(), next_level.error_correction(), next_level.residual(), + auto h_next_level_error_correction = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, next_level.error_correction()); + multigrid_V_Cycle(next_level.level_depth(), h_next_level_error_correction, next_level.residual(), next_level.solution()); /* -------------------------- */ @@ -74,12 +83,13 @@ void GMGPolar::extrapolated_multigri /* -------------------------- */ // Use 'residual' instead of 'level.error_correction()' as a temporary buffer. // Note: 'level.error_correction()' has size 0 at level depth = 0. + Kokkos::deep_copy(next_level.error_correction(), h_next_level_error_correction); extrapolatedProlongation(next_level.level_depth(), residual, next_level.error_correction()); /* ----------------------------------- */ /* Compute the corrected approximation */ /* ----------------------------------- */ - add(solution, HostConstVector(residual)); + add(solution, ConstVector(residual)); /* ------------- */ /* Postsmoothing */ @@ -94,6 +104,8 @@ void GMGPolar::extrapolated_multigri level.smoothing(solution, rhs, residual); } } + Kokkos::deep_copy(h_residual, residual); + Kokkos::deep_copy(h_solution, solution); auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now(); t_avg_MGC_postSmoothing_ += std::chrono::duration(end_MGC_postSmoothing - start_MGC_postSmoothing).count(); diff --git a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h index 67e3a1d0..249e6d7d 100644 --- a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h +++ b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h @@ -2,9 +2,9 @@ template void GMGPolar::extrapolated_multigrid_W_Cycle(int level_depth, - HostVector solution, - HostConstVector rhs, - HostVector residual) + HostVector h_solution, + HostConstVector h_rhs, + HostVector h_residual) { assert(level_depth == 0); assert(number_of_levels_ >= 2); @@ -25,6 +25,12 @@ void GMGPolar::extrapolated_multigri /* ------------ */ auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now(); + auto residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_residual); + auto solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_solution); + auto rhs = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_rhs); // const + auto next_level_residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.residual()); + auto next_level_solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.solution()); + for (int i = 0; i < pre_smoothing_steps_; i++) { if (level.level_depth() == 0 && !full_grid_smoothing_) { level.extrapolatedSmoothing(solution, rhs, residual); @@ -44,14 +50,14 @@ void GMGPolar::extrapolated_multigri // P_ex^T (f_l - A_l*u_l) level.computeResidual(residual, rhs, solution); - extrapolatedRestriction(level.level_depth(), next_level.residual(), residual); + extrapolatedRestriction(level.level_depth(), next_level_residual, residual); // f_{l-1} - A_{l-1}* Inject(u_l) - injection(level.level_depth(), next_level.solution(), solution); - next_level.computeResidual(next_level.error_correction(), next_level.rhs(), next_level.solution()); + injection(level.level_depth(), next_level_solution, solution); + next_level.computeResidual(next_level.error_correction(), next_level.rhs(), next_level_solution); // res_ex = 4/3 * P_ex^T (f_l - A_l*u_l) - 1/3 * (f_{l-1} - A_{l-1}* Inject(u_l)) - linear_combination(next_level.residual(), 4.0 / 3.0, HostConstVector(next_level.error_correction()), + linear_combination(next_level_residual, 4.0 / 3.0, ConstVector(next_level.error_correction()), -1.0 / 3.0); auto end_MGC_residual = std::chrono::high_resolution_clock::now(); @@ -66,14 +72,16 @@ void GMGPolar::extrapolated_multigri assign(next_level.error_correction(), 0.0); /* Solve for the error by recursively calling the multigrid cycle. */ - multigrid_W_Cycle(next_level.level_depth(), next_level.error_correction(), next_level.residual(), - next_level.solution()); + multigrid_W_Cycle(next_level.level_depth(), next_level.error_correction(), next_level_residual, + next_level_solution); /* Don't do a second recursion on the coarsest level since the DirectSolver is exact. */ if (next_level.level_depth() != number_of_levels_ - 1) { - multigrid_W_Cycle(next_level.level_depth(), next_level.error_correction(), next_level.residual(), - next_level.solution()); + multigrid_W_Cycle(next_level.level_depth(), next_level.error_correction(), next_level_residual, + next_level_solution); } + Kokkos::deep_copy(next_level.residual(), next_level_residual); + Kokkos::deep_copy(next_level.solution(), next_level_solution); /* -------------------------- */ /* Interpolate the correction */ @@ -85,7 +93,7 @@ void GMGPolar::extrapolated_multigri /* ----------------------------------- */ /* Compute the corrected approximation */ /* ----------------------------------- */ - add(solution, HostConstVector(residual)); + add(solution, ConstVector(residual)); /* ------------- */ /* Postsmoothing */ @@ -100,6 +108,8 @@ void GMGPolar::extrapolated_multigri level.smoothing(solution, rhs, residual); } } + Kokkos::deep_copy(h_residual, residual); + Kokkos::deep_copy(h_solution, solution); auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now(); t_avg_MGC_postSmoothing_ += std::chrono::duration(end_MGC_postSmoothing - start_MGC_postSmoothing).count(); diff --git a/include/GMGPolar/MultigridMethods/multigrid_F_Cycle.h b/include/GMGPolar/MultigridMethods/multigrid_F_Cycle.h index bd6493ef..72ccc7ce 100644 --- a/include/GMGPolar/MultigridMethods/multigrid_F_Cycle.h +++ b/include/GMGPolar/MultigridMethods/multigrid_F_Cycle.h @@ -2,10 +2,11 @@ template void GMGPolar::multigrid_F_Cycle(int level_depth, - HostVector solution, - HostConstVector rhs, - HostVector residual) + HostVector h_solution, + HostConstVector h_rhs, + HostVector h_residual) { + auto solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_solution); assert(0 <= level_depth && level_depth < number_of_levels_); std::chrono::high_resolution_clock::time_point start_MGC; @@ -15,14 +16,14 @@ void GMGPolar::multigrid_F_Cycle(int if (level_depth == number_of_levels_ - 1) { /* ---------------------------------------------------- */ - /* Coarsest level: solve A * x = rhs using DirectSolver */ + /* Coarsest level: solve A * x = h_rhs using DirectSolver */ /* ---------------------------------------------------- */ Level& coarsest_level = levels_[level_depth]; - /* Step 1: Copy rhs in solution */ - Kokkos::deep_copy(solution, rhs); + /* Step 1: Copy h_rhs in h_solution */ + Kokkos::deep_copy(h_solution, h_rhs); - /* Step 2: Solve for the solution in place */ + /* Step 2: Solve for the h_solution in place */ auto start_MGC_directSolver = std::chrono::high_resolution_clock::now(); coarsest_level.directSolveInPlace(solution); @@ -42,6 +43,9 @@ void GMGPolar::multigrid_F_Cycle(int /* ------------ */ auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now(); + auto residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_residual); + auto rhs = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_rhs); // const + auto next_level_residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.residual()); for (int i = 0; i < pre_smoothing_steps_; i++) { level.smoothing(solution, rhs, residual); } @@ -62,7 +66,8 @@ void GMGPolar::multigrid_F_Cycle(int /* --------------------- */ /* Restrict the residual */ /* --------------------- */ - restriction(level.level_depth(), next_level.residual(), residual); + restriction(level.level_depth(), next_level_residual, residual); + Kokkos::deep_copy(next_level.residual(), next_level_residual); /* ------------------------------------- */ /* Solve A * error = restricted residual */ @@ -72,12 +77,13 @@ void GMGPolar::multigrid_F_Cycle(int assign(next_level.error_correction(), 0.0); /* Solve for the error by recursively calling the multigrid cycle. */ - multigrid_F_Cycle(next_level.level_depth(), next_level.error_correction(), next_level.residual(), + auto h_next_level_error_correction = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), next_level.error_correction()); + multigrid_F_Cycle(next_level.level_depth(), h_next_level_error_correction, next_level.residual(), next_level.solution()); /* Don't do a second recursion on the coarsest level since the DirectSolver is exact. */ if (next_level.level_depth() != number_of_levels_ - 1) { - multigrid_V_Cycle(next_level.level_depth(), next_level.error_correction(), next_level.residual(), + multigrid_V_Cycle(next_level.level_depth(), h_next_level_error_correction, next_level.residual(), next_level.solution()); } @@ -86,12 +92,13 @@ void GMGPolar::multigrid_F_Cycle(int /* -------------------------- */ // Use 'residual' instead of 'level.error_correction()' as a temporary buffer. // Note: 'level.error_correction()' has size 0 at level depth = 0. + Kokkos::deep_copy(next_level.error_correction(), h_next_level_error_correction); prolongation(next_level.level_depth(), residual, next_level.error_correction()); /* ----------------------------------- */ /* Compute the corrected approximation */ /* ----------------------------------- */ - add(solution, HostConstVector(residual)); + add(solution, ConstVector(residual)); /* ------------- */ /* Postsmoothing */ @@ -101,6 +108,7 @@ void GMGPolar::multigrid_F_Cycle(int for (int i = 0; i < post_smoothing_steps_; i++) { level.smoothing(solution, rhs, residual); } + Kokkos::deep_copy(h_residual, residual); auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now(); t_avg_MGC_postSmoothing_ += @@ -111,4 +119,5 @@ void GMGPolar::multigrid_F_Cycle(int auto end_MGC = std::chrono::high_resolution_clock::now(); t_avg_MGC_total_ += std::chrono::duration(end_MGC - start_MGC).count(); } + Kokkos::deep_copy(h_solution, solution); } diff --git a/include/GMGPolar/MultigridMethods/multigrid_V_Cycle.h b/include/GMGPolar/MultigridMethods/multigrid_V_Cycle.h index 0bd95f1e..9d845429 100644 --- a/include/GMGPolar/MultigridMethods/multigrid_V_Cycle.h +++ b/include/GMGPolar/MultigridMethods/multigrid_V_Cycle.h @@ -2,10 +2,11 @@ template void GMGPolar::multigrid_V_Cycle(int level_depth, - HostVector solution, - HostConstVector rhs, - HostVector residual) + HostVector h_solution, + HostConstVector h_rhs, + HostVector h_residual) { + auto solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_solution); assert(0 <= level_depth && level_depth < number_of_levels_); std::chrono::high_resolution_clock::time_point start_MGC; @@ -15,12 +16,12 @@ void GMGPolar::multigrid_V_Cycle(int if (level_depth == number_of_levels_ - 1) { /* ---------------------------------------------------- */ - /* Coarsest level: solve A * x = rhs using DirectSolver */ + /* Coarsest level: solve A * x = h_rhs using DirectSolver */ /* ---------------------------------------------------- */ Level& coarsest_level = levels_[level_depth]; - /* Step 1: Copy rhs in solution */ - Kokkos::deep_copy(solution, rhs); + /* Step 1: Copy h_rhs in solution */ + Kokkos::deep_copy(solution, h_rhs); /* Step 2: Solve for the solution in place */ auto start_MGC_directSolver = std::chrono::high_resolution_clock::now(); @@ -42,6 +43,9 @@ void GMGPolar::multigrid_V_Cycle(int /* ------------ */ auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now(); + auto residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_residual); + auto rhs = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_rhs); // const + auto next_level_residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.residual()); for (int i = 0; i < pre_smoothing_steps_; i++) { level.smoothing(solution, rhs, residual); } @@ -62,7 +66,8 @@ void GMGPolar::multigrid_V_Cycle(int /* --------------------- */ /* Restrict the residual */ /* --------------------- */ - restriction(level.level_depth(), next_level.residual(), residual); + restriction(level.level_depth(), next_level_residual, residual); + Kokkos::deep_copy(next_level.residual(), next_level_residual); /* ------------------------------------- */ /* Solve A * error = restricted residual */ @@ -72,7 +77,8 @@ void GMGPolar::multigrid_V_Cycle(int assign(next_level.error_correction(), 0.0); /* Solve for the error by recursively calling the multigrid cycle. */ - multigrid_V_Cycle(next_level.level_depth(), next_level.error_correction(), next_level.residual(), + auto h_next_level_error_correction = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), next_level.error_correction()); + multigrid_V_Cycle(next_level.level_depth(), h_next_level_error_correction, next_level.residual(), next_level.solution()); /* -------------------------- */ @@ -80,12 +86,13 @@ void GMGPolar::multigrid_V_Cycle(int /* -------------------------- */ // Use 'residual' instead of 'level.error_correction()' as a temporary buffer. // Note: 'level.error_correction()' has size 0 at level depth = 0. + Kokkos::deep_copy(next_level.error_correction(), h_next_level_error_correction); prolongation(next_level.level_depth(), residual, next_level.error_correction()); /* ----------------------------------- */ /* Compute the corrected approximation */ /* ----------------------------------- */ - add(solution, HostConstVector(residual)); + add(solution, ConstVector(residual)); /* ------------- */ /* Postsmoothing */ @@ -95,6 +102,7 @@ void GMGPolar::multigrid_V_Cycle(int for (int i = 0; i < post_smoothing_steps_; i++) { level.smoothing(solution, rhs, residual); } + Kokkos::deep_copy(h_residual, residual); auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now(); t_avg_MGC_postSmoothing_ += @@ -105,4 +113,5 @@ void GMGPolar::multigrid_V_Cycle(int auto end_MGC = std::chrono::high_resolution_clock::now(); t_avg_MGC_total_ += std::chrono::duration(end_MGC - start_MGC).count(); } + Kokkos::deep_copy(h_solution, solution); } diff --git a/include/GMGPolar/MultigridMethods/multigrid_W_Cycle.h b/include/GMGPolar/MultigridMethods/multigrid_W_Cycle.h index 6c11283f..db9dc5a4 100644 --- a/include/GMGPolar/MultigridMethods/multigrid_W_Cycle.h +++ b/include/GMGPolar/MultigridMethods/multigrid_W_Cycle.h @@ -2,9 +2,9 @@ template void GMGPolar::multigrid_W_Cycle(int level_depth, - HostVector solution, - HostConstVector rhs, - HostVector residual) + Vector solution, + ConstVector rhs, + Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_); @@ -42,6 +42,7 @@ void GMGPolar::multigrid_W_Cycle(int /* ------------ */ auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now(); + auto next_level_residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.residual()); for (int i = 0; i < pre_smoothing_steps_; i++) { level.smoothing(solution, rhs, residual); } @@ -62,7 +63,7 @@ void GMGPolar::multigrid_W_Cycle(int /* --------------------- */ /* Restrict the residual */ /* --------------------- */ - restriction(level.level_depth(), next_level.residual(), residual); + restriction(level.level_depth(), next_level_residual, residual); /* ------------------------------------- */ /* Solve A * error = restricted residual */ @@ -72,14 +73,17 @@ void GMGPolar::multigrid_W_Cycle(int assign(next_level.error_correction(), 0.0); /* Solve for the error by recursively calling the multigrid cycle. */ - multigrid_W_Cycle(next_level.level_depth(), next_level.error_correction(), next_level.residual(), - next_level.solution()); + auto next_level_solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.solution()); + multigrid_W_Cycle(next_level.level_depth(), next_level.error_correction(), next_level_residual, + next_level_solution); /* Don't do a second recursion on the coarsest level since the DirectSolver is exact. */ if (next_level.level_depth() != number_of_levels_ - 1) { - multigrid_W_Cycle(next_level.level_depth(), next_level.error_correction(), next_level.residual(), - next_level.solution()); + multigrid_W_Cycle(next_level.level_depth(), next_level.error_correction(), next_level_residual, + next_level_solution); } + Kokkos::deep_copy(next_level.solution(), next_level_solution); + Kokkos::deep_copy(next_level.residual(), next_level_residual); /* -------------------------- */ /* Interpolate the correction */ @@ -91,7 +95,7 @@ void GMGPolar::multigrid_W_Cycle(int /* ----------------------------------- */ /* Compute the corrected approximation */ /* ----------------------------------- */ - add(solution, HostConstVector(residual)); + add(solution, ConstVector(residual)); /* ------------- */ /* Postsmoothing */ diff --git a/include/GMGPolar/gmgpolar.h b/include/GMGPolar/gmgpolar.h index 92c7a1e9..aa5af61a 100644 --- a/include/GMGPolar/gmgpolar.h +++ b/include/GMGPolar/gmgpolar.h @@ -154,7 +154,7 @@ class GMGPolar : public IGMGPolar /* --------------- */ /* Solve Functions */ - void applyExtrapolation(int current_level, HostVector fine_values, HostConstVector coarse_values); + void applyExtrapolation(int current_level, Vector fine_values, ConstVector coarse_values); private: /* --------------- */ @@ -194,8 +194,8 @@ class GMGPolar : public IGMGPolar /* Multigrid Functions */ void multigrid_V_Cycle(int level_depth, HostVector solution, HostConstVector rhs, HostVector residual); - void multigrid_W_Cycle(int level_depth, HostVector solution, HostConstVector rhs, - HostVector residual); + void multigrid_W_Cycle(int level_depth, Vector solution, ConstVector rhs, + Vector residual); void multigrid_F_Cycle(int level_depth, HostVector solution, HostConstVector rhs, HostVector residual); void extrapolated_multigrid_V_Cycle(int level_depth, HostVector solution, HostConstVector rhs, @@ -207,11 +207,11 @@ class GMGPolar : public IGMGPolar /* ----------------------- */ /* Interpolation functions */ - void prolongation(int current_level, HostVector result, HostConstVector x) const; - void restriction(int current_level, HostVector result, HostConstVector x) const; - void injection(int current_level, HostVector result, HostConstVector x) const; - void extrapolatedProlongation(int current_level, HostVector result, HostConstVector x) const; - void extrapolatedRestriction(int current_level, HostVector result, HostConstVector x) const; + void prolongation(int current_level, Vector result, ConstVector x) const; + void restriction(int current_level, Vector result, ConstVector x) const; + void injection(int current_level, Vector result, ConstVector x) const; + void extrapolatedProlongation(int current_level, Vector result, ConstVector x) const; + void extrapolatedRestriction(int current_level, Vector result, ConstVector x) const; void FMGInterpolation(int current_level, HostVector result, HostConstVector x) const; /* ------------- */ diff --git a/include/GMGPolar/solver.h b/include/GMGPolar/solver.h index 980887fb..37cbe35f 100644 --- a/include/GMGPolar/solver.h +++ b/include/GMGPolar/solver.h @@ -13,10 +13,8 @@ void GMGPolar::solve(const BoundaryC /* ------------------------------------- */ /* Build rhs_f on Level 0 (finest Level) */ /* ------------------------------------- */ - HostVector rhs_f_host = levels_[0].rhs(); - auto rhs_f = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), rhs_f_host); + Vector rhs_f = levels_[0].rhs(); build_rhs_f(levels_[0], rhs_f, boundary_conditions, source_term); - Kokkos::deep_copy(rhs_f_host, rhs_f); /* ---------------- */ /* Discretize rhs_f */ @@ -31,7 +29,9 @@ void GMGPolar::solve(const BoundaryC injection(level_depth, next_level.rhs(), current_level.rhs()); } // Discretize the rhs for the current level - discretize_rhs_f(current_level, current_level.rhs()); + auto h_rhs = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, current_level.rhs()); + discretize_rhs_f(current_level, h_rhs); + Kokkos::deep_copy(current_level.rhs(), h_rhs); } auto end_setup_rhs = std::chrono::high_resolution_clock::now(); @@ -207,7 +207,9 @@ void GMGPolar::fullMultigridApproxim // Solve directly on the coarsest level Kokkos::deep_copy(coarsest_level.solution(), coarsest_level.rhs()); - coarsest_level.directSolveInPlace(coarsest_level.solution()); + auto coarsest_level_solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), coarsest_level.solution()); + coarsest_level.directSolveInPlace(coarsest_level_solution); + Kokkos::deep_copy(coarsest_level.solution(), coarsest_level_solution); // Prolongate the solution from the coarsest level up to the finest, while applying Multigrid Cycles on each level for (int depth = coarsest_depth; depth > 0; --depth) { @@ -318,7 +320,9 @@ void GMGPolar::solvePCG(double& init // z = M^{-1} * r (preconditioned residual) if (PCG_FMG_) { - initRhsHierarchy(level.rhs()); + auto h_rhs = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, level.rhs()); + initRhsHierarchy(h_rhs); + Kokkos::deep_copy(level.rhs(), h_rhs); fullMultigridApproximation(PCG_FMG_cycle_, PCG_FMG_iterations_); } else { @@ -331,18 +335,26 @@ void GMGPolar::solvePCG(double& init Kokkos::deep_copy(pcg_search_direction_, level.solution()); // r^T * z - double r_z = dot_product(HostConstVector(level.rhs()), HostConstVector(level.solution())); + auto h_rhs = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, level.rhs()); + double r_z = dot_product(HostConstVector(h_rhs), HostConstVector(level.solution())); while (number_of_iterations_ < max_iterations_) { // A_p = A * p - level.applySystemOperator(level.residual(), pcg_search_direction_); + auto level_residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), level.residual()); + auto pcg_search_direction = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), pcg_search_direction_); + level.applySystemOperator(level_residual, pcg_search_direction); + Kokkos::deep_copy(level.residual(), level_residual); if (extrapolation_ != ExtrapolationType::NONE) { assert(number_of_levels_ > 1); Level& next_level = levels_[level.level_depth() + 1]; - injection(0, next_level.solution(), pcg_search_direction_); - next_level.applySystemOperator(next_level.residual(), next_level.solution()); - applyExtrapolation(0, level.residual(), next_level.residual()); + auto next_level_solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.solution()); + auto next_level_residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.residual()); + injection(0, next_level_solution, pcg_search_direction); + next_level.applySystemOperator(next_level_residual, next_level_solution); + applyExtrapolation(0, level_residual, next_level_residual); + Kokkos::deep_copy(level.residual(), level_residual); + Kokkos::deep_copy(next_level.residual(), next_level_residual); } // alpha = (r^T * z) / (p^T * A*p) @@ -353,14 +365,15 @@ void GMGPolar::solvePCG(double& init linear_combination(pcg_solution_, 1.0, HostConstVector(pcg_search_direction_), alpha); // r -= alpha * A*p - linear_combination(level.rhs(), 1.0, HostConstVector(level.residual()), -alpha); + Kokkos::deep_copy(h_rhs, level.rhs()); + linear_combination(h_rhs, 1.0, HostConstVector(level.residual()), -alpha); /* ---------------------------- */ /* Compute convergence criteria */ /* ---------------------------- */ auto start_check_convergence = std::chrono::high_resolution_clock::now(); - current_residual_norm = residualNorm(residual_norm_type_, level, level.rhs()); + current_residual_norm = residualNorm(residual_norm_type_, level, h_rhs); residual_norms_.push_back(current_residual_norm); current_relative_residual_norm = current_residual_norm / initial_residual_norm; @@ -392,17 +405,17 @@ void GMGPolar::solvePCG(double& init // z = M^{-1} * r (preconditioned residual) if (PCG_FMG_) { - initRhsHierarchy(level.rhs()); + initRhsHierarchy(h_rhs); fullMultigridApproximation(PCG_FMG_cycle_, PCG_FMG_iterations_); } else { // z = I^{-1} * r (no preconditioning) - Kokkos::deep_copy(level.solution(), level.rhs()); + Kokkos::deep_copy(level.solution(), h_rhs); } applyMultigridIterations(level, PCG_MG_cycle_, PCG_MG_iterations_); // r^T * z - double r_z_new = dot_product(HostConstVector(level.rhs()), HostConstVector(level.solution())); + double r_z_new = dot_product(HostConstVector(h_rhs), HostConstVector(level.solution())); // beta = (current r^T * z) / (previous r^T * z) double beta = r_z_new / r_z; // r_z = r^T * z for next iteration @@ -425,16 +438,21 @@ template ::applyMultigridIterations( Level& level, MultigridCycleType cycle, int iterations) { + auto h_rhs = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, level.rhs()); + auto residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace{}, level.residual()); + auto solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace{}, level.solution()); for (int i = 0; i < iterations; i++) { switch (cycle) { case MultigridCycleType::V_CYCLE: - multigrid_V_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); + multigrid_V_Cycle(level.level_depth(), level.solution(), h_rhs, level.residual()); break; case MultigridCycleType::W_CYCLE: - multigrid_W_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); + multigrid_W_Cycle(level.level_depth(), solution, level.rhs(), residual); + Kokkos::deep_copy(level.residual(), residual); + Kokkos::deep_copy(level.solution(), solution); break; case MultigridCycleType::F_CYCLE: - multigrid_F_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); + multigrid_F_Cycle(level.level_depth(), level.solution(), h_rhs, level.residual()); break; default: std::cerr << "Error: Unknown multigrid cycle type!" << std::endl; @@ -447,16 +465,17 @@ template ::applyExtrapolatedMultigridIterations( Level& level, MultigridCycleType cycle, int iterations) { + auto h_rhs = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, level.rhs()); for (int i = 0; i < iterations; i++) { switch (cycle) { case MultigridCycleType::V_CYCLE: - extrapolated_multigrid_V_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); + extrapolated_multigrid_V_Cycle(level.level_depth(), level.solution(), h_rhs, level.residual()); break; case MultigridCycleType::W_CYCLE: - extrapolated_multigrid_W_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); + extrapolated_multigrid_W_Cycle(level.level_depth(), level.solution(), h_rhs, level.residual()); break; case MultigridCycleType::F_CYCLE: - extrapolated_multigrid_F_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); + extrapolated_multigrid_F_Cycle(level.level_depth(), level.solution(), h_rhs, level.residual()); break; default: std::cerr << "Error: Unknown multigrid cycle type!" << std::endl; @@ -474,19 +493,28 @@ void GMGPolar::updateResidualNorms( Level& level, int iteration, double& initial_residual_norm, double& current_residual_norm, double& current_relative_residual_norm) { - level.computeResidual(level.residual(), level.rhs(), level.solution()); + auto solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), level.solution()); + auto residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), level.residual()); + level.computeResidual(residual, level.rhs(), solution); if (extrapolation_ != ExtrapolationType::NONE) { Level& next_level = levels_[level.level_depth() + 1]; - injection(level.level_depth(), next_level.solution(), level.solution()); - next_level.computeResidual(next_level.residual(), next_level.rhs(), next_level.solution()); - applyExtrapolation(level.level_depth(), level.residual(), next_level.residual()); + auto next_level_solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.solution()); + auto next_level_residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.residual()); + auto next_level_rhs = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.rhs()); + injection(level.level_depth(), next_level_solution, solution); + next_level.computeResidual(next_level_residual, next_level_rhs, next_level_solution); + applyExtrapolation(level.level_depth(), residual, next_level_residual); + Kokkos::deep_copy(next_level.solution(), next_level_solution); } + Kokkos::deep_copy(level.residual(), residual); + Kokkos::deep_copy(level.solution(), solution); current_residual_norm = residualNorm(residual_norm_type_, level, level.residual()); residual_norms_.push_back(current_residual_norm); if (number_of_iterations_ == 0) { - initial_residual_norm = !FMG_ ? current_residual_norm : residualNorm(residual_norm_type_, level, level.rhs()); + auto h_rhs = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, level.rhs()); + initial_residual_norm = !FMG_ ? current_residual_norm : residualNorm(residual_norm_type_, level, h_rhs); } current_relative_residual_norm = current_residual_norm / initial_residual_norm; @@ -523,11 +551,11 @@ double GMGPolar::residualNorm( template void GMGPolar::applyExtrapolation(int current_level, - HostVector fine_values, - HostConstVector coarse_values) + Vector fine_values, + ConstVector coarse_values) { - const PolarGrid fineGrid(levels_[current_level].grid()); - const PolarGrid coarseGrid(levels_[current_level + 1].grid()); + const PolarGrid& fineGrid(levels_[current_level].grid()); + const PolarGrid& coarseGrid(levels_[current_level + 1].grid()); assert(std::ssize(fine_values) == fineGrid.numberOfNodes()); assert(std::ssize(coarse_values) == coarseGrid.numberOfNodes()); @@ -539,7 +567,7 @@ void GMGPolar::applyExtrapolation(in /* For loop matches circular access pattern */ Kokkos::parallel_for( "Extrapolation: Apply Extrapolation (Circular)", - Kokkos::MDRangePolicy>( + Kokkos::MDRangePolicy>( {0, 0}, {fineGrid.numberSmootherCircles(), fineGrid.ntheta()}), KOKKOS_LAMBDA(const int i_r, const int i_theta) { const int fine_idx = fineGrid.index(i_r, i_theta); @@ -556,8 +584,8 @@ void GMGPolar::applyExtrapolation(in /* For loop matches radial access pattern */ Kokkos::parallel_for( "Extrapolation: Apply Extrapolation (Radial)", - Kokkos::MDRangePolicy>({0, fineGrid.numberSmootherCircles()}, - {fineGrid.ntheta(), fineGrid.nr()}), + Kokkos::MDRangePolicy>({0, fineGrid.numberSmootherCircles()}, + {fineGrid.ntheta(), fineGrid.nr()}), KOKKOS_LAMBDA(const int i_theta, const int i_r) { const int fine_idx = fineGrid.index(i_r, i_theta); if (i_r & 1 || i_theta & 1) { diff --git a/include/GMGPolar/utils.h b/include/GMGPolar/utils.h index 2804e51a..b7a8250c 100644 --- a/include/GMGPolar/utils.h +++ b/include/GMGPolar/utils.h @@ -4,78 +4,69 @@ /* Interpolation */ /* ---------------------------------------------------------------------- */ template -void GMGPolar::prolongation(int current_level, HostVector result, - HostConstVector x) const +void GMGPolar::prolongation(int current_level, Vector result, + ConstVector x) const { 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()); + const PolarGrid& current_grid = levels_[current_level].grid(); + const PolarGrid& previous_grid = levels_[current_level - 1].grid(); interpolation_->applyProlongation(current_grid, previous_grid, result, x); } template -void GMGPolar::restriction(int current_level, HostVector result, - HostConstVector x) const +void GMGPolar::restriction(int current_level, Vector result, + ConstVector x) const { 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()); + const PolarGrid& current_grid = levels_[current_level].grid(); + const PolarGrid& next_grid = levels_[current_level + 1].grid(); interpolation_->applyRestriction(current_grid, next_grid, result, x); } template -void GMGPolar::injection(int current_level, HostVector result, - HostConstVector x) const +void GMGPolar::injection(int current_level, Vector result, + ConstVector x) const { 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()); + const PolarGrid& current_grid = levels_[current_level].grid(); + const PolarGrid& next_grid = levels_[current_level + 1].grid(); interpolation_->applyInjection(current_grid, next_grid, result, x); } template void GMGPolar::extrapolatedProlongation( - int current_level, HostVector result_host, HostConstVector x_host) const + int current_level, Vector result, ConstVector x) 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(); + const PolarGrid& current_grid = levels_[current_level].grid(); + const 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_host, - HostConstVector x_host) const + Vector result, + ConstVector x) 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(); + const PolarGrid& current_grid = levels_[current_level].grid(); + const 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 7198458f..786dc4d5 100644 --- a/include/Interpolation/interpolation.h +++ b/include/Interpolation/interpolation.h @@ -19,22 +19,23 @@ class Interpolation explicit Interpolation(bool DirBC_Interior); /* Remark: This injection is not scaled. */ - void applyInjection(const PolarGrid& fine_grid, const PolarGrid& coarse_grid, - HostVector coarse_result, HostConstVector fine_values) const; + void applyInjection(const PolarGrid& fine_grid, + const PolarGrid& coarse_grid, Vector coarse_result, + ConstVector fine_values) const; /* Bilinear interpolation operator */ - void applyProlongation(const PolarGrid& coarse_grid, - const PolarGrid& fine_grid, HostVector fine_result, - HostConstVector coarse_values) const; + void applyProlongation(const PolarGrid& coarse_grid, + const PolarGrid& fine_grid, Vector fine_result, + ConstVector 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 applyRestriction(const PolarGrid& fine_grid, + const PolarGrid& coarse_grid, Vector coarse_result, + ConstVector fine_values) const; void applyExtrapolatedRestriction(const PolarGrid& fine_grid, const PolarGrid& coarse_grid, Vector coarse_result, diff --git a/include/Level/level.h b/include/Level/level.h index 363f4c4b..dd2b7db9 100644 --- a/include/Level/level.h +++ b/include/Level/level.h @@ -83,37 +83,37 @@ class Level const PolarGrid& grid() const; const LevelCacheType& levelCache() const; - HostVector rhs(); - HostConstVector rhs() const; + Vector rhs(); + ConstVector rhs() const; HostVector solution(); HostConstVector solution() const; HostVector residual(); HostConstVector residual() const; - HostVector error_correction(); - HostConstVector error_correction() const; + Vector error_correction(); + ConstVector error_correction() const; // -------------- // // Apply Residual // void initializeResidual(const bool DirBC_Interior, const StencilDistributionMethod stencil_distribution_method); - void computeResidual(HostVector result, HostConstVector rhs, HostConstVector x) const; - void applySystemOperator(HostVector result, HostConstVector x) const; + void computeResidual(Vector result, ConstVector rhs, ConstVector x) const; + void applySystemOperator(Vector result, ConstVector x) const; // ------------------- // // Solve coarse System // void initializeDirectSolver(const bool DirBC_Interior, const StencilDistributionMethod stencil_distribution_method); // Note: The rhs (right-hand side) vector gets overwritten by the solution. - void directSolveInPlace(HostVector x) const; + void directSolveInPlace(Vector x) const; // --------------- // // Apply Smoothing // void initializeSmoothing(const bool DirBC_Interior, const StencilDistributionMethod stencil_distribution_method); - void smoothing(HostVector x, HostConstVector rhs, HostVector temp) const; + void smoothing(Vector x, ConstVector rhs, Vector temp) const; // ---------------------------- // // Apply Extrapolated Smoothing // void initializeExtrapolatedSmoothing(const bool DirBC_Interior, const StencilDistributionMethod stencil_distribution_method); - void extrapolatedSmoothing(HostVector x, HostConstVector rhs, HostVector temp) const; + void extrapolatedSmoothing(Vector x, ConstVector rhs, Vector temp) const; private: const int level_depth_; @@ -125,10 +125,10 @@ class Level std::unique_ptr> op_smoother_; std::unique_ptr> op_extrapolated_smoother_; - HostVector rhs_; + Vector rhs_; HostVector solution_; HostVector residual_; - HostVector error_correction_; + Vector error_correction_; }; template diff --git a/include/Level/level.inl b/include/Level/level.inl index bb570e13..3ac1f3f7 100644 --- a/include/Level/level.inl +++ b/include/Level/level.inl @@ -41,13 +41,13 @@ Level::levelCache() const } template -HostVector Level::rhs() +Vector Level::rhs() { return rhs_; } template -HostConstVector Level::rhs() const +ConstVector Level::rhs() const { return rhs_; } @@ -77,13 +77,13 @@ HostConstVector Level::resid } template -HostVector Level::error_correction() +Vector Level::error_correction() { return error_correction_; } template -HostConstVector Level::error_correction() const +ConstVector Level::error_correction() const { return error_correction_; } @@ -105,17 +105,17 @@ void Level::initializeResidual( } template -void Level::computeResidual(HostVector result, - HostConstVector rhs, - HostConstVector x) const +void Level::computeResidual(Vector result, + ConstVector rhs, + ConstVector x) const { if (!op_residual_) throw std::runtime_error("Residual not initialized."); op_residual_->computeResidual(result, rhs, x); } template -void Level::applySystemOperator(HostVector result, - HostConstVector x) const +void Level::applySystemOperator(Vector result, + ConstVector x) const { if (!op_residual_) throw std::runtime_error("Residual not initialized."); @@ -140,7 +140,7 @@ void Level::initializeDirectSolver( } template -void Level::directSolveInPlace(HostVector x) const +void Level::directSolveInPlace(Vector x) const { if (!op_directSolver_) throw std::runtime_error("Coarse Solver not initialized."); @@ -164,8 +164,8 @@ void Level::initializeSmoothing( } template -void Level::smoothing(HostVector x, HostConstVector rhs, - HostVector temp) const +void Level::smoothing(Vector x, ConstVector rhs, + Vector temp) const { if (!op_smoother_) throw std::runtime_error("Smoother not initialized."); @@ -191,9 +191,8 @@ void Level::initializeExtrapolatedSm } template -void Level::extrapolatedSmoothing(HostVector x, - HostConstVector rhs, - HostVector temp) const +void Level::extrapolatedSmoothing(Vector x, ConstVector rhs, + Vector temp) const { if (!op_extrapolated_smoother_) throw std::runtime_error("Extrapolated Smoother not initialized."); diff --git a/include/LinearAlgebra/Vector/vector_operations.h b/include/LinearAlgebra/Vector/vector_operations.h index 1376b4fb..5b78e4a4 100644 --- a/include/LinearAlgebra/Vector/vector_operations.h +++ b/include/LinearAlgebra/Vector/vector_operations.h @@ -136,7 +136,7 @@ T l2_norm(ConstVector x) local_max = abs_val; } }, - Kokkos::Max(scale)); + Kokkos::Max(scale)); // 2) If the largest absolute value is zero, the norm is zero if (scale == T{0}) @@ -171,7 +171,7 @@ T infinity_norm(ConstVector x) local_max = abs_value; } }, - Kokkos::Max(result)); + Kokkos::Max(result)); return result; } diff --git a/include/Residual/ResidualGive/applyAGive.inl b/include/Residual/ResidualGive/applyAGive.inl index 0e2a4526..c5c0fe17 100644 --- a/include/Residual/ResidualGive/applyAGive.inl +++ b/include/Residual/ResidualGive/applyAGive.inl @@ -195,11 +195,8 @@ static KOKKOS_INLINE_FUNCTION void node_apply_a_give(int i_r, int i_theta, const } // namespace residual_give template -void ResidualGive::applySystemOperator(HostVector h_result, HostConstVector h_x) const +void ResidualGive::applySystemOperator(Vector result, ConstVector x) const { - auto x = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_x); - auto result = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_result); - assert(result.size() == x.size()); const PolarGrid& grid = Residual::grid_; @@ -264,6 +261,4 @@ void ResidualGive::applySystemOperator(HostVector h_resu }); Kokkos::fence(); } - - Kokkos::deep_copy(h_result, result); } diff --git a/include/Residual/ResidualGive/residualGive.h b/include/Residual/ResidualGive/residualGive.h index cdbaa0f5..781307c3 100644 --- a/include/Residual/ResidualGive/residualGive.h +++ b/include/Residual/ResidualGive/residualGive.h @@ -13,8 +13,8 @@ class ResidualGive : public Residual const bool DirBC_Interior); ~ResidualGive() override = default; - void applySystemOperator(HostVector result, HostConstVector x) const final; - void computeResidual(HostVector result, HostConstVector rhs, HostConstVector x) const final; + void applySystemOperator(Vector result, ConstVector x) const final; + void computeResidual(Vector result, ConstVector rhs, ConstVector x) const final; }; #include "residualGive.inl" diff --git a/include/Residual/ResidualGive/residualGive.inl b/include/Residual/ResidualGive/residualGive.inl index f6680d4e..fa4f5fc0 100644 --- a/include/Residual/ResidualGive/residualGive.inl +++ b/include/Residual/ResidualGive/residualGive.inl @@ -10,15 +10,12 @@ ResidualGive::ResidualGive(const PolarGrid& /* ------------------ */ /* result = rhs - A*x */ template -void ResidualGive::computeResidual(HostVector h_result, HostConstVector h_rhs, - HostConstVector h_x) const +void ResidualGive::computeResidual(Vector result, ConstVector rhs, + ConstVector x) const { - assert(h_result.size() == h_x.size()); + assert(result.size() == x.size()); - applySystemOperator(h_result, h_x); - - auto rhs = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_rhs); - auto result = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_result); + applySystemOperator(result, x); // Subtract A*x from rhs to get the residual. const int n = result.size(); @@ -28,6 +25,4 @@ void ResidualGive::computeResidual(HostVector h_result, KOKKOS_LAMBDA(const int i) { result[i] = rhs[i] - result[i]; }); Kokkos::fence(); - - Kokkos::deep_copy(h_result, result); } diff --git a/include/Residual/ResidualTake/applyATake.inl b/include/Residual/ResidualTake/applyATake.inl index 5e4d9665..9442515d 100644 --- a/include/Residual/ResidualTake/applyATake.inl +++ b/include/Residual/ResidualTake/applyATake.inl @@ -64,13 +64,10 @@ node_apply_a_take(const int i_r, const int i_theta, const PolarGrid -void ResidualTake::applySystemOperator(HostVector h_result, HostConstVector h_x) const +void ResidualTake::applySystemOperator(Vector result, ConstVector x) const { using residual_take::node_apply_a_take; - auto x = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_x); - auto result = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_result); - assert(result.size() == x.size()); const PolarGrid& grid = Residual::grid_; @@ -113,6 +110,4 @@ void ResidualTake::applySystemOperator(HostVector h_resu }); Kokkos::fence(); - - Kokkos::deep_copy(h_result, result); } diff --git a/include/Residual/ResidualTake/residualTake.h b/include/Residual/ResidualTake/residualTake.h index 8af3bd9f..403310ce 100644 --- a/include/Residual/ResidualTake/residualTake.h +++ b/include/Residual/ResidualTake/residualTake.h @@ -14,8 +14,8 @@ class ResidualTake : public Residual KOKKOS_DEFAULTED_FUNCTION ResidualTake(const ResidualTake&) = default; KOKKOS_DEFAULTED_FUNCTION ~ResidualTake() override = default; - void applySystemOperator(HostVector result, HostConstVector x) const final; - void computeResidual(HostVector result, HostConstVector rhs, HostConstVector x) const final; + void applySystemOperator(Vector result, ConstVector x) const final; + void computeResidual(Vector result, ConstVector rhs, ConstVector x) const final; }; #include "residualTake.inl" diff --git a/include/Residual/ResidualTake/residualTake.inl b/include/Residual/ResidualTake/residualTake.inl index 8f598d67..cf9996dd 100644 --- a/include/Residual/ResidualTake/residualTake.inl +++ b/include/Residual/ResidualTake/residualTake.inl @@ -10,15 +10,12 @@ ResidualTake::ResidualTake(const PolarGrid& /* ------------------ */ /* result = rhs - A*x */ template -void ResidualTake::computeResidual(HostVector h_result, HostConstVector h_rhs, - HostConstVector h_x) const +void ResidualTake::computeResidual(Vector result, ConstVector rhs, + ConstVector x) const { - assert(h_result.size() == h_x.size()); + assert(result.size() == x.size()); - applySystemOperator(h_result, h_x); - - auto rhs = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_rhs); - auto result = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_result); + applySystemOperator(result, x); // Subtract A*x from rhs to get the residual. const int n = result.size(); @@ -28,6 +25,4 @@ void ResidualTake::computeResidual(HostVector h_result, KOKKOS_LAMBDA(const int i) { result[i] = rhs[i] - result[i]; }); Kokkos::fence(); - - Kokkos::deep_copy(h_result, result); } diff --git a/include/Residual/residual.h b/include/Residual/residual.h index fcb04d00..b2815606 100644 --- a/include/Residual/residual.h +++ b/include/Residual/residual.h @@ -25,9 +25,9 @@ class Residual } virtual ~Residual() = default; - virtual void applySystemOperator(HostVector result, HostConstVector x) const = 0; - virtual void computeResidual(HostVector result, HostConstVector rhs, - HostConstVector x) const = 0; + virtual void applySystemOperator(Vector result, ConstVector x) const = 0; + virtual void computeResidual(Vector result, ConstVector rhs, + ConstVector x) const = 0; protected: /* ------------------- */ diff --git a/include/Smoother/SmootherGive/smootherGive.h b/include/Smoother/SmootherGive/smootherGive.h index 6237a14c..20ed82a0 100644 --- a/include/Smoother/SmootherGive/smootherGive.h +++ b/include/Smoother/SmootherGive/smootherGive.h @@ -61,7 +61,7 @@ class SmootherGive : public Smoother // Parallel implementation using OpenMP: // Scedule every 2nd/4th line in parallel to avoid race conditions arising from the A-Give distribution. // Sceduling every 3rd line in parallel would also be possible, but is less natural for the 2 coloring. - void smoothing(HostVector x, HostConstVector rhs, HostVector temp) final; + void smoothing(Vector x, ConstVector rhs, Vector temp) final; private: /* ------------------- */ diff --git a/include/Smoother/SmootherGive/smootherGive.inl b/include/Smoother/SmootherGive/smootherGive.inl index 6747cc95..76b718a8 100644 --- a/include/Smoother/SmootherGive/smootherGive.inl +++ b/include/Smoother/SmootherGive/smootherGive.inl @@ -41,13 +41,9 @@ SmootherGive::SmootherGive(const PolarGrid& // - The system is then solved in-place in temp, and the results // are copied back to x. template -void SmootherGive::smoothing(HostVector h_x, HostConstVector h_rhs, - HostVector h_temp) +void SmootherGive::smoothing(Vector x, ConstVector rhs, + Vector temp) { - auto x = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_x); - auto rhs = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_rhs); - auto temp = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_temp); - assert(x.size() == rhs.size()); assert(temp.size() == rhs.size()); @@ -64,7 +60,4 @@ void SmootherGive::smoothing(HostVector h_x, HostConstVe applyAscOrthoWhiteRadialSection(x, rhs, temp); solveWhiteRadialSection(x, temp); - - Kokkos::deep_copy(h_x, x); - Kokkos::deep_copy(h_temp, temp); } diff --git a/include/Smoother/SmootherTake/smootherTake.h b/include/Smoother/SmootherTake/smootherTake.h index 47ae4c1b..3321ae65 100644 --- a/include/Smoother/SmootherTake/smootherTake.h +++ b/include/Smoother/SmootherTake/smootherTake.h @@ -61,7 +61,7 @@ class SmootherTake : public Smoother // Performs one full coupled smoothing sweep: // BC -> WC -> BR -> WR // using temp as RHS workspace. - void smoothing(HostVector x, HostConstVector rhs, HostVector temp) final; + void smoothing(Vector x, ConstVector rhs, Vector temp) final; private: /* ------------------- */ diff --git a/include/Smoother/SmootherTake/smootherTake.inl b/include/Smoother/SmootherTake/smootherTake.inl index 349a2bf1..dcc99eaa 100644 --- a/include/Smoother/SmootherTake/smootherTake.inl +++ b/include/Smoother/SmootherTake/smootherTake.inl @@ -41,13 +41,9 @@ SmootherTake::SmootherTake(const PolarGrid& // - The system is then solved in-place in temp, and the results // are copied back to x. template -void SmootherTake::smoothing(HostVector h_x, HostConstVector h_rhs, - HostVector h_temp) +void SmootherTake::smoothing(Vector x, ConstVector rhs, + Vector temp) { - auto x = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_x); - auto rhs = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_rhs); - auto temp = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_temp); - assert(x.size() == rhs.size()); assert(temp.size() == rhs.size()); @@ -62,7 +58,4 @@ void SmootherTake::smoothing(HostVector h_x, HostConstVe applyAscOrthoWhiteRadialSection(x, rhs, temp); solveWhiteRadialSection(x, temp); - - Kokkos::deep_copy(h_x, x); - Kokkos::deep_copy(h_temp, temp); } diff --git a/include/Smoother/smoother.h b/include/Smoother/smoother.h index 41e2e017..092df381 100644 --- a/include/Smoother/smoother.h +++ b/include/Smoother/smoother.h @@ -36,7 +36,7 @@ class Smoother KOKKOS_DEFAULTED_FUNCTION Smoother(const Smoother&) = default; virtual ~Smoother() = default; - virtual void smoothing(HostVector x, HostConstVector rhs, HostVector temp) = 0; + virtual void smoothing(Vector x, ConstVector rhs, Vector temp) = 0; protected: const PolarGrid grid_; diff --git a/src/Interpolation/injection.cpp b/src/Interpolation/injection.cpp index 1895ee94..ce569f19 100644 --- a/src/Interpolation/injection.cpp +++ b/src/Interpolation/injection.cpp @@ -4,10 +4,9 @@ using namespace gmgpolar; /* Remark: This injection is not scaled. */ static KOKKOS_INLINE_FUNCTION void coarseNodeInjection(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_fine = i_r_coarse * 2; const int i_theta_fine = i_theta_coarse * 2; @@ -18,9 +17,9 @@ static KOKKOS_INLINE_FUNCTION void coarseNodeInjection(const int i_r_coarse, con coarse_result[coarse_index] = fine_values[fine_index]; } -void Interpolation::applyInjection(const PolarGrid& fine_grid, - const PolarGrid& coarse_grid, HostVector coarse_result, - HostConstVector fine_values) const +void Interpolation::applyInjection(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()); @@ -31,7 +30,7 @@ void Interpolation::applyInjection(const PolarGrid& fine_grid // The For loop matches circular access pattern */ Kokkos::parallel_for( "Interpolation: Injection (Circular)", - Kokkos::MDRangePolicy>( // 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 ), @@ -43,7 +42,7 @@ void Interpolation::applyInjection(const PolarGrid& fine_grid /* For loop matches radial access pattern */ Kokkos::parallel_for( "Interpolation: Injection (Radial)", - Kokkos::MDRangePolicy>( // 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/src/Interpolation/prolongation.cpp b/src/Interpolation/prolongation.cpp index b630a60c..c6a6846b 100644 --- a/src/Interpolation/prolongation.cpp +++ b/src/Interpolation/prolongation.cpp @@ -53,10 +53,10 @@ using namespace gmgpolar; */ static KOKKOS_INLINE_FUNCTION void fineNodeProlongation(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; @@ -109,9 +109,9 @@ static KOKKOS_INLINE_FUNCTION void fineNodeProlongation(const int i_r, const int } } -void Interpolation::applyProlongation(const PolarGrid& coarse_grid, - const PolarGrid& fine_grid, HostVector fine_result, - HostConstVector coarse_values) const +void Interpolation::applyProlongation(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()); @@ -122,7 +122,7 @@ void Interpolation::applyProlongation(const PolarGrid& coarse // The For loop matches circular access pattern */ Kokkos::parallel_for( "Interpolation: Prolongation (Circular)", - Kokkos::MDRangePolicy>( // 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 ), @@ -134,7 +134,7 @@ void Interpolation::applyProlongation(const PolarGrid& coarse /* For loop matches radial access pattern */ Kokkos::parallel_for( "Interpolation: Prolongation (Radial)", - Kokkos::MDRangePolicy>( // 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/restriction.cpp b/src/Interpolation/restriction.cpp index 98d1d739..409e3d98 100644 --- a/src/Interpolation/restriction.cpp +++ b/src/Interpolation/restriction.cpp @@ -24,10 +24,10 @@ using namespace gmgpolar; */ static KOKKOS_INLINE_FUNCTION void coarseNodeRestriction(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; @@ -69,9 +69,9 @@ static KOKKOS_INLINE_FUNCTION void coarseNodeRestriction(const int i_r_coarse, c coarse_result[coarse_grid.index(i_r_coarse, i_theta_coarse)] = value; } -void Interpolation::applyRestriction(const PolarGrid& fine_grid, - const PolarGrid& coarse_grid, HostVector coarse_result, - HostConstVector fine_values) const +void Interpolation::applyRestriction(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()); @@ -82,7 +82,7 @@ void Interpolation::applyRestriction(const PolarGrid& fine_gr // The For loop matches circular access pattern */ Kokkos::parallel_for( "Interpolation: Restriction (Circular)", - Kokkos::MDRangePolicy>( // 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 ), @@ -94,7 +94,7 @@ void Interpolation::applyRestriction(const PolarGrid& fine_gr /* For loop matches radial access pattern */ Kokkos::parallel_for( "Interpolation: Restriction (Radial)", - Kokkos::MDRangePolicy>( // 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/src/PolarGrid/polargrid.cpp b/src/PolarGrid/polargrid.cpp index 5337faa2..c458773d 100644 --- a/src/PolarGrid/polargrid.cpp +++ b/src/PolarGrid/polargrid.cpp @@ -214,14 +214,15 @@ void PolarGrid::initializeLineSplitting(std::optional split } } else { + auto h_radius = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, radii_); number_smoother_circles_ = 2; /* We assume numberSmootherCircles_ >= 2 in the further implementation */ for (int i_r = 2; i_r < nr() - 2; i_r++) { /* We assume lengthRadialSmoother_ >= 3 in the further implementation */ const double uniform_theta_k = (2 * M_PI) / ntheta(); - const double radial_dist_h = radius(i_r + 1) - radius(i_r); + const double radial_dist_h = h_radius[i_r + 1] - h_radius[i_r]; const double q = uniform_theta_k / radial_dist_h; - if (q * radius(i_r) > 1.0) { + if (q * h_radius[i_r] > 1.0) { number_smoother_circles_ = i_r; break; } diff --git a/tests/DirectSolver/directSolver.cpp b/tests/DirectSolver/directSolver.cpp index 34464e1b..4ebf8744 100644 --- a/tests/DirectSolver/directSolver.cpp +++ b/tests/DirectSolver/directSolver.cpp @@ -71,20 +71,23 @@ TEST(DirectSolverTest, directSolver_DirBC_Interior) DirectSolverTake directSolverGive_operator(level.grid(), level.levelCache(), DirBC_Interior); DirectSolverGive directSolverTake_operator(level.grid(), level.levelCache(), DirBC_Interior); - HostVector solution_Give = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector solution_Give = generate_random_sample_data(level.grid(), 69); directSolverGive_operator.solveInPlace(solution_Give); - HostVector solution_Take = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector solution_Take = generate_random_sample_data(level.grid(), 69); directSolverTake_operator.solveInPlace(solution_Take); + auto h_solution_Give = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), solution_Give); + auto h_solution_Take = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), solution_Take); + ASSERT_EQ(solution_Give.size(), solution_Take.size()); for (uint index = 0; index < solution_Give.size(); index++) { int i_r, i_theta; level.grid().multiIndex(index, i_r, i_theta); if (i_r == 0 && !DirBC_Interior) - ASSERT_NEAR(solution_Give(index), solution_Take(index), 1e-11); + ASSERT_NEAR(h_solution_Give(index), h_solution_Take(index), 1e-11); else - ASSERT_NEAR(solution_Give(index), solution_Take(index), 1e-11); + ASSERT_NEAR(h_solution_Give(index), h_solution_Take(index), 1e-11); } } @@ -120,20 +123,23 @@ TEST(DirectSolverTest, directSolver_AcrossOrigin) DirectSolverGive directSolverGive_operator(level.grid(), level.levelCache(), DirBC_Interior); DirectSolverTake directSolverTake_operator(level.grid(), level.levelCache(), DirBC_Interior); - HostVector solution_Give = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector solution_Give = generate_random_sample_data(level.grid(), 69); directSolverGive_operator.solveInPlace(solution_Give); - HostVector solution_Take = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector solution_Take = generate_random_sample_data(level.grid(), 69); directSolverTake_operator.solveInPlace(solution_Take); + auto h_solution_Give = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), solution_Give); + auto h_solution_Take = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), solution_Take); + ASSERT_EQ(solution_Give.size(), solution_Take.size()); for (uint index = 0; index < solution_Give.size(); index++) { int i_r, i_theta; level.grid().multiIndex(index, i_r, i_theta); if (i_r == 0 && !DirBC_Interior) - ASSERT_NEAR(solution_Give(index), solution_Take(index), 1e-8); + ASSERT_NEAR(h_solution_Give(index), h_solution_Take(index), 1e-8); else - ASSERT_NEAR(solution_Give(index), solution_Take(index), 1e-8); + ASSERT_NEAR(h_solution_Give(index), h_solution_Take(index), 1e-8); } } @@ -171,17 +177,17 @@ TEST(DirectSolverTest_CircularGeometry, DirectSolverDirBC_Interior_CircularGeome DirectSolverGive solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-12); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-12); } TEST(DirectSolverTest_CircularGeometry, DirectSolverAcrossOrigin_CircularGeometry) @@ -212,17 +218,17 @@ TEST(DirectSolverTest_CircularGeometry, DirectSolverAcrossOrigin_CircularGeometr DirectSolverGive solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-7); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-8); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-8); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-7); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-8); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-8); } /* --------- */ @@ -259,17 +265,17 @@ TEST(DirectSolverTest_ShafranovGeometry, DirectSolverDirBC_Interior_ShafranovGeo DirectSolverGive solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-12); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-12); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-12); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-12); } TEST(DirectSolverTest_ShafranovGeometry, DirectSolverAcrossOrigin_ShafranovGeometry) @@ -302,17 +308,17 @@ TEST(DirectSolverTest_ShafranovGeometry, DirectSolverAcrossOrigin_ShafranovGeome DirectSolverGive solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-7); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-8); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-8); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-7); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-8); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-8); } /* ------ */ @@ -349,17 +355,17 @@ TEST(DirectSolverTest_CzarnyGeometry, DirectSolverDirBC_Interior_CzarnyGeometry) DirectSolverGive solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-12); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-12); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-12); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-12); } TEST(DirectSolverTest_CzarnyGeometry, DirectSolverAcrossOrigin_CzarnyGeometry) @@ -392,17 +398,17 @@ TEST(DirectSolverTest_CzarnyGeometry, DirectSolverAcrossOrigin_CzarnyGeometry) DirectSolverGive solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-7); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-8); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-8); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-7); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-8); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-8); } /* ------ */ @@ -437,17 +443,17 @@ TEST(DirectSolverTest_CulhamGeometry, DirectSolverDirBC_Interior_CulhamGeometry) DirectSolverGive solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-12); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-12); } TEST(DirectSolverTest_CulhamGeometry, DirectSolverAcrossOrigin_CulhamGeometry) @@ -478,17 +484,17 @@ TEST(DirectSolverTest_CulhamGeometry, DirectSolverAcrossOrigin_CulhamGeometry) DirectSolverGive solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-7); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-7); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-8); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-7); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-7); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-8); } /* We adjust the PolarGrid to increase the precision */ @@ -531,17 +537,17 @@ TEST(DirectSolverTest_CircularGeometry, DirectSolverAcrossOriginHigherPrecision_ DirectSolverGive solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-9); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-10); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-10); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-9); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-10); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-10); } TEST(DirectSolverTest_CircularGeometry, DirectSolverAcrossOriginHigherPrecision2_CircularGeometry) @@ -572,17 +578,17 @@ TEST(DirectSolverTest_CircularGeometry, DirectSolverAcrossOriginHigherPrecision2 DirectSolverGive solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-12); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-12); } /* Same test now using Take */ @@ -614,17 +620,17 @@ TEST(DirectSolverTakeTest_CircularGeometry, DirectSolverDirBC_Interior_CircularG DirectSolverTake solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-12); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-12); } TEST(DirectSolverTakeTest_CircularGeometry, DirectSolverAcrossOrigin_CircularGeometry) @@ -655,17 +661,17 @@ TEST(DirectSolverTakeTest_CircularGeometry, DirectSolverAcrossOrigin_CircularGeo DirectSolverTake solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-7); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-8); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-8); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-7); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-8); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-8); } /* --------- */ @@ -702,17 +708,17 @@ TEST(DirectSolverTakeTest_ShafranovGeometry, DirectSolverDirBC_Interior_Shafrano DirectSolverTake solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-11); } TEST(DirectSolverTakeTest_ShafranovGeometry, DirectSolverAcrossOrigin_ShafranovGeometry) @@ -745,17 +751,17 @@ TEST(DirectSolverTakeTest_ShafranovGeometry, DirectSolverAcrossOrigin_ShafranovG DirectSolverTake solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-7); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-8); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-8); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-7); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-8); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-8); } /* ------ */ @@ -792,17 +798,17 @@ TEST(DirectSolverTakeTest_CzarnyGeometry, DirectSolverDirBC_Interior_CzarnyGeome DirectSolverTake solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-12); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-12); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-12); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-12); } TEST(DirectSolverTakeTest_CzarnyGeometry, DirectSolverAcrossOrigin_CzarnyGeometry) @@ -835,17 +841,17 @@ TEST(DirectSolverTakeTest_CzarnyGeometry, DirectSolverAcrossOrigin_CzarnyGeometr DirectSolverTake solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-7); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-8); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-8); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-7); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-8); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-8); } /* ------ */ @@ -880,17 +886,17 @@ TEST(DirectSolverTakeTest_CulhamGeometry, DirectSolverDirBC_Interior_CulhamGeome DirectSolverTake solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-11); } TEST(DirectSolverTakeTest_CulhamGeometry, DirectSolverAcrossOrigin_CulhamGeometry) @@ -921,17 +927,17 @@ TEST(DirectSolverTakeTest_CulhamGeometry, DirectSolverAcrossOrigin_CulhamGeometr DirectSolverTake solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-7); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-7); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-8); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-7); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-7); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-8); } TEST(DirectSolverTakeTest_CircularGeometry, DirectSolverAcrossOriginHigherPrecision_CircularGeometry) @@ -972,17 +978,17 @@ TEST(DirectSolverTakeTest_CircularGeometry, DirectSolverAcrossOriginHigherPrecis DirectSolverTake solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-10); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-10); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-10); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-10); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-10); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-10); } TEST(DirectSolverTakeTest_CircularGeometry, DirectSolverAcrossOriginHigherPrecision2_CircularGeometry) @@ -1013,15 +1019,15 @@ TEST(DirectSolverTakeTest_CircularGeometry, DirectSolverAcrossOriginHigherPrecis DirectSolverTake solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector solution("sol", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector solution("sol", rhs.size()); Kokkos::deep_copy(solution, rhs); solver_op.solveInPlace(solution); - HostVector residuum("residuum", level.grid().numberOfNodes()); + Vector residuum("residuum", level.grid().numberOfNodes()); residual_op.computeResidual(residuum, rhs, solution); - ASSERT_NEAR(l1_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(l2_norm(HostConstVector(residuum)), 0.0, 1e-11); - ASSERT_NEAR(infinity_norm(HostConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(l1_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(l2_norm(ConstVector(residuum)), 0.0, 1e-11); + ASSERT_NEAR(infinity_norm(ConstVector(residuum)), 0.0, 1e-11); } diff --git a/tests/ExtrapolatedSmoother/extrapolated_smoother.cpp b/tests/ExtrapolatedSmoother/extrapolated_smoother.cpp index 61e4f827..3644d0db 100644 --- a/tests/ExtrapolatedSmoother/extrapolated_smoother.cpp +++ b/tests/ExtrapolatedSmoother/extrapolated_smoother.cpp @@ -41,7 +41,7 @@ TEST(ExtrapolatedSmootherTest, extrapolatedSmoother_DirBC_Interior) using DomainGeometryType = CzarnyGeometry; DomainGeometryType domain_geometry(Rmax, kappa_eps, delta_e); - auto grid = std::make_unique>(radii, angles); + auto tmp_grid = std::make_unique>(radii, angles); double alpha_jump = 0.678 * Rmax; using DensityProfileCoefficientsType = ZoniShiftedCoefficients; @@ -54,34 +54,38 @@ TEST(ExtrapolatedSmootherTest, extrapolatedSmoother_DirBC_Interior) bool cache_domain_geometry = true; auto levelCache = std::make_unique>( - *grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); - Level level(0, std::move(grid), std::move(levelCache), + *tmp_grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); + Level level(0, std::move(tmp_grid), std::move(levelCache), ExtrapolationType::NONE, 0); + const PolarGrid& grid = level.grid(); ExtrapolatedSmootherGive smootherGive_operator(level.grid(), level.levelCache(), DirBC_Interior); ExtrapolatedSmootherTake smootherTake_operator(level.grid(), level.levelCache(), DirBC_Interior); - HostVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector rhs = generate_random_sample_data(level.grid(), 69); HostVector start = generate_random_sample_data(PolarGrid(level.grid()), 24); - HostVector temp = generate_random_sample_data(PolarGrid(level.grid()), 8); + Vector temp = generate_random_sample_data(level.grid(), 8); - HostVector solution_Give("solution_Give", start.size()); + Vector solution_Give("solution_Give", start.size()); Kokkos::deep_copy(solution_Give, start); smootherGive_operator.extrapolatedSmoothing(solution_Give, rhs, temp); - HostVector solution_Take("solution_Take", start.size()); + Vector solution_Take("solution_Take", start.size()); Kokkos::deep_copy(solution_Take, start); smootherTake_operator.extrapolatedSmoothing(solution_Take, rhs, temp); - ASSERT_EQ(solution_Give.size(), solution_Take.size()); - for (uint index = 0; index < solution_Give.size(); index++) { + auto h_solution_Give = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, solution_Give); + auto h_solution_Take = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, solution_Take); + + ASSERT_EQ(h_solution_Give.size(), h_solution_Take.size()); + for (uint index = 0; index < h_solution_Give.size(); index++) { int i_r, i_theta; level.grid().multiIndex(index, i_r, i_theta); if (i_r == 0 && !DirBC_Interior) - ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11); + ASSERT_NEAR(h_solution_Give[index], h_solution_Take[index], 1e-11); else - ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11); + ASSERT_NEAR(h_solution_Give[index], h_solution_Take[index], 1e-11); } } @@ -98,7 +102,7 @@ TEST(ExtrapolatedSmootherTest, extrapolatedSmoother_AcossOrigin) using DomainGeometryType = CzarnyGeometry; DomainGeometryType domain_geometry(Rmax, kappa_eps, delta_e); - auto grid = std::make_unique>(radii, angles); + auto tmp_grid = std::make_unique>(radii, angles); double alpha_jump = 0.678 * Rmax; using DensityProfileCoefficientsType = ZoniShiftedCoefficients; @@ -111,40 +115,43 @@ TEST(ExtrapolatedSmootherTest, extrapolatedSmoother_AcossOrigin) bool cache_domain_geometry = true; auto levelCache = std::make_unique>( - *grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); - Level level(0, std::move(grid), std::move(levelCache), + *tmp_grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); + Level level(0, std::move(tmp_grid), std::move(levelCache), ExtrapolationType::NONE, 0); ExtrapolatedSmootherGive smootherGive_operator(level.grid(), level.levelCache(), DirBC_Interior); ExtrapolatedSmootherTake smootherTake_operator(level.grid(), level.levelCache(), DirBC_Interior); - HostVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector rhs = generate_random_sample_data(level.grid(), 69); HostVector start = generate_random_sample_data(PolarGrid(level.grid()), 24); - HostVector temp = generate_random_sample_data(PolarGrid(level.grid()), 8); + Vector temp = generate_random_sample_data(level.grid(), 8); - HostVector solution_Give("solution_Give", start.size()); + Vector solution_Give("solution_Give", start.size()); Kokkos::deep_copy(solution_Give, start); smootherGive_operator.extrapolatedSmoothing(solution_Give, rhs, temp); - HostVector solution_Take("solution_Take", start.size()); + Vector solution_Take("solution_Take", start.size()); Kokkos::deep_copy(solution_Take, start); smootherTake_operator.extrapolatedSmoothing(solution_Take, rhs, temp); - ASSERT_EQ(solution_Give.size(), solution_Take.size()); - for (uint index = 0; index < solution_Give.size(); index++) { + auto h_solution_Give = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, solution_Give); + auto h_solution_Take = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, solution_Take); + + ASSERT_EQ(h_solution_Give.size(), h_solution_Take.size()); + for (uint index = 0; index < h_solution_Give.size(); index++) { int i_r, i_theta; level.grid().multiIndex(index, i_r, i_theta); if (i_r == 0 && !DirBC_Interior) - ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-8); + ASSERT_NEAR(h_solution_Give[index], h_solution_Take[index], 1e-8); else - ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11); + ASSERT_NEAR(h_solution_Give[index], h_solution_Take[index], 1e-11); } } /* Test 2/2: */ /* Does the smoother converge to the DirectSolverGive solution? */ -TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherDirBC_Interior) +void ExtrapolatedSmootherTest_ExtrapolatedSmootherDirBC_Interior() { std::vector radii = {1e-5, 0.2, 0.25, 0.5, 0.8, 0.9, 0.95, 1.2, 1.3}; std::vector angles = { @@ -157,7 +164,7 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherDirBC_Interior) using DomainGeometryType = CzarnyGeometry; DomainGeometryType domain_geometry(Rmax, kappa_eps, delta_e); - auto grid = std::make_unique>(radii, angles); + auto tmp_grid = std::make_unique>(radii, angles); double alpha_jump = 0.678 * Rmax; using DensityProfileCoefficientsType = ZoniShiftedCoefficients; @@ -169,32 +176,39 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherDirBC_Interior) bool cache_domain_geometry = false; auto levelCache = std::make_unique>( - *grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); - Level level(0, std::move(grid), std::move(levelCache), + *tmp_grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); + Level level(0, std::move(tmp_grid), std::move(levelCache), ExtrapolationType::IMPLICIT_EXTRAPOLATION, 0); + const PolarGrid& grid = level.grid(); DirectSolverGive solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); ExtrapolatedSmootherGive extrapolated_smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - HostVector error("error", level.grid().numberOfNodes()); - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); - for (int i_r = 0; i_r < level.grid().nr(); i_r += 2) { - for (int i_theta = 0; i_theta < level.grid().ntheta(); i_theta += 2) { - smoother_solution[level.grid().index(i_r, i_theta)] = discrete_solution[level.grid().index(i_r, i_theta)]; - } - } - - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Vector temp("temp", level.grid().numberOfNodes()); + Vector error("error", level.grid().numberOfNodes()); + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); + Kokkos::parallel_for( + "Fill smoother_solution", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {(grid.nr() + 1) / 2, (grid.ntheta() + 1) / 2} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int i_r_half, const int i_theta_half) { + const int i_r = 2 * i_r_half; + const int i_theta = 2 * i_theta_half; + smoother_solution[grid.index(i_r, i_theta)] = discrete_solution[grid.index(i_r, i_theta)]; + }); + + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); int iterations = 0; @@ -202,13 +216,12 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherDirBC_Interior) const int max_iterations = 10000; const double precision = 1e-12; - while (infinity_norm(HostConstVector(error)) > precision) { + while (infinity_norm(ConstVector(error)) > precision) { extrapolated_smoother_op.extrapolatedSmoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); iterations++; if (iterations >= max_iterations) { @@ -222,10 +235,14 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherDirBC_Interior) ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 300); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); +} +TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherDirBC_Interior) +{ + ExtrapolatedSmootherTest_ExtrapolatedSmootherDirBC_Interior(); } -TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherAcrossOrigin) +void ExtrapolatedSmootherTest_ExtrapolatedSmootherAcrossOrigin() { std::vector radii = {1e-5, 0.2, 0.25, 0.5, 0.8, 0.9, 0.95, 1.2, 1.3}; std::vector angles = { @@ -238,7 +255,7 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherAcrossOrigin) using DomainGeometryType = CzarnyGeometry; DomainGeometryType domain_geometry(Rmax, kappa_eps, delta_e); - auto grid = std::make_unique>(radii, angles); + auto tmp_grid = std::make_unique>(radii, angles); double alpha_jump = 0.678 * Rmax; using DensityProfileCoefficientsType = ZoniShiftedCoefficients; @@ -250,34 +267,41 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherAcrossOrigin) bool cache_domain_geometry = false; auto levelCache = std::make_unique>( - *grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); - Level level(0, std::move(grid), std::move(levelCache), + *tmp_grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); + Level level(0, std::move(tmp_grid), std::move(levelCache), ExtrapolationType::IMPLICIT_EXTRAPOLATION, 0); + const PolarGrid& grid = level.grid(); DirectSolverGive solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); ExtrapolatedSmootherGive extrapolated_smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - - HostVector error("error", level.grid().numberOfNodes()); - - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); - for (int i_r = 0; i_r < level.grid().nr(); i_r += 2) { - for (int i_theta = 0; i_theta < level.grid().ntheta(); i_theta += 2) { - smoother_solution[level.grid().index(i_r, i_theta)] = discrete_solution[level.grid().index(i_r, i_theta)]; - } - } - - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Vector temp("temp", level.grid().numberOfNodes()); + + Vector error("error", level.grid().numberOfNodes()); + + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); + Kokkos::parallel_for( + "Fill smoother_solution", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {(grid.nr() + 1) / 2, (grid.ntheta() + 1) / 2} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int i_r_half, const int i_theta_half) { + const int i_r = 2 * i_r_half; + const int i_theta = 2 * i_theta_half; + smoother_solution[grid.index(i_r, i_theta)] = discrete_solution[grid.index(i_r, i_theta)]; + }); + + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); int iterations = 0; @@ -285,13 +309,12 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherAcrossOrigin) const int max_iterations = 10000; const double precision = 1e-8; - while (infinity_norm(HostConstVector(error)) > precision) { + while (infinity_norm(ConstVector(error)) > precision) { extrapolated_smoother_op.extrapolatedSmoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); iterations++; if (iterations >= max_iterations) { @@ -305,10 +328,14 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherAcrossOrigin) ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 600); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); +} +TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherAcrossOrigin) +{ + ExtrapolatedSmootherTest_ExtrapolatedSmootherAcrossOrigin(); } -TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherDirBC_Interior_SmallestGrid) +void ExtrapolatedSmootherTest_ExtrapolatedSmootherDirBC_Interior_SmallestGrid() { std::vector radii = {1e-5, 0.2, 0.4, 0.45, 0.9, 1.2, 1.3}; std::vector angles = {0, M_PI / 8, M_PI, M_PI + M_PI / 8, M_PI + M_PI}; @@ -320,7 +347,7 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherDirBC_Interior_SmallestGrid) using DomainGeometryType = CzarnyGeometry; DomainGeometryType domain_geometry(Rmax, kappa_eps, delta_e); - auto grid = std::make_unique>(radii, angles); + auto tmp_grid = std::make_unique>(radii, angles); double alpha_jump = 0.678 * Rmax; using DensityProfileCoefficientsType = ZoniShiftedCoefficients; @@ -332,34 +359,41 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherDirBC_Interior_SmallestGrid) bool cache_domain_geometry = false; auto levelCache = std::make_unique>( - *grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); - Level level(0, std::move(grid), std::move(levelCache), + *tmp_grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); + Level level(0, std::move(tmp_grid), std::move(levelCache), ExtrapolationType::IMPLICIT_EXTRAPOLATION, 0); + const PolarGrid& grid = level.grid(); DirectSolverGive solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); ExtrapolatedSmootherGive extrapolated_smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - - HostVector error("error", level.grid().numberOfNodes()); - - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); - for (int i_r = 0; i_r < level.grid().nr(); i_r += 2) { - for (int i_theta = 0; i_theta < level.grid().ntheta(); i_theta += 2) { - smoother_solution[level.grid().index(i_r, i_theta)] = discrete_solution[level.grid().index(i_r, i_theta)]; - } - } - - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Vector temp("temp", level.grid().numberOfNodes()); + + Vector error("error", level.grid().numberOfNodes()); + + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); + Kokkos::parallel_for( + "Fill smoother_solution", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {(grid.nr() + 1) / 2, (grid.ntheta() + 1) / 2} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int i_r_half, const int i_theta_half) { + const int i_r = 2 * i_r_half; + const int i_theta = 2 * i_theta_half; + smoother_solution[grid.index(i_r, i_theta)] = discrete_solution[grid.index(i_r, i_theta)]; + }); + + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); int iterations = 0; @@ -367,13 +401,12 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherDirBC_Interior_SmallestGrid) const int max_iterations = 10000; double precision = 1e-12; - while (infinity_norm(HostConstVector(error)) > precision) { + while (infinity_norm(ConstVector(error)) > precision) { extrapolated_smoother_op.extrapolatedSmoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); iterations++; if (iterations >= max_iterations) { @@ -387,10 +420,14 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherDirBC_Interior_SmallestGrid) ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 200); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); +} +TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherDirBC_Interior_SmallestGrid) +{ + ExtrapolatedSmootherTest_ExtrapolatedSmootherDirBC_Interior_SmallestGrid(); } -TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherAcrossOrigin_SmallestGrid) +void ExtrapolatedSmootherTest_ExtrapolatedSmootherAcrossOrigin_SmallestGrid() { std::vector radii = {1e-5, 0.2, 0.4, 0.45, 0.9, 1.2, 1.3}; std::vector angles = {0, M_PI / 8, M_PI, M_PI + M_PI / 8, M_PI + M_PI}; @@ -402,7 +439,7 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherAcrossOrigin_SmallestGrid) using DomainGeometryType = CzarnyGeometry; DomainGeometryType domain_geometry(Rmax, kappa_eps, delta_e); - auto grid = std::make_unique>(radii, angles); + auto tmp_grid = std::make_unique>(radii, angles); double alpha_jump = 0.678 * Rmax; using DensityProfileCoefficientsType = ZoniShiftedCoefficients; @@ -414,32 +451,39 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherAcrossOrigin_SmallestGrid) bool cache_domain_geometry = false; auto levelCache = std::make_unique>( - *grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); - Level level(0, std::move(grid), std::move(levelCache), + *tmp_grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); + Level level(0, std::move(tmp_grid), std::move(levelCache), ExtrapolationType::IMPLICIT_EXTRAPOLATION, 0); + const PolarGrid& grid = level.grid(); DirectSolverGive solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); ExtrapolatedSmootherGive extrapolated_smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - HostVector error("error", level.grid().numberOfNodes()); - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); - for (int i_r = 0; i_r < level.grid().nr(); i_r += 2) { - for (int i_theta = 0; i_theta < level.grid().ntheta(); i_theta += 2) { - smoother_solution[level.grid().index(i_r, i_theta)] = discrete_solution[level.grid().index(i_r, i_theta)]; - } - } - - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Vector temp("temp", level.grid().numberOfNodes()); + Vector error("error", level.grid().numberOfNodes()); + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); + Kokkos::parallel_for( + "Fill smoother_solution", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {(grid.nr() + 1) / 2, (grid.ntheta() + 1) / 2} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int i_r_half, const int i_theta_half) { + const int i_r = 2 * i_r_half; + const int i_theta = 2 * i_theta_half; + smoother_solution[grid.index(i_r, i_theta)] = discrete_solution[grid.index(i_r, i_theta)]; + }); + + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); int iterations = 0; @@ -447,13 +491,12 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherAcrossOrigin_SmallestGrid) const int max_iterations = 10000; const double precision = 1e-8; - while (infinity_norm(HostConstVector(error)) > 1e-8) { + while (infinity_norm(ConstVector(error)) > 1e-8) { extrapolated_smoother_op.extrapolatedSmoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); iterations++; if (iterations >= max_iterations) { @@ -467,12 +510,16 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherAcrossOrigin_SmallestGrid) ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 150); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); +} +TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherAcrossOrigin_SmallestGrid) +{ + ExtrapolatedSmootherTest_ExtrapolatedSmootherAcrossOrigin_SmallestGrid(); } /* We now test using "Take" */ -TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeDirBC_Interior) +void ExtrapolatedSmootherTest_ExtrapolatedSmootherTakeDirBC_Interior() { std::vector radii = {1e-5, 0.2, 0.25, 0.5, 0.8, 0.9, 0.95, 1.2, 1.3}; std::vector angles = { @@ -485,7 +532,7 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeDirBC_Interior) using DomainGeometryType = CzarnyGeometry; DomainGeometryType domain_geometry(Rmax, kappa_eps, delta_e); - auto grid = std::make_unique>(radii, angles); + auto tmp_grid = std::make_unique>(radii, angles); double alpha_jump = 0.678 * Rmax; using DensityProfileCoefficientsType = ZoniShiftedCoefficients; @@ -497,32 +544,39 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeDirBC_Interior) bool cache_domain_geometry = true; auto levelCache = std::make_unique>( - *grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); - Level level(0, std::move(grid), std::move(levelCache), + *tmp_grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); + Level level(0, std::move(tmp_grid), std::move(levelCache), ExtrapolationType::IMPLICIT_EXTRAPOLATION, 0); + const PolarGrid& grid = level.grid(); DirectSolverTake solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualTake residual_op(level.grid(), level.levelCache(), DirBC_Interior); ExtrapolatedSmootherTake extrapolated_smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - HostVector error("error", level.grid().numberOfNodes()); - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); - for (int i_r = 0; i_r < level.grid().nr(); i_r += 2) { - for (int i_theta = 0; i_theta < level.grid().ntheta(); i_theta += 2) { - smoother_solution[level.grid().index(i_r, i_theta)] = discrete_solution[level.grid().index(i_r, i_theta)]; - } - } - - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Vector temp("temp", level.grid().numberOfNodes()); + Vector error("error", level.grid().numberOfNodes()); + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); + Kokkos::parallel_for( + "Fill smoother_solution", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {(grid.nr() + 1) / 2, (grid.ntheta() + 1) / 2} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int i_r_half, const int i_theta_half) { + const int i_r = 2 * i_r_half; + const int i_theta = 2 * i_theta_half; + smoother_solution[grid.index(i_r, i_theta)] = discrete_solution[grid.index(i_r, i_theta)]; + }); + + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); int iterations = 0; @@ -530,13 +584,12 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeDirBC_Interior) const int max_iterations = 10000; const double precision = 1e-12; - while (infinity_norm(HostConstVector(error)) > precision) { + while (infinity_norm(ConstVector(error)) > precision) { extrapolated_smoother_op.extrapolatedSmoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); iterations++; if (iterations >= max_iterations) { @@ -550,10 +603,14 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeDirBC_Interior) ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 300); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); +} +TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeDirBC_Interior) +{ + ExtrapolatedSmootherTest_ExtrapolatedSmootherTakeDirBC_Interior(); } -TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeAcrossOrigin) +void ExtrapolatedSmootherTest_ExtrapolatedSmootherTakeAcrossOrigin() { std::vector radii = {1e-5, 0.2, 0.25, 0.5, 0.8, 0.9, 0.95, 1.2, 1.3}; std::vector angles = { @@ -566,7 +623,7 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeAcrossOrigin) using DomainGeometryType = CzarnyGeometry; DomainGeometryType domain_geometry(Rmax, kappa_eps, delta_e); - auto grid = std::make_unique>(radii, angles); + auto tmp_grid = std::make_unique>(radii, angles); double alpha_jump = 0.678 * Rmax; using DensityProfileCoefficientsType = ZoniShiftedCoefficients; @@ -578,32 +635,39 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeAcrossOrigin) bool cache_domain_geometry = true; auto levelCache = std::make_unique>( - *grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); - Level level(0, std::move(grid), std::move(levelCache), + *tmp_grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); + Level level(0, std::move(tmp_grid), std::move(levelCache), ExtrapolationType::IMPLICIT_EXTRAPOLATION, 0); + const PolarGrid& grid = level.grid(); DirectSolverTake solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualTake residual_op(level.grid(), level.levelCache(), DirBC_Interior); ExtrapolatedSmootherTake extrapolated_smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - HostVector error("error", level.grid().numberOfNodes()); - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); - for (int i_r = 0; i_r < level.grid().nr(); i_r += 2) { - for (int i_theta = 0; i_theta < level.grid().ntheta(); i_theta += 2) { - smoother_solution[level.grid().index(i_r, i_theta)] = discrete_solution[level.grid().index(i_r, i_theta)]; - } - } - - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Vector temp("temp", level.grid().numberOfNodes()); + Vector error("error", level.grid().numberOfNodes()); + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); + Kokkos::parallel_for( + "Fill smoother_solution", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {(grid.nr() + 1) / 2, (grid.ntheta() + 1) / 2} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int i_r_half, const int i_theta_half) { + const int i_r = 2 * i_r_half; + const int i_theta = 2 * i_theta_half; + smoother_solution[grid.index(i_r, i_theta)] = discrete_solution[grid.index(i_r, i_theta)]; + }); + + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); int iterations = 0; @@ -611,13 +675,12 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeAcrossOrigin) const int max_iterations = 10000; const double precision = 1e-8; - while (infinity_norm(HostConstVector(error)) > precision) { + while (infinity_norm(ConstVector(error)) > precision) { extrapolated_smoother_op.extrapolatedSmoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); iterations++; if (iterations >= max_iterations) { @@ -631,10 +694,14 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeAcrossOrigin) ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 600); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); +} +TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeAcrossOrigin) +{ + ExtrapolatedSmootherTest_ExtrapolatedSmootherTakeAcrossOrigin(); } -TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeDirBC_Interior_SmallestGrid) +void ExtrapolatedSmootherTest_ExtrapolatedSmootherTakeDirBC_Interior_SmallestGrid() { std::vector radii = {1e-5, 0.2, 0.4, 0.45, 0.9, 1.2, 1.3}; std::vector angles = {0, M_PI / 8, M_PI, M_PI + M_PI / 8, M_PI + M_PI}; @@ -646,7 +713,7 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeDirBC_Interior_SmallestGr using DomainGeometryType = CzarnyGeometry; DomainGeometryType domain_geometry(Rmax, kappa_eps, delta_e); - auto grid = std::make_unique>(radii, angles); + auto tmp_grid = std::make_unique>(radii, angles); double alpha_jump = 0.678 * Rmax; using DensityProfileCoefficientsType = ZoniShiftedCoefficients; @@ -658,32 +725,39 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeDirBC_Interior_SmallestGr bool cache_domain_geometry = true; auto levelCache = std::make_unique>( - *grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); - Level level(0, std::move(grid), std::move(levelCache), + *tmp_grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); + Level level(0, std::move(tmp_grid), std::move(levelCache), ExtrapolationType::IMPLICIT_EXTRAPOLATION, 0); + const PolarGrid& grid = level.grid(); DirectSolverTake solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualTake residual_op(level.grid(), level.levelCache(), DirBC_Interior); ExtrapolatedSmootherTake extrapolated_smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - HostVector error("error", level.grid().numberOfNodes()); - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); - for (int i_r = 0; i_r < level.grid().nr(); i_r += 2) { - for (int i_theta = 0; i_theta < level.grid().ntheta(); i_theta += 2) { - smoother_solution[level.grid().index(i_r, i_theta)] = discrete_solution[level.grid().index(i_r, i_theta)]; - } - } - - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Vector temp("temp", level.grid().numberOfNodes()); + Vector error("error", level.grid().numberOfNodes()); + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); + Kokkos::parallel_for( + "Fill smoother_solution", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {(grid.nr() + 1) / 2, (grid.ntheta() + 1) / 2} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int i_r_half, const int i_theta_half) { + const int i_r = 2 * i_r_half; + const int i_theta = 2 * i_theta_half; + smoother_solution[grid.index(i_r, i_theta)] = discrete_solution[grid.index(i_r, i_theta)]; + }); + + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); int iterations = 0; @@ -691,13 +765,12 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeDirBC_Interior_SmallestGr const int max_iterations = 10000; double precision = 1e-12; - while (infinity_norm(HostConstVector(error)) > precision) { + while (infinity_norm(ConstVector(error)) > precision) { extrapolated_smoother_op.extrapolatedSmoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); iterations++; if (iterations >= max_iterations) { @@ -711,10 +784,14 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeDirBC_Interior_SmallestGr ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 200); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); +} +TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeDirBC_Interior_SmallestGrid) +{ + ExtrapolatedSmootherTest_ExtrapolatedSmootherTakeDirBC_Interior_SmallestGrid(); } -TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeAcrossOrigin_SmallestGrid) +void ExtrapolatedSmootherTest_ExtrapolatedSmootherTakeAcrossOrigin_SmallestGrid() { std::vector radii = {1e-5, 0.2, 0.4, 0.45, 0.9, 1.2, 1.3}; std::vector angles = {0, M_PI / 8, M_PI, M_PI + M_PI / 8, M_PI + M_PI}; @@ -726,7 +803,7 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeAcrossOrigin_SmallestGrid using DomainGeometryType = CzarnyGeometry; DomainGeometryType domain_geometry(Rmax, kappa_eps, delta_e); - auto grid = std::make_unique>(radii, angles); + auto tmp_grid = std::make_unique>(radii, angles); double alpha_jump = 0.678 * Rmax; using DensityProfileCoefficientsType = ZoniShiftedCoefficients; @@ -738,32 +815,39 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeAcrossOrigin_SmallestGrid bool cache_domain_geometry = true; auto levelCache = std::make_unique>( - *grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); - Level level(0, std::move(grid), std::move(levelCache), + *tmp_grid, coefficients, domain_geometry, cache_density_rpofile_coefficients, cache_domain_geometry); + Level level(0, std::move(tmp_grid), std::move(levelCache), ExtrapolationType::IMPLICIT_EXTRAPOLATION, 0); + const PolarGrid& grid = level.grid(); DirectSolverTake solver_op(level.grid(), level.levelCache(), DirBC_Interior); ResidualTake residual_op(level.grid(), level.levelCache(), DirBC_Interior); ExtrapolatedSmootherTake extrapolated_smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - HostVector error("error", level.grid().numberOfNodes()); - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); - for (int i_r = 0; i_r < level.grid().nr(); i_r += 2) { - for (int i_theta = 0; i_theta < level.grid().ntheta(); i_theta += 2) { - smoother_solution[level.grid().index(i_r, i_theta)] = discrete_solution[level.grid().index(i_r, i_theta)]; - } - } - - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Vector temp("temp", level.grid().numberOfNodes()); + Vector error("error", level.grid().numberOfNodes()); + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); + Kokkos::parallel_for( + "Fill smoother_solution", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {(grid.nr() + 1) / 2, (grid.ntheta() + 1) / 2} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int i_r_half, const int i_theta_half) { + const int i_r = 2 * i_r_half; + const int i_theta = 2 * i_theta_half; + smoother_solution[grid.index(i_r, i_theta)] = discrete_solution[grid.index(i_r, i_theta)]; + }); + + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); int iterations = 0; @@ -771,13 +855,12 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeAcrossOrigin_SmallestGrid const int max_iterations = 10000; const double precision = 1e-8; - while (infinity_norm(HostConstVector(error)) > 1e-8) { + while (infinity_norm(ConstVector(error)) > 1e-8) { extrapolated_smoother_op.extrapolatedSmoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { - error[i] = discrete_solution[i] - smoother_solution[i]; - }); + Kokkos::parallel_for( + "get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error[i] = discrete_solution[i] - smoother_solution[i]; }); Kokkos::fence(); iterations++; if (iterations >= max_iterations) { @@ -791,5 +874,9 @@ TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeAcrossOrigin_SmallestGrid ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 150); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); +} +TEST(ExtrapolatedSmootherTest, ExtrapolatedSmootherTakeAcrossOrigin_SmallestGrid) +{ + ExtrapolatedSmootherTest_ExtrapolatedSmootherTakeAcrossOrigin_SmallestGrid(); } diff --git a/tests/Interpolation/prolongation.cpp b/tests/Interpolation/prolongation.cpp index a9018f1b..5b10f186 100644 --- a/tests/Interpolation/prolongation.cpp +++ b/tests/Interpolation/prolongation.cpp @@ -61,20 +61,26 @@ TEST(ProlongationTest, ProlongationMatchesStencil) 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); Interpolation I(/*DirBC*/ true); - HostVector coarse_values = generate_random_sample_data(coarse_grid, 1234); - HostVector fine_result("fine_result", fine_grid.numberOfNodes()); + Vector coarse_values = generate_random_sample_data(coarse_grid, 1234); + Vector fine_result("fine_result", fine_grid.numberOfNodes()); I.applyProlongation(coarse_grid, fine_grid, fine_result, coarse_values); + PolarGrid h_fine_grid(fine_grid); + PolarGrid h_coarse_grid(coarse_grid); + + auto h_coarse_values = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, 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_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_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/restriction.cpp b/tests/Interpolation/restriction.cpp index 29821ff1..8bcaf97a 100644 --- a/tests/Interpolation/restriction.cpp +++ b/tests/Interpolation/restriction.cpp @@ -57,21 +57,26 @@ TEST(RestrictionTest, RestrictionMatchesStencil) 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); Interpolation I(/*DirBC*/ true); - HostVector fine_values = generate_random_sample_data(fine_grid, 5678, 0.0, 1.0); - HostVector coarse_result("coarse_result", coarse_grid.numberOfNodes()); + Vector fine_values = generate_random_sample_data(fine_grid, 5678, 0.0, 1.0); + Vector coarse_result("coarse_result", coarse_grid.numberOfNodes()); I.applyRestriction(fine_grid, coarse_grid, coarse_result, fine_values); + PolarGrid h_fine_grid(fine_grid); + PolarGrid h_coarse_grid(coarse_grid); + auto h_fine_values = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, 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_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)]; + expected_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 << ")"; } } diff --git a/tests/Residual/residual.cpp b/tests/Residual/residual.cpp index c74dc5d7..0c902c86 100644 --- a/tests/Residual/residual.cpp +++ b/tests/Residual/residual.cpp @@ -55,23 +55,26 @@ TEST(OperatorATest, applyA_DirBC_Interior) ResidualGive residualGive_operator(level.grid(), level.levelCache(), DirBC_Interior); ResidualTake residualTake_operator(level.grid(), level.levelCache(), DirBC_Interior); - HostVector x = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector x = generate_random_sample_data(level.grid(), 42); + Vector rhs = generate_random_sample_data(level.grid(), 69); - HostVector result_Give("result_Give", level.grid().numberOfNodes()); + Vector result_Give("result_Give", level.grid().numberOfNodes()); residualGive_operator.computeResidual(result_Give, rhs, x); - HostVector result_Take("result_Take", level.grid().numberOfNodes()); + Vector result_Take("result_Take", level.grid().numberOfNodes()); residualTake_operator.computeResidual(result_Take, rhs, x); + auto h_result_Give = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, result_Give); + auto h_result_Take = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, result_Take); + ASSERT_EQ(result_Give.size(), result_Take.size()); for (uint index = 0; index < result_Give.size(); index++) { int i_r, i_theta; level.grid().multiIndex(index, i_r, i_theta); if (i_r == 0 && !DirBC_Interior) - ASSERT_NEAR(result_Give[index], result_Take[index], 1e-8); + ASSERT_NEAR(h_result_Give[index], h_result_Take[index], 1e-8); else - ASSERT_NEAR(result_Give[index], result_Take[index], 1e-11); + ASSERT_NEAR(h_result_Give[index], h_result_Take[index], 1e-11); } } @@ -107,22 +110,25 @@ TEST(OperatorATest, applyA_AcrossOrigin) ResidualGive residualGive_operator(level.grid(), level.levelCache(), DirBC_Interior); ResidualTake residualTake_operator(level.grid(), level.levelCache(), DirBC_Interior); - HostVector x = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector x = generate_random_sample_data(level.grid(), 42); + Vector rhs = generate_random_sample_data(level.grid(), 69); - HostVector result_Give("result_Give", level.grid().numberOfNodes()); + Vector result_Give("result_Give", level.grid().numberOfNodes()); residualGive_operator.computeResidual(result_Give, rhs, x); - HostVector result_Take("result_Take", level.grid().numberOfNodes()); + Vector result_Take("result_Take", level.grid().numberOfNodes()); residualTake_operator.computeResidual(result_Take, rhs, x); + auto h_result_Give = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, result_Give); + auto h_result_Take = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, result_Take); + ASSERT_EQ(result_Give.size(), result_Take.size()); for (uint index = 0; index < result_Give.size(); index++) { int i_r, i_theta; level.grid().multiIndex(index, i_r, i_theta); if (i_r == 0 && !DirBC_Interior) - ASSERT_NEAR(result_Give[index], result_Take[index], 1e-8); + ASSERT_NEAR(h_result_Give[index], h_result_Take[index], 1e-8); else - ASSERT_NEAR(result_Give[index], result_Take[index], 1e-11); + ASSERT_NEAR(h_result_Give[index], h_result_Take[index], 1e-11); } } diff --git a/tests/Smoother/smoother.cpp b/tests/Smoother/smoother.cpp index 28f015ea..e66a6499 100644 --- a/tests/Smoother/smoother.cpp +++ b/tests/Smoother/smoother.cpp @@ -27,7 +27,7 @@ using namespace gmgpolar; /* Test 1/2: */ /* Does the Take and Give Implementation match up? */ -TEST(SmootherTest, smoother_DirBC_Interior) +void SmootherTest_smoother_DirBC_Interior() { std::vector radii = {1e-5, 0.2, 0.25, 0.5, 0.8, 0.9, 0.95, 1.2, 1.3}; std::vector angles = { @@ -61,30 +61,35 @@ TEST(SmootherTest, smoother_DirBC_Interior) SmootherGive smootherGive_operator(level.grid(), level.levelCache(), DirBC_Interior); SmootherTake smootherTake_operator(level.grid(), level.levelCache(), DirBC_Interior); - HostVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 69); - HostVector start = generate_random_sample_data(PolarGrid(level.grid()), 24); - HostVector temp = generate_random_sample_data(PolarGrid(level.grid()), 8); + Vector rhs = generate_random_sample_data(level.grid(), 69); + Vector start = generate_random_sample_data(level.grid(), 24); + Vector temp = generate_random_sample_data(level.grid(), 8); - HostVector solution_Give("solution_Give", start.size()); + Vector solution_Give("solution_Give", start.size()); Kokkos::deep_copy(solution_Give, start); smootherGive_operator.smoothing(solution_Give, rhs, temp); - HostVector solution_Take("solution_Take", start.size()); + Vector solution_Take("solution_Take", start.size()); Kokkos::deep_copy(solution_Take, start); smootherTake_operator.smoothing(solution_Take, rhs, temp); + auto h_solution_Give = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), solution_Give); + auto h_solution_Take = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), solution_Take); + ASSERT_EQ(solution_Give.size(), solution_Take.size()); - for (uint index = 0; index < solution_Give.size(); index++) { + for (uint index = 0; index < h_solution_Give.size(); index++) { int i_r, i_theta; level.grid().multiIndex(index, i_r, i_theta); if (i_r == 0 && !DirBC_Interior) - ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11); + ASSERT_NEAR(h_solution_Give[index], h_solution_Take[index], 1e-11); else - ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11); + ASSERT_NEAR(h_solution_Give[index], h_solution_Take[index], 1e-11); } } +TEST(SmootherTest, smoother_DirBC_Interior) +{ SmootherTest_smoother_DirBC_Interior();} -TEST(SmootherTest, smoother_AcrossOrigin) +void SmootherTest_smoother_AcrossOrigin() { std::vector radii = {1e-5, 0.2, 0.25, 0.5, 0.8, 0.9, 0.95, 1.2, 1.3}; std::vector angles = { @@ -118,33 +123,38 @@ TEST(SmootherTest, smoother_AcrossOrigin) SmootherGive smootherGive_operator(level.grid(), level.levelCache(), DirBC_Interior); SmootherTake smootherTake_operator(level.grid(), level.levelCache(), DirBC_Interior); - HostVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 69); - HostVector start = generate_random_sample_data(PolarGrid(level.grid()), 24); - HostVector temp = generate_random_sample_data(PolarGrid(level.grid()), 8); + Vector rhs = generate_random_sample_data(level.grid(), 69); + Vector start = generate_random_sample_data(level.grid(), 24); + Vector temp = generate_random_sample_data(level.grid(), 8); - HostVector solution_Give("solution_Give", start.size()); + Vector solution_Give("solution_Give", start.size()); Kokkos::deep_copy(solution_Give, start); smootherGive_operator.smoothing(solution_Give, rhs, temp); - HostVector solution_Take("solution_Take", start.size()); + Vector solution_Take("solution_Take", start.size()); Kokkos::deep_copy(solution_Take, start); smootherTake_operator.smoothing(solution_Take, rhs, temp); + auto h_solution_Give = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), solution_Give); + auto h_solution_Take = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), solution_Take); + ASSERT_EQ(solution_Give.size(), solution_Take.size()); for (uint index = 0; index < solution_Give.size(); index++) { int i_r, i_theta; level.grid().multiIndex(index, i_r, i_theta); if (i_r == 0 && !DirBC_Interior) - ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-8); + ASSERT_NEAR(h_solution_Give[index], h_solution_Take[index], 1e-8); else - ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-10); + ASSERT_NEAR(h_solution_Give[index], h_solution_Take[index], 1e-10); } } +TEST(SmootherTest, smoother_AcrossOrigin) +{ SmootherTest_smoother_AcrossOrigin(); } /* Test 2/2: */ /* Does the smoother converge to the directSolver solution? */ -TEST(SmootherTest, SmootherDirBC_Interior) +void SmootherTest_SmootherDirBC_Interior() { std::vector radii = {1e-5, 0.2, 0.25, 0.5, 0.8, 0.9, 0.95, 1.2, 1.3}; std::vector angles = { @@ -176,17 +186,17 @@ TEST(SmootherTest, SmootherDirBC_Interior) ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); SmootherGive smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - HostVector error("error", level.grid().numberOfNodes()); - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector temp("temp", level.grid().numberOfNodes()); + Vector error("error", level.grid().numberOfNodes()); + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -196,11 +206,11 @@ TEST(SmootherTest, SmootherDirBC_Interior) const int max_iterations = 10000; const double precision = 1e-12; - while (infinity_norm(HostConstVector(error)) > precision) { + while (infinity_norm(ConstVector(error)) > precision) { smoother_op.smoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -216,10 +226,12 @@ TEST(SmootherTest, SmootherDirBC_Interior) ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 300); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); } +TEST(SmootherTest, SmootherDirBC_Interior) +{SmootherTest_SmootherDirBC_Interior(); } -TEST(SmootherTest, SmootherAcrossOrigin) +void SmootherTest_SmootherAcrossOrigin() { std::vector radii = {1e-5, 0.2, 0.25, 0.5, 0.8, 0.9, 0.95, 1.2, 1.3}; std::vector angles = { @@ -251,17 +263,17 @@ TEST(SmootherTest, SmootherAcrossOrigin) ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); SmootherGive smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - HostVector error("error", level.grid().numberOfNodes()); - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector temp("temp", level.grid().numberOfNodes()); + Vector error("error", level.grid().numberOfNodes()); + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -271,11 +283,11 @@ TEST(SmootherTest, SmootherAcrossOrigin) const int max_iterations = 10000; const double precision = 1e-8; - while (infinity_norm(HostConstVector(error)) > precision) { + while (infinity_norm(ConstVector(error)) > precision) { smoother_op.smoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -291,10 +303,14 @@ TEST(SmootherTest, SmootherAcrossOrigin) ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 600); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); +} +TEST(SmootherTest, SmootherAcrossOrigin) +{ + SmootherTest_SmootherAcrossOrigin(); } -TEST(SmootherTest, SmootherDirBC_Interior_SmallestGrid) +void SmootherTest_SmootherDirBC_Interior_SmallestGrid() { std::vector radii = {1e-5, 0.2, 0.9, 1.2, 1.3}; std::vector angles = {0, M_PI / 8, M_PI, M_PI + M_PI / 8, M_PI + M_PI}; @@ -325,17 +341,17 @@ TEST(SmootherTest, SmootherDirBC_Interior_SmallestGrid) ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); SmootherGive smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - HostVector error("error", level.grid().numberOfNodes()); - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector temp("temp", level.grid().numberOfNodes()); + Vector error("error", level.grid().numberOfNodes()); + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -345,11 +361,11 @@ TEST(SmootherTest, SmootherDirBC_Interior_SmallestGrid) const int max_iterations = 10000; double precision = 1e-12; - while (infinity_norm(HostConstVector(error)) > precision) { + while (infinity_norm(ConstVector(error)) > precision) { smoother_op.smoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -365,10 +381,12 @@ TEST(SmootherTest, SmootherDirBC_Interior_SmallestGrid) ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 30); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); } +TEST(SmootherTest, SmootherDirBC_Interior_SmallestGrid) +{SmootherTest_SmootherDirBC_Interior_SmallestGrid(); } -TEST(SmootherTest, SmootherAcrossOrigin_SmallestGrid) +void SmootherTest_SmootherAcrossOrigin_SmallestGrid() { std::vector radii = {1e-5, 0.2, 0.9, 1.2, 1.3}; std::vector angles = {0, M_PI / 8, M_PI, M_PI + M_PI / 8, M_PI + M_PI}; @@ -399,17 +417,17 @@ TEST(SmootherTest, SmootherAcrossOrigin_SmallestGrid) ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); SmootherGive smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - HostVector error("error", level.grid().numberOfNodes()); - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector temp("temp", level.grid().numberOfNodes()); + Vector error("error", level.grid().numberOfNodes()); + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -419,11 +437,11 @@ TEST(SmootherTest, SmootherAcrossOrigin_SmallestGrid) const int max_iterations = 10000; const double precision = 1e-8; - while (infinity_norm(HostConstVector(error)) > 1e-8) { + while (infinity_norm(ConstVector(error)) > 1e-8) { smoother_op.smoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -439,12 +457,14 @@ TEST(SmootherTest, SmootherAcrossOrigin_SmallestGrid) ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 80); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); } +TEST(SmootherTest, SmootherAcrossOrigin_SmallestGrid) +{ SmootherTest_SmootherAcrossOrigin_SmallestGrid(); } /* Using "Take" */ -TEST(SmootherTest, SmootherTakeDirBC_Interior) +void SmootherTest_SmootherTakeDirBC_Interior() { std::vector radii = {1e-5, 0.2, 0.25, 0.5, 0.8, 0.9, 0.95, 1.2, 1.3}; std::vector angles = { @@ -476,17 +496,17 @@ TEST(SmootherTest, SmootherTakeDirBC_Interior) ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); SmootherTake smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - HostVector error("error", level.grid().numberOfNodes()); - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector temp("temp", level.grid().numberOfNodes()); + Vector error("error", level.grid().numberOfNodes()); + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -496,11 +516,11 @@ TEST(SmootherTest, SmootherTakeDirBC_Interior) const int max_iterations = 10000; const double precision = 1e-12; - while (infinity_norm(HostConstVector(error)) > precision) { + while (infinity_norm(ConstVector(error)) > precision) { smoother_op.smoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -516,10 +536,12 @@ TEST(SmootherTest, SmootherTakeDirBC_Interior) ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 300); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); } +TEST(SmootherTest, SmootherTakeDirBC_Interior) +{ SmootherTest_SmootherTakeDirBC_Interior(); } -TEST(SmootherTest, SmootherTakeAcrossOrigin) +void SmootherTest_SmootherTakeAcrossOrigin() { std::vector radii = {1e-5, 0.2, 0.25, 0.5, 0.8, 0.9, 0.95, 1.2, 1.3}; std::vector angles = { @@ -551,17 +573,17 @@ TEST(SmootherTest, SmootherTakeAcrossOrigin) ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); SmootherTake smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - HostVector error("error", level.grid().numberOfNodes()); - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector temp("temp", level.grid().numberOfNodes()); + Vector error("error", level.grid().numberOfNodes()); + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -571,11 +593,11 @@ TEST(SmootherTest, SmootherTakeAcrossOrigin) const int max_iterations = 10000; const double precision = 1e-8; - while (infinity_norm(HostConstVector(error)) > precision) { + while (infinity_norm(ConstVector(error)) > precision) { smoother_op.smoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -591,10 +613,12 @@ TEST(SmootherTest, SmootherTakeAcrossOrigin) ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 600); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); } +TEST(SmootherTest, SmootherTakeAcrossOrigin) +{SmootherTest_SmootherTakeAcrossOrigin(); } -TEST(SmootherTest, SmootherTakeDirBC_Interior_SmallestGrid) +void SmootherTest_SmootherTakeDirBC_Interior_SmallestGrid() { std::vector radii = {1e-5, 0.2, 0.9, 1.2, 1.3}; std::vector angles = {0, M_PI / 8, M_PI, M_PI + M_PI / 8, M_PI + M_PI}; @@ -625,17 +649,17 @@ TEST(SmootherTest, SmootherTakeDirBC_Interior_SmallestGrid) ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); SmootherTake smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - HostVector error("error", level.grid().numberOfNodes()); - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector temp("temp", level.grid().numberOfNodes()); + Vector error("error", level.grid().numberOfNodes()); + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -645,11 +669,11 @@ TEST(SmootherTest, SmootherTakeDirBC_Interior_SmallestGrid) const int max_iterations = 10000; double precision = 1e-12; - while (infinity_norm(HostConstVector(error)) > precision) { + while (infinity_norm(ConstVector(error)) > precision) { smoother_op.smoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -665,10 +689,12 @@ TEST(SmootherTest, SmootherTakeDirBC_Interior_SmallestGrid) ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 30); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); } +TEST(SmootherTest, SmootherTakeDirBC_Interior_SmallestGrid) +{ SmootherTest_SmootherTakeDirBC_Interior_SmallestGrid(); } -TEST(SmootherTest, SmootherTakeAcrossOrigin_SmallestGrid) +void SmootherTest_SmootherTakeAcrossOrigin_SmallestGrid() { std::vector radii = {1e-5, 0.2, 0.9, 1.2, 1.3}; std::vector angles = {0, M_PI / 8, M_PI, M_PI + M_PI / 8, M_PI + M_PI}; @@ -699,17 +725,17 @@ TEST(SmootherTest, SmootherTakeAcrossOrigin_SmallestGrid) ResidualGive residual_op(level.grid(), level.levelCache(), DirBC_Interior); SmootherTake smoother_op(level.grid(), level.levelCache(), DirBC_Interior); - HostConstVector rhs = generate_random_sample_data(PolarGrid(level.grid()), 42); - HostVector discrete_solution("discrete_solution", rhs.size()); + ConstVector rhs = generate_random_sample_data(level.grid(), 42); + Vector discrete_solution("discrete_solution", rhs.size()); Kokkos::deep_copy(discrete_solution, rhs); solver_op.solveInPlace(discrete_solution); - HostVector temp("temp", level.grid().numberOfNodes()); - HostVector error("error", level.grid().numberOfNodes()); - HostVector smoother_solution = generate_random_sample_data(PolarGrid(level.grid()), 69); + Vector temp("temp", level.grid().numberOfNodes()); + Vector error("error", level.grid().numberOfNodes()); + Vector smoother_solution = generate_random_sample_data(level.grid(), 69); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -719,11 +745,11 @@ TEST(SmootherTest, SmootherTakeAcrossOrigin_SmallestGrid) const int max_iterations = 10000; const double precision = 1e-8; - while (infinity_norm(HostConstVector(error)) > 1e-8) { + while (infinity_norm(ConstVector(error)) > 1e-8) { smoother_op.smoothing(smoother_solution, rhs, temp); - Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), - [&](uint i) { + Kokkos::parallel_for("get error", Kokkos::RangePolicy(0, error.size()), + KOKKOS_LAMBDA(uint i) { error(i) = discrete_solution(i) - smoother_solution(i); }); Kokkos::fence(); @@ -739,5 +765,7 @@ TEST(SmootherTest, SmootherTakeAcrossOrigin_SmallestGrid) ASSERT_TRUE(!max_iterations_reached); ASSERT_LT(iterations, 80); - ASSERT_NEAR(infinity_norm(HostConstVector(error)), 0.0, precision); + ASSERT_NEAR(infinity_norm(ConstVector(error)), 0.0, precision); } +TEST(SmootherTest, SmootherTakeAcrossOrigin_SmallestGrid) +{ SmootherTest_SmootherTakeAcrossOrigin_SmallestGrid(); } diff --git a/tests/test_tools.h b/tests/test_tools.h index b4371332..c68792b6 100644 --- a/tests/test_tools.h +++ b/tests/test_tools.h @@ -6,7 +6,8 @@ using namespace gmgpolar; -inline HostVector generate_random_sample_data(const PolarGrid& grid, unsigned int seed, +template +inline Vector generate_random_sample_data(const PolarGrid& grid, unsigned int seed, double min_val = -100.0, double max_val = 100.0) { HostVector x("x", grid.numberOfNodes()); @@ -15,5 +16,6 @@ inline HostVector generate_random_sample_data(const PolarGrid