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/include/Level/level.h b/include/Level/level.h index fdfb265c..dd4ffab0 100644 --- a/include/Level/level.h +++ b/include/Level/level.h @@ -114,7 +114,6 @@ class LevelCache explicit LevelCache(const PolarGrid& grid, const DensityProfileCoefficients& density_profile_coefficients, const DomainGeometry& domain_geometry, const bool cache_density_profile_coefficients, const bool cache_domain_geometry); - explicit LevelCache(const Level& previous_level, const PolarGrid& current_grid); const DomainGeometry& domainGeometry() const; const DensityProfileCoefficients& densityProfileCoefficients() const; @@ -132,18 +131,8 @@ class LevelCache inline void obtainValues(const int i_r, const int i_theta, const int global_index, double r, double theta, double& coeff_beta, double& arr, double& att, double& art, double& detDF) const { - if (cache_density_profile_coefficients_) - coeff_beta = coeff_beta_[global_index]; - else - coeff_beta = density_profile_coefficients_.beta(r, theta); - - double coeff_alpha; - if (!cache_domain_geometry_) { - if (cache_density_profile_coefficients_) - coeff_alpha = coeff_alpha_[global_index]; - else - coeff_alpha = density_profile_coefficients_.alpha(r, theta); - } + coeff_beta = cache_density_profile_coefficients_ ? coeff_beta_[global_index] + : density_profile_coefficients_.beta(r, theta); if (cache_domain_geometry_) { arr = arr_[global_index]; @@ -152,6 +141,9 @@ class LevelCache detDF = detDF_[global_index]; } else { + double coeff_alpha = cache_density_profile_coefficients_ ? coeff_alpha_[global_index] + : density_profile_coefficients_.alpha(r, theta); + compute_jacobian_elements(domain_geometry_, r, theta, coeff_alpha, arr, att, art, detDF); } } @@ -160,7 +152,7 @@ class LevelCache const DomainGeometry& domain_geometry_; const DensityProfileCoefficients& density_profile_coefficients_; - bool cache_density_profile_coefficients_; // cache alpha(r_i), beta(r_i) + bool cache_density_profile_coefficients_; // cache alpha(r, theta), beta(r, theta) Vector coeff_alpha_; Vector coeff_beta_; diff --git a/include/PolarGrid/multiindex.h b/include/PolarGrid/multiindex.h deleted file mode 100755 index 32508f0b..00000000 --- a/include/PolarGrid/multiindex.h +++ /dev/null @@ -1,46 +0,0 @@ -#pragma once - -#include -#include -#include - -#include "../common/space_dimension.h" - -class MultiIndex -{ -public: - static_assert(space_dimension > 0 && space_dimension <= 3, "Invalid space dimension"); - - //! Creates a multi index with undefined coordinate values. - MultiIndex() = default; - - //! Initializes all coordinate values of the multi index with `value`. - explicit MultiIndex(int value); - //! Initializes 2-dimensional multi index coordinate values with `i` and `j`. - //! Works only if `space_dimension` is 2. - explicit MultiIndex(int i, int j); - //! Initializes 3-dimensional multi index coordinate values with `i`, `j` and `k`. - //! Works only if `space_dimension` is 3. - explicit MultiIndex(int i, int j, int k); - - //! Returns the size of the multi index. - //! This is equal to `space_dimension`. - int size() const; - - //! Tests if two multi indices are equal - bool operator==(const MultiIndex& other) const; - //! Tests if two multi indices are inequal - bool operator!=(const MultiIndex& other) const; - - //! Returns the `i`th coordinate value of the multi index. - const int& operator[](int i) const; - - //! Returns a mutable reference to the `i`th coordinate value of the multi index. - int& operator[](int i); - - int* data(); - const int* data() const; - -private: - int data_[space_dimension]; -}; \ No newline at end of file diff --git a/include/PolarGrid/point.h b/include/PolarGrid/point.h deleted file mode 100755 index 2a288e9c..00000000 --- a/include/PolarGrid/point.h +++ /dev/null @@ -1,45 +0,0 @@ -#pragma once - -#include -#include -#include - -#include "../common/equals.h" -#include "../common/space_dimension.h" - -class Point -{ -public: - static_assert(space_dimension > 0 && space_dimension <= 3, "Invalid space dimension"); - - //! Creates a point with undefined coordinate values. - Point() = default; - - //! Initializes all coordinate values of the point with `value`. - explicit Point(double value); - //! Initializes 2-dimensional point coordinate values with `x` and `y`. - //! Works only if `space_dimension` is 2. - explicit Point(double x, double y); - //! Initializes 3-dimensional point coordinate values with `x`, `y` and `z`. - //! Works only if `space_dimension` is 3. - explicit Point(double x, double y, double z); - - //! Returns the size of the point. - //! This is equal to `space_dimension`. - int size() const; - - //! Returns the `i`th coordinate value of the point. - double operator[](int i) const; - - //! Returns a mutable reference to the `i`th coordinate value of the point. - double& operator[](int i); - -private: - double data_[space_dimension]; -}; - -bool equals(const Point& lhs, const Point& rhs); -double norm(const Point& point); -void add(Point& result, const Point& lhs, const Point& rhs); -void subtract(Point& result, const Point& lhs, const Point& rhs); -void multiply(Point& result, double scalar, const Point& lhs); \ No newline at end of file diff --git a/include/PolarGrid/polargrid.h b/include/PolarGrid/polargrid.h index 15b9e05d..a05c44d9 100644 --- a/include/PolarGrid/polargrid.h +++ b/include/PolarGrid/polargrid.h @@ -19,13 +19,6 @@ #include "../common/equals.h" -#include "../PolarGrid/multiindex.h" -#include "../PolarGrid/point.h" - -// The PolarGrid class implements a donut-shaped 2D grid. -// It is designed to handle polar coordinates, providing functionalities -// for storing data points and performing operations on them. - class PolarGrid { public: @@ -44,12 +37,11 @@ class PolarGrid // Optimized, inlined indexing. int wrapThetaIndex(const int unwrapped_theta_index) const; int index(const int r_index, const int unwrapped_theta_index) const; - int fastIndex(const int r_index, const int theta_index) const; void multiIndex(const int node_index, int& r_index, int& theta_index) const; + // Grid Parameters // Number of grid nodes int numberOfNodes() const; - // Grid Parameters // Get the number of grid points in radial direction int nr() const; // Get the number of angular divisions @@ -58,9 +50,6 @@ class PolarGrid double radius(const int r_index) const; // Get the angle at a specific angular index double theta(const int theta_index) const; - // Get all radii and angles available which define the grid - const std::vector& radii() const; - const std::vector& angles() const; // Grid distances // Get the radial distance to the next consecutive radial node at a specified radial index. @@ -80,48 +69,8 @@ class PolarGrid // Get the number of nodes in radial smoother. int numberRadialSmootherNodes() const; - // ------------------------------------------- // - // Unoptimized Indexing and Neighbor Retrieval // - // ------------------------------------------- // - // Node Indexing (based on the combined circle-radial smoother) - // Get the index of a node based on its position. - int index(const MultiIndex& position) const; - // Get the position of a node based on its index. - MultiIndex multiIndex(const int node_index) const; - // Get the polar coordinates of a node based on its position. - Point polarCoordinates(const MultiIndex& position) const; - // Get adjacent neighbors of a node. - // If the neighbor index is -1, then there is no neighboring node in that direction. - // - The first entry (neighbors[0]) represents the radial direction (r): - // - First: inward neighbor (r - 1) - // - Second: outward neighbor (r + 1) - // - The second entry (neighbors[1]) represents the angular direction (theta): - // - First: counterclockwise neighbor (theta - 1) - // - Second: clockwise neighbor (theta + 1) - void adjacentNeighborsOf(const MultiIndex& position, - std::array, space_dimension>& neighbors) const; - // Get diagonal neighbors of a node. - // If the neighbor index is -1, then there is no neighboring node in that direction. - // - The first entry (neighbors[0]) represents: - // - First: bottom left neighbor (r - 1, theta - 1) - // - Second: bottom right neighbor (r + 1, theta - 1) - // - The second entry (neighbors[1]) represents: - // - First: top left neighbor (r - 1, theta + 1) - // - Second: top right neighbor (r + 1, theta + 1) - void diagonalNeighborsOf(const MultiIndex& position, - std::array, space_dimension>& neighbors) const; - // Neighbor distances - // Get distances to adjacent neighbors of a node. - // If there is no neighboring node in that direction, then the neighbor distance is 0. - void adjacentNeighborDistances(const MultiIndex& position, - std::array, space_dimension>& neighbor_distances) const; - private: - // --------------- // - // Private members // - // --------------- // - - // We will use the convention: + // We use the convention: // radii = [R0, ..., R], angles = [0, ..., 2*pi] // Note that ntheta will be one less than the size of angles since 0 and 2pi are the same point. int nr_; // number of nodes in radial direction @@ -155,15 +104,11 @@ class PolarGrid * - numberCircularSmootherNodes + numberRadialSmootherNodes = number_of_nodes() */ - // ------------------------ // - // Private Helper Functions // - // ------------------------ // - // Check parameter validity void checkParameters(const std::vector& radii, const std::vector& angles) const; + // Initialize radial_spacings_, angular_spacings_ void initializeDistances(); - // Initializes line splitting parameters for Circle/radial indexing. // splitting_radius: The radius value used for dividing the smoother into a circular and radial section. // If std::nullopt, automatic line-splitting is enabled. @@ -192,4 +137,4 @@ class PolarGrid // ---------------------------------------------------- // PolarGrid coarseningGrid(const PolarGrid& grid); -#include "polargrid.inl" // Include the inline function definitions \ No newline at end of file +#include "polargrid.inl" // Include the inline function definitions diff --git a/include/PolarGrid/polargrid.inl b/include/PolarGrid/polargrid.inl index 512343dc..9d9df122 100755 --- a/include/PolarGrid/polargrid.inl +++ b/include/PolarGrid/polargrid.inl @@ -28,6 +28,11 @@ inline double PolarGrid::theta(const int theta_index) const return angles_[theta_index]; } +inline double PolarGrid::smootherSplittingRadius() const +{ + return smoother_splitting_radius_; +} + // Get the number of circles in the circular smoother. inline int PolarGrid::numberSmootherCircles() const { @@ -60,14 +65,9 @@ inline double PolarGrid::angularSpacing(const int unwrapped_theta_index) const { // unwrapped_theta_index may be negative or larger than ntheta() to allow for periodicity. const int theta_index = wrapThetaIndex(unwrapped_theta_index); - assert(theta_index >= 0 && theta_index < ntheta()); return angular_spacings_[theta_index]; } -// ------------------ // -// Optimized indexing // -// ------------------ // - inline int PolarGrid::wrapThetaIndex(const int unwrapped_theta_index) const { // The unwrapped_theta_index may be negative or exceed the number of theta steps (ntheta()), @@ -80,40 +80,34 @@ inline int PolarGrid::wrapThetaIndex(const int unwrapped_theta_index) const // This effectively computes unwrapped_theta_index % ntheta(), because it discards all higher bits. // // If ntheta is not a power of two, we use the standard modulo approach to handle wrapping. - return is_ntheta_PowerOfTwo_ ? unwrapped_theta_index & (ntheta() - 1) : (unwrapped_theta_index % ntheta() + ntheta()) % ntheta(); + int theta_index = is_ntheta_PowerOfTwo_ ? unwrapped_theta_index & (ntheta() - 1) + : (unwrapped_theta_index % ntheta() + ntheta()) % ntheta(); + assert(0 <= theta_index && theta_index < ntheta()); + return theta_index; } inline int PolarGrid::index(const int r_index, const int unwrapped_theta_index) const { // unwrapped_theta_index may be negative or larger than ntheta() to allow for periodicity. assert(0 <= r_index && r_index < nr()); - const int theta_index = wrapThetaIndex(unwrapped_theta_index); - assert(0 <= theta_index && theta_index < ntheta()); - return r_index < numberSmootherCircles() - ? theta_index + ntheta() * r_index - : numberCircularSmootherNodes() + r_index - numberSmootherCircles() + lengthSmootherRadial() * theta_index; -} - -inline int PolarGrid::fastIndex(const int r_index, const int theta_index) const -{ - assert(0 <= r_index && r_index < nr()); - assert(0 <= theta_index && theta_index < ntheta()); - return r_index < numberSmootherCircles() - ? theta_index + ntheta() * r_index - : numberCircularSmootherNodes() + r_index - numberSmootherCircles() + lengthSmootherRadial() * theta_index; + int theta_index = wrapThetaIndex(unwrapped_theta_index); + int global_index = + r_index < numberSmootherCircles() + ? theta_index + ntheta() * r_index + : numberCircularSmootherNodes() + r_index - numberSmootherCircles() + lengthSmootherRadial() * theta_index; + assert(0 <= global_index && global_index < numberOfNodes()); + return global_index; } inline void PolarGrid::multiIndex(const int node_index, int& r_index, int& theta_index) const { assert(0 <= node_index && node_index < numberOfNodes()); - if (node_index < numberCircularSmootherNodes()) - { - r_index = node_index / ntheta(); - theta_index = is_ntheta_PowerOfTwo_ ? node_index & (ntheta() - 1) : node_index % ntheta(); + if (node_index < numberCircularSmootherNodes()) { + r_index = node_index / ntheta(); + theta_index = wrapThetaIndex(node_index); } - else - { + else { theta_index = (node_index - numberCircularSmootherNodes()) / lengthSmootherRadial(); - r_index = numberSmootherCircles() + (node_index - numberCircularSmootherNodes()) % lengthSmootherRadial(); + r_index = numberSmootherCircles() + (node_index - numberCircularSmootherNodes()) % lengthSmootherRadial(); } } diff --git a/include/common/space_dimension.h b/include/common/space_dimension.h deleted file mode 100644 index 9ba369e1..00000000 --- a/include/common/space_dimension.h +++ /dev/null @@ -1,3 +0,0 @@ -#pragma once - -constexpr inline int space_dimension = 2; \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1b7ea608..eb06eece 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -6,10 +6,8 @@ add_subdirectory(InputFunctions) # Gather all source files # file(GLOB_RECURSE POLAR_GRID_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/PolarGrid/*.cpp) set(POLAR_GRID_SOURCES - ${CMAKE_CURRENT_SOURCE_DIR}/PolarGrid/anisotropic_division.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/PolarGrid/multiindex.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/PolarGrid/point.cpp ${CMAKE_CURRENT_SOURCE_DIR}/PolarGrid/polargrid.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/PolarGrid/anisotropic_division.cpp ) # file(GLOB_RECURSE GMG_POLAR_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/*.cpp) diff --git a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/smootherSolver.cpp b/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/smootherSolver.cpp index 034c8ee5..8124c66f 100644 --- a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/smootherSolver.cpp +++ b/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/smootherSolver.cpp @@ -826,30 +826,8 @@ void ExtrapolatedSmootherGive::applyAscOrthoCircleSection(const int i_r, const S const double theta = grid_.theta(i_theta); const int index = grid_.index(i_r, i_theta); - double coeff_beta; - if (level_cache_.cacheDensityProfileCoefficients()) - coeff_beta = level_cache_.coeff_beta()[index]; - else - coeff_beta = density_profile_coefficients_.beta(r, theta); - - double coeff_alpha; - if (!level_cache_.cacheDomainGeometry()) { - if (level_cache_.cacheDensityProfileCoefficients()) - coeff_alpha = level_cache_.coeff_alpha()[index]; - else - coeff_alpha = density_profile_coefficients_.alpha(r, theta); - } - - double arr, att, art, detDF; - if (level_cache_.cacheDomainGeometry()) { - arr = level_cache_.arr()[index]; - att = level_cache_.att()[index]; - art = level_cache_.art()[index]; - detDF = level_cache_.detDF()[index]; - } - else { - compute_jacobian_elements(domain_geometry_, r, theta, coeff_alpha, arr, att, art, detDF); - } + double coeff_beta, arr, att, art, detDF; + level_cache_.obtainValues(i_r, i_theta, index, r, theta, coeff_beta, arr, att, art, detDF); // Apply Asc Ortho at the current node NODE_APPLY_ASC_ORTHO_CIRCLE_GIVE(i_r, i_theta, r, theta, grid_, DirBC_Interior_, smoother_color, x, rhs, temp, @@ -868,30 +846,8 @@ void ExtrapolatedSmootherGive::applyAscOrthoRadialSection(const int i_theta, con const double r = grid_.radius(i_r); const int index = grid_.index(i_r, i_theta); - double coeff_beta; - if (level_cache_.cacheDensityProfileCoefficients()) - coeff_beta = level_cache_.coeff_beta()[index]; - else - coeff_beta = density_profile_coefficients_.beta(r, theta); - - double coeff_alpha; - if (!level_cache_.cacheDomainGeometry()) { - if (level_cache_.cacheDensityProfileCoefficients()) - coeff_alpha = level_cache_.coeff_alpha()[index]; - else - coeff_alpha = density_profile_coefficients_.alpha(r, theta); - } - - double arr, att, art, detDF; - if (level_cache_.cacheDomainGeometry()) { - arr = level_cache_.arr()[index]; - att = level_cache_.att()[index]; - art = level_cache_.art()[index]; - detDF = level_cache_.detDF()[index]; - } - else { - compute_jacobian_elements(domain_geometry_, r, theta, coeff_alpha, arr, att, art, detDF); - } + double coeff_beta, arr, att, art, detDF; + level_cache_.obtainValues(i_r, i_theta, index, r, theta, coeff_beta, arr, att, art, detDF); // Apply Asc Ortho at the current node NODE_APPLY_ASC_ORTHO_RADIAL_GIVE(i_r, i_theta, r, theta, grid_, DirBC_Interior_, smoother_color, x, rhs, temp, diff --git a/src/GMGPolar/setup.cpp b/src/GMGPolar/setup.cpp index fa5dd591..4dc67886 100644 --- a/src/GMGPolar/setup.cpp +++ b/src/GMGPolar/setup.cpp @@ -61,8 +61,10 @@ void GMGPolar::setup() levels_.emplace_back(level_depth, std::move(finest_grid), std::move(finest_levelCache), extrapolation_, FMG_); for (level_depth = 1; level_depth < number_of_levels_; level_depth++) { - auto current_grid = std::make_unique(coarseningGrid(levels_[level_depth - 1].grid())); - auto current_levelCache = std::make_unique(levels_[level_depth - 1], *current_grid); + auto current_grid = std::make_unique(coarseningGrid(levels_[level_depth - 1].grid())); + auto current_levelCache = + std::make_unique(*current_grid, density_profile_coefficients_, domain_geometry_, + cache_density_profile_coefficients_, cache_domain_geometry_); levels_.emplace_back(level_depth, std::move(current_grid), std::move(current_levelCache), extrapolation_, FMG_); } @@ -271,7 +273,8 @@ void GMGPolar::printSettings() const const PolarGrid& finest_grid = levels_.front().grid(); const PolarGrid& coarsest_grid = levels_.back().grid(); - std::cout << "r ∈ [" << finest_grid.radii().front() << ", " << finest_grid.radii().back() << "], θ ∈ [0, 2π]\n"; + std::cout << "r ∈ [" << finest_grid.radius(0) << ", " << finest_grid.radius(finest_grid.nr() - 1) + << "], θ ∈ [0, 2π]\n"; std::cout << "(nr × nθ) = (" << finest_grid.nr() << " × " << finest_grid.ntheta() << ") → (" << coarsest_grid.nr() << " × " << coarsest_grid.ntheta() << ")\n"; std::cout << "Smoother: " << finest_grid.numberSmootherCircles() << " circles\n"; 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/src/Level/levelCache.cpp b/src/Level/levelCache.cpp index c8d31243..8271e293 100644 --- a/src/Level/levelCache.cpp +++ b/src/Level/levelCache.cpp @@ -16,6 +16,8 @@ LevelCache::LevelCache(const PolarGrid& grid, const DensityProfileCoefficients& , art_("art", cache_domain_geometry ? grid.numberOfNodes() : 0) , detDF_("detDF", cache_domain_geometry ? grid.numberOfNodes() : 0) { + // Pre-compute and store alpha/beta coefficients at all grid nodes to avoid + // repeated expensive evaluations during runtime computations if (cache_density_profile_coefficients_) { #pragma omp parallel for for (int i_r = 0; i_r < grid.nr(); i_r++) { @@ -31,15 +33,20 @@ LevelCache::LevelCache(const PolarGrid& grid, const DensityProfileCoefficients& } } + // Pre-compute and store Jacobian matrix elements (arr, att, art, detDF) at all grid nodes + // to avoid repeated coordinate transformation calculations during domain operations if (cache_domain_geometry_) { + // We split the loops into two regions to better respect the + // access patterns of the smoother and improve cache locality + + // Circular Indexing Section #pragma omp parallel for for (int i_r = 0; i_r < grid.numberSmootherCircles(); i_r++) { const double r = grid.radius(i_r); for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { - const double theta = grid.theta(i_theta); - const int index = grid.index(i_r, i_theta); - - double coeff_alpha = density_profile_coefficients.alpha(r, theta); + const double theta = grid.theta(i_theta); + const int index = grid.index(i_r, i_theta); + const double coeff_alpha = density_profile_coefficients.alpha(r, theta); double arr, att, art, detDF; compute_jacobian_elements(domain_geometry_, r, theta, coeff_alpha, arr, att, art, detDF); @@ -49,21 +56,14 @@ LevelCache::LevelCache(const PolarGrid& grid, const DensityProfileCoefficients& art_(index) = art; } } - + // Radial Indexing Section #pragma omp parallel for for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { const double theta = grid.theta(i_theta); for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { - const double r = grid.radius(i_r); - const int index = grid.index(i_r, i_theta); - - double coeff_alpha; - if (cache_density_profile_coefficients_ && !cache_domain_geometry_) { - coeff_alpha = coeff_alpha_(index); - } - else { - coeff_alpha = density_profile_coefficients.alpha(r, theta); - } + const double r = grid.radius(i_r); + const int index = grid.index(i_r, i_theta); + const double coeff_alpha = density_profile_coefficients.alpha(r, theta); double arr, att, art, detDF; compute_jacobian_elements(domain_geometry_, r, theta, coeff_alpha, arr, att, art, detDF); @@ -76,62 +76,6 @@ LevelCache::LevelCache(const PolarGrid& grid, const DensityProfileCoefficients& } } -LevelCache::LevelCache(const Level& previous_level, const PolarGrid& current_grid) - : domain_geometry_(previous_level.levelCache().domainGeometry()) - , density_profile_coefficients_(previous_level.levelCache().densityProfileCoefficients()) - , cache_density_profile_coefficients_(previous_level.levelCache().cacheDensityProfileCoefficients()) - , coeff_alpha_("coeff_alpha", - previous_level.levelCache().coeff_alpha().size() > 0 ? current_grid.numberOfNodes() : 0) - , coeff_beta_("coeff_beta", previous_level.levelCache().coeff_beta().size() > 0 ? current_grid.numberOfNodes() : 0) - , cache_domain_geometry_(previous_level.levelCache().cacheDomainGeometry()) - , arr_("arr", previous_level.levelCache().arr().size() > 0 ? current_grid.numberOfNodes() : 0) - , att_("att", previous_level.levelCache().att().size() > 0 ? current_grid.numberOfNodes() : 0) - , art_("art", previous_level.levelCache().art().size() > 0 ? current_grid.numberOfNodes() : 0) - , detDF_("detDF", previous_level.levelCache().detDF().size() > 0 ? current_grid.numberOfNodes() : 0) -{ - const auto& previous_level_cache = previous_level.levelCache(); - - if (previous_level_cache.cacheDensityProfileCoefficients()) { -#pragma omp parallel for - for (int i_r = 0; i_r < current_grid.nr(); i_r++) { - for (int i_theta = 0; i_theta < current_grid.ntheta(); i_theta++) { - const int current_index = current_grid.index(i_r, i_theta); - const int previous_index = previous_level.grid().index(2 * i_r, 2 * i_theta); - - if (!previous_level_cache.cacheDomainGeometry()) { - coeff_alpha_[current_index] = previous_level_cache.coeff_alpha()[previous_index]; - } - coeff_beta_[current_index] = previous_level_cache.coeff_beta()[previous_index]; - } - } - } - - if (previous_level_cache.cacheDomainGeometry()) { -#pragma omp parallel for - for (int i_r = 0; i_r < current_grid.numberSmootherCircles(); i_r++) { - for (int i_theta = 0; i_theta < current_grid.ntheta(); i_theta++) { - const int current_index = current_grid.index(i_r, i_theta); - const int previous_index = previous_level.grid().index(2 * i_r, 2 * i_theta); - arr_[current_index] = previous_level_cache.arr()[previous_index]; - att_[current_index] = previous_level_cache.att()[previous_index]; - art_[current_index] = previous_level_cache.art()[previous_index]; - detDF_[current_index] = previous_level_cache.detDF()[previous_index]; - } - } -#pragma omp parallel for - for (int i_theta = 0; i_theta < current_grid.ntheta(); i_theta++) { - for (int i_r = current_grid.numberSmootherCircles(); i_r < current_grid.nr(); i_r++) { - const int current_index = current_grid.index(i_r, i_theta); - const int previous_index = previous_level.grid().index(2 * i_r, 2 * i_theta); - arr_[current_index] = previous_level_cache.arr()[previous_index]; - att_[current_index] = previous_level_cache.att()[previous_index]; - art_[current_index] = previous_level_cache.art()[previous_index]; - detDF_[current_index] = previous_level_cache.detDF()[previous_index]; - } - } - } -} - const DensityProfileCoefficients& LevelCache::densityProfileCoefficients() const { return density_profile_coefficients_; diff --git a/src/PolarGrid/multiindex.cpp b/src/PolarGrid/multiindex.cpp deleted file mode 100644 index ec3c16a3..00000000 --- a/src/PolarGrid/multiindex.cpp +++ /dev/null @@ -1,71 +0,0 @@ -#include "../../include/PolarGrid/multiindex.h" - -MultiIndex::MultiIndex(int value) -{ - assert(space_dimension >= 0); - std::fill(std::begin(data_), std::end(data_), value); -} - -MultiIndex::MultiIndex(int i, int j) -{ - assert(space_dimension == 2); - data_[0] = i; - data_[1] = j; -} - -MultiIndex::MultiIndex(int i, int j, int k) -{ - assert(space_dimension == 3); - data_[0] = i; - data_[1] = j; - data_[2] = k; -} - -int MultiIndex::size() const -{ - return space_dimension; -} - -bool MultiIndex::operator==(const MultiIndex& other) const -{ - for (int d = 0; d < space_dimension; d++) { - if (data_[d] != other.data_[d]) { - return false; - } - } - return true; -} - -bool MultiIndex::operator!=(const MultiIndex& other) const -{ - for (int d = 0; d < space_dimension; d++) { - if (data_[d] != other.data_[d]) { - return true; - } - } - return false; -} - -const int& MultiIndex::operator[](int i) const -{ - assert(i >= 0); - assert(i < space_dimension); - return data_[i]; -} - -int& MultiIndex::operator[](int i) -{ - assert(i >= 0); - assert(i < space_dimension); - return data_[i]; -} - -int* MultiIndex::data() -{ - return data_; -} - -const int* MultiIndex::data() const -{ - return data_; -} \ No newline at end of file diff --git a/src/PolarGrid/point.cpp b/src/PolarGrid/point.cpp deleted file mode 100644 index cb0bdf75..00000000 --- a/src/PolarGrid/point.cpp +++ /dev/null @@ -1,85 +0,0 @@ -#include "../../include/PolarGrid/point.h" - -Point::Point(double value) -{ - assert(space_dimension >= 0); - std::fill(std::begin(data_), std::end(data_), value); -} - -Point::Point(double x, double y) -{ - assert(space_dimension == 2); - data_[0] = x; - data_[1] = y; -} - -Point::Point(double x, double y, double z) -{ - assert(space_dimension == 3); - data_[0] = x; - data_[1] = y; - data_[2] = z; -} - -int Point::size() const -{ - return space_dimension; -} - -double Point::operator[](int i) const -{ - assert(i >= 0); - assert(i < space_dimension); - return data_[i]; -} - -double& Point::operator[](int i) -{ - assert(i >= 0); - assert(i < space_dimension); - return data_[i]; -} - -bool equals(const Point& lhs, const Point& rhs) -{ - for (int i = 0; i < space_dimension; ++i) { - if (!equals(lhs[i], rhs[i])) - return false; - } - return true; -} - -double norm(const Point& point) -{ - double result = 0; - for (int i = 0; i < space_dimension; ++i) { - result += point[i] * point[i]; - } - return sqrt(result); -} - -void add(Point& result, const Point& lhs, const Point& rhs) -{ - assert(lhs.size() == rhs.size()); - assert(result.size() == lhs.size()); - for (int i = 0; i < space_dimension; i++) { - result[i] = lhs[i] + rhs[i]; - } -} - -void subtract(Point& result, const Point& lhs, const Point& rhs) -{ - assert(lhs.size() == rhs.size()); - assert(result.size() == lhs.size()); - for (int i = 0; i < space_dimension; i++) { - result[i] = lhs[i] - rhs[i]; - } -} - -void multiply(Point& result, double scalar, const Point& lhs) -{ - assert(result.size() == lhs.size()); - for (int i = 0; i < space_dimension; i++) { - result[i] = scalar * lhs[i]; - } -} \ No newline at end of file diff --git a/src/PolarGrid/polargrid.cpp b/src/PolarGrid/polargrid.cpp index b2a1f732..3eb07423 100644 --- a/src/PolarGrid/polargrid.cpp +++ b/src/PolarGrid/polargrid.cpp @@ -232,143 +232,9 @@ PolarGrid coarseningGrid(const PolarGrid& fineGrid) } } -// ---------------- // -// Getter Functions // -// ---------------- // - -const std::vector& PolarGrid::radii() const -{ - return radii_; -} -const std::vector& PolarGrid::angles() const -{ - return angles_; -} - -// Get the radius at which the grid is split into circular and radial smoothing -double PolarGrid::smootherSplittingRadius() const -{ - return smoother_splitting_radius_; -} - -// ------------------------------------- // -// Definition of node indexing. // -// Based on the circular-radial smoother // -// ------------------------------------- // - -/* OPTIMIZED INDEXING IS DEFINED IN include/PolarGrid/polargrid.inl */ - -int PolarGrid::index(const MultiIndex& position) const -{ - assert(position[0] >= 0 && position[0] < nr()); - assert(position[1] >= 0 && position[1] < ntheta()); - if (position[0] < numberSmootherCircles()) { - return position[1] + ntheta() * position[0]; - } - else { - return numberCircularSmootherNodes() + - (position[0] - numberSmootherCircles() + lengthSmootherRadial() * position[1]); - } -} - -MultiIndex PolarGrid::multiIndex(const int node_index) const -{ - assert(0 <= node_index && node_index < numberOfNodes()); - if (node_index < numberCircularSmootherNodes()) { - auto result = std::div(node_index, ntheta()); - return MultiIndex(result.quot, result.rem); - } - else { - auto result = std::div(node_index - numberCircularSmootherNodes(), lengthSmootherRadial()); - return MultiIndex(numberSmootherCircles() + result.rem, result.quot); - } -} - -Point PolarGrid::polarCoordinates(const MultiIndex& position) const -{ - assert(position[0] >= 0 && position[0] < nr()); - assert(position[1] >= 0 && position[1] < ntheta()); - return Point(radii_[position[0]], angles_[(position[1])]); -} - -void PolarGrid::adjacentNeighborDistances( - const MultiIndex& position, std::array, space_dimension>& neighbor_distance) const -{ - assert(position[0] >= 0 && position[0] < nr()); - assert(position[1] >= 0 && position[1] < ntheta()); - - neighbor_distance[0].first = (position[0] <= 0) ? 0.0 : radialSpacing(position[0] - 1); - neighbor_distance[0].second = (position[0] >= nr() - 1) ? 0.0 : radialSpacing(position[0]); - - neighbor_distance[1].first = angularSpacing(position[1] - 1); - neighbor_distance[1].second = angularSpacing(position[1]); -} - -void PolarGrid::adjacentNeighborsOf(const MultiIndex& position, - std::array, space_dimension>& neighbors) const -{ - assert(position[0] >= 0 && position[0] < nr()); - assert(position[1] >= 0 && position[1] < ntheta()); - - MultiIndex neigbor_position = position; - neigbor_position[0] -= 1; - neighbors[0].first = (neigbor_position[0] < 0) ? -1 : index(neigbor_position); - - neigbor_position = position; - neigbor_position[0] += 1; - neighbors[0].second = (neigbor_position[0] >= nr()) ? -1 : index(neigbor_position); - - neigbor_position = position; - neigbor_position[1] -= 1; - if (neigbor_position[1] < 0) - neigbor_position[1] += ntheta(); - neighbors[1].first = index(neigbor_position); - - neigbor_position = position; - neigbor_position[1] += 1; - if (neigbor_position[1] >= ntheta()) - neigbor_position[1] -= ntheta(); - neighbors[1].second = index(neigbor_position); -} - -void PolarGrid::diagonalNeighborsOf(const MultiIndex& position, - std::array, space_dimension>& neighbors) const -{ - assert(position[0] >= 0 && position[0] < nr()); - assert(position[1] >= 0 && position[1] < ntheta()); - - MultiIndex neigbor_position = position; - neigbor_position[0] -= 1; - neigbor_position[1] -= 1; - if (neigbor_position[1] < 0) - neigbor_position[1] += ntheta(); - neighbors[0].first = (neigbor_position[0] < 0) ? -1 : index(neigbor_position); - - neigbor_position = position; - neigbor_position[0] += 1; - neigbor_position[1] -= 1; - if (neigbor_position[1] < 0) - neigbor_position[1] += ntheta(); - neighbors[0].second = (neigbor_position[0] >= nr()) ? -1 : index(neigbor_position); - - neigbor_position = position; - neigbor_position[0] -= 1; - neigbor_position[1] += 1; - if (neigbor_position[1] >= ntheta()) - neigbor_position[1] -= ntheta(); - neighbors[1].first = (neigbor_position[0] < 0) ? -1 : index(neigbor_position); - - neigbor_position = position; - neigbor_position[0] += 1; - neigbor_position[1] += 1; - if (neigbor_position[1] >= ntheta()) - neigbor_position[1] -= ntheta(); - neighbors[1].second = (neigbor_position[0] >= nr()) ? -1 : index(neigbor_position); -} - // ------------------------ // // Check parameter validity // -// ---------------------..- // +// ------------------------ // void PolarGrid::checkParameters(const std::vector& radii, const std::vector& angles) const { diff --git a/src/Smoother/SmootherGive/smootherSolver.cpp b/src/Smoother/SmootherGive/smootherSolver.cpp index 86e4474f..28843757 100644 --- a/src/Smoother/SmootherGive/smootherSolver.cpp +++ b/src/Smoother/SmootherGive/smootherSolver.cpp @@ -334,30 +334,8 @@ void SmootherGive::applyAscOrthoCircleSection(const int i_r, const SmootherColor const double theta = grid_.theta(i_theta); const int index = grid_.index(i_r, i_theta); - double coeff_beta; - if (level_cache_.cacheDensityProfileCoefficients()) - coeff_beta = level_cache_.coeff_beta()[index]; - else - coeff_beta = density_profile_coefficients_.beta(r, theta); - - double coeff_alpha; - if (!level_cache_.cacheDomainGeometry()) { - if (level_cache_.cacheDensityProfileCoefficients()) - coeff_alpha = level_cache_.coeff_alpha()[index]; - else - coeff_alpha = density_profile_coefficients_.alpha(r, theta); - } - - double arr, att, art, detDF; - if (level_cache_.cacheDomainGeometry()) { - arr = level_cache_.arr()[index]; - att = level_cache_.att()[index]; - art = level_cache_.art()[index]; - detDF = level_cache_.detDF()[index]; - } - else { - compute_jacobian_elements(domain_geometry_, r, theta, coeff_alpha, arr, att, art, detDF); - } + double coeff_beta, arr, att, art, detDF; + level_cache_.obtainValues(i_r, i_theta, index, r, theta, coeff_beta, arr, att, art, detDF); // Apply Asc Ortho at the current node NODE_APPLY_ASC_ORTHO_CIRCLE_GIVE(i_r, i_theta, r, theta, grid_, DirBC_Interior_, smoother_color, x, rhs, temp, @@ -375,30 +353,9 @@ void SmootherGive::applyAscOrthoRadialSection(const int i_theta, const SmootherC const double r = grid_.radius(i_r); const int index = grid_.index(i_r, i_theta); - double coeff_beta; - if (level_cache_.cacheDensityProfileCoefficients()) - coeff_beta = level_cache_.coeff_beta()[index]; - else - coeff_beta = density_profile_coefficients_.beta(r, theta); - - double coeff_alpha; - if (!level_cache_.cacheDomainGeometry()) { - if (level_cache_.cacheDensityProfileCoefficients()) - coeff_alpha = level_cache_.coeff_alpha()[index]; - else - coeff_alpha = density_profile_coefficients_.alpha(r, theta); - } + double coeff_beta, arr, att, art, detDF; + level_cache_.obtainValues(i_r, i_theta, index, r, theta, coeff_beta, arr, att, art, detDF); - double arr, att, art, detDF; - if (level_cache_.cacheDomainGeometry()) { - arr = level_cache_.arr()[index]; - att = level_cache_.att()[index]; - art = level_cache_.art()[index]; - detDF = level_cache_.detDF()[index]; - } - else { - compute_jacobian_elements(domain_geometry_, r, theta, coeff_alpha, arr, att, art, detDF); - } // Apply Asc Ortho at the current node NODE_APPLY_ASC_ORTHO_RADIAL_GIVE(i_r, i_theta, r, theta, grid_, DirBC_Interior_, smoother_color, x, rhs, temp, arr, att, art, detDF, coeff_beta); diff --git a/tests/DirectSolver/directSolver.cpp b/tests/DirectSolver/directSolver.cpp index 84e8c5a3..c9896f18 100644 --- a/tests/DirectSolver/directSolver.cpp +++ b/tests/DirectSolver/directSolver.cpp @@ -112,8 +112,9 @@ TEST(DirectSolverTest, directSolver_DirBC_Interior) ASSERT_EQ(solution_Give.size(), solution_Take.size()); for (uint index = 0; index < solution_Give.size(); index++) { - MultiIndex alpha = level.grid().multiIndex(index); - if (alpha[0] == 0 && !DirBC_Interior) + int i_r, i_theta; + level.grid().multiIndex(index, i_r, i_theta); + if (i_r == 0 && !DirBC_Interior) ASSERT_NEAR(solution_Give(index), solution_Take(index), 1e-11); else ASSERT_NEAR(solution_Give(index), solution_Take(index), 1e-11); @@ -167,8 +168,9 @@ TEST(DirectSolverTest, directSolver_AcrossOrigin) ASSERT_EQ(solution_Give.size(), solution_Take.size()); for (uint index = 0; index < solution_Give.size(); index++) { - MultiIndex alpha = level.grid().multiIndex(index); - if (alpha[0] == 0 && !DirBC_Interior) + int i_r, i_theta; + level.grid().multiIndex(index, i_r, i_theta); + if (i_r == 0 && !DirBC_Interior) ASSERT_NEAR(solution_Give(index), solution_Take(index), 1e-8); else ASSERT_NEAR(solution_Give(index), solution_Take(index), 1e-8); diff --git a/tests/DirectSolver/directSolverNoMumps.cpp b/tests/DirectSolver/directSolverNoMumps.cpp index 49207283..89066c85 100644 --- a/tests/DirectSolver/directSolverNoMumps.cpp +++ b/tests/DirectSolver/directSolverNoMumps.cpp @@ -110,8 +110,9 @@ TEST(DirectSolverTestNoMumps, directSolver_DirBC_Interior) ASSERT_EQ(solution_Give.size(), solution_Take.size()); for (uint index = 0; index < solution_Give.size(); index++) { - MultiIndex alpha = level.grid().multiIndex(index); - if (alpha[0] == 0 && !DirBC_Interior) + int i_r, i_theta; + level.grid().multiIndex(index, i_r, i_theta); + if (i_r == 0 && !DirBC_Interior) ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11); else ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11); @@ -165,8 +166,9 @@ TEST(DirectSolverTestNoMumps, directSolver_AcrossOrigin) ASSERT_EQ(solution_Give.size(), solution_Take.size()); for (uint index = 0; index < solution_Give.size(); index++) { - MultiIndex alpha = level.grid().multiIndex(index); - if (alpha[0] == 0 && !DirBC_Interior) + int i_r, i_theta; + level.grid().multiIndex(index, i_r, i_theta); + if (i_r == 0 && !DirBC_Interior) ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-8); else ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-9); diff --git a/tests/ExtrapolatedSmoother/extrapolated_smoother.cpp b/tests/ExtrapolatedSmoother/extrapolated_smoother.cpp index 90729ec6..a8bddbbf 100644 --- a/tests/ExtrapolatedSmoother/extrapolated_smoother.cpp +++ b/tests/ExtrapolatedSmoother/extrapolated_smoother.cpp @@ -97,8 +97,9 @@ TEST(ExtrapolatedSmootherTest, extrapolatedSmoother_DirBC_Interior) ASSERT_EQ(solution_Give.size(), solution_Take.size()); for (uint index = 0; index < solution_Give.size(); index++) { - MultiIndex alpha = level.grid().multiIndex(index); - if (alpha[0] == 0 && !DirBC_Interior) + int i_r, i_theta; + level.grid().multiIndex(index, i_r, i_theta); + if (i_r == 0 && !DirBC_Interior) ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11); else ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11); @@ -156,8 +157,9 @@ TEST(ExtrapolatedSmootherTest, extrapolatedSmoother_AcossOrigin) ASSERT_EQ(solution_Give.size(), solution_Take.size()); for (uint index = 0; index < solution_Give.size(); index++) { - MultiIndex alpha = level.grid().multiIndex(index); - if (alpha[0] == 0 && !DirBC_Interior) + int i_r, i_theta; + level.grid().multiIndex(index, i_r, i_theta); + if (i_r == 0 && !DirBC_Interior) ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-8); else ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11); 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 << ")"; + } } } diff --git a/tests/PolarGrid/polargrid.cpp b/tests/PolarGrid/polargrid.cpp index c3915acc..31a87f2a 100644 --- a/tests/PolarGrid/polargrid.cpp +++ b/tests/PolarGrid/polargrid.cpp @@ -58,17 +58,18 @@ TEST(PolarGridTest, IndexingTest) for (int i = 0; i < grid.nr(); i++) { for (int j = 0; j < grid.ntheta(); j++) { - MultiIndex alpha(i, j); - int node_index = grid.index(alpha); - MultiIndex beta = grid.multiIndex(node_index); - ASSERT_EQ(alpha[0], beta[0]); - ASSERT_EQ(alpha[1], beta[1]); + int node_index = grid.index(i, j); + int r_out, theta_out; + grid.multiIndex(node_index, r_out, theta_out); + ASSERT_EQ(i, r_out); + ASSERT_EQ(j, theta_out); } } for (int i = 0; i < grid.numberOfNodes(); i++) { - MultiIndex alpha = grid.multiIndex(i); - int node_index = grid.index(alpha); + int r_out, theta_out; + grid.multiIndex(i, r_out, theta_out); + int node_index = grid.index(r_out, theta_out); ASSERT_EQ(node_index, i); } } @@ -82,39 +83,39 @@ TEST(PolarGridTest, IndexingValuesTest) PolarGrid grid(radii, angles, splitting_radius); { - MultiIndex alpha(2, 6); - int node_index = grid.index(alpha); + int node_index = grid.index(2, 6); ASSERT_EQ(node_index, 22); - MultiIndex beta = grid.multiIndex(node_index); - ASSERT_EQ(alpha[0], beta[0]); - ASSERT_EQ(alpha[1], beta[1]); + int r_out, theta_out; + grid.multiIndex(node_index, r_out, theta_out); + ASSERT_EQ(2, r_out); + ASSERT_EQ(6, theta_out); } { - MultiIndex alpha(3, 2); - int node_index = grid.index(alpha); + int node_index = grid.index(3, 2); ASSERT_EQ(node_index, 26); - MultiIndex beta = grid.multiIndex(node_index); - ASSERT_EQ(alpha[0], beta[0]); - ASSERT_EQ(alpha[1], beta[1]); + int r_out, theta_out; + grid.multiIndex(node_index, r_out, theta_out); + ASSERT_EQ(3, r_out); + ASSERT_EQ(2, theta_out); } { - MultiIndex alpha(6, 4); - int node_index = grid.index(alpha); + int node_index = grid.index(6, 4); ASSERT_EQ(node_index, 54); - MultiIndex beta = grid.multiIndex(node_index); - ASSERT_EQ(alpha[0], beta[0]); - ASSERT_EQ(alpha[1], beta[1]); + int r_out, theta_out; + grid.multiIndex(node_index, r_out, theta_out); + ASSERT_EQ(6, r_out); + ASSERT_EQ(4, theta_out); } { - MultiIndex alpha(4, 7); - int node_index = grid.index(alpha); + int node_index = grid.index(4, 7); ASSERT_EQ(node_index, 67); - MultiIndex beta = grid.multiIndex(node_index); - ASSERT_EQ(alpha[0], beta[0]); - ASSERT_EQ(alpha[1], beta[1]); + int r_out, theta_out; + grid.multiIndex(node_index, r_out, theta_out); + ASSERT_EQ(4, r_out); + ASSERT_EQ(7, theta_out); } } @@ -126,129 +127,27 @@ TEST(PolarGridTest, CoordinatesTest) double splitting_radius = 0.6; PolarGrid grid(radii, angles, splitting_radius); - MultiIndex alpha(3, 2); - Point P1 = grid.polarCoordinates(alpha); - ASSERT_DOUBLE_EQ(P1[0], 0.5); - ASSERT_DOUBLE_EQ(P1[1], M_PI / 8); + ASSERT_DOUBLE_EQ(grid.radius(3), 0.5); + ASSERT_DOUBLE_EQ(grid.theta(2), M_PI / 8); - alpha[0] += 4; - alpha[1] -= 1; - P1 = grid.polarCoordinates(alpha); - ASSERT_DOUBLE_EQ(P1[0], 1.4); - ASSERT_DOUBLE_EQ(P1[1], M_PI / 16); + ASSERT_DOUBLE_EQ(grid.radius(7), 1.4); + ASSERT_DOUBLE_EQ(grid.theta(1), M_PI / 16); } -TEST(PolarGridTest, NeighborDistanceTest) +TEST(PolarGridTest, SpacingTest) { std::vector radii = {0.1, 0.2, 0.25, 0.5, 0.8, 0.9, 1.3, 1.4, 2.0}; std::vector 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}; - double splitting_radius = 0.6; - PolarGrid grid(radii, angles, splitting_radius); - - std::array, space_dimension> neighbor_distance; - { - MultiIndex alpha(3, 2); - grid.adjacentNeighborDistances(alpha, neighbor_distance); - ASSERT_DOUBLE_EQ(neighbor_distance[0].first, 0.5 - 0.25); - ASSERT_DOUBLE_EQ(neighbor_distance[0].second, 0.8 - 0.5); - ASSERT_DOUBLE_EQ(neighbor_distance[1].first, M_PI / 8 - M_PI / 16); - ASSERT_DOUBLE_EQ(neighbor_distance[1].second, M_PI / 2 - M_PI / 8); - } - - { - MultiIndex alpha(8, 7); - grid.adjacentNeighborDistances(alpha, neighbor_distance); - ASSERT_DOUBLE_EQ(neighbor_distance[0].first, 2.0 - 1.4); - ASSERT_DOUBLE_EQ(neighbor_distance[0].second, 0.0); - ASSERT_DOUBLE_EQ(neighbor_distance[1].first, M_PI / 2 - M_PI / 8); - ASSERT_DOUBLE_EQ(neighbor_distance[1].second, M_PI / 2); - } - - { - MultiIndex alpha(4, 0); - grid.adjacentNeighborDistances(alpha, neighbor_distance); - ASSERT_DOUBLE_EQ(neighbor_distance[0].first, 0.8 - 0.5); - ASSERT_DOUBLE_EQ(neighbor_distance[0].second, 0.9 - 0.8); - ASSERT_DOUBLE_EQ(neighbor_distance[1].first, M_PI / 2); - ASSERT_DOUBLE_EQ(neighbor_distance[1].second, M_PI / 16); - } - - { - MultiIndex alpha(0, 5); - grid.adjacentNeighborDistances(alpha, neighbor_distance); - ASSERT_DOUBLE_EQ(neighbor_distance[0].first, 0.0); - ASSERT_DOUBLE_EQ(neighbor_distance[0].second, 0.2 - 0.1); - ASSERT_NEAR(neighbor_distance[1].first, M_PI / 16, 1e-15); - ASSERT_NEAR(neighbor_distance[1].second, M_PI / 8 - M_PI / 16, 1e-15); - ASSERT_EQ(equals(neighbor_distance[1].first, M_PI / 16), true); - ASSERT_EQ(equals(neighbor_distance[1].second, M_PI / 8 - M_PI / 16), true); - } -} - -TEST(PolarGridTest, NeighborsTest) -{ - std::vector radii = {0.1, 0.2, 0.25, 0.5, 0.8, 0.9, 1.3, 1.4, 2.0}; - std::vector 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}; - double splitting_radius = 0.6; - PolarGrid grid(radii, angles, splitting_radius); - - std::array, space_dimension> neighbors; - - { - MultiIndex alpha(3, 2); - grid.adjacentNeighborsOf(alpha, neighbors); - ASSERT_EQ(neighbors[0].first, 18); - ASSERT_EQ(neighbors[0].second, 42); - ASSERT_EQ(neighbors[1].first, 25); - ASSERT_EQ(neighbors[1].second, 27); - grid.diagonalNeighborsOf(alpha, neighbors); - ASSERT_EQ(neighbors[0].first, 17); - ASSERT_EQ(neighbors[0].second, 37); - ASSERT_EQ(neighbors[1].first, 19); - ASSERT_EQ(neighbors[1].second, 47); - } - - { - MultiIndex alpha(8, 7); - grid.adjacentNeighborsOf(alpha, neighbors); - ASSERT_EQ(neighbors[0].first, 70); - ASSERT_EQ(neighbors[0].second, -1); - ASSERT_EQ(neighbors[1].first, 66); - ASSERT_EQ(neighbors[1].second, 36); - grid.diagonalNeighborsOf(alpha, neighbors); - ASSERT_EQ(neighbors[0].first, 65); - ASSERT_EQ(neighbors[0].second, -1); - ASSERT_EQ(neighbors[1].first, 35); - ASSERT_EQ(neighbors[1].second, -1); - } + PolarGrid grid(radii, angles); - { - MultiIndex alpha(4, 0); - grid.adjacentNeighborsOf(alpha, neighbors); - ASSERT_EQ(neighbors[0].first, 24); - ASSERT_EQ(neighbors[0].second, 33); - ASSERT_EQ(neighbors[1].first, 67); - ASSERT_EQ(neighbors[1].second, 37); - grid.diagonalNeighborsOf(alpha, neighbors); - ASSERT_EQ(neighbors[0].first, 31); - ASSERT_EQ(neighbors[0].second, 68); - ASSERT_EQ(neighbors[1].first, 25); - ASSERT_EQ(neighbors[1].second, 38); - } + // Test radial spacings + ASSERT_DOUBLE_EQ(grid.radialSpacing(0), 0.2 - 0.1); + ASSERT_DOUBLE_EQ(grid.radialSpacing(1), 0.25 - 0.2); + ASSERT_DOUBLE_EQ(grid.radialSpacing(2), 0.5 - 0.25); - { - MultiIndex alpha(0, 5); - grid.adjacentNeighborsOf(alpha, neighbors); - ASSERT_EQ(neighbors[0].first, -1); - ASSERT_EQ(neighbors[0].second, 13); - ASSERT_EQ(neighbors[1].first, 4); - ASSERT_EQ(neighbors[1].second, 6); - grid.diagonalNeighborsOf(alpha, neighbors); - ASSERT_EQ(neighbors[0].first, -1); - ASSERT_EQ(neighbors[0].second, 12); - ASSERT_EQ(neighbors[1].first, -1); - ASSERT_EQ(neighbors[1].second, 14); - } + // Test angular spacings + ASSERT_DOUBLE_EQ(grid.angularSpacing(0), M_PI / 16 - 0); + ASSERT_DOUBLE_EQ(grid.angularSpacing(1), M_PI / 8 - M_PI / 16); + ASSERT_DOUBLE_EQ(grid.angularSpacing(2), M_PI / 2 - M_PI / 8); } diff --git a/tests/Residual/residual.cpp b/tests/Residual/residual.cpp index 7522f8c3..08f44302 100644 --- a/tests/Residual/residual.cpp +++ b/tests/Residual/residual.cpp @@ -87,8 +87,9 @@ TEST(OperatorATest, applyA_DirBC_Interior) ASSERT_EQ(result_Give.size(), result_Take.size()); for (uint index = 0; index < result_Give.size(); index++) { - MultiIndex alpha = level.grid().multiIndex(index); - if (alpha[0] == 0 && !DirBC_Interior) + int i_r, i_theta; + level.grid().multiIndex(index, i_r, i_theta); + if (i_r == 0 && !DirBC_Interior) ASSERT_NEAR(result_Give[index], result_Take[index], 1e-8); else ASSERT_NEAR(result_Give[index], result_Take[index], 1e-11); @@ -143,8 +144,9 @@ TEST(OperatorATest, applyA_AcrossOrigin) ASSERT_EQ(result_Give.size(), result_Take.size()); for (uint index = 0; index < result_Give.size(); index++) { - MultiIndex alpha = level.grid().multiIndex(index); - if (alpha[0] == 0 && !DirBC_Interior) + int i_r, i_theta; + level.grid().multiIndex(index, i_r, i_theta); + if (i_r == 0 && !DirBC_Interior) ASSERT_NEAR(result_Give[index], result_Take[index], 1e-8); else ASSERT_NEAR(result_Give[index], result_Take[index], 1e-11); diff --git a/tests/Smoother/smoother.cpp b/tests/Smoother/smoother.cpp index 40ea3003..0392e887 100644 --- a/tests/Smoother/smoother.cpp +++ b/tests/Smoother/smoother.cpp @@ -94,8 +94,9 @@ TEST(SmootherTest, smoother_DirBC_Interior) ASSERT_EQ(solution_Give.size(), solution_Take.size()); for (uint index = 0; index < solution_Give.size(); index++) { - MultiIndex alpha = level.grid().multiIndex(index); - if (alpha[0] == 0 && !DirBC_Interior) + int i_r, i_theta; + level.grid().multiIndex(index, i_r, i_theta); + if (i_r == 0 && !DirBC_Interior) ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11); else ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11); @@ -153,8 +154,9 @@ TEST(SmootherTest, smoother_AcrossOrigin) ASSERT_EQ(solution_Give.size(), solution_Take.size()); for (uint index = 0; index < solution_Give.size(); index++) { - MultiIndex alpha = level.grid().multiIndex(index); - if (alpha[0] == 0 && !DirBC_Interior) + int i_r, i_theta; + level.grid().multiIndex(index, i_r, i_theta); + if (i_r == 0 && !DirBC_Interior) ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-8); else ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-10);