diff --git a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h index b36e98a4..246ec4da 100644 --- a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h +++ b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h @@ -44,7 +44,7 @@ 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); + restriction(level.level_depth(), next_level.residual(), residual); // f_{l-1} - A_{l-1}* Inject(u_l) injection(level.level_depth(), next_level.solution(), solution); @@ -80,7 +80,7 @@ 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()); + prolongation(next_level.level_depth(), residual, next_level.error_correction()); /* ----------------------------------- */ /* Compute the corrected approximation */ diff --git a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h index 36e6c67e..e2838367 100644 --- a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h +++ b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h @@ -44,7 +44,7 @@ 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); + restriction(level.level_depth(), next_level.residual(), residual); // f_{l-1} - A_{l-1}* Inject(u_l) injection(level.level_depth(), next_level.solution(), solution); @@ -74,7 +74,7 @@ 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()); + prolongation(next_level.level_depth(), residual, next_level.error_correction()); /* ----------------------------------- */ /* Compute the corrected approximation */ diff --git a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h index 0b157735..a38187dc 100644 --- a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h +++ b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h @@ -44,7 +44,7 @@ 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); + restriction(level.level_depth(), next_level.residual(), residual); // f_{l-1} - A_{l-1}* Inject(u_l) injection(level.level_depth(), next_level.solution(), solution); @@ -80,7 +80,7 @@ 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()); + prolongation(next_level.level_depth(), residual, next_level.error_correction()); /* ----------------------------------- */ /* Compute the corrected approximation */ diff --git a/include/GMGPolar/gmgpolar.h b/include/GMGPolar/gmgpolar.h index e874028b..9e156bf4 100644 --- a/include/GMGPolar/gmgpolar.h +++ b/include/GMGPolar/gmgpolar.h @@ -189,8 +189,6 @@ class GMGPolar : public IGMGPolar 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, Vector result, ConstVector x) const; /* ------------- */ diff --git a/include/GMGPolar/utils.h b/include/GMGPolar/utils.h index 77b60675..e3ddf0f2 100644 --- a/include/GMGPolar/utils.h +++ b/include/GMGPolar/utils.h @@ -34,32 +34,6 @@ void GMGPolar::injection(int current interpolation_->applyInjection(levels_[current_level].grid(), levels_[current_level + 1].grid(), result, x); } -template -void GMGPolar::extrapolatedProlongation(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); -} - -template -void GMGPolar::extrapolatedRestriction(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); -} - template void GMGPolar::FMGInterpolation(int current_level, Vector result, ConstVector x) const diff --git a/include/Interpolation/interpolation.h b/include/Interpolation/interpolation.h index 11a0ca91..057936da 100644 --- a/include/Interpolation/interpolation.h +++ b/include/Interpolation/interpolation.h @@ -23,16 +23,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 b0d90afd..51338e6a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -35,8 +35,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/Interpolation/extrapolated_prolongation.cpp b/src/Interpolation/extrapolated_prolongation.cpp deleted file mode 100644 index 9145fc9a..00000000 --- a/src/Interpolation/extrapolated_prolongation.cpp +++ /dev/null @@ -1,114 +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 - * - */ - -static inline void fineNodeExtrapolatedProlongation(int i_r, int i_theta, int i_r_coarse, int i_theta_coarse, - const PolarGrid& coarse_grid, const PolarGrid& fine_grid, - Vector& fine_result, ConstVector& coarse_values) -{ - 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 */ - } - } -} - -void Interpolation::applyExtrapolatedProlongation(const PolarGrid& coarse_grid, const PolarGrid& fine_grid, - Vector fine_result, ConstVector coarse_values) const -{ - assert(std::ssize(coarse_values) == coarse_grid.numberOfNodes()); - assert(std::ssize(fine_result) == fine_grid.numberOfNodes()); - - /* 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; - fineNodeExtrapolatedProlongation(i_r, i_theta, i_r_coarse, i_theta_coarse, coarse_grid, fine_grid, - fine_result, coarse_values); - } - } - - /* 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; - fineNodeExtrapolatedProlongation(i_r, i_theta, i_r_coarse, i_theta_coarse, coarse_grid, fine_grid, - fine_result, coarse_values); - } - } - } -} diff --git a/src/Interpolation/extrapolated_restriction.cpp b/src/Interpolation/extrapolated_restriction.cpp deleted file mode 100644 index a248f940..00000000 --- a/src/Interpolation/extrapolated_restriction.cpp +++ /dev/null @@ -1,79 +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 - */ - -static inline void coarseNodeExtrapolatedRestriction(int i_r_coarse, int i_theta_coarse, const PolarGrid& fine_grid, - const PolarGrid& coarse_grid, Vector& coarse_result, - ConstVector& fine_values) -{ - 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; -} - -void Interpolation::applyExtrapolatedRestriction(const PolarGrid& fine_grid, const PolarGrid& coarse_grid, - Vector coarse_result, ConstVector fine_values) const -{ - assert(std::ssize(fine_values) == fine_grid.numberOfNodes()); - assert(std::ssize(coarse_result) == coarse_grid.numberOfNodes()); - - /* 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++) { - coarseNodeExtrapolatedRestriction(i_r_coarse, i_theta_coarse, fine_grid, coarse_grid, coarse_result, - fine_values); - } - } - - /* 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++) { - coarseNodeExtrapolatedRestriction(i_r_coarse, i_theta_coarse, fine_grid, coarse_grid, coarse_result, - fine_values); - } - } - } -} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 038278b1..27559661 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 << ")"; - } - } -}