From ebaacb60a6fba95ef4827812f4e1794438bbcbee Mon Sep 17 00:00:00 2001 From: julianlitz Date: Sat, 20 Dec 2025 16:48:10 +0000 Subject: [PATCH] Remove unoptimized interpolation --- include/Interpolation/interpolation.h | 10 +- .../extrapolated_prolongation.cpp | 148 ++++----- .../extrapolated_restriction.cpp | 196 +++-------- src/Interpolation/prolongation.cpp | 221 +++++-------- src/Interpolation/restriction.cpp | 303 ++++-------------- .../extrapolated_prolongation.cpp | 90 +++--- .../extrapolated_restriction.cpp | 159 +++------ tests/Interpolation/prolongation.cpp | 100 ++++-- tests/Interpolation/restriction.cpp | 180 ++++------- 9 files changed, 510 insertions(+), 897 deletions(-) diff --git a/include/Interpolation/interpolation.h b/include/Interpolation/interpolation.h index 10e093de..92e83f75 100644 --- a/include/Interpolation/interpolation.h +++ b/include/Interpolation/interpolation.h @@ -23,22 +23,16 @@ class Interpolation ConstVector fine_values) const; /* Bilinear interpolation operator */ - void applyProlongation0(const PolarGrid& coarse_grid, const PolarGrid& fine_grid, Vector fine_result, - ConstVector coarse_values) const; void applyProlongation(const PolarGrid& coarse_grid, const PolarGrid& fine_grid, Vector fine_result, ConstVector coarse_values) const; - void applyExtrapolatedProlongation0(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 applyRestriction0(const PolarGrid& fine_grid, const PolarGrid& coarse_grid, Vector coarse_result, - ConstVector fine_values) const; void applyRestriction(const PolarGrid& fine_grid, const PolarGrid& coarse_grid, Vector coarse_result, ConstVector fine_values) const; - void applyExtrapolatedRestriction0(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; diff --git a/src/Interpolation/extrapolated_prolongation.cpp b/src/Interpolation/extrapolated_prolongation.cpp index 2c406099..343770a1 100644 --- a/src/Interpolation/extrapolated_prolongation.cpp +++ b/src/Interpolation/extrapolated_prolongation.cpp @@ -1,101 +1,82 @@ #include "../../include/Interpolation/interpolation.h" -void Interpolation::applyExtrapolatedProlongation0(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())); - -#pragma omp parallel for num_threads(max_omp_threads_) - for (int index = 0; index < fine_grid.numberOfNodes(); index++) { - std::array, space_dimension> neighbor_distance; - - MultiIndex fine_node = fine_grid.multiIndex(index); - MultiIndex coarse_node(fine_node[0] / 2, fine_node[1] / 2); // Nearest lower left coarse node in the fine grid. - - if (fine_node[0] % 2 == 0 && fine_node[1] % 2 == 0) { - // Fine node appears in coarse grid - fine_result[index] = coarse_values[coarse_grid.index(coarse_node)]; - } - - if (fine_node[0] % 2 == 0 && fine_node[1] % 2 == 1) { - // Fine node between two coarse nodes in theta direction - // X - // | - // O - // | - // X - MultiIndex bottomNeighbor(coarse_node[0], coarse_node[1]); - MultiIndex topNeighbor(coarse_node[0], (coarse_node[1] + 1) % coarse_grid.ntheta()); - fine_result[index] = 0.5 * (coarse_values[coarse_grid.index(bottomNeighbor)] + - coarse_values[coarse_grid.index(topNeighbor)]); - } - - if (fine_node[0] % 2 == 1 && fine_node[1] % 2 == 0) { - // Fine node between two coarse nodes in radial direction - // X -- O -- X - MultiIndex leftNeighbor(coarse_node[0], coarse_node[1]); - MultiIndex rightNeighbor(coarse_node[0] + 1, coarse_node[1]); - fine_result[index] = 0.5 * (coarse_values[coarse_grid.index(leftNeighbor)] + - coarse_values[coarse_grid.index(rightNeighbor)]); - } - - if (fine_node[0] % 2 == 1 && fine_node[1] % 2 == 1) { - // Interpolates a fine node value based on two neighboring coarse nodes. - // Fine node lies in the center of four coarse nodes forming a cross shape: - // - // X - /* \ */ - // O <-- Fine Node (i_r, i_theta) - /* \ */ - // X - // - MultiIndex bottom_right_neighbor(coarse_node[0] + 1, coarse_node[1]); - MultiIndex top_left_neighbor(coarse_node[0], (coarse_node[1] + 1) % coarse_grid.ntheta()); - fine_result[index] = 0.5 * (coarse_values[coarse_grid.index(bottom_right_neighbor)] + - coarse_values[coarse_grid.index(top_left_neighbor)]); - } - } -} - -// --------------------------------------- // -// Optimized version of applyProlongation0 // -// --------------------------------------- // +/* + * 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) { \ - /* i_r % 2 == 1, i_theta % 2 == 1 */ \ - /* Fine node in the center of four coarse nodes */ \ - fine_result[fine_grid.index(i_r, i_theta)] = \ + /* (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 { \ - /* i_r % 2 == 1, i_theta % 2 == 0 */ \ - /* Fine node between coarse nodes in radial direction */ \ - fine_result[fine_grid.index(i_r, i_theta)] = \ - 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 */ \ - ); \ + /* (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) { \ - /* i_r % 2 == 0, i_theta % 2 == 1 */ \ - /* Fine node between coarse nodes in theta direction */ \ - fine_result[fine_grid.index(i_r, i_theta)] = \ - 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 */ \ - ); \ + /* (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 { \ - /* i_r % 2 == 0, i_theta % 2 == 0 */ \ - /* Fine node appears in coarse grid */ \ + /* (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 */ \ + coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)]; /* Center */ \ } \ } \ } while (0) @@ -106,10 +87,12 @@ void Interpolation::applyExtrapolatedProlongation(const PolarGrid& coarse_grid, 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_) { -/* Circluar Indexing Section */ -/* For loop matches circular access pattern */ + /* 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; @@ -119,8 +102,7 @@ void Interpolation::applyExtrapolatedProlongation(const PolarGrid& coarse_grid, } } -/* Radial Indexing Section */ -/* For loop matches radial access pattern */ + /* 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; diff --git a/src/Interpolation/extrapolated_restriction.cpp b/src/Interpolation/extrapolated_restriction.cpp index f44ea56d..ae692536 100644 --- a/src/Interpolation/extrapolated_restriction.cpp +++ b/src/Interpolation/extrapolated_restriction.cpp @@ -1,64 +1,50 @@ #include "../../include/Interpolation/interpolation.h" -/* For the restriction we use R_ex = P_ex^T */ - -void Interpolation::applyExtrapolatedRestriction0(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())); - -#pragma omp parallel for num_threads(max_omp_threads_) - for (int index = 0; index < coarse_grid.numberOfNodes(); index++) { - MultiIndex coarse_node = coarse_grid.multiIndex(index); - MultiIndex fine_node(2 * coarse_node[0], 2 * coarse_node[1]); - - std::array, space_dimension> neighbors; - - // Center - double value = fine_values[fine_grid.index(fine_node)]; - - fine_grid.adjacentNeighborsOf(fine_node, neighbors); - - // Left - if (neighbors[0].first != -1) { - value += 0.5 * fine_values[neighbors[0].first]; - } - - // Right - if (neighbors[0].second != -1) { - value += 0.5 * fine_values[neighbors[0].second]; - } - - // Bottom - if (neighbors[1].first != -1) { - value += 0.5 * fine_values[neighbors[1].first]; - } - - // Top - if (neighbors[1].second != -1) { - value += 0.5 * fine_values[neighbors[1].second]; - } - - fine_grid.diagonalNeighborsOf(fine_node, neighbors); - - // Bottom Right - if (neighbors[0].second != -1) { - value += 0.5 * fine_values[neighbors[0].second]; - } - - // Top Left - if (neighbors[1].first != -1) { - value += 0.5 * fine_values[neighbors[1].first]; - } - - coarse_result[index] = value; - } -} - -// -------------------------------------- // -// Optimized version of applyRestriction0 // -// -------------------------------------- // +/* + * 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 @@ -66,100 +52,24 @@ void Interpolation::applyExtrapolatedRestriction(const PolarGrid& fine_grid, con assert(fine_values.size() == static_cast(fine_grid.numberOfNodes())); assert(coarse_result.size() == static_cast(coarse_grid.numberOfNodes())); - const int coarseNumberSmootherCircles = coarse_grid.numberSmootherCircles(); + /* 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_) { -/* For loop matches circular access pattern */ + /* Circular Indexing Section */ #pragma omp for nowait - for (int i_r_coarse = 0; i_r_coarse < coarseNumberSmootherCircles; i_r_coarse++) { - int i_r = i_r_coarse * 2; + 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++) { - int i_theta = i_theta_coarse * 2; - - if (0 < i_r_coarse && i_r_coarse < coarseNumberSmootherCircles - 1) { - int i_theta_M1 = fine_grid.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = fine_grid.wrapThetaIndex(i_theta + 1); - - coarse_result[coarse_grid.index(i_r_coarse, i_theta_coarse)] = - // Center - fine_values[fine_grid.index(i_r, i_theta)] + - // Left, Right, Bottom, Top - 0.5 * fine_values[fine_grid.index(i_r - 1, i_theta)] + - 0.5 * fine_values[fine_grid.index(i_r + 1, i_theta)] + - 0.5 * fine_values[fine_grid.index(i_r, i_theta_M1)] + - 0.5 * fine_values[fine_grid.index(i_r, i_theta_P1)] + - // Bottom Right, Top Left - 0.5 * fine_values[fine_grid.index(i_r + 1, i_theta_M1)] + - 0.5 * fine_values[fine_grid.index(i_r - 1, i_theta_P1)]; - } - else { - /* First and Last Circle have to be checked for domain boundary */ - int i_theta_M1 = fine_grid.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = fine_grid.wrapThetaIndex(i_theta + 1); - // Center, Bottom, Top - double value = fine_values[fine_grid.index(i_r, i_theta)] + - 0.5 * fine_values[fine_grid.index(i_r, i_theta_M1)] + - 0.5 * fine_values[fine_grid.index(i_r, i_theta_P1)]; - - if (i_r_coarse > 0) { - // Left, Top Left - 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_P1)]; - } - if (i_r_coarse < coarse_grid.nr() - 1) { - // Right, Bottom Right - 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_M1)]; - } - coarse_result[coarse_grid.index(i_r_coarse, i_theta_coarse)] = value; - } + COARSE_NODE_EXTRAPOLATED_RESTRICTION(); } } -/* For loop matches radial access pattern */ + /* Radial Indexing Section */ #pragma omp for nowait for (int i_theta_coarse = 0; i_theta_coarse < coarse_grid.ntheta(); i_theta_coarse++) { - int i_theta = i_theta_coarse * 2; - for (int i_r_coarse = coarseNumberSmootherCircles; i_r_coarse < coarse_grid.nr(); i_r_coarse++) { - int i_r = i_r_coarse * 2; - - if (coarse_grid.numberSmootherCircles() < i_r_coarse && i_r_coarse < coarse_grid.nr() - 1) { - int i_theta_M1 = fine_grid.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = fine_grid.wrapThetaIndex(i_theta + 1); - - coarse_result[coarse_grid.index(i_r_coarse, i_theta_coarse)] = - // Center - fine_values[fine_grid.index(i_r, i_theta)] + - // Left, Right, Bottom, Top - 0.5 * fine_values[fine_grid.index(i_r - 1, i_theta)] + - 0.5 * fine_values[fine_grid.index(i_r + 1, i_theta)] + - 0.5 * fine_values[fine_grid.index(i_r, i_theta_M1)] + - 0.5 * fine_values[fine_grid.index(i_r, i_theta_P1)] + - // Bottom Right, Top Left - 0.5 * fine_values[fine_grid.index(i_r + 1, i_theta_M1)] + - 0.5 * fine_values[fine_grid.index(i_r - 1, i_theta_P1)]; - } - else { - /* First and Last radial nodes have to be checked for domain boundary */ - int i_theta_M1 = fine_grid.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = fine_grid.wrapThetaIndex(i_theta + 1); - // Center, Bottom, Top - double value = fine_values[fine_grid.index(i_r, i_theta)] + - 0.5 * fine_values[fine_grid.index(i_r, i_theta_M1)] + - 0.5 * fine_values[fine_grid.index(i_r, i_theta_P1)]; - if (i_r_coarse > 0) { - // Left, Top Left - 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_P1)]; - } - if (i_r_coarse < coarse_grid.nr() - 1) { - // Right, Bottom Right - 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_M1)]; - } - coarse_result[coarse_grid.index(i_r_coarse, i_theta_coarse)] = value; - } + for (int i_r_coarse = coarse_grid.numberSmootherCircles(); i_r_coarse < coarse_grid.nr(); i_r_coarse++) { + COARSE_NODE_EXTRAPOLATED_RESTRICTION(); } } } diff --git a/src/Interpolation/prolongation.cpp b/src/Interpolation/prolongation.cpp index 79ebca10..d0cf221f 100644 --- a/src/Interpolation/prolongation.cpp +++ b/src/Interpolation/prolongation.cpp @@ -1,150 +1,106 @@ #include "../../include/Interpolation/interpolation.h" -// We use the anisotropic bilinear interpolation stencil. -// For an isotropic mesh, this stencil reduces to -// |1 2 1| -// P = 1/4 * |2 4 2| -// |1 2 1| - -void Interpolation::applyProlongation0(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())); - -#pragma omp parallel for num_threads(max_omp_threads_) - for (int index = 0; index < fine_grid.numberOfNodes(); index++) { - std::array, space_dimension> neighbor_distance; - - MultiIndex fine_node = fine_grid.multiIndex(index); - MultiIndex coarse_node(fine_node[0] / 2, fine_node[1] / 2); // Nearest lower left coarse node in the fine grid. - - if (fine_node[0] % 2 == 0 && fine_node[1] % 2 == 0) { - // Fine node appears in coarse grid - fine_result[index] = coarse_values[coarse_grid.index(coarse_node)]; - } - - if (fine_node[0] % 2 == 0 && fine_node[1] % 2 == 1) { - // Fine node between two coarse nodes in theta direction - // X - // | - // O - // | - // X - fine_grid.adjacentNeighborDistances(fine_node, neighbor_distance); - - double k1 = neighbor_distance[1].first; - double k2 = neighbor_distance[1].second; - - MultiIndex bottomNeighbor(coarse_node[0], coarse_node[1]); - MultiIndex topNeighbor(coarse_node[0], (coarse_node[1] + 1) % coarse_grid.ntheta()); - - fine_result[index] = (k1 * coarse_values[coarse_grid.index(bottomNeighbor)] + - k2 * coarse_values[coarse_grid.index(topNeighbor)]) / - (k1 + k2); - } - - if (fine_node[0] % 2 == 1 && fine_node[1] % 2 == 0) { - // Fine node between two coarse nodes in radial direction - // X -- O -- X - fine_grid.adjacentNeighborDistances(fine_node, neighbor_distance); - - double h1 = neighbor_distance[0].first; - double h2 = neighbor_distance[0].second; - - MultiIndex leftNeighbor(coarse_node[0], coarse_node[1]); - MultiIndex rightNeighbor(coarse_node[0] + 1, coarse_node[1]); - - fine_result[index] = (h1 * coarse_values[coarse_grid.index(leftNeighbor)] + - h2 * coarse_values[coarse_grid.index(rightNeighbor)]) / - (h1 + h2); - } - - if (fine_node[0] % 2 == 1 && fine_node[1] % 2 == 1) { - // Interpolates a fine node value based on four neighboring coarse nodes. - // Fine node lies in the center of four coarse nodes forming a cross shape: - // - // X X - /* \ / */ - // O <-- Fine Node (i_r, i_theta) - /* / \ */ - // X X - // - fine_grid.adjacentNeighborDistances(fine_node, neighbor_distance); - - double h1 = neighbor_distance[0].first; - double h2 = neighbor_distance[0].second; - double k1 = neighbor_distance[1].first; - double k2 = neighbor_distance[1].second; - - MultiIndex bottom_left_neighbor(coarse_node[0], coarse_node[1]); - MultiIndex bottom_right_neighbor(coarse_node[0] + 1, coarse_node[1]); - MultiIndex top_left_neighbor(coarse_node[0], (coarse_node[1] + 1) % coarse_grid.ntheta()); - MultiIndex top_right_neighbor(coarse_node[0] + 1, (coarse_node[1] + 1) % coarse_grid.ntheta()); - - fine_result[index] = (h1 * k1 * coarse_values[coarse_grid.index(bottom_left_neighbor)] + - h2 * k1 * coarse_values[coarse_grid.index(bottom_right_neighbor)] + - h1 * k2 * coarse_values[coarse_grid.index(top_left_neighbor)] + - h2 * k2 * coarse_values[coarse_grid.index(top_right_neighbor)]) / - ((h1 + h2) * (k1 + k2)); - } - } -} - -// --------------------------------------- // -// Optimized version of applyProlongation0 // -// --------------------------------------- // +/* + * Prolongation Operator + * --------------------- + * + * We use an anisotropic bilinear interpolation stencil. + * For an isotropic (uniform) mesh, this stencil reduces to: + * + * |1 2 1| + * P = 1/4 |2 4 2| + * |1 2 1| + * + * A fine node is classified by the parity of its (i_r, i_theta) indices: + * + * 1) (even, even) + * The fine node coincides with a coarse node: + * -> Copy value + * + * 2) (odd, even) + * Node lies between two coarse nodes in radial direction: + * + * X ---- O ---- X + * + * -> 1D radial interpolation with anisotropic weights + * + * 3) (even, odd) + * Node lies between two coarse nodes in angular direction: + * + * X + * | + * O + * | + * X + * + * -> 1D angular interpolation with anisotropic weights + * + * 4) (odd, odd) + * Node lies inside a coarse cell, surrounded by 4 coarse nodes: + * + * X X + * \ / + * O + * / \ + * X X + * + * -> full anisotropic bilinear interpolation + * + * All weights are defined via local mesh spacings: + * - h1, h2 in radial direction + * - k1, k2 in angular direction + */ #define FINE_NODE_PROLONGATION() \ do { \ if (i_r & 1) { \ if (i_theta & 1) { \ - /* i_r % 2 == 1, i_theta % 2 == 1 */ \ - /* Fine node in the center of four coarse nodes */ \ - double h1 = fine_grid.radialSpacing(i_r - 1); \ - double h2 = fine_grid.radialSpacing(i_r); \ - double k1 = fine_grid.angularSpacing(i_theta - 1); \ - double k2 = fine_grid.angularSpacing(i_theta); \ - int i_theta_coarse_P1 = coarse_grid.wrapThetaIndex(i_theta_coarse + 1); \ - double divisor = (h1 + h2) * (k1 + k2); \ + /* (odd, odd) -> fine node in center of coarse cell */ \ + double h1 = fine_grid.radialSpacing(i_r - 1); \ + double h2 = fine_grid.radialSpacing(i_r); \ + double k1 = fine_grid.angularSpacing(i_theta - 1); \ + double k2 = fine_grid.angularSpacing(i_theta); \ + \ double value = \ (h1 * k1 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* Bottom left */ \ h2 * k1 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse)] + /* Bottom right */ \ - h1 * k2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_P1)] + /* Top left */ \ - h2 * k2 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse_P1)] /* Top right */ \ - ); \ - fine_result[fine_grid.index(i_r, i_theta)] = value / divisor; \ + h1 * k2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 1)] + /* Top left */ \ + h2 * k2 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse + 1)] /* Top right */ \ + ) / \ + ((h1 + h2) * (k1 + k2)); \ + \ + fine_result[fine_grid.index(i_r, i_theta)] = value; \ } \ else { \ - /* i_r % 2 == 1, i_theta % 2 == 0 */ \ - /* Fine node between coarse nodes in radial direction */ \ - double h1 = fine_grid.radialSpacing(i_r - 1); \ - double h2 = fine_grid.radialSpacing(i_r); \ - double divisor = (h1 + h2); \ - double value = (h1 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* left */ \ - h2 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse)] /* right */ \ - ); \ - fine_result[fine_grid.index(i_r, i_theta)] = value / divisor; \ + /* (odd, even) -> between coarse nodes in radial direction */ \ + double h1 = fine_grid.radialSpacing(i_r - 1); \ + double h2 = fine_grid.radialSpacing(i_r); \ + \ + double value = (h1 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* Left */ \ + h2 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse)] /* Right */ \ + ) / \ + (h1 + h2); \ + \ + fine_result[fine_grid.index(i_r, i_theta)] = value; \ } \ } \ else { \ if (i_theta & 1) { \ - /* i_r % 2 == 0, i_theta % 2 == 1 */ \ - /* Fine node between coarse nodes in theta direction */ \ - double k1 = fine_grid.angularSpacing(i_theta - 1); \ - double k2 = fine_grid.angularSpacing(i_theta); \ - int i_theta_coarse_P1 = coarse_grid.wrapThetaIndex(i_theta_coarse + 1); \ - double divisor = (k1 + k2); \ - double value = (k1 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* bottom */ \ - k2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_P1)] /* top */ \ - ); \ - fine_result[fine_grid.index(i_r, i_theta)] = value / divisor; \ + /* (even, odd) -> between coarse nodes in angular direction */ \ + double k1 = fine_grid.angularSpacing(i_theta - 1); \ + double k2 = fine_grid.angularSpacing(i_theta); \ + \ + double value = (k1 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* Bottom */ \ + k2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 1)] /* Top */ \ + ) / \ + (k1 + k2); \ + \ + fine_result[fine_grid.index(i_r, i_theta)] = value; \ } \ else { \ - /* i_r % 2 == 0, i_theta % 2 == 0 */ \ - /* Fine node appears in coarse grid */ \ + /* (even, even) -> node lies on coarse grid */ \ fine_result[fine_grid.index(i_r, i_theta)] = \ - coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)]; /* center */ \ + coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)]; /* Center */ \ } \ } \ } while (0) @@ -155,10 +111,12 @@ void Interpolation::applyProlongation(const PolarGrid& coarse_grid, const PolarG 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_) { -/* Circluar Indexing Section */ -/* For loop matches circular access pattern */ + /* 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; @@ -168,8 +126,7 @@ void Interpolation::applyProlongation(const PolarGrid& coarse_grid, const PolarG } } -/* Radial Indexing Section */ -/* For loop matches radial access pattern */ + /* 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; diff --git a/src/Interpolation/restriction.cpp b/src/Interpolation/restriction.cpp index 24d9ff9d..ca99cace 100644 --- a/src/Interpolation/restriction.cpp +++ b/src/Interpolation/restriction.cpp @@ -1,99 +1,68 @@ #include "../../include/Interpolation/interpolation.h" -// For the restriction we use R = P^T. -// The restriction for an siotropic mesh reduces to -// |1 2 1| -// R = 1/4 * |2 4 2| = P^T -// |1 2 1| - -void Interpolation::applyRestriction0(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())); - -#pragma omp parallel for num_threads(max_omp_threads_) - for (int index = 0; index < coarse_grid.numberOfNodes(); index++) { - MultiIndex coarse_node = coarse_grid.multiIndex(index); - MultiIndex fine_node(2 * coarse_node[0], 2 * coarse_node[1]); - - std::array, space_dimension> neighbor_distance; - std::array, space_dimension> neighbors; - - // Center - double value = fine_values[fine_grid.index(fine_node)]; - - fine_grid.adjacentNeighborsOf(fine_node, neighbors); - - // Left - if (neighbors[0].first != -1) { - fine_grid.adjacentNeighborDistances(fine_grid.multiIndex(neighbors[0].first), neighbor_distance); - value += neighbor_distance[0].second * fine_values[neighbors[0].first] / - (neighbor_distance[0].first + neighbor_distance[0].second); - } - - // Right - if (neighbors[0].second != -1) { - fine_grid.adjacentNeighborDistances(fine_grid.multiIndex(neighbors[0].second), neighbor_distance); - value += neighbor_distance[0].first * fine_values[neighbors[0].second] / - (neighbor_distance[0].first + neighbor_distance[0].second); - } - - // Bottom - if (neighbors[1].first != -1) { - fine_grid.adjacentNeighborDistances(fine_grid.multiIndex(neighbors[1].first), neighbor_distance); - value += neighbor_distance[1].second * fine_values[neighbors[1].first] / - (neighbor_distance[1].first + neighbor_distance[1].second); - } - - // Top - if (neighbors[1].second != -1) { - fine_grid.adjacentNeighborDistances(fine_grid.multiIndex(neighbors[1].second), neighbor_distance); - value += neighbor_distance[1].first * fine_values[neighbors[1].second] / - (neighbor_distance[1].first + neighbor_distance[1].second); - } - - fine_grid.diagonalNeighborsOf(fine_node, neighbors); - - // Bottom Left - if (neighbors[0].first != -1) { - fine_grid.adjacentNeighborDistances(fine_grid.multiIndex(neighbors[0].first), neighbor_distance); - value += neighbor_distance[0].second * neighbor_distance[1].second * fine_values[neighbors[0].first] / - ((neighbor_distance[0].first + neighbor_distance[0].second) * - (neighbor_distance[1].first + neighbor_distance[1].second)); - } - - // Bottom Right - if (neighbors[0].second != -1) { - fine_grid.adjacentNeighborDistances(fine_grid.multiIndex(neighbors[0].second), neighbor_distance); - value += neighbor_distance[0].first * neighbor_distance[1].second * fine_values[neighbors[0].second] / - ((neighbor_distance[0].first + neighbor_distance[0].second) * - (neighbor_distance[1].first + neighbor_distance[1].second)); - } - - // Top Left - if (neighbors[1].first != -1) { - fine_grid.adjacentNeighborDistances(fine_grid.multiIndex(neighbors[1].first), neighbor_distance); - value += neighbor_distance[0].second * neighbor_distance[1].first * fine_values[neighbors[1].first] / - ((neighbor_distance[0].first + neighbor_distance[0].second) * - (neighbor_distance[1].first + neighbor_distance[1].second)); - } - - // Top Right - if (neighbors[1].second != -1) { - fine_grid.adjacentNeighborDistances(fine_grid.multiIndex(neighbors[1].second), neighbor_distance); - value += neighbor_distance[0].first * neighbor_distance[1].first * fine_values[neighbors[1].second] / - ((neighbor_distance[0].first + neighbor_distance[0].second) * - (neighbor_distance[1].first + neighbor_distance[1].second)); - } - - coarse_result[index] = value; - } -} - -// -------------------------------------- // -// Optimized version of applyRestriction0 // -// -------------------------------------- // +/* + * Restriction Operator + * -------------------- + * + * We use the transpose of the anisotropic bilinear interpolation stencil: R = P^T + * For an isotropic (uniform) mesh, this stencil reduces to: + * + * |1 2 1| + * R = 1/4 |2 4 2| + * |1 2 1| + * + * Each coarse node accumulates contributions from its corresponding 9 surrounding fine nodes. + * + * Weights are determined by the anisotropic mesh spacings: + * - h1, h2, h3, h4 in radial direction + * - k1, k2, k3, k4 in angular direction + * + * Boundary handling: + * - Angular direction: periodic (wrapping) + * - Radial direction: check domain boundaries (inner/outer radius) + */ + +#define COARSE_NODE_RESTRICTION() \ + do { \ + int i_r = i_r_coarse * 2; \ + int i_theta = i_theta_coarse * 2; \ + \ + /* Angular indices with periodic wrapping */ \ + int i_theta_M2 = fine_grid.wrapThetaIndex(i_theta - 2); \ + int i_theta_M1 = fine_grid.wrapThetaIndex(i_theta - 1); \ + int i_theta_P1 = fine_grid.wrapThetaIndex(i_theta + 1); \ + \ + /* Angular spacings */ \ + double k1 = fine_grid.angularSpacing(i_theta_M2); \ + double k2 = fine_grid.angularSpacing(i_theta_M1); \ + double k3 = fine_grid.angularSpacing(i_theta); \ + double k4 = fine_grid.angularSpacing(i_theta_P1); \ + \ + /* Center + Angular contributions (always present) */ \ + double value = fine_values[fine_grid.index(i_r, i_theta)] + \ + k2 / (k1 + k2) * fine_values[fine_grid.index(i_r, i_theta_M1)] + \ + k3 / (k3 + k4) * fine_values[fine_grid.index(i_r, i_theta_P1)]; \ + \ + /* Left contributions (if not at inner boundary) */ \ + if (i_r_coarse > 0) { \ + double h1 = fine_grid.radialSpacing(i_r - 2); \ + double h2 = fine_grid.radialSpacing(i_r - 1); \ + value += h2 / (h1 + h2) * fine_values[fine_grid.index(i_r - 1, i_theta)] + \ + h2 * k2 / ((h1 + h2) * (k1 + k2)) * fine_values[fine_grid.index(i_r - 1, i_theta_M1)] + \ + h2 * k3 / ((h1 + h2) * (k3 + k4)) * fine_values[fine_grid.index(i_r - 1, i_theta_P1)]; \ + } \ + \ + /* Right contributions (if not at outer boundary) */ \ + if (i_r_coarse < coarse_grid.nr() - 1) { \ + double h3 = fine_grid.radialSpacing(i_r); \ + double h4 = fine_grid.radialSpacing(i_r + 1); \ + value += h3 / (h3 + h4) * fine_values[fine_grid.index(i_r + 1, i_theta)] + \ + h3 * k2 / ((h3 + h4) * (k1 + k2)) * fine_values[fine_grid.index(i_r + 1, i_theta_M1)] + \ + h3 * k3 / ((h3 + h4) * (k3 + k4)) * fine_values[fine_grid.index(i_r + 1, i_theta_P1)]; \ + } \ + \ + coarse_result[coarse_grid.index(i_r_coarse, i_theta_coarse)] = value; \ + } while (0) void Interpolation::applyRestriction(const PolarGrid& fine_grid, const PolarGrid& coarse_grid, Vector coarse_result, ConstVector fine_values) const @@ -101,155 +70,25 @@ void Interpolation::applyRestriction(const PolarGrid& fine_grid, const PolarGrid assert(fine_values.size() == static_cast(fine_grid.numberOfNodes())); assert(coarse_result.size() == static_cast(coarse_grid.numberOfNodes())); - const int coarseNumberSmootherCircles = coarse_grid.numberSmootherCircles(); + /* 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_) { -/* For loop matches circular access pattern */ + /* Circular Indexing Section */ #pragma omp for nowait - for (int i_r_coarse = 0; i_r_coarse < coarseNumberSmootherCircles; i_r_coarse++) { - int i_r = i_r_coarse * 2; + 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++) { - int i_theta = i_theta_coarse * 2; - - if (0 < i_r_coarse && i_r_coarse < coarseNumberSmootherCircles - 1) { - double h1 = fine_grid.radialSpacing(i_r - 2); - double h2 = fine_grid.radialSpacing(i_r - 1); - double h3 = fine_grid.radialSpacing(i_r); - double h4 = fine_grid.radialSpacing(i_r + 1); - - int i_theta_M2 = fine_grid.wrapThetaIndex(i_theta - 2); - int i_theta_M1 = fine_grid.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = fine_grid.wrapThetaIndex(i_theta + 1); - - double k1 = fine_grid.angularSpacing(i_theta_M2); - double k2 = fine_grid.angularSpacing(i_theta_M1); - double k3 = fine_grid.angularSpacing(i_theta); - double k4 = fine_grid.angularSpacing(i_theta_P1); - - coarse_result[coarse_grid.index(i_r_coarse, i_theta_coarse)] = - // Center - fine_values[fine_grid.index(i_r, i_theta)] + - // Left, Right, Bottom, Top - h2 * fine_values[fine_grid.index(i_r - 1, i_theta)] / (h1 + h2) + - h3 * fine_values[fine_grid.index(i_r + 1, i_theta)] / (h3 + h4) + - k2 * fine_values[fine_grid.index(i_r, i_theta_M1)] / (k1 + k2) + - k3 * fine_values[fine_grid.index(i_r, i_theta_P1)] / (k3 + k4) + - // Bottom Left, Bottom Right, Top Left, Top Right - h2 * k2 * fine_values[fine_grid.index(i_r - 1, i_theta_M1)] / ((h1 + h2) * (k1 + k2)) + - h3 * k2 * fine_values[fine_grid.index(i_r + 1, i_theta_M1)] / ((h3 + h4) * (k1 + k2)) + - h2 * k3 * fine_values[fine_grid.index(i_r - 1, i_theta_P1)] / ((h1 + h2) * (k3 + k4)) + - h3 * k3 * fine_values[fine_grid.index(i_r + 1, i_theta_P1)] / ((h3 + h4) * (k3 + k4)); - } - else { - /* First and Last Circle have to be checked for domain boundary */ - // Middle Part - int i_theta_M2 = fine_grid.wrapThetaIndex(i_theta - 2); - int i_theta_M1 = fine_grid.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = fine_grid.wrapThetaIndex(i_theta + 1); - double k1 = fine_grid.angularSpacing(i_theta_M2); - double k2 = fine_grid.angularSpacing(i_theta_M1); - double k3 = fine_grid.angularSpacing(i_theta); - double k4 = fine_grid.angularSpacing(i_theta_P1); - // Center, Bottom, Top - double value = fine_values[fine_grid.index(i_r, i_theta)] + - k2 * fine_values[fine_grid.index(i_r, i_theta_M1)] / (k1 + k2) + - k3 * fine_values[fine_grid.index(i_r, i_theta_P1)] / (k3 + k4); - - if (i_r_coarse > 0) { - // Left Part - double h1 = fine_grid.radialSpacing(i_r - 2); - double h2 = fine_grid.radialSpacing(i_r - 1); - // Left, Bottom Left, Top Left - value += h2 * fine_values[fine_grid.index(i_r - 1, i_theta)] / (h1 + h2) + - h2 * k2 * fine_values[fine_grid.index(i_r - 1, i_theta_M1)] / ((h1 + h2) * (k1 + k2)) + - h2 * k3 * fine_values[fine_grid.index(i_r - 1, i_theta_P1)] / ((h1 + h2) * (k3 + k4)); - } - if (i_r_coarse < coarse_grid.nr() - 1) { - // Right Part - double h3 = fine_grid.radialSpacing(i_r); - double h4 = fine_grid.radialSpacing(i_r + 1); - // Right, Bottom Right, Top Right - value += h3 * fine_values[fine_grid.index(i_r + 1, i_theta)] / (h3 + h4) + - h3 * k2 * fine_values[fine_grid.index(i_r + 1, i_theta_M1)] / ((h3 + h4) * (k1 + k2)) + - h3 * k3 * fine_values[fine_grid.index(i_r + 1, i_theta_P1)] / ((h3 + h4) * (k3 + k4)); - } - coarse_result[coarse_grid.index(i_r_coarse, i_theta_coarse)] = value; - } + COARSE_NODE_RESTRICTION(); } } -/* For loop matches radial access pattern */ + /* Radial Indexing Section */ #pragma omp for nowait for (int i_theta_coarse = 0; i_theta_coarse < coarse_grid.ntheta(); i_theta_coarse++) { - int i_theta = i_theta_coarse * 2; - for (int i_r_coarse = coarseNumberSmootherCircles; i_r_coarse < coarse_grid.nr(); i_r_coarse++) { - int i_r = i_r_coarse * 2; - - if (coarse_grid.numberSmootherCircles() < i_r_coarse && i_r_coarse < coarse_grid.nr() - 1) { - double h1 = fine_grid.radialSpacing(i_r - 2); - double h2 = fine_grid.radialSpacing(i_r - 1); - double h3 = fine_grid.radialSpacing(i_r); - double h4 = fine_grid.radialSpacing(i_r + 1); - - int i_theta_M2 = fine_grid.wrapThetaIndex(i_theta - 2); - int i_theta_M1 = fine_grid.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = fine_grid.wrapThetaIndex(i_theta + 1); - double k1 = fine_grid.angularSpacing(i_theta_M2); - double k2 = fine_grid.angularSpacing(i_theta_M1); - double k3 = fine_grid.angularSpacing(i_theta); - double k4 = fine_grid.angularSpacing(i_theta_P1); - - coarse_result[coarse_grid.index(i_r_coarse, i_theta_coarse)] = - // Center - fine_values[fine_grid.index(i_r, i_theta)] + - // Left, Right, Bottom, Top - h2 * fine_values[fine_grid.index(i_r - 1, i_theta)] / (h1 + h2) + - h3 * fine_values[fine_grid.index(i_r + 1, i_theta)] / (h3 + h4) + - k2 * fine_values[fine_grid.index(i_r, i_theta_M1)] / (k1 + k2) + - k3 * fine_values[fine_grid.index(i_r, i_theta_P1)] / (k3 + k4) + - // Bottom Left, Bottom Right, Top Left, Top Right - h2 * k2 * fine_values[fine_grid.index(i_r - 1, i_theta_M1)] / ((h1 + h2) * (k1 + k2)) + - h3 * k2 * fine_values[fine_grid.index(i_r + 1, i_theta_M1)] / ((h3 + h4) * (k1 + k2)) + - h2 * k3 * fine_values[fine_grid.index(i_r - 1, i_theta_P1)] / ((h1 + h2) * (k3 + k4)) + - h3 * k3 * fine_values[fine_grid.index(i_r + 1, i_theta_P1)] / ((h3 + h4) * (k3 + k4)); - } - else { - /* First and Last radial nodes have to be checked for domain boundary */ - // Middle Part - int i_theta_M2 = fine_grid.wrapThetaIndex(i_theta - 2); - int i_theta_M1 = fine_grid.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = fine_grid.wrapThetaIndex(i_theta + 1); - double k1 = fine_grid.angularSpacing(i_theta_M2); - double k2 = fine_grid.angularSpacing(i_theta_M1); - double k3 = fine_grid.angularSpacing(i_theta); - double k4 = fine_grid.angularSpacing(i_theta_P1); - // Center, Bottom, Top - double value = fine_values[fine_grid.index(i_r, i_theta)] + - k2 * fine_values[fine_grid.index(i_r, i_theta_M1)] / (k1 + k2) + - k3 * fine_values[fine_grid.index(i_r, i_theta_P1)] / (k3 + k4); - if (i_r_coarse > 0) { - // Left Part - double h1 = fine_grid.radialSpacing(i_r - 2); - double h2 = fine_grid.radialSpacing(i_r - 1); - // Left, Bottom Left, Top Left - value += h2 * fine_values[fine_grid.index(i_r - 1, i_theta)] / (h1 + h2) + - h2 * k2 * fine_values[fine_grid.index(i_r - 1, i_theta_M1)] / ((h1 + h2) * (k1 + k2)) + - h2 * k3 * fine_values[fine_grid.index(i_r - 1, i_theta_P1)] / ((h1 + h2) * (k3 + k4)); - } - if (i_r_coarse < coarse_grid.nr() - 1) { - // Right Part - double h3 = fine_grid.radialSpacing(i_r); - double h4 = fine_grid.radialSpacing(i_r + 1); - // Right, Bottom Right, Top Right - value += h3 * fine_values[fine_grid.index(i_r + 1, i_theta)] / (h3 + h4) + - h3 * k2 * fine_values[fine_grid.index(i_r + 1, i_theta_M1)] / ((h3 + h4) * (k1 + k2)) + - h3 * k3 * fine_values[fine_grid.index(i_r + 1, i_theta_P1)] / ((h3 + h4) * (k3 + k4)); - } - coarse_result[coarse_grid.index(i_r_coarse, i_theta_coarse)] = value; - } + for (int i_r_coarse = coarse_grid.numberSmootherCircles(); i_r_coarse < coarse_grid.nr(); i_r_coarse++) { + COARSE_NODE_RESTRICTION(); } } } } -//clang-format on diff --git a/tests/Interpolation/extrapolated_prolongation.cpp b/tests/Interpolation/extrapolated_prolongation.cpp index 0e749047..07225b8f 100644 --- a/tests/Interpolation/extrapolated_prolongation.cpp +++ b/tests/Interpolation/extrapolated_prolongation.cpp @@ -1,58 +1,74 @@ #include - #include -#include "../../include/GMGPolar/gmgpolar.h" #include "../../include/Interpolation/interpolation.h" -#include "../../include/InputFunctions/domainGeometry.h" -#include "../../include/InputFunctions/densityProfileCoefficients.h" - -#include "../../include/InputFunctions/DomainGeometry/circularGeometry.h" #include "../../include/InputFunctions/DensityProfileCoefficients/poissonCoefficients.h" -namespace ExtrapolatedProlongationTest +static Vector generate_random_sample_data(const PolarGrid& grid, unsigned int seed) { -// Function to generate sample data for vector x using random values with seed -Vector generate_random_sample_data(const PolarGrid& grid, unsigned int seed) + Vector vec("vec", grid.numberOfNodes()); + std::mt19937 gen(seed); + std::uniform_real_distribution dist(-20.0, 20.0); + + for (uint i = 0; i < vec.size(); ++i) + vec[i] = dist(gen); + + return vec; +} + +// 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) { - Vector vector("vector", grid.numberOfNodes()); - std::mt19937 gen(seed); // Standard mersenne_twister_engine seeded with seed - std::uniform_real_distribution dist(0.0, 1.0); // Generate random double between 0 and 1 - for (uint i = 0; i < vector.size(); ++i) { - vector[i] = dist(gen); + 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)]); } - return vector; -} -} // namespace ExtrapolatedProlongationTest -using namespace ExtrapolatedProlongationTest; + // 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, ExtrapolatedProlongationSmoothingRadius) +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, M_PI + M_PI}; - - int maxOpenMPThreads = 16; - bool DirBC_Interior = true; - - PolarGrid finest_grid(fine_radii, fine_angles); - PolarGrid coarse_grid = coarseningGrid(finest_grid); + 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}; - Interpolation interpolation_operator(maxOpenMPThreads, DirBC_Interior); + PolarGrid fine_grid(fine_radii, fine_angles); + PolarGrid coarse_grid = coarseningGrid(fine_grid); - unsigned int seed = 42; - Vector x = generate_random_sample_data(coarse_grid, seed); + Interpolation I(/*threads*/ 16, /*DirBC*/ true); - // Apply prolongation to both functions - Vector result1("result1", finest_grid.numberOfNodes()); - Vector result2("result2", finest_grid.numberOfNodes()); + Vector coarse_values = generate_random_sample_data(coarse_grid, 1234); + Vector fine_result("fine_result", fine_grid.numberOfNodes()); - interpolation_operator.applyExtrapolatedProlongation0(coarse_grid, finest_grid, result1, x); - interpolation_operator.applyExtrapolatedProlongation(coarse_grid, finest_grid, result2, x); + I.applyExtrapolatedProlongation(coarse_grid, fine_grid, fine_result, coarse_values); - ASSERT_EQ(result1.size(), result2.size()); - for (uint i = 0; i < result1.size(); ++i) { - ASSERT_DOUBLE_EQ(result1[i], result2[i]); + 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 index 6d619a7b..552c83a1 100644 --- a/tests/Interpolation/extrapolated_restriction.cpp +++ b/tests/Interpolation/extrapolated_restriction.cpp @@ -1,140 +1,73 @@ #include - #include -#include "../../include/GMGPolar/gmgpolar.h" #include "../../include/Interpolation/interpolation.h" -#include "../../include/InputFunctions/domainGeometry.h" -#include "../../include/InputFunctions/densityProfileCoefficients.h" - -#include "../../include/InputFunctions/DomainGeometry/circularGeometry.h" #include "../../include/InputFunctions/DensityProfileCoefficients/poissonCoefficients.h" -namespace ExtrapolatedRestrictionTest -{ -// Function to generate sample data for vector x using random values with seed -Vector generate_random_sample_data(const PolarGrid& grid, unsigned int seed) -{ - Vector vector("vector", grid.numberOfNodes()); - std::mt19937 gen(seed); // Standard mersenne_twister_engine seeded with seed - std::uniform_real_distribution dist(0.0, 1.0); // Generate random double between 0 and 1 - for (uint i = 0; i < vector.size(); ++i) { - vector[i] = dist(gen); - } - return vector; -} - -/* In src/Interpolation/restriction.cpp the Restriction Operator is implemented with "Take". */ -/* Here we test against the "Give" version. */ - -void applyExtrapolatedRestrictionGive0(const PolarGrid& fine_grid, const PolarGrid& coarse_grid, - Vector coarse_result, ConstVector fine_values) +static Vector generate_random_sample_data(const PolarGrid& grid, unsigned int seed) { - assert(fine_values.size() == static_cast(fine_grid.numberOfNodes())); - assert(coarse_result.size() == static_cast(coarse_grid.numberOfNodes())); - - assign(coarse_result, 0.0); - - for (int index = 0; index < fine_grid.numberOfNodes(); index++) { - std::array, space_dimension> neighbor_distance; - std::array, space_dimension> neighbors; - - MultiIndex fine_node = fine_grid.multiIndex(index); - - // Fine node appears in coarse grid - if (fine_node[0] % 2 == 0 && fine_node[1] % 2 == 0) { - // Input x needs a fine grid index: fine_values[FINE_INDEX] - // Result needs a coarse grid index: coarse_result[COARSE_INDEX] - MultiIndex coarse_node(fine_node[0] / 2, fine_node[1] / 2); - coarse_result[coarse_grid.index(coarse_node)] += fine_values[index]; - } - - // Fine node between coarse nodes in theta direction - if (fine_node[0] % 2 == 0 && fine_node[1] % 2 == 1) { - fine_grid.adjacentNeighborDistances(fine_node, neighbor_distance); - double k1 = neighbor_distance[1].first; - double k2 = neighbor_distance[1].second; - - fine_grid.adjacentNeighborsOf(fine_node, neighbors); - - MultiIndex bottom_coarse_node(fine_node[0] / 2, fine_node[1] / 2); - MultiIndex top_coarse_node(fine_node[0] / 2, (fine_node[1] / 2 + 1) % coarse_grid.ntheta()); + Vector vec("vec", grid.numberOfNodes()); + std::mt19937 gen(seed); + std::uniform_real_distribution dist(-20.0, 20.0); - coarse_result[coarse_grid.index(bottom_coarse_node)] += fine_values[index] / 2.0; - coarse_result[coarse_grid.index(top_coarse_node)] += fine_values[index] / 2.0; - } - - // Fine node between coarse nodes in radial direction - if (fine_node[0] % 2 == 1 && fine_node[1] % 2 == 0) { - fine_grid.adjacentNeighborDistances(fine_node, neighbor_distance); - double h1 = neighbor_distance[0].first; - double h2 = neighbor_distance[0].second; - - fine_grid.adjacentNeighborsOf(fine_node, neighbors); - - MultiIndex left_coarse_node(fine_node[0] / 2, fine_node[1] / 2); - MultiIndex right_coarse_node(fine_node[0] / 2 + 1, fine_node[1] / 2); + for (uint i = 0; i < vec.size(); ++i) + vec[i] = dist(gen); - coarse_result[coarse_grid.index(left_coarse_node)] += fine_values[index] / 2.0; - coarse_result[coarse_grid.index(right_coarse_node)] += fine_values[index] / 2.0; - } + return vec; +} - //Fine node in the center of four coarse nodes - if (fine_node[0] % 2 == 1 && fine_node[1] % 2 == 1) { +// 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; - fine_grid.adjacentNeighborDistances(fine_node, neighbor_distance); - double h1 = neighbor_distance[0].first; - double h2 = neighbor_distance[0].second; - double k1 = neighbor_distance[1].first; - double k2 = neighbor_distance[1].second; + // Angular indices with periodic wrapping + int i_theta_M1 = fine.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = fine.wrapThetaIndex(i_theta + 1); - fine_grid.adjacentNeighborsOf(fine_node, neighbors); + // 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)]; - MultiIndex bottom_right_coarse_node(fine_node[0] / 2 + 1, fine_node[1] / 2); - MultiIndex top_left_coarse_node(fine_node[0] / 2, (fine_node[1] / 2 + 1) % coarse_grid.ntheta()); + // 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 + } - coarse_result[coarse_grid.index(bottom_right_coarse_node)] += fine_values[index] / 2.0; - coarse_result[coarse_grid.index(top_left_coarse_node)] += fine_values[index] / 2.0; - } + // 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 } -} -} // namespace ExtrapolatedRestrictionTest -using namespace ExtrapolatedRestrictionTest; + return value; +} -TEST(ExtrapolatedRestrictionTest, applyExtrapolatedRestriction) +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, M_PI + M_PI}; - - int maxOpenMPThreads = 16; - bool DirBC_Interior = true; + 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 finest_grid(fine_radii, fine_angles); - PolarGrid coarse_grid = coarseningGrid(finest_grid); + PolarGrid fine_grid(fine_radii, fine_angles); + PolarGrid coarse_grid = coarseningGrid(fine_grid); - Interpolation interpolation_operator(maxOpenMPThreads, DirBC_Interior); + Interpolation I(/*threads*/ 16, /*DirBC*/ true); - unsigned int seed = 42; - Vector x = generate_random_sample_data(finest_grid, seed); + Vector fine_values = generate_random_sample_data(fine_grid, 9012); + Vector coarse_result("coarse_result", coarse_grid.numberOfNodes()); - // Apply prolongation to both functions - Vector result1("result1", coarse_grid.numberOfNodes()); - Vector result2("result2", coarse_grid.numberOfNodes()); - Vector result3("result3", coarse_grid.numberOfNodes()); + I.applyExtrapolatedRestriction(fine_grid, coarse_grid, coarse_result, fine_values); - interpolation_operator.applyExtrapolatedRestriction0(finest_grid, coarse_grid, result1, x); - interpolation_operator.applyExtrapolatedRestriction(finest_grid, coarse_grid, result2, x); - - applyExtrapolatedRestrictionGive0(finest_grid, coarse_grid, result3, x); - - ASSERT_EQ(result1.size(), result2.size()); - for (uint i = 0; i < result1.size(); ++i) { - ASSERT_DOUBLE_EQ(result1[i], result2[i]); - } - ASSERT_EQ(result2.size(), result3.size()); - for (uint i = 0; i < result2.size(); ++i) { - ASSERT_DOUBLE_EQ(result2[i], result3[i]); + 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 << ")"; + } } } diff --git a/tests/Interpolation/prolongation.cpp b/tests/Interpolation/prolongation.cpp index 01655b11..207634a4 100644 --- a/tests/Interpolation/prolongation.cpp +++ b/tests/Interpolation/prolongation.cpp @@ -1,52 +1,90 @@ #include - #include -#include "../../include/GMGPolar/gmgpolar.h" #include "../../include/Interpolation/interpolation.h" #include "../../include/InputFunctions/DensityProfileCoefficients/poissonCoefficients.h" -namespace ProlongationTest -{ -Vector generate_random_sample_data(const PolarGrid& grid, unsigned int seed) +static Vector generate_random_sample_data(const PolarGrid& grid, unsigned int seed) { - Vector vector("vector", grid.numberOfNodes()); + Vector vec("vec", grid.numberOfNodes()); std::mt19937 gen(seed); - std::uniform_real_distribution dist(-100.0, 100.0); - for (uint i = 0; i < vector.size(); ++i) { - vector[i] = dist(gen); - } - return vector; + std::uniform_real_distribution dist(-20.0, 20.0); + + for (uint i = 0; i < vec.size(); ++i) + vec[i] = dist(gen); + + return vec; } -} // namespace ProlongationTest -using namespace ProlongationTest; +// Helper that computes the mathematically expected prolongation value +static double expected_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 + double h1 = fine.radialSpacing(i_r - 1); + double h2 = fine.radialSpacing(i_r); -TEST(ProlongationTest, ProlongationTest) + return (h1 * coarse_vals[coarse.index(i_r_coarse, i_theta_coarse)] + + h2 * coarse_vals[coarse.index(i_r_coarse + 1, i_theta_coarse)]) / + (h1 + h2); + } + + if (r_even && !t_even) { + // Angular midpoint + double k1 = fine.angularSpacing(i_theta - 1); + double k2 = fine.angularSpacing(i_theta); + + return (k1 * coarse_vals[coarse.index(i_r_coarse, t_even ? i_theta_coarse : i_theta_coarse)] + + k2 * coarse_vals[coarse.index(i_r_coarse, i_theta_coarse + 1)]) / + (k1 + k2); + } + + // Center of coarse cell + double h1 = fine.radialSpacing(i_r - 1); + double h2 = fine.radialSpacing(i_r); + double k1 = fine.angularSpacing(i_theta - 1); + double k2 = fine.angularSpacing(i_theta); + + return (h1 * k1 * coarse_vals[coarse.index(i_r_coarse, i_theta_coarse)] + + h2 * k1 * coarse_vals[coarse.index(i_r_coarse + 1, i_theta_coarse)] + + h1 * k2 * coarse_vals[coarse.index(i_r_coarse, i_theta_coarse + 1)] + + h2 * k2 * coarse_vals[coarse.index(i_r_coarse + 1, i_theta_coarse + 1)]) / + ((h1 + h2) * (k1 + k2)); +} + +TEST(ProlongationTest, ProlongationMatchesStencil) { 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, M_PI + M_PI}; - - int maxOpenMPThreads = 16; - bool DirBC_Interior = true; - - PolarGrid finest_grid(fine_radii, fine_angles); - PolarGrid coarse_grid = coarseningGrid(finest_grid); + 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}; - Interpolation interpolation_operator(maxOpenMPThreads, DirBC_Interior); + PolarGrid fine_grid(fine_radii, fine_angles); + PolarGrid coarse_grid = coarseningGrid(fine_grid); - Vector x = generate_random_sample_data(coarse_grid, 42); + Interpolation I(/*threads*/ 16, /*DirBC*/ true); - // Apply prolongation to both functions - Vector result1("result1", finest_grid.numberOfNodes()); - Vector result2("result2", finest_grid.numberOfNodes()); + Vector coarse_values = generate_random_sample_data(coarse_grid, 1234); + Vector fine_result("fine_result", fine_grid.numberOfNodes()); - interpolation_operator.applyProlongation0(coarse_grid, finest_grid, result1, x); - interpolation_operator.applyProlongation(coarse_grid, finest_grid, result2, x); + I.applyProlongation(coarse_grid, fine_grid, fine_result, coarse_values); - ASSERT_EQ(result1.size(), result2.size()); - for (uint i = 0; i < result1.size(); ++i) { - ASSERT_NEAR(result1[i], result2[i], 1e-10); + 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)]; + 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 006f8c75..6445af78 100644 --- a/tests/Interpolation/restriction.cpp +++ b/tests/Interpolation/restriction.cpp @@ -1,142 +1,86 @@ #include - #include -#include "../../include/GMGPolar/gmgpolar.h" #include "../../include/Interpolation/interpolation.h" #include "../../include/InputFunctions/DensityProfileCoefficients/poissonCoefficients.h" -namespace RestrictionTest -{ -// Function to generate sample data for vector x using random values with seed -Vector generate_random_sample_data(const PolarGrid& grid, unsigned int seed) -{ - Vector vector("vector", grid.numberOfNodes()); - std::mt19937 gen(seed); // Standard mersenne_twister_engine seeded with seed - std::uniform_real_distribution dist(0.0, 1.0); // Generate random double between 0 and 1 - for (uint i = 0; i < vector.size(); ++i) { - vector[i] = dist(gen); - } - return vector; -} - -/* In src/Interpolation/restriction.cpp the Restriction Operator is implemented with "Take". */ -/* Here we test against the "Give" version. */ -void applyRestrictionGive0(const PolarGrid& fine_grid, const PolarGrid& coarse_grid, Vector coarse_result, - const Vector fine_values) +static Vector generate_random_sample_data(const PolarGrid& grid, unsigned int seed) { - assert(fine_values.size() == static_cast(fine_grid.numberOfNodes())); - assert(coarse_result.size() == static_cast(coarse_grid.numberOfNodes())); - - assign(coarse_result, 0.0); - - for (int index = 0; index < fine_grid.numberOfNodes(); index++) { - std::array, space_dimension> neighbor_distance; - std::array, space_dimension> neighbors; - - MultiIndex fine_node = fine_grid.multiIndex(index); - - // Fine node appears in coarse grid - if (fine_node[0] % 2 == 0 && fine_node[1] % 2 == 0) { - // Input x needs a fine grid index: fine_values[FINE_INDEX] - // Result needs a coarse grid index: coarse_result[COARSE_INDEX] - MultiIndex coarse_node(fine_node[0] / 2, fine_node[1] / 2); - coarse_result[coarse_grid.index(coarse_node)] += fine_values[index]; - } - - // Fine node between coarse nodes in theta direction - if (fine_node[0] % 2 == 0 && fine_node[1] % 2 == 1) { - fine_grid.adjacentNeighborDistances(fine_node, neighbor_distance); - double k1 = neighbor_distance[1].first; - double k2 = neighbor_distance[1].second; - - fine_grid.adjacentNeighborsOf(fine_node, neighbors); - - MultiIndex bottom_coarse_node(fine_node[0] / 2, fine_node[1] / 2); - MultiIndex top_coarse_node(fine_node[0] / 2, (fine_node[1] / 2 + 1) % coarse_grid.ntheta()); - - coarse_result[coarse_grid.index(bottom_coarse_node)] += k1 * fine_values[index] / (k1 + k2); - coarse_result[coarse_grid.index(top_coarse_node)] += k2 * fine_values[index] / (k1 + k2); - } - - // Fine node between coarse nodes in radial direction - if (fine_node[0] % 2 == 1 && fine_node[1] % 2 == 0) { - fine_grid.adjacentNeighborDistances(fine_node, neighbor_distance); - double h1 = neighbor_distance[0].first; - double h2 = neighbor_distance[0].second; + Vector vec("vec", grid.numberOfNodes()); + std::mt19937 gen(seed); + std::uniform_real_distribution dist(-20.0, 20.0); - fine_grid.adjacentNeighborsOf(fine_node, neighbors); + for (uint i = 0; i < vec.size(); ++i) + vec[i] = dist(gen); - MultiIndex left_coarse_node(fine_node[0] / 2, fine_node[1] / 2); - MultiIndex right_coarse_node(fine_node[0] / 2 + 1, fine_node[1] / 2); + return vec; +} - coarse_result[coarse_grid.index(left_coarse_node)] += (h1 * fine_values[index]) / (h1 + h2); - coarse_result[coarse_grid.index(right_coarse_node)] += (h2 * fine_values[index]) / (h1 + h2); - } +// Helper that computes the mathematically expected restriction value +static double expected_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_M2 = fine.wrapThetaIndex(i_theta - 2); + int i_theta_M1 = fine.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = fine.wrapThetaIndex(i_theta + 1); + + // Angular spacings + double k1 = fine.angularSpacing(i_theta_M2); + double k2 = fine.angularSpacing(i_theta_M1); + double k3 = fine.angularSpacing(i_theta); + double k4 = fine.angularSpacing(i_theta_P1); + + // Center + Angular contributions (always present) + double value = fine_vals[fine.index(i_r, i_theta)] + k2 / (k1 + k2) * fine_vals[fine.index(i_r, i_theta_M1)] + + k3 / (k3 + k4) * fine_vals[fine.index(i_r, i_theta_P1)]; + + // Left contributions (if not at inner boundary) + if (i_r_coarse > 0) { + double h1 = fine.radialSpacing(i_r - 2); + double h2 = fine.radialSpacing(i_r - 1); + value += h2 / (h1 + h2) * fine_vals[fine.index(i_r - 1, i_theta)] + + h2 * k2 / ((h1 + h2) * (k1 + k2)) * fine_vals[fine.index(i_r - 1, i_theta_M1)] + + h2 * k3 / ((h1 + h2) * (k3 + k4)) * fine_vals[fine.index(i_r - 1, i_theta_P1)]; + } - //Fine node in the center of four coarse nodes - if (fine_node[0] % 2 == 1 && fine_node[1] % 2 == 1) { - - fine_grid.adjacentNeighborDistances(fine_node, neighbor_distance); - double h1 = neighbor_distance[0].first; - double h2 = neighbor_distance[0].second; - double k1 = neighbor_distance[1].first; - double k2 = neighbor_distance[1].second; - - fine_grid.adjacentNeighborsOf(fine_node, neighbors); - - MultiIndex bottom_left_coarse_node(fine_node[0] / 2, fine_node[1] / 2); - MultiIndex bottom_right_coarse_node(fine_node[0] / 2 + 1, fine_node[1] / 2); - MultiIndex top_left_coarse_node(fine_node[0] / 2, (fine_node[1] / 2 + 1) % coarse_grid.ntheta()); - MultiIndex top_right_node(fine_node[0] / 2 + 1, (fine_node[1] / 2 + 1) % coarse_grid.ntheta()); - - coarse_result[coarse_grid.index(bottom_left_coarse_node)] += - h1 * k1 * fine_values[index] / ((h1 + h2) * (k1 + k2)); - coarse_result[coarse_grid.index(bottom_right_coarse_node)] += - h2 * k1 * fine_values[index] / ((h1 + h2) * (k1 + k2)); - coarse_result[coarse_grid.index(top_left_coarse_node)] += - h1 * k2 * fine_values[index] / ((h1 + h2) * (k1 + k2)); - coarse_result[coarse_grid.index(top_right_node)] += h2 * k2 * fine_values[index] / ((h1 + h2) * (k1 + k2)); - } + // Right contributions (if not at outer boundary) + if (i_r_coarse < coarse.nr() - 1) { + double h3 = fine.radialSpacing(i_r); + double h4 = fine.radialSpacing(i_r + 1); + value += h3 / (h3 + h4) * fine_vals[fine.index(i_r + 1, i_theta)] + + h3 * k2 / ((h3 + h4) * (k1 + k2)) * fine_vals[fine.index(i_r + 1, i_theta_M1)] + + h3 * k3 / ((h3 + h4) * (k3 + k4)) * fine_vals[fine.index(i_r + 1, i_theta_P1)]; } -} -} // namespace RestrictionTest -using namespace RestrictionTest; + return value; +} -TEST(RestrictionTest, applyRestriction) +TEST(RestrictionTest, RestrictionMatchesStencil) { 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, M_PI + M_PI}; - - int maxOpenMPThreads = 16; - bool DirBC_Interior = true; - - PolarGrid finest_grid(fine_radii, fine_angles); - PolarGrid coarse_grid = coarseningGrid(finest_grid); - - Interpolation interpolation_operator(maxOpenMPThreads, DirBC_Interior); + 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}; - Vector x = generate_random_sample_data(finest_grid, 42); + PolarGrid fine_grid(fine_radii, fine_angles); + PolarGrid coarse_grid = coarseningGrid(fine_grid); - // Apply prolongation to both functions - Vector result1("result1", coarse_grid.numberOfNodes()); - Vector result2("result2", coarse_grid.numberOfNodes()); - Vector result3("result3", coarse_grid.numberOfNodes()); + Interpolation I(/*threads*/ 16, /*DirBC*/ true); - interpolation_operator.applyRestriction0(finest_grid, coarse_grid, result1, x); - interpolation_operator.applyRestriction(finest_grid, coarse_grid, result2, x); + Vector fine_values = generate_random_sample_data(fine_grid, 5678); + Vector coarse_result("coarse_result", coarse_grid.numberOfNodes()); - applyRestrictionGive0(finest_grid, coarse_grid, result3, x); + I.applyRestriction(fine_grid, coarse_grid, coarse_result, fine_values); - ASSERT_EQ(result1.size(), result2.size()); - for (uint i = 0; i < result1.size(); ++i) { - ASSERT_DOUBLE_EQ(result1[i], result2[i]); - } - - ASSERT_EQ(result2.size(), result3.size()); - for (uint i = 0; i < result2.size(); ++i) { - ASSERT_DOUBLE_EQ(result2[i], result3[i]); + 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)]; + ASSERT_NEAR(expected, got, 1e-10) << "Mismatch at (" << i_r_coarse << ", " << i_theta_coarse << ")"; + } } }