From 51fb816d72a05497bf8315f29bed3c1e2895366c Mon Sep 17 00:00:00 2001 From: julianlitz Date: Sat, 16 May 2026 19:26:39 +0200 Subject: [PATCH 1/9] Use Kokkos::Parallel in Solver Utils --- include/GMGPolar/solver.h | 136 ++++++++++++++++++++------------------ 1 file changed, 73 insertions(+), 63 deletions(-) diff --git a/include/GMGPolar/solver.h b/include/GMGPolar/solver.h index 5ff24ea52..cc1219935 100644 --- a/include/GMGPolar/solver.h +++ b/include/GMGPolar/solver.h @@ -511,46 +511,48 @@ void GMGPolar::applyExtrapolation(in assert(std::ssize(fine_values) == fineGrid.numberOfNodes()); assert(std::ssize(coarse_values) == coarseGrid.numberOfNodes()); -#pragma omp parallel num_threads(max_omp_threads_) - { -/* Circluar Indexing Section */ -/* For loop matches circular access pattern */ -#pragma omp for nowait - for (int i_r = 0; i_r < fineGrid.numberSmootherCircles(); i_r++) { - int i_r_coarse = i_r / 2; - for (int i_theta = 0; i_theta < fineGrid.ntheta(); i_theta++) { - int i_theta_coarse = i_theta / 2; - - if (i_r & 1 || i_theta & 1) { - fine_values[fineGrid.index(i_r, i_theta)] *= 4.0 / 3.0; - } - else { - int fine_idx = fineGrid.index(i_r, i_theta); - int coarse_idx = coarseGrid.index(i_r_coarse, i_theta_coarse); - fine_values[fine_idx] = (4.0 * fine_values[fine_idx] - coarse_values[coarse_idx]) / 3.0; - } + /* We split the loops into two regions to better respect the */ + /* access patterns of the smoother and improve cache locality. */ + + /* Circular Indexing Section */ + /* For loop matches circular access pattern */ + Kokkos::parallel_for( + "Extrapolation: Apply Extrapolation (Circular)", + Kokkos::MDRangePolicy>( + {0, 0}, {fineGrid.numberSmootherCircles(), fineGrid.ntheta()}), + KOKKOS_CLASS_LAMBDA(const int i_r, const int i_theta) { + const int i_r_coarse = i_r / 2; + const int i_theta_coarse = i_theta / 2; + if (i_r & 1 || i_theta & 1) { + fine_values[fineGrid.index(i_r, i_theta)] *= 4.0 / 3.0; } - } - -/* Radial Indexing Section */ -/* For loop matches radial access pattern */ -#pragma omp for nowait - for (int i_theta = 0; i_theta < fineGrid.ntheta(); i_theta++) { - int i_theta_coarse = i_theta / 2; - for (int i_r = fineGrid.numberSmootherCircles(); i_r < fineGrid.nr(); i_r++) { - int i_r_coarse = i_r / 2; - - if (i_r & 1 || i_theta & 1) { - fine_values[fineGrid.index(i_r, i_theta)] *= 4.0 / 3.0; - } - else { - int fine_idx = fineGrid.index(i_r, i_theta); - int coarse_idx = coarseGrid.index(i_r_coarse, i_theta_coarse); - fine_values[fine_idx] = (4.0 * fine_values[fine_idx] - coarse_values[coarse_idx]) / 3.0; - } + else { + const int fine_idx = fineGrid.index(i_r, i_theta); + const int coarse_idx = coarseGrid.index(i_r_coarse, i_theta_coarse); + fine_values[fine_idx] = (4.0 * fine_values[fine_idx] - coarse_values[coarse_idx]) / 3.0; } - } - } + }); + + /* Radial Indexing Section */ + /* For loop matches radial access pattern */ + Kokkos::parallel_for( + "Extrapolation: Apply Extrapolation (Radial)", + Kokkos::MDRangePolicy>({0, fineGrid.numberSmootherCircles()}, + {fineGrid.ntheta(), fineGrid.nr()}), + KOKKOS_CLASS_LAMBDA(const int i_theta, const int i_r) { + const int i_r_coarse = i_r / 2; + const int i_theta_coarse = i_theta / 2; + if (i_r & 1 || i_theta & 1) { + fine_values[fineGrid.index(i_r, i_theta)] *= 4.0 / 3.0; + } + else { + const int fine_idx = fineGrid.index(i_r, i_theta); + const int coarse_idx = coarseGrid.index(i_r_coarse, i_theta_coarse); + fine_values[fine_idx] = (4.0 * fine_values[fine_idx] - coarse_values[coarse_idx]) / 3.0; + } + }); + + Kokkos::fence(); } // ============================================================================= @@ -592,31 +594,39 @@ std::pair GMGPolar:: assert(solution.size() == error.size()); assert(std::ssize(solution) == grid.numberOfNodes()); -#pragma omp parallel num_threads(max_omp_threads_) - { -#pragma omp for nowait - for (int i_r = 0; i_r < grid.numberSmootherCircles(); i_r++) { - double r = grid.radius(i_r); - for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { - double theta = grid.theta(i_theta); - error[grid.index(i_r, i_theta)] = - exact_solution.exact_solution(r, theta) - solution[grid.index(i_r, i_theta)]; - } - } -#pragma omp for nowait - for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { - double theta = grid.theta(i_theta); - for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { - double r = grid.radius(i_r); - error[grid.index(i_r, i_theta)] = - exact_solution.exact_solution(r, theta) - solution[grid.index(i_r, i_theta)]; - } - } - } - ConstVector c_error = error; - double weighted_euclidean_error = l2_norm(c_error) / std::sqrt(grid.numberOfNodes()); - double infinity_error = infinity_norm(c_error); - + /* We split the loops into two regions to better respect the */ + /* access patterns of the smoother and improve cache locality. */ + + /* Circular Indexing Section */ + /* For loop matches circular access pattern */ + Kokkos::parallel_for( + "Exact Error: Compute Error (Circular)", + Kokkos::MDRangePolicy>( + {0, 0}, {grid.numberSmootherCircles(), grid.ntheta()}), + KOKKOS_CLASS_LAMBDA(const int i_r, const int i_theta) { + const double radius = grid.radius(i_r); + const double theta = grid.theta(i_theta); + error[grid.index(i_r, i_theta)] = + exact_solution.exact_solution(radius, theta) - solution[grid.index(i_r, i_theta)]; + }); + + /* Radial Indexing Section */ + /* For loop matches radial access pattern */ + Kokkos::parallel_for( + "Exact Error: Compute Error (Radial)", + Kokkos::MDRangePolicy>({0, grid.numberSmootherCircles()}, + {grid.ntheta(), grid.nr()}), + KOKKOS_CLASS_LAMBDA(const int i_theta, const int i_r) { + const double radius = grid.radius(i_r); + const double theta = grid.theta(i_theta); + error[grid.index(i_r, i_theta)] = + exact_solution.exact_solution(radius, theta) - solution[grid.index(i_r, i_theta)]; + }); + + Kokkos::fence(); + + const double weighted_euclidean_error = l2_norm(ConstVector(error)) / std::sqrt(grid.numberOfNodes()); + const double infinity_error = infinity_norm(ConstVector(error)); return std::make_pair(weighted_euclidean_error, infinity_error); } From 3b62a4b2c61e54671a822a9074c1c15c2210008b Mon Sep 17 00:00:00 2001 From: julianlitz Date: Sat, 16 May 2026 19:42:02 +0200 Subject: [PATCH 2/9] Yeet class lambda --- include/GMGPolar/solver.h | 8 ++++---- include/Residual/ResidualTake/applyATake.inl | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/include/GMGPolar/solver.h b/include/GMGPolar/solver.h index cc1219935..a7acc9ae7 100644 --- a/include/GMGPolar/solver.h +++ b/include/GMGPolar/solver.h @@ -520,7 +520,7 @@ void GMGPolar::applyExtrapolation(in "Extrapolation: Apply Extrapolation (Circular)", Kokkos::MDRangePolicy>( {0, 0}, {fineGrid.numberSmootherCircles(), fineGrid.ntheta()}), - KOKKOS_CLASS_LAMBDA(const int i_r, const int i_theta) { + KOKKOS_LAMBDA(const int i_r, const int i_theta) { const int i_r_coarse = i_r / 2; const int i_theta_coarse = i_theta / 2; if (i_r & 1 || i_theta & 1) { @@ -539,7 +539,7 @@ void GMGPolar::applyExtrapolation(in "Extrapolation: Apply Extrapolation (Radial)", Kokkos::MDRangePolicy>({0, fineGrid.numberSmootherCircles()}, {fineGrid.ntheta(), fineGrid.nr()}), - KOKKOS_CLASS_LAMBDA(const int i_theta, const int i_r) { + KOKKOS_LAMBDA(const int i_theta, const int i_r) { const int i_r_coarse = i_r / 2; const int i_theta_coarse = i_theta / 2; if (i_r & 1 || i_theta & 1) { @@ -603,7 +603,7 @@ std::pair GMGPolar:: "Exact Error: Compute Error (Circular)", Kokkos::MDRangePolicy>( {0, 0}, {grid.numberSmootherCircles(), grid.ntheta()}), - KOKKOS_CLASS_LAMBDA(const int i_r, const int i_theta) { + KOKKOS_LAMBDA(const int i_r, const int i_theta) { const double radius = grid.radius(i_r); const double theta = grid.theta(i_theta); error[grid.index(i_r, i_theta)] = @@ -616,7 +616,7 @@ std::pair GMGPolar:: "Exact Error: Compute Error (Radial)", Kokkos::MDRangePolicy>({0, grid.numberSmootherCircles()}, {grid.ntheta(), grid.nr()}), - KOKKOS_CLASS_LAMBDA(const int i_theta, const int i_r) { + KOKKOS_LAMBDA(const int i_theta, const int i_r) { const double radius = grid.radius(i_r); const double theta = grid.theta(i_theta); error[grid.index(i_r, i_theta)] = diff --git a/include/Residual/ResidualTake/applyATake.inl b/include/Residual/ResidualTake/applyATake.inl index 062d9c317..3cd44fc32 100644 --- a/include/Residual/ResidualTake/applyATake.inl +++ b/include/Residual/ResidualTake/applyATake.inl @@ -94,7 +94,7 @@ void ResidualTake::applySystemOperator(Vector result, Co {grid.numberSmootherCircles(), grid.ntheta()} // Ending point of the index space ), // Kokkos lambda function to execute for each point in the index space - KOKKOS_CLASS_LAMBDA(const int i_r, const int i_theta) { + KOKKOS_LAMBDA(const int i_r, const int i_theta) { node_apply_a_take(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF, coeff_beta); }); @@ -106,7 +106,7 @@ void ResidualTake::applySystemOperator(Vector result, Co {grid.ntheta(), grid.nr()} // Ending point of the index space ), // Kokkos lambda function to execute for each point in the index space - KOKKOS_CLASS_LAMBDA(const int i_theta, const int i_r) { + KOKKOS_LAMBDA(const int i_theta, const int i_r) { node_apply_a_take(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF, coeff_beta); }); From cc30b20e6273545a004e6a6fd90332bb6b635f53 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Wed, 20 May 2026 23:13:26 +0200 Subject: [PATCH 3/9] Refactor exact_solution checks and improve comments --- include/GMGPolar/solver.h | 67 +++++++++++++++++++++------------------ 1 file changed, 37 insertions(+), 30 deletions(-) diff --git a/include/GMGPolar/solver.h b/include/GMGPolar/solver.h index 3916c83de..0d69605dc 100644 --- a/include/GMGPolar/solver.h +++ b/include/GMGPolar/solver.h @@ -96,8 +96,9 @@ void GMGPolar::solve(const BoundaryC LIKWID_STOP("Solver"); auto start_check_exact_error = std::chrono::high_resolution_clock::now(); - if (exact_solution_ != nullptr) + if (exact_solution_) { evaluateExactError(level, *exact_solution_); + } auto end_check_exact_error = std::chrono::high_resolution_clock::now(); t_check_exact_error_ += std::chrono::duration(end_check_exact_error - start_check_exact_error).count(); @@ -169,7 +170,7 @@ void GMGPolar::solve(const BoundaryC if (paraview_) { writeToVTK("output_solution", level, level.solution()); - if (exact_solution_ != nullptr) { + if (exact_solution_) { computeExactError(level, level.solution(), level.residual(), *exact_solution_); writeToVTK("output_error", level, level.residual()); } @@ -244,8 +245,9 @@ void GMGPolar::solveMultigrid(double LIKWID_STOP("Solver"); auto start_check_exact_error = std::chrono::high_resolution_clock::now(); - if (exact_solution_ != nullptr) + if (exact_solution_) { evaluateExactError(level, *exact_solution_); + } auto end_check_exact_error = std::chrono::high_resolution_clock::now(); t_check_exact_error_ += std::chrono::duration(end_check_exact_error - start_check_exact_error).count(); @@ -593,37 +595,42 @@ std::pair GMGPolar:: assert(solution.size() == error.size()); assert(std::ssize(solution) == grid.numberOfNodes()); - /* We split the loops into two regions to better respect the */ - /* access patterns of the smoother and improve cache locality. */ - - /* Circular Indexing Section */ - /* For loop matches circular access pattern */ + /* Fill exact solution on the host. */ + // Note: We must compute the exact solution values on the host because the + // ExactSolution object is not designed to be used on the device. + // As it is planned that the solution and error vector will use device memory during GPU porting, + // we need to transfer the exact solution values from host to device using the Kokkos::deep_copy operation. Kokkos::parallel_for( - "Exact Error: Compute Error (Circular)", - Kokkos::MDRangePolicy>( - {0, 0}, {grid.numberSmootherCircles(), grid.ntheta()}), + "Exact Error: Compute Exact Solution Values (Circular)", + Kokkos::MDRangePolicy>( // Host parallel loop + {0, 0}, // Starting index + {grid.numberSmootherCircles(), grid.ntheta()}), // Ending index KOKKOS_LAMBDA(const int i_r, const int i_theta) { - const double radius = grid.radius(i_r); - const double theta = grid.theta(i_theta); - error[grid.index(i_r, i_theta)] = - exact_solution.exact_solution(radius, theta) - solution[grid.index(i_r, i_theta)]; + const double radius = grid.radius(i_r); + const double theta = grid.theta(i_theta); + const int index = grid.index(i_r, i_theta); + exact_solution_values_host[index] = exact_solution.exact_solution(radius, theta); }); - - /* Radial Indexing Section */ - /* For loop matches radial access pattern */ Kokkos::parallel_for( - "Exact Error: Compute Error (Radial)", - Kokkos::MDRangePolicy>({0, grid.numberSmootherCircles()}, - {grid.ntheta(), grid.nr()}), + "Exact Error: Compute Exact Solution Values (Radial)", + Kokkos::MDRangePolicy>( // Host parallel loop + {0, grid.numberSmootherCircles()}, // Starting index + {grid.ntheta(), grid.nr()}), // Ending index KOKKOS_LAMBDA(const int i_theta, const int i_r) { - const double radius = grid.radius(i_r); - const double theta = grid.theta(i_theta); - error[grid.index(i_r, i_theta)] = - exact_solution.exact_solution(radius, theta) - solution[grid.index(i_r, i_theta)]; + const double radius = grid.radius(i_r); + const double theta = grid.theta(i_theta); + const int index = grid.index(i_r, i_theta); + exact_solution_values_host[index] = exact_solution.exact_solution(radius, theta); }); - Kokkos::fence(); + // Transfer the exact solution values from host to device memory for error computation. + Kokkos::deep_copy(error, exact_solution_values_host); + + // Compute the error as the difference between the exact and numerical solutions. + subtract(error, HostConstVector(solution)); + + // Compute the weighted L2 norm and infinity norm of the error. const double weighted_euclidean_error = l2_norm(HostConstVector(error)) / std::sqrt(grid.numberOfNodes()); const double infinity_error = infinity_norm(HostConstVector(error)); return std::make_pair(weighted_euclidean_error, infinity_error); @@ -655,7 +662,7 @@ bool GMGPolar::converged(double resi // ============================================================================= template -void GMGPolar::printIterationHeader(const ExactSolution* exact_solution) +void GMGPolar::printIterationHeader(const bool is_exact_solution_provided) { if (verbose_ <= 0) return; @@ -673,7 +680,7 @@ void GMGPolar::printIterationHeader( else std::cout << std::setw(15 + table_spacing) << "||r_k||/||rhs||"; } - if (exact_solution != nullptr) { + if (is_exact_solution_provided) { std::cout << std::setw(12 + table_spacing) << "||u-u_k||_l2"; std::cout << std::setw(13 + table_spacing) << "||u-u_k||_inf"; } @@ -685,7 +692,7 @@ template ::printIterationInfo(int iteration, double current_residual_norm, double current_relative_residual_norm, - const ExactSolution* exact_solution) + const bool is_exact_solution_provided) { if (verbose_ <= 0) return; @@ -697,7 +704,7 @@ void GMGPolar::printIterationInfo(in std::cout << std::setw(9 + table_spacing + 2) << current_residual_norm; std::cout << std::setw(15 + table_spacing) << current_relative_residual_norm; } - if (exact_solution != nullptr) { + if (is_exact_solution_provided) { std::cout << std::setw(12 + table_spacing) << exact_errors_.back().first; std::cout << std::setw(13 + table_spacing) << exact_errors_.back().second; } From cba304a101ec93ff87f35a361377461dcb0057d8 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Wed, 20 May 2026 23:13:54 +0200 Subject: [PATCH 4/9] Change printIterationHeader parameter type --- include/GMGPolar/gmgpolar.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/GMGPolar/gmgpolar.h b/include/GMGPolar/gmgpolar.h index 1f46bfe39..3061679ee 100644 --- a/include/GMGPolar/gmgpolar.h +++ b/include/GMGPolar/gmgpolar.h @@ -174,9 +174,9 @@ class GMGPolar : public IGMGPolar /* ----------------- */ /* Print information */ void printSettings(const PolarGrid& finest_grid, const PolarGrid& coarsest_grid) const; - void printIterationHeader(const ExactSolution* exact_solution); + void printIterationHeader(const bool is_exact_solution_provided); void printIterationInfo(int iteration, double current_residual_norm, double current_relative_residual_norm, - const ExactSolution* exact_solution); + const bool is_exact_solution_provided); /* ------------------- */ /* Multigrid Functions */ From eedcebe7cbddcc16bdc32824d746f4dc0ad1bbf1 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Wed, 20 May 2026 23:18:30 +0200 Subject: [PATCH 5/9] Introduce host vector for exact solution values Added a host vector for exact solution values in solver.h. --- include/GMGPolar/solver.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/GMGPolar/solver.h b/include/GMGPolar/solver.h index 0d69605dc..56dbe6b8b 100644 --- a/include/GMGPolar/solver.h +++ b/include/GMGPolar/solver.h @@ -600,6 +600,8 @@ std::pair GMGPolar:: // ExactSolution object is not designed to be used on the device. // As it is planned that the solution and error vector will use device memory during GPU porting, // we need to transfer the exact solution values from host to device using the Kokkos::deep_copy operation. + HostVector exact_solution_values_host("Exact Solution Values Host", grid.numberOfNodes()); + Kokkos::parallel_for( "Exact Error: Compute Exact Solution Values (Circular)", Kokkos::MDRangePolicy>( // Host parallel loop From 3337bc5ed0478d18b778ba2697cfff4815e56838 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Wed, 20 May 2026 23:20:58 +0200 Subject: [PATCH 6/9] Fix formatting issues in solver.h comments --- include/GMGPolar/solver.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/GMGPolar/solver.h b/include/GMGPolar/solver.h index 56dbe6b8b..131d2e3e9 100644 --- a/include/GMGPolar/solver.h +++ b/include/GMGPolar/solver.h @@ -601,7 +601,7 @@ std::pair GMGPolar:: // As it is planned that the solution and error vector will use device memory during GPU porting, // we need to transfer the exact solution values from host to device using the Kokkos::deep_copy operation. HostVector exact_solution_values_host("Exact Solution Values Host", grid.numberOfNodes()); - + Kokkos::parallel_for( "Exact Error: Compute Exact Solution Values (Circular)", Kokkos::MDRangePolicy>( // Host parallel loop From dedc1450ee4dba1cb952c45cbc41d72d4cfe3935 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Wed, 20 May 2026 23:28:51 +0200 Subject: [PATCH 7/9] Refactor lambda expressions to use capture by reference --- include/GMGPolar/solver.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/GMGPolar/solver.h b/include/GMGPolar/solver.h index 131d2e3e9..35c35628b 100644 --- a/include/GMGPolar/solver.h +++ b/include/GMGPolar/solver.h @@ -607,7 +607,7 @@ std::pair GMGPolar:: Kokkos::MDRangePolicy>( // Host parallel loop {0, 0}, // Starting index {grid.numberSmootherCircles(), grid.ntheta()}), // Ending index - KOKKOS_LAMBDA(const int i_r, const int i_theta) { + [&](const int i_r, const int i_theta) { const double radius = grid.radius(i_r); const double theta = grid.theta(i_theta); const int index = grid.index(i_r, i_theta); @@ -618,7 +618,7 @@ std::pair GMGPolar:: Kokkos::MDRangePolicy>( // Host parallel loop {0, grid.numberSmootherCircles()}, // Starting index {grid.ntheta(), grid.nr()}), // Ending index - KOKKOS_LAMBDA(const int i_theta, const int i_r) { + [&](const int i_theta, const int i_r) { const double radius = grid.radius(i_r); const double theta = grid.theta(i_theta); const int index = grid.index(i_r, i_theta); From b35ce04184b8a7b68683d9374641abd7437d1b31 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Wed, 20 May 2026 23:42:20 +0200 Subject: [PATCH 8/9] Update gmgpolar.h --- include/GMGPolar/gmgpolar.h | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/include/GMGPolar/gmgpolar.h b/include/GMGPolar/gmgpolar.h index 3061679ee..a5a18e847 100644 --- a/include/GMGPolar/gmgpolar.h +++ b/include/GMGPolar/gmgpolar.h @@ -128,6 +128,7 @@ class GMGPolar : public IGMGPolar double mean_residual_reduction_factor_; bool converged(double current_residual_norm, double first_residual_norm); +public: // Public due to cuda restrictions /* ---------------------------------------------------- */ /* Compute exact error if an exact solution is provided */ // The results are stored as a pair: (weighted L2 error, infinity error). @@ -138,15 +139,13 @@ class GMGPolar : public IGMGPolar /* --------------- */ /* Setup Functions */ - int chooseNumberOfLevels(const PolarGrid& finest_grid); - // Public due to cuda restrictions -public: template void build_rhs_f(const Level& level, HostVector rhs_f, const BoundaryConditions& boundary_conditions, const SourceTerm& source_term); void discretize_rhs_f(const Level& level, HostVector rhs_f); private: + int chooseNumberOfLevels(const PolarGrid& finest_grid); bool checkUniformRefinement(const PolarGrid& grid, double tolerance) const; /* --------------- */ From 1847b8754d5a0f064158ee38647253f38f79f4a7 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Wed, 20 May 2026 23:44:31 +0200 Subject: [PATCH 9/9] Update gmgpolar.h --- include/GMGPolar/gmgpolar.h | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/include/GMGPolar/gmgpolar.h b/include/GMGPolar/gmgpolar.h index a5a18e847..c1d5478d3 100644 --- a/include/GMGPolar/gmgpolar.h +++ b/include/GMGPolar/gmgpolar.h @@ -144,7 +144,13 @@ class GMGPolar : public IGMGPolar const BoundaryConditions& boundary_conditions, const SourceTerm& source_term); void discretize_rhs_f(const Level& level, HostVector rhs_f); + /* --------------- */ + /* Solve Functions */ + void applyExtrapolation(int current_level, HostVector fine_values, HostConstVector coarse_values); + private: + /* --------------- */ + /* Setup Functions */ int chooseNumberOfLevels(const PolarGrid& finest_grid); bool checkUniformRefinement(const PolarGrid& grid, double tolerance) const; @@ -167,8 +173,6 @@ class GMGPolar : public IGMGPolar int iterations); void applyExtrapolatedMultigridIterations(Level& level, MultigridCycleType cycle, int iterations); - // Compute the extrapolated values: u_ex = 4/3 u_fine - 1/3 u_coarse - void applyExtrapolation(int current_level, HostVector fine_values, HostConstVector coarse_values); /* ----------------- */ /* Print information */ @@ -218,4 +222,5 @@ class GMGPolar : public IGMGPolar #include "MultigridMethods/extrapolated_multigrid_V_Cycle.h" #include "MultigridMethods/multigrid_F_Cycle.h" #include "MultigridMethods/multigrid_W_Cycle.h" + } // namespace gmgpolar