diff --git a/include/GMGPolar/gmgpolar.h b/include/GMGPolar/gmgpolar.h index c2d90a29..1853346c 100644 --- a/include/GMGPolar/gmgpolar.h +++ b/include/GMGPolar/gmgpolar.h @@ -296,8 +296,6 @@ class GMGPolar void prolongation(const int current_level, Vector result, ConstVector x) const; void restriction(const int current_level, Vector result, ConstVector x) const; void injection(const int current_level, Vector result, ConstVector x) const; - void extrapolatedProlongation(const int current_level, Vector result, ConstVector x) const; - void extrapolatedRestriction(const int current_level, Vector result, ConstVector x) const; void FMGInterpolation(const int current_level, Vector result, ConstVector x) const; /* ------------- */ diff --git a/include/Interpolation/interpolation.h b/include/Interpolation/interpolation.h index 92e83f75..fccea283 100644 --- a/include/Interpolation/interpolation.h +++ b/include/Interpolation/interpolation.h @@ -26,16 +26,10 @@ class Interpolation 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, Vector coarse_result, ConstVector fine_values) const; - void applyExtrapolatedRestriction(const PolarGrid& fine_grid, const PolarGrid& coarse_grid, - Vector coarse_result, ConstVector fine_values) const; - /* Bicubic FMG interpolator 1/16 * [-1, 9, 9, -1] */ void applyFMGInterpolation(const PolarGrid& coarse_grid, const PolarGrid& fine_grid, Vector fine_result, ConstVector coarse_values) const; diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index eb06eece..f0750574 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -43,8 +43,6 @@ set(STENCIL_SOURCES # file(GLOB_RECURSE INTERPOLATION_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/Interpolation/*.cpp) set(INTERPOLATION_SOURCES - ${CMAKE_CURRENT_SOURCE_DIR}/Interpolation/extrapolated_prolongation.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/Interpolation/extrapolated_restriction.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Interpolation/fmg_interpolation.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Interpolation/injection.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Interpolation/interpolation.cpp diff --git a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp index 63b00848..458369d7 100644 --- a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp @@ -45,7 +45,7 @@ void GMGPolar::implicitlyExtrapolatedMultigrid_F_Cycle(const int level_depth, Ve // P_ex^T (f_l - A_l*u_l) level.computeResidual(residual, rhs, solution); - extrapolatedRestriction(level_depth, next_level.residual(), residual); + restriction(level_depth, next_level.residual(), residual); // f_{l-1} - A_{l-1}* Inject(u_l) injection(level_depth, next_level.solution(), solution); @@ -76,7 +76,7 @@ void GMGPolar::implicitlyExtrapolatedMultigrid_F_Cycle(const int level_depth, Ve // P_ex^T (f_l - A_l*u_l) level.computeResidual(residual, rhs, solution); - extrapolatedRestriction(level_depth, next_level.error_correction(), residual); + restriction(level_depth, next_level.error_correction(), residual); // f_{l-1} - A_{l-1}* Inject(u_l) injection(level_depth, next_level.solution(), solution); @@ -98,7 +98,7 @@ void GMGPolar::implicitlyExtrapolatedMultigrid_F_Cycle(const int level_depth, Ve } /* Interpolate the correction */ - extrapolatedProlongation(level_depth + 1, residual, next_level.residual()); + prolongation(level_depth + 1, residual, next_level.residual()); /* Compute the corrected approximation: u = u + error */ add(solution, ConstVector(residual)); diff --git a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp index 8ef84850..13131545 100644 --- a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp @@ -45,7 +45,7 @@ void GMGPolar::implicitlyExtrapolatedMultigrid_V_Cycle(const int level_depth, Ve // P_ex^T (f_l - A_l*u_l) level.computeResidual(residual, rhs, solution); - extrapolatedRestriction(level_depth, next_level.residual(), residual); + restriction(level_depth, next_level.residual(), residual); // f_{l-1} - A_{l-1}* Inject(u_l) injection(level_depth, next_level.solution(), solution); @@ -76,7 +76,7 @@ void GMGPolar::implicitlyExtrapolatedMultigrid_V_Cycle(const int level_depth, Ve // P_ex^T (f_l - A_l*u_l) level.computeResidual(residual, rhs, solution); - extrapolatedRestriction(level_depth, next_level.error_correction(), residual); + restriction(level_depth, next_level.error_correction(), residual); // f_{l-1} - A_{l-1}* Inject(u_l) injection(level_depth, next_level.solution(), solution); @@ -97,7 +97,7 @@ void GMGPolar::implicitlyExtrapolatedMultigrid_V_Cycle(const int level_depth, Ve } /* Interpolate the correction */ - extrapolatedProlongation(level_depth + 1, residual, next_level.residual()); + prolongation(level_depth + 1, residual, next_level.residual()); /* Compute the corrected approximation: u = u + error */ add(solution, ConstVector(residual)); diff --git a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp index cffcccc7..60ed0ff4 100644 --- a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp @@ -45,7 +45,7 @@ void GMGPolar::implicitlyExtrapolatedMultigrid_W_Cycle(const int level_depth, Ve // P_ex^T (f_l - A_l*u_l) level.computeResidual(residual, rhs, solution); - extrapolatedRestriction(level_depth, next_level.residual(), residual); + restriction(level_depth, next_level.residual(), residual); // f_{l-1} - A_{l-1}* Inject(u_l) injection(level_depth, next_level.solution(), solution); @@ -76,7 +76,7 @@ void GMGPolar::implicitlyExtrapolatedMultigrid_W_Cycle(const int level_depth, Ve // P_ex^T (f_l - A_l*u_l) level.computeResidual(residual, rhs, solution); - extrapolatedRestriction(level_depth, next_level.error_correction(), residual); + restriction(level_depth, next_level.error_correction(), residual); // f_{l-1} - A_{l-1}* Inject(u_l) injection(level_depth, next_level.solution(), solution); @@ -98,7 +98,7 @@ void GMGPolar::implicitlyExtrapolatedMultigrid_W_Cycle(const int level_depth, Ve } /* Interpolate the correction */ - extrapolatedProlongation(level_depth + 1, residual, next_level.residual()); + prolongation(level_depth + 1, residual, next_level.residual()); /* Compute the corrected approximation: u = u + error */ add(solution, ConstVector(residual)); diff --git a/src/GMGPolar/level_interpolation.cpp b/src/GMGPolar/level_interpolation.cpp index 58b311ca..fdf544ae 100644 --- a/src/GMGPolar/level_interpolation.cpp +++ b/src/GMGPolar/level_interpolation.cpp @@ -27,26 +27,6 @@ void GMGPolar::injection(const int current_level, Vector result, ConstVe interpolation_->applyInjection(levels_[current_level].grid(), levels_[current_level + 1].grid(), result, x); } -void GMGPolar::extrapolatedProlongation(const 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."); - - interpolation_->applyExtrapolatedProlongation(levels_[current_level].grid(), levels_[current_level - 1].grid(), - result, x); -} - -void GMGPolar::extrapolatedRestriction(const 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."); - - interpolation_->applyExtrapolatedRestriction(levels_[current_level].grid(), levels_[current_level + 1].grid(), - result, x); -} - void GMGPolar::FMGInterpolation(const int current_level, Vector result, ConstVector x) const { assert(current_level < number_of_levels_ && 1 <= current_level); diff --git a/src/Interpolation/extrapolated_prolongation.cpp b/src/Interpolation/extrapolated_prolongation.cpp deleted file mode 100644 index 343770a1..00000000 --- a/src/Interpolation/extrapolated_prolongation.cpp +++ /dev/null @@ -1,115 +0,0 @@ -#include "../../include/Interpolation/interpolation.h" - -/* - * Extrapolated Prolongation Operator - * ---------------------------------- - * - * Extrapolated prolongation is used between the finest most grids in the multigrid hierarchy. - * It assumes the fine grid comes from a uniform refinement of the coarse grid, so all spacings are equal. - * Thus fine values can be computed simply by averaging the neighboring coarse nodes. - * - * A fine node is classified by the parity of its (i_r, i_theta) indices: - * - * 1) (even, even) - * Fine node coincides with a coarse node - * -> copy value - * - * 2) (odd, even) - * Node lies between two coarse nodes in radial direction - * - * X ---- O ---- X - * - * -> arithmetic mean of left + right coarse node - * - * 3) (even, odd) - * Node lies between two coarse nodes in angular direction - * - * X - * | - * O - * | - * X - * - * -> arithmetic mean of bottom + top coarse node - * - * 4) (odd, odd) - * Node lies inside a coarse cell - * We extrapolate/average across the diagonal: - * - * X - * \ - * O - * \ - * X - * - * -> arithmetic mean of bottom right + top left coarse node - * - */ - -#define FINE_NODE_EXTRAPOLATED_PROLONGATION() \ - do { \ - if (i_r & 1) { \ - if (i_theta & 1) { \ - /* (odd, odd) -> node in center of coarse cell */ \ - double value = \ - 0.5 * (coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse)] + /* Bottom right */ \ - coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 1)] /* Top left */ \ - ); \ - fine_result[fine_grid.index(i_r, i_theta)] = value; \ - } \ - else { \ - /* (odd, even) -> between coarse nodes in radial direction */ \ - double value = 0.5 * (coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* Left */ \ - coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse)] /* Right */ \ - ); \ - fine_result[fine_grid.index(i_r, i_theta)] = value; \ - } \ - } \ - else { \ - if (i_theta & 1) { \ - /* (even, odd) -> between coarse nodes in angular direction */ \ - double value = 0.5 * (coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* Bottom */ \ - coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 1)] /* Top */ \ - ); \ - fine_result[fine_grid.index(i_r, i_theta)] = value; \ - } \ - else { \ - /* (even, even) -> node lies exactly on coarse grid */ \ - fine_result[fine_grid.index(i_r, i_theta)] = \ - coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)]; /* Center */ \ - } \ - } \ - } while (0) - -void Interpolation::applyExtrapolatedProlongation(const PolarGrid& coarse_grid, const PolarGrid& fine_grid, - Vector fine_result, ConstVector coarse_values) const -{ - assert(coarse_values.size() == static_cast(coarse_grid.numberOfNodes())); - assert(fine_result.size() == static_cast(fine_grid.numberOfNodes())); - - /* We split the loops into two regions to better respect the */ - /* access patterns of the smoother and improve cache locality. */ - -#pragma omp parallel num_threads(max_omp_threads_) - { - /* Circular Indexing Section */ -#pragma omp for nowait - for (int i_r = 0; i_r < fine_grid.numberSmootherCircles(); i_r++) { - int i_r_coarse = i_r / 2; - for (int i_theta = 0; i_theta < fine_grid.ntheta(); i_theta++) { - int i_theta_coarse = i_theta / 2; - FINE_NODE_EXTRAPOLATED_PROLONGATION(); - } - } - - /* Radial Indexing Section */ -#pragma omp for nowait - for (int i_theta = 0; i_theta < fine_grid.ntheta(); i_theta++) { - int i_theta_coarse = i_theta / 2; - for (int i_r = fine_grid.numberSmootherCircles(); i_r < fine_grid.nr(); i_r++) { - int i_r_coarse = i_r / 2; - FINE_NODE_EXTRAPOLATED_PROLONGATION(); - } - } - } -} diff --git a/src/Interpolation/extrapolated_restriction.cpp b/src/Interpolation/extrapolated_restriction.cpp deleted file mode 100644 index ae692536..00000000 --- a/src/Interpolation/extrapolated_restriction.cpp +++ /dev/null @@ -1,76 +0,0 @@ -#include "../../include/Interpolation/interpolation.h" - -/* - * Extrapolated Restriction Operator - * ---------------------------------- - * - * This is the transpose of the extrapolated prolongation operator: R = P^T - * Used between the finest grids in the multigrid hierarchy where uniform refinement is assumed. - * All spacings are equal, so weights are simple factors of 0.5. - * Each coarse node accumulates contributions from its corresponding 9 surrounding fine nodes. - * - * The stencil accumulates: - * - Center (weight 1.0) - * - Left, Right, Bottom, Top (weight 0.5 each) - * - Bottom-Right and Top-Left diagonal (weight 0.5 each) - * - * Note: Bottom-Left and Top-Right are NOT included (consistent with diagonal averaging in prolongation) - * - * Boundary handling: - * - Angular direction: periodic (wrapping) - * - Radial direction: check domain boundaries - */ - -#define COARSE_NODE_EXTRAPOLATED_RESTRICTION() \ - do { \ - int i_r = i_r_coarse * 2; \ - int i_theta = i_theta_coarse * 2; \ - \ - /* Center + Angular contributions (always present) */ \ - double value = fine_values[fine_grid.index(i_r, i_theta)] + \ - 0.5 * fine_values[fine_grid.index(i_r, i_theta - 1)] + \ - 0.5 * fine_values[fine_grid.index(i_r, i_theta + 1)]; \ - \ - /* Left contributions (if not at inner boundary) */ \ - if (i_r_coarse > 0) { \ - value += 0.5 * fine_values[fine_grid.index(i_r - 1, i_theta)] + \ - 0.5 * fine_values[fine_grid.index(i_r - 1, i_theta + 1)]; /* Top-Left diagonal */ \ - } \ - \ - /* Right contributions (if not at outer boundary) */ \ - if (i_r_coarse < coarse_grid.nr() - 1) { \ - value += 0.5 * fine_values[fine_grid.index(i_r + 1, i_theta)] + \ - 0.5 * fine_values[fine_grid.index(i_r + 1, i_theta - 1)]; /* Bottom-Right diagonal */ \ - } \ - \ - coarse_result[coarse_grid.index(i_r_coarse, i_theta_coarse)] = value; \ - } while (0) - -void Interpolation::applyExtrapolatedRestriction(const PolarGrid& fine_grid, const PolarGrid& coarse_grid, - Vector coarse_result, ConstVector fine_values) const -{ - assert(fine_values.size() == static_cast(fine_grid.numberOfNodes())); - assert(coarse_result.size() == static_cast(coarse_grid.numberOfNodes())); - - /* We split the loops into two regions to better respect the */ - /* access patterns of the smoother and improve cache locality. */ - -#pragma omp parallel num_threads(max_omp_threads_) - { - /* Circular Indexing Section */ -#pragma omp for nowait - for (int i_r_coarse = 0; i_r_coarse < coarse_grid.numberSmootherCircles(); i_r_coarse++) { - for (int i_theta_coarse = 0; i_theta_coarse < coarse_grid.ntheta(); i_theta_coarse++) { - COARSE_NODE_EXTRAPOLATED_RESTRICTION(); - } - } - - /* Radial Indexing Section */ -#pragma omp for nowait - for (int i_theta_coarse = 0; i_theta_coarse < coarse_grid.ntheta(); i_theta_coarse++) { - for (int i_r_coarse = coarse_grid.numberSmootherCircles(); i_r_coarse < coarse_grid.nr(); i_r_coarse++) { - COARSE_NODE_EXTRAPOLATED_RESTRICTION(); - } - } - } -} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index fad28d0a..828398ce 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -17,8 +17,6 @@ add_executable(gmgpolar_tests PolarGrid/polargrid.cpp Interpolation/prolongation.cpp Interpolation/restriction.cpp - Interpolation/extrapolated_prolongation.cpp - Interpolation/extrapolated_restriction.cpp Residual/residual.cpp DirectSolver/directSolver.cpp DirectSolver/directSolverNoMumps.cpp diff --git a/tests/Interpolation/extrapolated_prolongation.cpp b/tests/Interpolation/extrapolated_prolongation.cpp deleted file mode 100644 index 4650f7fb..00000000 --- a/tests/Interpolation/extrapolated_prolongation.cpp +++ /dev/null @@ -1,65 +0,0 @@ -#include -#include - -#include "../test_tools.h" - -#include "../../include/GMGPolar/gmgpolar.h" -#include "../../include/Interpolation/interpolation.h" -#include "../../include/InputFunctions/DensityProfileCoefficients/poissonCoefficients.h" - -// Helper that computes the mathematically expected extrapolated prolongation value -static double expected_extrapolated_value(const PolarGrid& coarse, const PolarGrid& fine, - ConstVector coarse_vals, int i_r, int i_theta) -{ - int i_r_coarse = i_r / 2; - int i_theta_coarse = i_theta / 2; - - bool r_even = (i_r % 2 == 0); - bool t_even = (i_theta % 2 == 0); - - if (r_even && t_even) { - // Node coincides with a coarse node - return coarse_vals[coarse.index(i_r_coarse, i_theta_coarse)]; - } - - if (!r_even && t_even) { - // Radial midpoint - arithmetic mean of left and right - return 0.5 * (coarse_vals[coarse.index(i_r_coarse, i_theta_coarse)] + - coarse_vals[coarse.index(i_r_coarse + 1, i_theta_coarse)]); - } - - if (r_even && !t_even) { - // Angular midpoint - arithmetic mean of bottom and top - return 0.5 * (coarse_vals[coarse.index(i_r_coarse, i_theta_coarse)] + - coarse_vals[coarse.index(i_r_coarse, i_theta_coarse + 1)]); - } - - // Center of coarse cell - arithmetic mean of diagonal nodes (bottom-right + top-left) - return 0.5 * (coarse_vals[coarse.index(i_r_coarse + 1, i_theta_coarse)] + - coarse_vals[coarse.index(i_r_coarse, i_theta_coarse + 1)]); -} - -TEST(ExtrapolatedProlongationTest, ExtrapolatedProlongationMatchesStencil) -{ - std::vector fine_radii = {0.1, 0.2, 0.25, 0.5, 0.8, 0.9, 1.3, 1.4, 2.0}; - 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); - - Interpolation I(/*threads*/ 16, /*DirBC*/ true); - - Vector coarse_values = generate_random_sample_data(coarse_grid, 1234, 0.0, 1.0); - Vector fine_result("fine_result", fine_grid.numberOfNodes()); - - I.applyExtrapolatedProlongation(coarse_grid, fine_grid, fine_result, coarse_values); - - for (int i_r = 0; i_r < fine_grid.nr(); ++i_r) { - for (int i_theta = 0; i_theta < fine_grid.ntheta(); ++i_theta) { - double expected = expected_extrapolated_value(coarse_grid, fine_grid, coarse_values, i_r, i_theta); - double got = fine_result[fine_grid.index(i_r, i_theta)]; - ASSERT_NEAR(expected, got, 1e-10) << "Mismatch at (" << i_r << ", " << i_theta << ")"; - } - } -} diff --git a/tests/Interpolation/extrapolated_restriction.cpp b/tests/Interpolation/extrapolated_restriction.cpp deleted file mode 100644 index 2af24650..00000000 --- a/tests/Interpolation/extrapolated_restriction.cpp +++ /dev/null @@ -1,64 +0,0 @@ -#include -#include - -#include "../test_tools.h" - -#include "../../include/GMGPolar/gmgpolar.h" -#include "../../include/Interpolation/interpolation.h" -#include "../../include/InputFunctions/DensityProfileCoefficients/poissonCoefficients.h" - -// Helper that computes the mathematically expected extrapolated restriction value -static double expected_extrapolated_restriction_value(const PolarGrid& fine, const PolarGrid& coarse, - ConstVector fine_vals, int i_r_coarse, int i_theta_coarse) -{ - int i_r = i_r_coarse * 2; - int i_theta = i_theta_coarse * 2; - - // Angular indices with periodic wrapping - int i_theta_M1 = fine.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = fine.wrapThetaIndex(i_theta + 1); - - // Center + Angular contributions (always present) - double value = fine_vals[fine.index(i_r, i_theta)] + 0.5 * fine_vals[fine.index(i_r, i_theta_M1)] + - 0.5 * fine_vals[fine.index(i_r, i_theta_P1)]; - - // Left contributions (if not at inner boundary) - if (i_r_coarse > 0) { - value += 0.5 * fine_vals[fine.index(i_r - 1, i_theta)] + - 0.5 * fine_vals[fine.index(i_r - 1, i_theta_P1)]; // Top-Left diagonal - } - - // Right contributions (if not at outer boundary) - if (i_r_coarse < coarse.nr() - 1) { - value += 0.5 * fine_vals[fine.index(i_r + 1, i_theta)] + - 0.5 * fine_vals[fine.index(i_r + 1, i_theta_M1)]; // Bottom-Right diagonal - } - - return value; -} - -TEST(ExtrapolatedRestrictionTest, ExtrapolatedRestrictionMatchesStencil) -{ - std::vector fine_radii = {0.1, 0.2, 0.25, 0.5, 0.8, 0.9, 1.3, 1.4, 2.0}; - 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); - - Interpolation I(/*threads*/ 16, /*DirBC*/ true); - - Vector fine_values = generate_random_sample_data(fine_grid, 9012, 0.0, 1.0); - Vector coarse_result("coarse_result", coarse_grid.numberOfNodes()); - - I.applyExtrapolatedRestriction(fine_grid, coarse_grid, coarse_result, fine_values); - - for (int i_r_coarse = 0; i_r_coarse < coarse_grid.nr(); ++i_r_coarse) { - for (int i_theta_coarse = 0; i_theta_coarse < coarse_grid.ntheta(); ++i_theta_coarse) { - double expected = expected_extrapolated_restriction_value(fine_grid, coarse_grid, fine_values, i_r_coarse, - i_theta_coarse); - double got = coarse_result[coarse_grid.index(i_r_coarse, i_theta_coarse)]; - ASSERT_NEAR(expected, got, 1e-10) << "Mismatch at (" << i_r_coarse << ", " << i_theta_coarse << ")"; - } - } -}