From 8317d86c444005707f811d10c53a437fc645d48e Mon Sep 17 00:00:00 2001 From: julianlitz Date: Fri, 27 Feb 2026 22:11:03 +0100 Subject: [PATCH 01/10] Working PCG --- include/ConfigParser/config_parser.h | 18 +- include/GMGPolar/igmgpolar.h | 43 +++++ include/GMGPolar/setup.h | 21 ++- include/GMGPolar/solver.h | 107 +++++------- include/Level/level.h | 5 +- scripts/compile.sh | 2 +- scripts/tutorial/run.sh | 34 +++- src/CMakeLists.txt | 2 +- src/ConfigParser/config_parser.cpp | 60 +++++++ src/GMGPolar/gmgpolar.cpp | 56 +++++++ src/GMGPolar/solver.cpp | 240 ++++++++++++++++++++++++--- src/Level/level.cpp | 11 +- src/main.cpp | 8 + 13 files changed, 501 insertions(+), 106 deletions(-) diff --git a/include/ConfigParser/config_parser.h b/include/ConfigParser/config_parser.h index 00f65a3c..60c18f89 100644 --- a/include/ConfigParser/config_parser.h +++ b/include/ConfigParser/config_parser.h @@ -42,10 +42,20 @@ class ConfigParser const PolarGrid& grid() const; - // Multigrid Parameters + // Full Multigrid Method bool FMG() const; int FMG_iterations() const; MultigridCycleType FMG_cycle() const; + + // Preconditioned Conjugate Gradient Method + bool PCG() const; + bool PCG_FMG() const; + int PCG_FMG_iterations() const; + MultigridCycleType PCG_FMG_cycle() const; + int PCG_MG_iterations() const; + MultigridCycleType PCG_MG_cycle() const; + + // Multigrid Parameters ExtrapolationType extrapolation() const; int maxLevels() const; int preSmoothingSteps() const; @@ -86,6 +96,12 @@ class ConfigParser bool FMG_; int FMG_iterations_; MultigridCycleType FMG_cycle_; + bool PCG_; + bool PCG_FMG_; + int PCG_FMG_iterations_; + MultigridCycleType PCG_FMG_cycle_; + int PCG_MG_iterations_; + MultigridCycleType PCG_MG_cycle_; // Iterative solver controls int max_iterations_; ResidualNormType residual_norm_type_; diff --git a/include/GMGPolar/igmgpolar.h b/include/GMGPolar/igmgpolar.h index 52c317ef..9770aacc 100644 --- a/include/GMGPolar/igmgpolar.h +++ b/include/GMGPolar/igmgpolar.h @@ -105,6 +105,20 @@ class IGMGPolar MultigridCycleType FMG_cycle() const; void FMG_cycle(MultigridCycleType FMG_cycle); + // Preconditioned Conjugate Gradient (PCG) control. + bool PCG() const; + void PCG(bool PCG); + bool PCG_FMG() const; + void PCG_FMG(bool PCG_FMG); + int PCG_FMG_iterations() const; + void PCG_FMG_iterations(int PCG_FMG_iterations); + MultigridCycleType PCG_FMG_cycle() const; + void PCG_FMG_cycle(MultigridCycleType PCG_FMG_cycle); + int PCG_MG_iterations() const; + void PCG_MG_iterations(int PCG_MG_iterations); + MultigridCycleType PCG_MG_cycle() const; + void PCG_MG_cycle(MultigridCycleType PCG_MG_cycle); + /* ---------------------------------------------------------------------- */ /* Iterative solver termination */ /* ---------------------------------------------------------------------- */ @@ -215,6 +229,13 @@ class IGMGPolar bool FMG_; int FMG_iterations_; MultigridCycleType FMG_cycle_; + // PCG settings + bool PCG_; + bool PCG_FMG_; + int PCG_FMG_iterations_; + MultigridCycleType PCG_FMG_cycle_; + int PCG_MG_iterations_; + MultigridCycleType PCG_MG_cycle_; // Convergence settings int max_iterations_; ResidualNormType residual_norm_type_; @@ -234,6 +255,22 @@ class IGMGPolar /* Chooses if full grid smoothing is active on level 0 for extrapolation > 0 */ bool full_grid_smoothing_ = false; + /* -------------------------------------------------- */ + /* Vectors for PCG (Preconditioned Conjugate Gradient) + * https://en.wikipedia.org/wiki/Conjugate_gradient_method#The_preconditioned_conjugate_gradient_method + * + * Dedicated vectors: + * x (solution) -> pcg_solution_ + * p (search direction) -> pcg_search_direction_ + * + * Reused vectors (to avoid extra allocations): + * r (residual) -> level.rhs() + * z (preconditioned residual) -> level.solution() + * A*p (matrix applied to search dir.) -> level.residual() + */ + AllocatableVector pcg_solution_; // x (solution) + AllocatableVector pcg_search_direction_; // p (search direction) + /* -------------------- */ /* Convergence criteria */ int number_of_iterations_; @@ -263,10 +300,16 @@ class IGMGPolar /* --------------- */ /* Solve Functions */ void fullMultigridApproximation(MultigridCycleType FMG_cycle, int FMG_iterations); + void solveMultigrid(double& initial_residual_norm, double& current_residual_norm, + double& current_relative_residual_norm); + void solvePCG(double& initial_residual_norm, double& current_residual_norm, double& current_relative_residual_norm); double residualNorm(const ResidualNormType& norm_type, const Level& level, ConstVector residual) const; void evaluateExactError(Level& level, const ExactSolution& exact_solution); void updateResidualNorms(Level& level, int iteration, double& initial_residual_norm, double& current_residual_norm, double& current_relative_residual_norm); + void initRhsHierarchy(Vector rhs); + void applyMultigridIterations(Level& level, MultigridCycleType cycle, int iterations); + void applyExtrapolatedMultigridIterations(Level& level, MultigridCycleType cycle, int iterations); /* ----------------- */ /* Print information */ diff --git a/include/GMGPolar/setup.h b/include/GMGPolar/setup.h index f85ffbf6..de07e019 100644 --- a/include/GMGPolar/setup.h +++ b/include/GMGPolar/setup.h @@ -56,14 +56,15 @@ void GMGPolar::setup() auto finest_levelCache = std::make_unique(*finest_grid, density_profile_coefficients_, domain_geometry_, cache_density_profile_coefficients_, cache_domain_geometry_); - levels_.emplace_back(0, std::move(finest_grid), std::move(finest_levelCache), extrapolation_, FMG_); + levels_.emplace_back(0, std::move(finest_grid), std::move(finest_levelCache), extrapolation_, FMG_, PCG_FMG_); for (int 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(*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_); + levels_.emplace_back(level_depth, std::move(current_grid), std::move(current_levelCache), extrapolation_, FMG_, + PCG_FMG_); } auto end_setup_createLevels = std::chrono::high_resolution_clock::now(); @@ -80,6 +81,19 @@ void GMGPolar::setup() if (verbose_ > 0) printSettings(); + // ------------------------------- // + // PCG-specific vector allocations // + // ------------------------------- // + if (PCG_) { + const int grid_size = levels_[0].grid().numberOfNodes(); + if (std::ssize(pcg_solution_) != grid_size) { + pcg_solution_ = Vector("pcg_solution", grid_size); + } + if (std::ssize(pcg_search_direction_) != grid_size) { + pcg_search_direction_ = Vector("pcg_search_direction", grid_size); + } + } + // ------------------------------------------------ // // Define residual operator on all multigrid levels // // ------------------------------------------------ // @@ -129,7 +143,8 @@ void GMGPolar::setup() full_grid_smoothing_ = do_full_grid_smoothing; if (number_of_levels_ > 1) { - if (do_full_grid_smoothing) { + // PCG uses non-extrapolated smoothing on level 0, so we need to initialize it if PCG is enabled. + if (do_full_grid_smoothing || (PCG_ && PCG_MG_iterations_ > 0)) { levels_[0].initializeSmoothing(domain_geometry_, density_profile_coefficients_, DirBC_Interior_, max_omp_threads_, stencil_distribution_method_); } diff --git a/include/GMGPolar/solver.h b/include/GMGPolar/solver.h index 66bcaee9..dcbe51a0 100644 --- a/include/GMGPolar/solver.h +++ b/include/GMGPolar/solver.h @@ -86,77 +86,46 @@ void IGMGPolar::solve(const BoundaryConditions& boundary_conditions, const Sourc printIterationHeader(exact_solution_); - while (number_of_iterations_ < max_iterations_) { - /* ---------------------------------------------- */ - /* Test solution against exact solution if given. */ - /* ---------------------------------------------- */ - LIKWID_STOP("Solver"); - auto start_check_exact_error = std::chrono::high_resolution_clock::now(); - - if (exact_solution_ != nullptr) - evaluateExactError(level, *exact_solution_); - - auto end_check_exact_error = std::chrono::high_resolution_clock::now(); - t_check_exact_error_ += std::chrono::duration(end_check_exact_error - start_check_exact_error).count(); - LIKWID_START("Solver"); - - /* ---------------------------- */ - /* Compute convergence criteria */ - /* ---------------------------- */ - auto start_check_convergence = std::chrono::high_resolution_clock::now(); - - if (absolute_tolerance_.has_value() || relative_tolerance_.has_value()) { - updateResidualNorms(level, number_of_iterations_, initial_residual_norm, current_residual_norm, - current_relative_residual_norm); - } + /* ---------------------------------------------- */ + /* Test solution against exact solution if given. */ + /* ---------------------------------------------- */ + LIKWID_STOP("Solver"); + auto start_check_exact_error = std::chrono::high_resolution_clock::now(); - auto end_check_convergence = std::chrono::high_resolution_clock::now(); - t_check_convergence_ += std::chrono::duration(end_check_convergence - start_check_convergence).count(); - - printIterationInfo(number_of_iterations_, current_residual_norm, current_relative_residual_norm, - exact_solution_); - - if (converged(current_residual_norm, current_relative_residual_norm)) - break; - - /* ----------------------- */ - /* Perform Multigrid Cycle */ - /* ----------------------- */ - auto start_solve_multigrid_iterations = std::chrono::high_resolution_clock::now(); - - switch (multigrid_cycle_) { - case MultigridCycleType::V_CYCLE: - if (extrapolation_ == ExtrapolationType::NONE) { - multigrid_V_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); - } - else { - extrapolated_multigrid_V_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); - } - break; - case MultigridCycleType::W_CYCLE: - if (extrapolation_ == ExtrapolationType::NONE) { - multigrid_W_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); - } - else { - extrapolated_multigrid_W_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); - } - break; - case MultigridCycleType::F_CYCLE: - if (extrapolation_ == ExtrapolationType::NONE) { - multigrid_F_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); - } - else { - extrapolated_multigrid_F_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); - } - break; - default: - throw std::invalid_argument("Unknown MultigridCycleType"); - } - number_of_iterations_++; + if (exact_solution_ != nullptr) + evaluateExactError(level, *exact_solution_); + + auto end_check_exact_error = std::chrono::high_resolution_clock::now(); + t_check_exact_error_ += std::chrono::duration(end_check_exact_error - start_check_exact_error).count(); + LIKWID_START("Solver"); + + /* ---------------------------- */ + /* Compute convergence criteria */ + /* ---------------------------- */ + auto start_check_convergence = std::chrono::high_resolution_clock::now(); - auto end_solve_multigrid_iterations = std::chrono::high_resolution_clock::now(); - t_solve_multigrid_iterations_ += - std::chrono::duration(end_solve_multigrid_iterations - start_solve_multigrid_iterations).count(); + if (absolute_tolerance_.has_value() || relative_tolerance_.has_value()) { + updateResidualNorms(level, number_of_iterations_, initial_residual_norm, current_residual_norm, + current_relative_residual_norm); + } + + auto end_check_convergence = std::chrono::high_resolution_clock::now(); + t_check_convergence_ += std::chrono::duration(end_check_convergence - start_check_convergence).count(); + + printIterationInfo(number_of_iterations_, current_residual_norm, current_relative_residual_norm, exact_solution_); + + if (!converged(current_residual_norm, current_relative_residual_norm)) { + if (!PCG_) { + // Solve A*x = b directly using multigrid cycles (V/W/F-cycle) + // until convergence or max_iterations_ is reached. + solveMultigrid(initial_residual_norm, current_residual_norm, current_relative_residual_norm); + } + else { + // Solve A*x = b using Preconditioned Conjugate Gradient (PCG), + // with multigrid cycles as the preconditioner (i.e., one multigrid + // cycle approximates the action of A^{-1} at each PCG iteration). + solvePCG(initial_residual_norm, current_residual_norm, current_relative_residual_norm); + } } /* ---------------------- */ diff --git a/include/Level/level.h b/include/Level/level.h index 0a83b52c..657bea64 100644 --- a/include/Level/level.h +++ b/include/Level/level.h @@ -43,8 +43,8 @@ class Level // ----------- // // Constructor // explicit Level(const int level_depth, std::unique_ptr grid, - std::unique_ptr level_cache, const ExtrapolationType extrapolation, - const bool FMG); + std::unique_ptr level_cache, const ExtrapolationType extrapolation, const bool FMG, + const bool PCG_FMG = false); // ---------------- // // Getter Functions // @@ -67,6 +67,7 @@ class Level const DensityProfileCoefficients& density_profile_coefficients, const bool DirBC_Interior, const int num_omp_threads, const StencilDistributionMethod stencil_distribution_method); void computeResidual(Vector result, ConstVector rhs, ConstVector x) const; + void applySystemOperator(Vector result, ConstVector x) const; // ------------------- // // Solve coarse System // diff --git a/scripts/compile.sh b/scripts/compile.sh index 77122d2e..80dac594 100755 --- a/scripts/compile.sh +++ b/scripts/compile.sh @@ -59,4 +59,4 @@ if [ -n "$build_type" ]; then -DCMAKE_BUILD_TYPE="$build_type" .. || { echo "CMake configuration failed"; exit 1; } fi -cmake --build ${PWD}/../build -j 16 +cmake --build ${PWD}/../build -j 2 diff --git a/scripts/tutorial/run.sh b/scripts/tutorial/run.sh index 59310dfb..a282e08f 100755 --- a/scripts/tutorial/run.sh +++ b/scripts/tutorial/run.sh @@ -16,17 +16,17 @@ paraview=0 # OpenMP settings: # Maximum number of threads OpenMP can use for parallel execution -maxOpenMPThreads=32 +maxOpenMPThreads=1 # Stencil distribution method: # 0 - CPU "Take": Each node independently applies the stencil # 1 - CPU "Give": The stencil operation is distributed across adjacent neighboring nodes -stencilDistributionMethod=1 +stencilDistributionMethod=0 # Caching behavior: # 0 - Recompute values on each iteration: Uses less memory but results in slower execution. # 1 - Reuse cached values: Consumes more memory but significantly improves performance. cacheDensityProfileCoefficients=1 -cacheDomainGeometry=0 +cacheDomainGeometry=1 # Note: In the "Take" approach (stencilDistributionMethod=0), # caching is required for optimal performance, # so both density profile coefficients and domain geometry need to be cached. @@ -62,10 +62,24 @@ beta_coeff=1 # Zero(0), Gyro - Alpha Inverse(1) # Full Multigrid Method: # 0: Initial approximation is set to zero # 1: Initial approximation obtained by nested iteration (recommended) -FMG=0 -FMG_iterations=3 +FMG=1 +FMG_iterations=2 FMG_cycle=2 # V-Cycle(0), W-Cycle(1), F-Cycle(2) +# Preconditioned Conjugate Gradient Method: +# 0: GMGPolar as iterative solver +# 1: GMGPolar solver as preconditioner for Conjugate Gradient (recommended) +PCG=1 +# Initial approximation for PCG: +# 0: Initial approximation is set to residual -> no preconditioning +# 1: FMG-approximation as initial guess (recommended) +PCG_FMG=1 +PCG_FMG_iterations=1 +PCG_FMG_cycle=0 # V-Cycle(0), W-Cycle(1), F-Cycle(2) +# Additional multigrid iterations after initial approximation to solve the linear system in each PCG iteration +PCG_MG_iterations=1 +PCG_MG_cycle=0 # V-Cycle(0), W-Cycle(1), F-Cycle(2) + # Extrapolation Method: # 0: No extrapolation # 1: Implicit extrapolation (recommended) @@ -86,8 +100,8 @@ multigridCycle=0 # Convergence criteria: maxIterations=150 residualNormType=0 # L2-Norm(0) = 0, Weighted L2-Norm(1), Infinity-Norm(2) -absoluteTolerance=1e-8 -relativeTolerance=1e-8 +absoluteTolerance=1e-12 +relativeTolerance=1e-12 # Define additional geometry parameters kappa_eps=0.0 @@ -141,6 +155,12 @@ export OMP_NUM_THREADS=$maxOpenMPThreads --FMG $FMG \ --FMG_iterations $FMG_iterations \ --FMG_cycle $FMG_cycle \ + --PCG $PCG \ + --PCG_FMG $PCG_FMG \ + --PCG_FMG_iterations $PCG_FMG_iterations \ + --PCG_FMG_cycle $PCG_FMG_cycle \ + --PCG_MG_iterations $PCG_MG_iterations \ + --PCG_MG_cycle $PCG_MG_cycle \ --extrapolation $extrapolation \ --maxLevels $maxLevels \ --preSmoothingSteps $preSmoothingSteps \ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1427f30b..b950d960 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -195,7 +195,7 @@ find_package(OpenMP REQUIRED) if(OpenMP_CXX_FOUND) target_link_libraries(GMGPolarLib PUBLIC OpenMP::OpenMP_CXX) endif() -find_package(Kokkos 4.4.1...<5 QUIET REQUIRED) +find_package(Kokkos 4.4.1...<5.1 QUIET REQUIRED) target_link_libraries(GMGPolarLib PUBLIC Kokkos::kokkos) diff --git a/src/ConfigParser/config_parser.cpp b/src/ConfigParser/config_parser.cpp index e4e51551..63b2b8bd 100644 --- a/src/ConfigParser/config_parser.cpp +++ b/src/ConfigParser/config_parser.cpp @@ -28,10 +28,21 @@ ConfigParser::ConfigParser() parser_.add("maxLevels", 'l', "Max multigrid levels.", OPTIONAL, -1); parser_.add("preSmoothingSteps", '\0', "Pre-smoothing steps.", OPTIONAL, 1); parser_.add("postSmoothingSteps", '\0', "Post-smoothing steps.", OPTIONAL, 1); + parser_.add("multigridCycle", '\0', "Cycle type (0=V,1=W,2=F).", OPTIONAL, 0, cmdline::oneof(0, 1, 2)); parser_.add("FMG", '\0', "Use Full Multigrid (0/1).", OPTIONAL, 0, cmdline::oneof(0, 1)); parser_.add("FMG_iterations", '\0', "FMG iterations.", OPTIONAL, 2); parser_.add("FMG_cycle", '\0', "FMG cycle type (0=V,1=W,2=F).", OPTIONAL, 0, cmdline::oneof(0, 1, 2)); + + parser_.add("PCG", '\0', "Use Preconditioned Conjugate Gradient (0/1).", OPTIONAL, 0, cmdline::oneof(0, 1)); + parser_.add("PCG_FMG", '\0', "Use FMG as preconditioner for PCG (0/1).", OPTIONAL, 0, cmdline::oneof(0, 1)); + parser_.add("PCG_FMG_iterations", '\0', "FMG iterations for PCG preconditioner.", OPTIONAL, 2); + parser_.add("PCG_FMG_cycle", '\0', "FMG cycle type for PCG preconditioner (0=V,1=W,2=F).", OPTIONAL, 0, + cmdline::oneof(0, 1, 2)); + parser_.add("PCG_MG_iterations", '\0', "Multigrid iterations for PCG preconditioner.", OPTIONAL, 1); + parser_.add("PCG_MG_cycle", '\0', "Multigrid cycle type for PCG preconditioner (0=V,1=W,2=F).", OPTIONAL, 0, + cmdline::oneof(0, 1, 2)); + parser_.add("maxIterations", '\0', "Max solver iterations.", OPTIONAL, 150); parser_.add("residualNormType", '\0', "Residual norm (0=Euclidean,1=Weighted,2=Infinity)", OPTIONAL, 0, cmdline::oneof(0, 1, 2)); @@ -106,6 +117,30 @@ bool ConfigParser::parse(int argc, char* argv[]) else { throw std::runtime_error("Invalid FMG cycle type."); } + + PCG_ = parser_.get("PCG") != 0; + PCG_FMG_ = parser_.get("PCG_FMG") != 0; + PCG_FMG_iterations_ = parser_.get("PCG_FMG_iterations"); + const int PCG_FMG_cycleValue = parser_.get("PCG_FMG_cycle"); + if (PCG_FMG_cycleValue == static_cast(MultigridCycleType::V_CYCLE) || + PCG_FMG_cycleValue == static_cast(MultigridCycleType::W_CYCLE) || + PCG_FMG_cycleValue == static_cast(MultigridCycleType::F_CYCLE)) { + PCG_FMG_cycle_ = static_cast(PCG_FMG_cycleValue); + } + else { + throw std::runtime_error("Invalid PCG FMG cycle type."); + } + PCG_MG_iterations_ = parser_.get("PCG_MG_iterations"); + const int PCG_MG_cycleValue = parser_.get("PCG_MG_cycle"); + if (PCG_MG_cycleValue == static_cast(MultigridCycleType::V_CYCLE) || + PCG_MG_cycleValue == static_cast(MultigridCycleType::W_CYCLE) || + PCG_MG_cycleValue == static_cast(MultigridCycleType::F_CYCLE)) { + PCG_MG_cycle_ = static_cast(PCG_MG_cycleValue); + } + else { + throw std::runtime_error("Invalid PCG MG cycle type."); + } + const int extrapolationValue = parser_.get("extrapolation"); if (extrapolationValue == static_cast(ExtrapolationType::NONE) || extrapolationValue == static_cast(ExtrapolationType::IMPLICIT_EXTRAPOLATION) || @@ -330,6 +365,31 @@ MultigridCycleType ConfigParser::FMG_cycle() const { return FMG_cycle_; } + +bool ConfigParser::PCG() const +{ + return PCG_; +} +bool ConfigParser::PCG_FMG() const +{ + return PCG_FMG_; +} +int ConfigParser::PCG_FMG_iterations() const +{ + return PCG_FMG_iterations_; +} +MultigridCycleType ConfigParser::PCG_FMG_cycle() const +{ + return PCG_FMG_cycle_; +} +int ConfigParser::PCG_MG_iterations() const +{ + return PCG_MG_iterations_; +} +MultigridCycleType ConfigParser::PCG_MG_cycle() const +{ + return PCG_MG_cycle_; +} ExtrapolationType ConfigParser::extrapolation() const { return extrapolation_; diff --git a/src/GMGPolar/gmgpolar.cpp b/src/GMGPolar/gmgpolar.cpp index ee630f70..037cf23b 100644 --- a/src/GMGPolar/gmgpolar.cpp +++ b/src/GMGPolar/gmgpolar.cpp @@ -26,6 +26,13 @@ IGMGPolar::IGMGPolar(const PolarGrid& grid) , FMG_(false) , FMG_iterations_(3) , FMG_cycle_(MultigridCycleType::F_CYCLE) + // PCG settings + , PCG_(false) + , PCG_FMG_(true) + , PCG_FMG_iterations_(1) + , PCG_FMG_cycle_(MultigridCycleType::V_CYCLE) + , PCG_MG_iterations_(1) + , PCG_MG_cycle_(MultigridCycleType::V_CYCLE) // Convergence settings , max_iterations_(300) , residual_norm_type_(ResidualNormType::WEIGHTED_EUCLIDEAN) @@ -195,6 +202,55 @@ void IGMGPolar::FMG_cycle(MultigridCycleType FMG_cycle) FMG_cycle_ = FMG_cycle; } +bool IGMGPolar::PCG() const +{ + return PCG_; +} +void IGMGPolar::PCG(bool PCG) +{ + PCG_ = PCG; +} +bool IGMGPolar::PCG_FMG() const +{ + return PCG_FMG_; +} +void IGMGPolar::PCG_FMG(bool PCG_FMG) +{ + PCG_FMG_ = PCG_FMG; +} +int IGMGPolar::PCG_FMG_iterations() const +{ + return PCG_FMG_iterations_; +} +void IGMGPolar::PCG_FMG_iterations(int PCG_FMG_iterations) +{ + PCG_FMG_iterations_ = PCG_FMG_iterations; +} +MultigridCycleType IGMGPolar::PCG_FMG_cycle() const +{ + return PCG_FMG_cycle_; +} +void IGMGPolar::PCG_FMG_cycle(MultigridCycleType PCG_FMG_cycle) +{ + PCG_FMG_cycle_ = PCG_FMG_cycle; +} +int IGMGPolar::PCG_MG_iterations() const +{ + return PCG_MG_iterations_; +} +void IGMGPolar::PCG_MG_iterations(int PCG_MG_iterations) +{ + PCG_MG_iterations_ = PCG_MG_iterations; +} +MultigridCycleType IGMGPolar::PCG_MG_cycle() const +{ + return PCG_MG_cycle_; +} +void IGMGPolar::PCG_MG_cycle(MultigridCycleType PCG_MG_cycle) +{ + PCG_MG_cycle_ = PCG_MG_cycle; +} + /* ---------------------------------------------------------------------- */ /* Iterative solver termination */ /* ---------------------------------------------------------------------- */ diff --git a/src/GMGPolar/solver.cpp b/src/GMGPolar/solver.cpp index b56ae2d5..3604ad61 100644 --- a/src/GMGPolar/solver.cpp +++ b/src/GMGPolar/solver.cpp @@ -25,26 +25,7 @@ void IGMGPolar::fullMultigridApproximation(MultigridCycleType FMG_cycle, int FMG // Apply some FMG iterations, except on the finest level, // where the interpolated solution is sufficintly accurate as an initial guess if (fine_level.level_depth() > 0) { - for (int i = 0; i < FMG_iterations; i++) { - switch (FMG_cycle) { - case MultigridCycleType::V_CYCLE: - multigrid_V_Cycle(fine_level.level_depth(), fine_level.solution(), fine_level.rhs(), - fine_level.residual()); - break; - case MultigridCycleType::W_CYCLE: - multigrid_W_Cycle(fine_level.level_depth(), fine_level.solution(), fine_level.rhs(), - fine_level.residual()); - break; - case MultigridCycleType::F_CYCLE: - multigrid_F_Cycle(fine_level.level_depth(), fine_level.solution(), fine_level.rhs(), - fine_level.residual()); - break; - default: - std::cerr << "Error: Unknown multigrid cycle type!" << std::endl; - throw std::runtime_error("Invalid multigrid cycle type encountered."); - break; - } - } + applyMultigridIterations(fine_level, FMG_cycle, FMG_iterations); } } } @@ -265,3 +246,222 @@ void IGMGPolar::printIterationInfo(int iteration, double current_residual_norm, std::cout << "\n"; std::cout << std::right << std::defaultfloat << std::setprecision(6) << std::setfill(' '); } + +void IGMGPolar::solveMultigrid(double& initial_residual_norm, double& current_residual_norm, + double& current_relative_residual_norm) +{ + Level& level = levels_[0]; + + while (number_of_iterations_ < max_iterations_) { + /* ----------------------- */ + /* Perform Multigrid Cycle */ + /* ----------------------- */ + auto start_solve_multigrid_iterations = std::chrono::high_resolution_clock::now(); + + if (extrapolation_ == ExtrapolationType::NONE) { + applyMultigridIterations(level, multigrid_cycle_, 1); + } + else { + applyExtrapolatedMultigridIterations(level, multigrid_cycle_, 1); + } + + number_of_iterations_++; + + auto end_solve_multigrid_iterations = std::chrono::high_resolution_clock::now(); + t_solve_multigrid_iterations_ += + std::chrono::duration(end_solve_multigrid_iterations - start_solve_multigrid_iterations).count(); + + /* ---------------------------------------------- */ + /* Test solution against exact solution if given. */ + /* ---------------------------------------------- */ + LIKWID_STOP("Solver"); + auto start_check_exact_error = std::chrono::high_resolution_clock::now(); + + if (exact_solution_ != nullptr) + evaluateExactError(level, *exact_solution_); + + auto end_check_exact_error = std::chrono::high_resolution_clock::now(); + t_check_exact_error_ += std::chrono::duration(end_check_exact_error - start_check_exact_error).count(); + LIKWID_START("Solver"); + + /* ---------------------------- */ + /* Compute convergence criteria */ + /* ---------------------------- */ + auto start_check_convergence = std::chrono::high_resolution_clock::now(); + + if (absolute_tolerance_.has_value() || relative_tolerance_.has_value()) { + updateResidualNorms(level, number_of_iterations_, initial_residual_norm, current_residual_norm, + current_relative_residual_norm); + } + + auto end_check_convergence = std::chrono::high_resolution_clock::now(); + t_check_convergence_ += std::chrono::duration(end_check_convergence - start_check_convergence).count(); + + printIterationInfo(number_of_iterations_, current_residual_norm, current_relative_residual_norm, + exact_solution_); + + if (converged(current_residual_norm, current_relative_residual_norm)) + break; + } +} + +/* -------------------------------------------------- */ +/* Vectors for PCG (Preconditioned Conjugate Gradient) + * https://en.wikipedia.org/wiki/Conjugate_gradient_method#The_preconditioned_conjugate_gradient_method + * + * Dedicated vectors: + * x (solution) -> pcg_solution_ + * p (search direction) -> pcg_search_direction_ + * + * Reused vectors (to avoid extra allocations): + * r (residual) -> level.rhs() + * z (preconditioned residual) -> level.solution() + * A*p (matrix applied to search dir.) -> level.residual() + */ +void IGMGPolar::solvePCG(double& initial_residual_norm, double& current_residual_norm, + double& current_relative_residual_norm) +{ + Level& level = levels_[0]; + + // x = initial guess + Kokkos::deep_copy(pcg_solution_, level.solution()); + // r = residual of initial guess + Kokkos::deep_copy(level.rhs(), level.residual()); + + // z = M^{-1} * r (preconditioned residual) + if (PCG_FMG_) { + initRhsHierarchy(level.rhs()); + fullMultigridApproximation(PCG_FMG_cycle_, PCG_FMG_iterations_); + applyMultigridIterations(level, PCG_MG_cycle_, PCG_MG_iterations_); + } + else { + // z = I^{-1} * r (no preconditioning) + Kokkos::deep_copy(level.solution(), level.rhs()); + } + + // p = z + Kokkos::deep_copy(pcg_search_direction_, level.solution()); + + // r^T * z + double r_z = dot_product(ConstVector(level.rhs()), ConstVector(level.solution())); + + while (number_of_iterations_ < max_iterations_) { + + // A_p = A * p + level.applySystemOperator(level.residual(), pcg_search_direction_); + if (extrapolation_ != ExtrapolationType::NONE) { + assert(number_of_levels_ > 1); + Level& next_level = levels_[level.level_depth() + 1]; + injection(0, next_level.solution(), pcg_search_direction_); + next_level.applySystemOperator(next_level.residual(), next_level.solution()); + extrapolatedResidual(0, level.residual(), next_level.residual()); + } + + // alpha = (r^T * z) / (p^T * A*p) + double alpha = + r_z / dot_product(ConstVector(pcg_search_direction_), ConstVector(level.residual())); + + // x += alpha * p + linear_combination(pcg_solution_, 1.0, ConstVector(pcg_search_direction_), alpha); + + // r -= alpha * A*p + linear_combination(level.rhs(), 1.0, ConstVector(level.residual()), -alpha); + + current_residual_norm = residualNorm(residual_norm_type_, level, level.rhs()); + residual_norms_.push_back(current_residual_norm); + current_relative_residual_norm = current_residual_norm / initial_residual_norm; + + // --- Check exact error, excluded from timings + auto start_check_exact_error = std::chrono::high_resolution_clock::now(); + if (exact_solution_ != nullptr) + exact_errors_.push_back(computeExactError(level, pcg_solution_, level.residual(), *exact_solution_)); + auto end_check_exact_error = std::chrono::high_resolution_clock::now(); + t_check_exact_error_ += std::chrono::duration(end_check_exact_error - start_check_exact_error).count(); + + number_of_iterations_++; + + printIterationInfo(number_of_iterations_, current_residual_norm, current_relative_residual_norm, + exact_solution_); + + if (converged(current_residual_norm, current_relative_residual_norm)) + break; + + // z = M^{-1} * r (preconditioned residual) + if (PCG_FMG_) { + initRhsHierarchy(level.rhs()); + fullMultigridApproximation(PCG_FMG_cycle_, PCG_FMG_iterations_); + applyMultigridIterations(level, PCG_MG_cycle_, PCG_MG_iterations_); + } + else { + // z = I^{-1} * r (no preconditioning) + Kokkos::deep_copy(level.solution(), level.rhs()); + } + + // r^T * z + double r_z_new = dot_product(ConstVector(level.rhs()), ConstVector(level.solution())); + // beta = (current r^T * z) / (previous r^T * z) + double beta = r_z_new / r_z; + // r_z = r^T * z for next iteration + r_z = r_z_new; + + // p *= beta + multiply(pcg_search_direction_, beta); + // p += z + add(pcg_search_direction_, ConstVector(level.solution())); + } + // level.solution() = x + Kokkos::deep_copy(level.solution(), pcg_solution_); +} + +void IGMGPolar::initRhsHierarchy(Vector rhs) +{ + Level& level = levels_[0]; + Kokkos::deep_copy(level.rhs(), rhs); + for (int level_depth = 0; level_depth < number_of_levels_; ++level_depth) { + Level& current_level = levels_[level_depth]; + if (level_depth + 1 < number_of_levels_) { + Level& next_level = levels_[level_depth + 1]; + restriction(level_depth, next_level.rhs(), current_level.rhs()); + } + } +} + +void IGMGPolar::applyMultigridIterations(Level& level, MultigridCycleType cycle, int iterations) +{ + for (int i = 0; i < iterations; i++) { + switch (cycle) { + case MultigridCycleType::V_CYCLE: + multigrid_V_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); + break; + case MultigridCycleType::W_CYCLE: + multigrid_W_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); + break; + case MultigridCycleType::F_CYCLE: + multigrid_F_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); + break; + default: + std::cerr << "Error: Unknown multigrid cycle type!" << std::endl; + break; + } + } +} + +void IGMGPolar::applyExtrapolatedMultigridIterations(Level& level, MultigridCycleType cycle, int iterations) +{ + for (int i = 0; i < iterations; i++) { + switch (cycle) { + case MultigridCycleType::V_CYCLE: + extrapolated_multigrid_V_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); + break; + case MultigridCycleType::W_CYCLE: + extrapolated_multigrid_W_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); + break; + case MultigridCycleType::F_CYCLE: + extrapolated_multigrid_F_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual()); + break; + default: + std::cerr << "Error: Unknown multigrid cycle type!" << std::endl; + break; + } + } +} \ No newline at end of file diff --git a/src/Level/level.cpp b/src/Level/level.cpp index e5e4dac8..389b785e 100644 --- a/src/Level/level.cpp +++ b/src/Level/level.cpp @@ -17,11 +17,12 @@ // ----------- // // Constructor // Level::Level(const int level_depth, std::unique_ptr grid, - std::unique_ptr level_cache, const ExtrapolationType extrapolation, const bool FMG) + std::unique_ptr level_cache, const ExtrapolationType extrapolation, const bool FMG, + const bool PCG_FMG) : level_depth_(level_depth) , grid_(std::move(grid)) , level_cache_(std::move(level_cache)) - , rhs_("rhs", (FMG || level_depth == 0 || (level_depth == 1 && extrapolation != ExtrapolationType::NONE)) + , rhs_("rhs", (FMG || PCG_FMG || level_depth == 0 || (level_depth == 1 && extrapolation != ExtrapolationType::NONE)) ? grid_->numberOfNodes() : 0) , solution_("solution", grid_->numberOfNodes()) @@ -104,6 +105,12 @@ void Level::computeResidual(Vector result, ConstVector rhs, Cons throw std::runtime_error("Residual not initialized."); op_residual_->computeResidual(result, rhs, x); } +void Level::applySystemOperator(Vector result, ConstVector x) const +{ + if (!op_residual_) + throw std::runtime_error("Residual not initialized."); + op_residual_->applySystemOperator(result, x); +} // ------------------- // // Solve coarse System // diff --git a/src/main.cpp b/src/main.cpp index 36b3fbaa..b1a573ba 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -41,6 +41,14 @@ int main(int argc, char* argv[]) solver->FMG_iterations(parser.FMG_iterations()); // FMG iteration count solver->FMG_cycle(parser.FMG_cycle()); // FMG cycle type + // --- Preconditioned Conjugate Gradient settings --- // + solver->PCG(parser.PCG()); // Preconditioned Conjugate Gradient mode on/off + solver->PCG_FMG(parser.PCG_FMG()); // Use FMG as preconditioner for PCG + solver->PCG_FMG_iterations(parser.PCG_FMG_iterations()); // FMG iterations for PCG preconditioner + solver->PCG_FMG_cycle(parser.PCG_FMG_cycle()); // FMG cycle type for PCG preconditioner + solver->PCG_MG_iterations(parser.PCG_MG_iterations()); // Multigrid iterations for PCG preconditioner + solver->PCG_MG_cycle(parser.PCG_MG_cycle()); // Multigrid cycle type for PCG preconditioner + // --- Iterative solver controls --- // solver->maxIterations(parser.maxIterations()); // Max number of iterations solver->residualNormType(parser.residualNormType()); // Residual norm type (L2, weighted-L2, L∞) From 6d53e03954d9e46b119375a5ca8e52a4ebf6e925 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Fri, 27 Feb 2026 22:24:01 +0100 Subject: [PATCH 02/10] Add timings --- include/GMGPolar/igmgpolar.h | 1 + include/GMGPolar/solver.h | 6 ++++++ scripts/tutorial/run.sh | 4 ++-- src/GMGPolar/gmgpolar.cpp | 32 +++++++++++++++++++++----------- src/GMGPolar/solver.cpp | 18 ++++++++++++++++-- 5 files changed, 46 insertions(+), 15 deletions(-) diff --git a/include/GMGPolar/igmgpolar.h b/include/GMGPolar/igmgpolar.h index 9770aacc..3d06bec8 100644 --- a/include/GMGPolar/igmgpolar.h +++ b/include/GMGPolar/igmgpolar.h @@ -362,6 +362,7 @@ class IGMGPolar double t_solve_multigrid_iterations_; double t_check_convergence_; double t_check_exact_error_; + double t_conjugate_gradient_; void resetAvgMultigridCycleTimings(); double t_avg_MGC_total_; diff --git a/include/GMGPolar/solver.h b/include/GMGPolar/solver.h index dcbe51a0..2b750ff2 100644 --- a/include/GMGPolar/solver.h +++ b/include/GMGPolar/solver.h @@ -124,7 +124,13 @@ void IGMGPolar::solve(const BoundaryConditions& boundary_conditions, const Sourc // Solve A*x = b using Preconditioned Conjugate Gradient (PCG), // with multigrid cycles as the preconditioner (i.e., one multigrid // cycle approximates the action of A^{-1} at each PCG iteration). + auto start_conjugate_gradient = std::chrono::high_resolution_clock::now(); + solvePCG(initial_residual_norm, current_residual_norm, current_relative_residual_norm); + + auto end_conjugate_gradient = std::chrono::high_resolution_clock::now(); + t_conjugate_gradient_ += + std::chrono::duration(end_conjugate_gradient - start_conjugate_gradient).count(); } } diff --git a/scripts/tutorial/run.sh b/scripts/tutorial/run.sh index a282e08f..45271a30 100755 --- a/scripts/tutorial/run.sh +++ b/scripts/tutorial/run.sh @@ -100,8 +100,8 @@ multigridCycle=0 # Convergence criteria: maxIterations=150 residualNormType=0 # L2-Norm(0) = 0, Weighted L2-Norm(1), Infinity-Norm(2) -absoluteTolerance=1e-12 -relativeTolerance=1e-12 +absoluteTolerance=1e-10 +relativeTolerance=1e-10 # Define additional geometry parameters kappa_eps=0.0 diff --git a/src/GMGPolar/gmgpolar.cpp b/src/GMGPolar/gmgpolar.cpp index 037cf23b..58e791d8 100644 --- a/src/GMGPolar/gmgpolar.cpp +++ b/src/GMGPolar/gmgpolar.cpp @@ -409,6 +409,7 @@ void IGMGPolar::resetSolvePhaseTimings() t_solve_multigrid_iterations_ = 0.0; t_check_convergence_ = 0.0; t_check_exact_error_ = 0.0; + t_conjugate_gradient_ = 0.0; } void IGMGPolar::resetAvgMultigridCycleTimings() @@ -437,19 +438,28 @@ void IGMGPolar::printTimings() const std::cout << " (Build rhs: " << t_setup_rhs_ << " seconds)" << std::endl; std::cout << "\nSolve Time: " << t_solve_total_ << " seconds" << std::endl; std::cout << " Initial Approximation: " << t_solve_initial_approximation_ << " seconds" << std::endl; - std::cout << " Multigrid Iteration: " << t_solve_multigrid_iterations_ << " seconds" << std::endl; + if (!PCG_) { + std::cout << " Multigrid Iterations: " << t_solve_multigrid_iterations_ << " seconds" << std::endl; + } + else { + std::cout << " Preconditioned Conjugate Gradient: " + << t_conjugate_gradient_ - t_check_convergence_ - t_check_exact_error_ << " seconds" << std::endl; + } std::cout << " Check Convergence: " << t_check_convergence_ << " seconds" << std::endl; std::cout << " (Check Exact Error: " << t_check_exact_error_ << " seconds)" << std::endl; - std::cout << "\nAverage Multigrid Iteration: " << t_avg_MGC_total_ << " seconds" << std::endl; - std::cout << " PreSmoothing: " << t_avg_MGC_preSmoothing_ << " seconds" << std::endl; - std::cout << " PostSmoothing: " << t_avg_MGC_postSmoothing_ << " seconds" << std::endl; - std::cout << " Residual: " << t_avg_MGC_residual_ << " seconds" << std::endl; - std::cout << " DirectSolve: " << t_avg_MGC_directSolver_ << " seconds" << std::endl; - std::cout << " Other Computations: " - << std::max(t_avg_MGC_total_ - t_avg_MGC_preSmoothing_ - t_avg_MGC_postSmoothing_ - t_avg_MGC_residual_ - - t_avg_MGC_directSolver_, - 0.0) - << " seconds" << std::endl; + + if (!PCG_) { + std::cout << "\nAverage Multigrid Iteration: " << t_avg_MGC_total_ << " seconds" << std::endl; + std::cout << " PreSmoothing: " << t_avg_MGC_preSmoothing_ << " seconds" << std::endl; + std::cout << " PostSmoothing: " << t_avg_MGC_postSmoothing_ << " seconds" << std::endl; + std::cout << " Residual: " << t_avg_MGC_residual_ << " seconds" << std::endl; + std::cout << " DirectSolve: " << t_avg_MGC_directSolver_ << " seconds" << std::endl; + std::cout << " Other Computations: " + << std::max(t_avg_MGC_total_ - t_avg_MGC_preSmoothing_ - t_avg_MGC_postSmoothing_ - + t_avg_MGC_residual_ - t_avg_MGC_directSolver_, + 0.0) + << " seconds" << std::endl; + } } // Number of iterations taken by last solve. diff --git a/src/GMGPolar/solver.cpp b/src/GMGPolar/solver.cpp index 3604ad61..58107d13 100644 --- a/src/GMGPolar/solver.cpp +++ b/src/GMGPolar/solver.cpp @@ -367,16 +367,30 @@ void IGMGPolar::solvePCG(double& initial_residual_norm, double& current_residual // r -= alpha * A*p linear_combination(level.rhs(), 1.0, ConstVector(level.residual()), -alpha); + /* ---------------------------- */ + /* Compute convergence criteria */ + /* ---------------------------- */ + auto start_check_convergence = std::chrono::high_resolution_clock::now(); + current_residual_norm = residualNorm(residual_norm_type_, level, level.rhs()); residual_norms_.push_back(current_residual_norm); current_relative_residual_norm = current_residual_norm / initial_residual_norm; - // --- Check exact error, excluded from timings + auto end_check_convergence = std::chrono::high_resolution_clock::now(); + t_check_convergence_ += std::chrono::duration(end_check_convergence - start_check_convergence).count(); + + /* ---------------------------------------------- */ + /* Test solution against exact solution if given. */ + /* ---------------------------------------------- */ + LIKWID_STOP("Solver"); auto start_check_exact_error = std::chrono::high_resolution_clock::now(); + if (exact_solution_ != nullptr) - exact_errors_.push_back(computeExactError(level, pcg_solution_, level.residual(), *exact_solution_)); + evaluateExactError(level, *exact_solution_); + auto end_check_exact_error = std::chrono::high_resolution_clock::now(); t_check_exact_error_ += std::chrono::duration(end_check_exact_error - start_check_exact_error).count(); + LIKWID_START("Solver"); number_of_iterations_++; From 8de99382fc84b8fecbcdb46d4d68fd4fe2b727c1 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Fri, 27 Feb 2026 23:36:46 +0100 Subject: [PATCH 03/10] Add PCG google test --- include/GMGPolar/igmgpolar.h | 1 + include/GMGPolar/setup.h | 3 +- src/GMGPolar/gmgpolar.cpp | 5 + src/GMGPolar/setup.cpp | 45 ++ src/GMGPolar/solver.cpp | 2 +- tests/CMakeLists.txt | 1 + tests/GMGPolar/pcg_tests.cpp | 783 +++++++++++++++++++++++++++++++++++ 7 files changed, 838 insertions(+), 2 deletions(-) create mode 100644 tests/GMGPolar/pcg_tests.cpp diff --git a/include/GMGPolar/igmgpolar.h b/include/GMGPolar/igmgpolar.h index 3d06bec8..0ba8861f 100644 --- a/include/GMGPolar/igmgpolar.h +++ b/include/GMGPolar/igmgpolar.h @@ -192,6 +192,7 @@ class IGMGPolar double timeSolveMultigridIterations() const; double timeCheckConvergence() const; double timeCheckExactError() const; + double timeConjugateGradient() const; double timeAvgMGCTotal() const; double timeAvgMGCPreSmoothing() const; diff --git a/include/GMGPolar/setup.h b/include/GMGPolar/setup.h index de07e019..1306bee4 100644 --- a/include/GMGPolar/setup.h +++ b/include/GMGPolar/setup.h @@ -148,7 +148,8 @@ void GMGPolar::setup() levels_[0].initializeSmoothing(domain_geometry_, density_profile_coefficients_, DirBC_Interior_, max_omp_threads_, stencil_distribution_method_); } - if (do_extrapolated_smoothing) { + // PCG doesn't use extrapolated smoothing, so we only initialize it if PCG is disabled. + if (do_extrapolated_smoothing && !PCG_) { levels_[0].initializeExtrapolatedSmoothing(domain_geometry_, density_profile_coefficients_, DirBC_Interior_, max_omp_threads_, stencil_distribution_method_); } diff --git a/src/GMGPolar/gmgpolar.cpp b/src/GMGPolar/gmgpolar.cpp index 58e791d8..52ad76d4 100644 --- a/src/GMGPolar/gmgpolar.cpp +++ b/src/GMGPolar/gmgpolar.cpp @@ -358,6 +358,11 @@ double IGMGPolar::timeCheckExactError() const return t_check_exact_error_; } +double IGMGPolar::timeConjugateGradient() const +{ + return t_conjugate_gradient_; +} + /* ---------------------------------------------------------------------- */ /* Average Multigrid Cycle timings */ /* ---------------------------------------------------------------------- */ diff --git a/src/GMGPolar/setup.cpp b/src/GMGPolar/setup.cpp index ff896ba1..188cbfdf 100644 --- a/src/GMGPolar/setup.cpp +++ b/src/GMGPolar/setup.cpp @@ -198,4 +198,49 @@ void IGMGPolar::printSettings() const else { std::cout << "Full-Multigrid: Disabled\n"; } + + if (PCG_) { + std::cout << "Preconditioned Conjugate Gradient: Enabled\n"; + if (PCG_FMG_) { + std::cout << "- PCG Full-Multigrid: " << PCG_FMG_iterations_ << "x "; + switch (PCG_FMG_cycle_) { + case MultigridCycleType::V_CYCLE: + std::cout << "V(" << pre_smoothing_steps_ << "," << post_smoothing_steps_ << ")-Cycle\n"; + break; + case MultigridCycleType::W_CYCLE: + std::cout << "W(" << pre_smoothing_steps_ << "," << post_smoothing_steps_ << ")-Cycle\n"; + break; + case MultigridCycleType::F_CYCLE: + std::cout << "F(" << pre_smoothing_steps_ << "," << post_smoothing_steps_ << ")-Cycle\n"; + break; + default: + std::cout << "Unknown Configuration\n"; + break; + } + } + else { + std::cout << "PCG Full-Multigrid: Disabled\n"; + } + + std::cout << "- PCG Multigrid Iteration: " << PCG_MG_iterations_ << "x "; + switch (PCG_MG_cycle_) { + case MultigridCycleType::V_CYCLE: + std::cout << "V(" << pre_smoothing_steps_ << "," << post_smoothing_steps_ << ")-Cycle\n"; + break; + case MultigridCycleType::W_CYCLE: + std::cout << "W(" << pre_smoothing_steps_ << "," << post_smoothing_steps_ << ")-Cycle\n"; + break; + case MultigridCycleType::F_CYCLE: + std::cout << "F(" << pre_smoothing_steps_ << "," << post_smoothing_steps_ << ")-Cycle\n"; + break; + default: + std::cout << "Unknown Configuration\n"; + break; + } + } + else { + std::cout << "Preconditioned Conjugate Gradient: Disabled\n"; + } + + std::cout << "------------------------------\n"; } diff --git a/src/GMGPolar/solver.cpp b/src/GMGPolar/solver.cpp index 58107d13..39d168c8 100644 --- a/src/GMGPolar/solver.cpp +++ b/src/GMGPolar/solver.cpp @@ -386,7 +386,7 @@ void IGMGPolar::solvePCG(double& initial_residual_norm, double& current_residual auto start_check_exact_error = std::chrono::high_resolution_clock::now(); if (exact_solution_ != nullptr) - evaluateExactError(level, *exact_solution_); + exact_errors_.push_back(computeExactError(level, pcg_solution_, level.residual(), *exact_solution_)); auto end_check_exact_error = std::chrono::high_resolution_clock::now(); t_check_exact_error_ += std::chrono::duration(end_check_exact_error - start_check_exact_error).count(); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f5699a4e..e02b7577 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -25,6 +25,7 @@ add_executable(gmgpolar_tests ExtrapolatedSmoother/extrapolated_smoother.cpp ConfigParser/config_parser.cpp GMGPolar/solve_tests.cpp + GMGPolar/pcg_tests.cpp GMGPolar/convergence_order.cpp ) diff --git a/tests/GMGPolar/pcg_tests.cpp b/tests/GMGPolar/pcg_tests.cpp new file mode 100644 index 00000000..e1e799fb --- /dev/null +++ b/tests/GMGPolar/pcg_tests.cpp @@ -0,0 +1,783 @@ +#include +#include +#include +#include +#include + +// Including the necessary header from the project +#include "../../include/GMGPolar/gmgpolar.h" + +template +class PCGTestCase; + +// clang-format off +template < + // Parameters + class DomainGeometryType, + class DensityProfileCoefficientsType, + class BoundaryConditionsType, + class SourceTermType, + class ExactSolutionType, + double R0_, + double Rmax_, + int nrExp_, + int nthetaExp_, + double refinementRadius_, + int anisotropicFactor_, + int divideBy2_, + int verbose_, + int maxOpenMPThreads_, + bool DirBC_Interior_, + StencilDistributionMethod stencilDistributionMethod_, + bool cacheDensityProfileCoefficients_, + bool cacheDomainGeometry_, + ExtrapolationType extrapolation_, + int maxLevels_, + MultigridCycleType multigridCycle_, + bool FMG_, + int FMG_iterations_, + MultigridCycleType FMG_cycle_, + bool PCG_, + bool PCG_FMG_, + int PCG_FMG_iterations_, + MultigridCycleType PCG_FMG_cycle_, + int PCG_MG_iterations_, + MultigridCycleType PCG_MG_cycle_, + int maxIterations_, + ResidualNormType residualNormType_, + double absoluteTolerance_, + double relativeTolerance_, + // Results + int expected_iterations_, + double expected_l2_error_, + double expected_inf_error_, + double expected_residual_reduction_ +> +class PCGTestCase< + std::tuple< + DomainGeometryType, + DensityProfileCoefficientsType, + BoundaryConditionsType, + SourceTermType, + ExactSolutionType, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant, + std::integral_constant + > +> : public testing::Test { // clang-format on +public: + // Renaming static constexpr variables to avoid shadowing the template parameters + using DomainGeometry = DomainGeometryType; + using DensityProfileCoefficients = DensityProfileCoefficientsType; + using BoundaryConditions = BoundaryConditionsType; + using SourceTerm = SourceTermType; + using ExactSolution = ExactSolutionType; + static constexpr double R0 = R0_; + static constexpr double Rmax = Rmax_; + static constexpr int nrExp = nrExp_; + static constexpr int nthetaExp = nthetaExp_; + static constexpr double refinementRadius = refinementRadius_; + static constexpr int anisotropicFactor = anisotropicFactor_; + static constexpr int divideBy2 = divideBy2_; + static constexpr int verbose = verbose_; + static constexpr int maxOpenMPThreads = maxOpenMPThreads_; + static constexpr bool DirBC_Interior = DirBC_Interior_; + static constexpr StencilDistributionMethod stencilDistributionMethod = stencilDistributionMethod_; + static constexpr bool cacheDensityProfileCoefficients = cacheDensityProfileCoefficients_; + static constexpr bool cacheDomainGeometry = cacheDomainGeometry_; + static constexpr ExtrapolationType extrapolation = extrapolation_; + static constexpr int maxLevels = maxLevels_; + static constexpr MultigridCycleType multigridCycle = multigridCycle_; + static constexpr bool FMG = FMG_; + static constexpr int FMG_iterations = FMG_iterations_; + static constexpr MultigridCycleType FMG_cycle = FMG_cycle_; + static constexpr bool PCG = PCG_; + static constexpr bool PCG_FMG = PCG_FMG_; + static constexpr int PCG_FMG_iterations = PCG_FMG_iterations_; + static constexpr MultigridCycleType PCG_FMG_cycle = PCG_FMG_cycle_; + static constexpr int PCG_MG_iterations = PCG_MG_iterations_; + static constexpr MultigridCycleType PCG_MG_cycle = PCG_MG_cycle_; + static constexpr int maxIterations = maxIterations_; + static constexpr ResidualNormType residualNormType = residualNormType_; + static constexpr double absoluteTolerance = absoluteTolerance_; + static constexpr double relativeTolerance = relativeTolerance_; + static constexpr int expected_iterations = expected_iterations_; + static constexpr double expected_l2_error = expected_l2_error_; + static constexpr double expected_inf_error = expected_inf_error_; + static constexpr double expected_residual_reduction = expected_residual_reduction_; +}; + +// clang-format off +using gmgpolar_test_suite = testing::Types< + /* --------------- */ + /* 1. Verbose Test */ + /* --------------- */ + std::tuple< + CzarnyGeometry, + ZoniShiftedGyroCoefficients, + PolarR6_Boundary_CzarnyGeometry, + PolarR6_ZoniShiftedGyro_CzarnyGeometry, + PolarR6_CzarnyGeometry, + std::integral_constant, // R0 + std::integral_constant, // Rmax + std::integral_constant, // nrExp + std::integral_constant, // nthetaExp + std::integral_constant, // refinementRadius + std::integral_constant, // anisotropicFactor + std::integral_constant, // divideBy2 + std::integral_constant, // verbose + std::integral_constant, // maxOpenMPThreads + std::integral_constant, // DirBC_Interior + std::integral_constant, // StencilDistributionMethod + std::integral_constant, // cacheDensityProfileCoefficient + std::integral_constant, // cacheDomainGeometry + std::integral_constant, // extrapolation + std::integral_constant, // maxLevels + std::integral_constant, // multigridCycle + std::integral_constant, // FMG + std::integral_constant, // FMG_iterations + std::integral_constant, // FMG_cycle + std::integral_constant, // PCG + std::integral_constant, // PCG_FMG + std::integral_constant, // PCG_FMG_iterations + std::integral_constant, // PCG_FMG_cycle + std::integral_constant, // PCG_MG_iterations + std::integral_constant, // PCG_MG_cycle + std::integral_constant, // maxIterations + std::integral_constant, // residualNormType + std::integral_constant, // absoluteTolerance + std::integral_constant, // relativeTolerance + std::integral_constant, // expected_iterations + std::integral_constant, // expected_l2_error + std::integral_constant, // expected_inf_error + std::integral_constant // expected_residual_reduction + >, + /* ---------------------------------------- */ + /* 2.1 Stencil Test: Give (Without caching) */ + /* ---------------------------------------- */ + std::tuple< + CzarnyGeometry, + ZoniShiftedGyroCoefficients, + PolarR6_Boundary_CzarnyGeometry, + PolarR6_ZoniShiftedGyro_CzarnyGeometry, + PolarR6_CzarnyGeometry, + std::integral_constant, // R0 + std::integral_constant, // Rmax + std::integral_constant, // nrExp + std::integral_constant, // nthetaExp + std::integral_constant, // refinementRadius + std::integral_constant, // anisotropicFactor + std::integral_constant, // divideBy2 + std::integral_constant, // verbose + std::integral_constant, // maxOpenMPThreads + std::integral_constant, // DirBC_Interior + std::integral_constant, // StencilDistributionMethod + std::integral_constant, // cacheDensityProfileCoefficient + std::integral_constant, // cacheDomainGeometry + std::integral_constant, // extrapolation + std::integral_constant, // maxLevels + std::integral_constant, // multigridCycle + std::integral_constant, // FMG + std::integral_constant, // FMG_iterations + std::integral_constant, // FMG_cycle + std::integral_constant, // PCG + std::integral_constant, // PCG_FMG + std::integral_constant, // PCG_FMG_iterations + std::integral_constant, // PCG_FMG_cycle + std::integral_constant, // PCG_MG_iterations + std::integral_constant, // PCG_MG_cycle + std::integral_constant, // maxIterations + std::integral_constant, // residualNormType + std::integral_constant, // absoluteTolerance + std::integral_constant, // relativeTolerance + std::integral_constant, // expected_iterations + std::integral_constant, // expected_l2_error + std::integral_constant, // expected_inf_error + std::integral_constant // expected_residual_reduction + >, + /* ------------------------------------- */ + /* 2.2 Stencil Test: Give (With caching) */ + /* ------------------------------------- */ + std::tuple< + CzarnyGeometry, + ZoniShiftedGyroCoefficients, + PolarR6_Boundary_CzarnyGeometry, + PolarR6_ZoniShiftedGyro_CzarnyGeometry, + PolarR6_CzarnyGeometry, + std::integral_constant, // R0 + std::integral_constant, // Rmax + std::integral_constant, // nrExp + std::integral_constant, // nthetaExp + std::integral_constant, // refinementRadius + std::integral_constant, // anisotropicFactor + std::integral_constant, // divideBy2 + std::integral_constant, // verbose + std::integral_constant, // maxOpenMPThreads + std::integral_constant, // DirBC_Interior + std::integral_constant, // StencilDistributionMethod + std::integral_constant, // cacheDensityProfileCoefficient + std::integral_constant, // cacheDomainGeometry + std::integral_constant, // extrapolation + std::integral_constant, // maxLevels + std::integral_constant, // multigridCycle + std::integral_constant, // FMG + std::integral_constant, // FMG_iterations + std::integral_constant, // FMG_cycle + std::integral_constant, // PCG + std::integral_constant, // PCG_FMG + std::integral_constant, // PCG_FMG_iterations + std::integral_constant, // PCG_FMG_cycle + std::integral_constant, // PCG_MG_iterations + std::integral_constant, // PCG_MG_cycle + std::integral_constant, // maxIterations + std::integral_constant, // residualNormType + std::integral_constant, // absoluteTolerance + std::integral_constant, // relativeTolerance + std::integral_constant, // expected_iterations + std::integral_constant, // expected_l2_error + std::integral_constant, // expected_inf_error + std::integral_constant // expected_residual_reduction + >, + /* ------------------------------------- */ + /* 2.3 Stencil Test: Take (With caching) */ + /* ------------------------------------- */ + std::tuple< + CzarnyGeometry, + ZoniShiftedGyroCoefficients, + PolarR6_Boundary_CzarnyGeometry, + PolarR6_ZoniShiftedGyro_CzarnyGeometry, + PolarR6_CzarnyGeometry, + std::integral_constant, // R0 + std::integral_constant, // Rmax + std::integral_constant, // nrExp + std::integral_constant, // nthetaExp + std::integral_constant, // refinementRadius + std::integral_constant, // anisotropicFactor + std::integral_constant, // divideBy2 + std::integral_constant, // verbose + std::integral_constant, // maxOpenMPThreads + std::integral_constant, // DirBC_Interior + std::integral_constant, // StencilDistributionMethod + std::integral_constant, // cacheDensityProfileCoefficient + std::integral_constant, // cacheDomainGeometry + std::integral_constant, // extrapolation + std::integral_constant, // maxLevels + std::integral_constant, // multigridCycle + std::integral_constant, // FMG + std::integral_constant, // FMG_iterations + std::integral_constant, // FMG_cycle + std::integral_constant, // PCG + std::integral_constant, // PCG_FMG + std::integral_constant, // PCG_FMG_iterations + std::integral_constant, // PCG_FMG_cycle + std::integral_constant, // PCG_MG_iterations + std::integral_constant, // PCG_MG_cycle + std::integral_constant, // maxIterations + std::integral_constant, // residualNormType + std::integral_constant, // absoluteTolerance + std::integral_constant, // relativeTolerance + std::integral_constant, // expected_iterations + std::integral_constant, // expected_l2_error + std::integral_constant, // expected_inf_error + std::integral_constant // expected_residual_reduction + >, + /* No Extrapolation */ + std::tuple< + CzarnyGeometry, + ZoniShiftedGyroCoefficients, + PolarR6_Boundary_CzarnyGeometry, + PolarR6_ZoniShiftedGyro_CzarnyGeometry, + PolarR6_CzarnyGeometry, + std::integral_constant, // R0 + std::integral_constant, // Rmax + std::integral_constant, // nrExp + std::integral_constant, // nthetaExp + std::integral_constant, // refinementRadius + std::integral_constant, // anisotropicFactor + std::integral_constant, // divideBy2 + std::integral_constant, // verbose + std::integral_constant, // maxOpenMPThreads + std::integral_constant, // DirBC_Interior + std::integral_constant, // StencilDistributionMethod + std::integral_constant, // cacheDensityProfileCoefficient + std::integral_constant, // cacheDomainGeometry + std::integral_constant, // extrapolation + std::integral_constant, // maxLevels + std::integral_constant, // multigridCycle + std::integral_constant, // FMG + std::integral_constant, // FMG_iterations + std::integral_constant, // FMG_cycle + std::integral_constant, // PCG + std::integral_constant, // PCG_FMG + std::integral_constant, // PCG_FMG_iterations + std::integral_constant, // PCG_FMG_cycle + std::integral_constant, // PCG_MG_iterations + std::integral_constant, // PCG_MG_cycle + std::integral_constant, // maxIterations + std::integral_constant, // residualNormType + std::integral_constant, // absoluteTolerance + std::integral_constant, // relativeTolerance + std::integral_constant, // expected_iterations + std::integral_constant, // expected_l2_error + std::integral_constant, // expected_inf_error + std::integral_constant // expected_residual_reduction + >, + /* Implicit Extrapolation */ + std::tuple< + CzarnyGeometry, + ZoniShiftedGyroCoefficients, + PolarR6_Boundary_CzarnyGeometry, + PolarR6_ZoniShiftedGyro_CzarnyGeometry, + PolarR6_CzarnyGeometry, + std::integral_constant, // R0 + std::integral_constant, // Rmax + std::integral_constant, // nrExp + std::integral_constant, // nthetaExp + std::integral_constant, // refinementRadius + std::integral_constant, // anisotropicFactor + std::integral_constant, // divideBy2 + std::integral_constant, // verbose + std::integral_constant, // maxOpenMPThreads + std::integral_constant, // DirBC_Interior + std::integral_constant, // StencilDistributionMethod + std::integral_constant, // cacheDensityProfileCoefficient + std::integral_constant, // cacheDomainGeometry + std::integral_constant, // extrapolation + std::integral_constant, // maxLevels + std::integral_constant, // multigridCycle + std::integral_constant, // FMG + std::integral_constant, // FMG_iterations + std::integral_constant, // FMG_cycle + std::integral_constant, // PCG + std::integral_constant, // PCG_FMG + std::integral_constant, // PCG_FMG_iterations + std::integral_constant, // PCG_FMG_cycle + std::integral_constant, // PCG_MG_iterations + std::integral_constant, // PCG_MG_cycle + std::integral_constant, // maxIterations + std::integral_constant, // residualNormType + std::integral_constant, // absoluteTolerance + std::integral_constant, // relativeTolerance + std::integral_constant, // expected_iterations + std::integral_constant, // expected_l2_error + std::integral_constant, // expected_inf_error + std::integral_constant // expected_residual_reduction + >, + /* Combined Extrapolation */ + std::tuple< + CzarnyGeometry, + ZoniShiftedGyroCoefficients, + PolarR6_Boundary_CzarnyGeometry, + PolarR6_ZoniShiftedGyro_CzarnyGeometry, + PolarR6_CzarnyGeometry, + std::integral_constant, // R0 + std::integral_constant, // Rmax + std::integral_constant, // nrExp + std::integral_constant, // nthetaExp + std::integral_constant, // refinementRadius + std::integral_constant, // anisotropicFactor + std::integral_constant, // divideBy2 + std::integral_constant, // verbose + std::integral_constant, // maxOpenMPThreads + std::integral_constant, // DirBC_Interior + std::integral_constant, // StencilDistributionMethod + std::integral_constant, // cacheDensityProfileCoefficient + std::integral_constant, // cacheDomainGeometry + std::integral_constant, // extrapolation + std::integral_constant, // maxLevels + std::integral_constant, // multigridCycle + std::integral_constant, // FMG + std::integral_constant, // FMG_iterations + std::integral_constant, // FMG_cycle + std::integral_constant, // PCG + std::integral_constant, // PCG_FMG + std::integral_constant, // PCG_FMG_iterations + std::integral_constant, // PCG_FMG_cycle + std::integral_constant, // PCG_MG_iterations + std::integral_constant, // PCG_MG_cycle + std::integral_constant, // maxIterations + std::integral_constant, // residualNormType + std::integral_constant, // absoluteTolerance + std::integral_constant, // relativeTolerance + std::integral_constant, // expected_iterations + std::integral_constant, // expected_l2_error + std::integral_constant, // expected_inf_error + std::integral_constant // expected_residual_reduction + >, + /* Implicit Extrapolation with full grid smoothing */ + std::tuple< + CzarnyGeometry, + ZoniShiftedGyroCoefficients, + PolarR6_Boundary_CzarnyGeometry, + PolarR6_ZoniShiftedGyro_CzarnyGeometry, + PolarR6_CzarnyGeometry, + std::integral_constant, // R0 + std::integral_constant, // Rmax + std::integral_constant, // nrExp + std::integral_constant, // nthetaExp + std::integral_constant, // refinementRadius + std::integral_constant, // anisotropicFactor + std::integral_constant, // divideBy2 + std::integral_constant, // verbose + std::integral_constant, // maxOpenMPThreads + std::integral_constant, // DirBC_Interior + std::integral_constant, // StencilDistributionMethod + std::integral_constant, // cacheDensityProfileCoefficient + std::integral_constant, // cacheDomainGeometry + std::integral_constant, // extrapolation + std::integral_constant, // maxLevels + std::integral_constant, // multigridCycle + std::integral_constant, // FMG + std::integral_constant, // FMG_iterations + std::integral_constant, // FMG_cycle + std::integral_constant, // PCG + std::integral_constant, // PCG_FMG + std::integral_constant, // PCG_FMG_iterations + std::integral_constant, // PCG_FMG_cycle + std::integral_constant, // PCG_MG_iterations + std::integral_constant, // PCG_MG_cycle + std::integral_constant, // maxIterations + std::integral_constant, // residualNormType + std::integral_constant, // absoluteTolerance + std::integral_constant, // relativeTolerance + std::integral_constant, // expected_iterations + std::integral_constant, // expected_l2_error + std::integral_constant, // expected_inf_error + std::integral_constant // expected_residual_reduction + >, + std::tuple< + CzarnyGeometry, + ZoniShiftedGyroCoefficients, + PolarR6_Boundary_CzarnyGeometry, + PolarR6_ZoniShiftedGyro_CzarnyGeometry, + PolarR6_CzarnyGeometry, + std::integral_constant, // R0 + std::integral_constant, // Rmax + std::integral_constant, // nrExp + std::integral_constant, // nthetaExp + std::integral_constant, // refinementRadius + std::integral_constant, // anisotropicFactor + std::integral_constant, // divideBy2 + std::integral_constant, // verbose + std::integral_constant, // maxOpenMPThreads + std::integral_constant, // DirBC_Interior + std::integral_constant, // StencilDistributionMethod + std::integral_constant, // cacheDensityProfileCoefficient + std::integral_constant, // cacheDomainGeometry + std::integral_constant, // extrapolation + std::integral_constant, // maxLevels + std::integral_constant, // multigridCycle + std::integral_constant, // FMG + std::integral_constant, // FMG_iterations + std::integral_constant, // FMG_cycle + std::integral_constant, // PCG + std::integral_constant, // PCG_FMG + std::integral_constant, // PCG_FMG_iterations + std::integral_constant, // PCG_FMG_cycle + std::integral_constant, // PCG_MG_iterations + std::integral_constant, // PCG_MG_cycle + std::integral_constant, // maxIterations + std::integral_constant, // residualNormType + std::integral_constant, // absoluteTolerance + std::integral_constant, // relativeTolerance + std::integral_constant, // expected_iterations + std::integral_constant, // expected_l2_error + std::integral_constant, // expected_inf_error + std::integral_constant // expected_residual_reduction + >, + std::tuple< + CzarnyGeometry, + ZoniShiftedGyroCoefficients, + PolarR6_Boundary_CzarnyGeometry, + PolarR6_ZoniShiftedGyro_CzarnyGeometry, + PolarR6_CzarnyGeometry, + std::integral_constant, // R0 + std::integral_constant, // Rmax + std::integral_constant, // nrExp + std::integral_constant, // nthetaExp + std::integral_constant, // refinementRadius + std::integral_constant, // anisotropicFactor + std::integral_constant, // divideBy2 + std::integral_constant, // verbose + std::integral_constant, // maxOpenMPThreads + std::integral_constant, // DirBC_Interior + std::integral_constant, // StencilDistributionMethod + std::integral_constant, // cacheDensityProfileCoefficient + std::integral_constant, // cacheDomainGeometry + std::integral_constant, // extrapolation + std::integral_constant, // maxLevels + std::integral_constant, // multigridCycle + std::integral_constant, // FMG + std::integral_constant, // FMG_iterations + std::integral_constant, // FMG_cycle + std::integral_constant, // PCG + std::integral_constant, // PCG_FMG + std::integral_constant, // PCG_FMG_iterations + std::integral_constant, // PCG_FMG_cycle + std::integral_constant, // PCG_MG_iterations + std::integral_constant, // PCG_MG_cycle + std::integral_constant, // maxIterations + std::integral_constant, // residualNormType + std::integral_constant, // absoluteTolerance + std::integral_constant, // relativeTolerance + std::integral_constant, // expected_iterations + std::integral_constant, // expected_l2_error + std::integral_constant, // expected_inf_error + std::integral_constant // expected_residual_reduction + >, + std::tuple< + CzarnyGeometry, + ZoniShiftedGyroCoefficients, + PolarR6_Boundary_CzarnyGeometry, + PolarR6_ZoniShiftedGyro_CzarnyGeometry, + PolarR6_CzarnyGeometry, + std::integral_constant, // R0 + std::integral_constant, // Rmax + std::integral_constant, // nrExp + std::integral_constant, // nthetaExp + std::integral_constant, // refinementRadius + std::integral_constant, // anisotropicFactor + std::integral_constant, // divideBy2 + std::integral_constant, // verbose + std::integral_constant, // maxOpenMPThreads + std::integral_constant, // DirBC_Interior + std::integral_constant, // StencilDistributionMethod + std::integral_constant, // cacheDensityProfileCoefficient + std::integral_constant, // cacheDomainGeometry + std::integral_constant, // extrapolation + std::integral_constant, // maxLevels + std::integral_constant, // multigridCycle + std::integral_constant, // FMG + std::integral_constant, // FMG_iterations + std::integral_constant, // FMG_cycle + std::integral_constant, // PCG + std::integral_constant, // PCG_FMG + std::integral_constant, // PCG_FMG_iterations + std::integral_constant, // PCG_FMG_cycle + std::integral_constant, // PCG_MG_iterations + std::integral_constant, // PCG_MG_cycle + std::integral_constant, // maxIterations + std::integral_constant, // residualNormType + std::integral_constant, // absoluteTolerance + std::integral_constant, // relativeTolerance + std::integral_constant, // expected_iterations + std::integral_constant, // expected_l2_error + std::integral_constant, // expected_inf_error + std::integral_constant // expected_residual_reduction + >, + std::tuple< + CzarnyGeometry, + SonnendruckerGyroCoefficients, + PolarR6_Boundary_CzarnyGeometry, + PolarR6_Sonnendrucker_CzarnyGeometry, + PolarR6_CzarnyGeometry, + std::integral_constant, // R0 + std::integral_constant, // Rmax + std::integral_constant, // nrExp + std::integral_constant, // nthetaExp + std::integral_constant, // refinementRadius + std::integral_constant, // anisotropicFactor + std::integral_constant, // divideBy2 + std::integral_constant, // verbose + std::integral_constant, // maxOpenMPThreads + std::integral_constant, // DirBC_Interior + std::integral_constant, // StencilDistributionMethod + std::integral_constant, // cacheDensityProfileCoefficient + std::integral_constant, // cacheDomainGeometry + std::integral_constant, // extrapolation + std::integral_constant, // maxLevels + std::integral_constant, // multigridCycle + std::integral_constant, // FMG + std::integral_constant, // FMG_iterations + std::integral_constant, // FMG_cycle + std::integral_constant, // PCG + std::integral_constant, // PCG_FMG + std::integral_constant, // PCG_FMG_iterations + std::integral_constant, // PCG_FMG_cycle + std::integral_constant, // PCG_MG_iterations + std::integral_constant, // PCG_MG_cycle + std::integral_constant, // maxIterations + std::integral_constant, // residualNormType + std::integral_constant, // absoluteTolerance + std::integral_constant, // relativeTolerance + std::integral_constant, // expected_iterations + std::integral_constant, // expected_l2_error + std::integral_constant, // expected_inf_error + std::integral_constant // expected_residual_reduction + > +>; +// clang-format on + +TYPED_TEST_SUITE(PCGTestCase, gmgpolar_test_suite); + +template +void run_gmgpolar() +{ + PolarGrid grid(TestFixture::R0, TestFixture::Rmax, TestFixture::nrExp, TestFixture::nthetaExp, + TestFixture::refinementRadius, TestFixture::anisotropicFactor, TestFixture::divideBy2); + + const double inverse_aspect_ratio_epsilon = 0.3; + const double ellipticity_e = 1.4; + + typename TestFixture::DomainGeometry domain(TestFixture::Rmax, inverse_aspect_ratio_epsilon, ellipticity_e); + typename TestFixture::DensityProfileCoefficients profile_coefficients(TestFixture::Rmax, 0.0); + typename TestFixture::BoundaryConditions boundary_conditions(TestFixture::Rmax, inverse_aspect_ratio_epsilon, + ellipticity_e); + typename TestFixture::SourceTerm source_term(grid, TestFixture::Rmax, inverse_aspect_ratio_epsilon, ellipticity_e); + typename TestFixture::ExactSolution exact_solution(TestFixture::Rmax, inverse_aspect_ratio_epsilon, ellipticity_e); + + GMGPolar solver(grid, domain, profile_coefficients); + + bool paraview = false; + int preSmoothingSteps = 1; + int postSmoothingSteps = 1; + + // --- General solver output and visualization settings --- // + solver.verbose(TestFixture::verbose); + solver.paraview(paraview); + // --- Parallelization and threading settings --- // + solver.maxOpenMPThreads(TestFixture::maxOpenMPThreads); + omp_set_num_threads(TestFixture::maxOpenMPThreads); + // --- Numerical method setup --- // + solver.DirBC_Interior(TestFixture::DirBC_Interior); + solver.stencilDistributionMethod(TestFixture::stencilDistributionMethod); + solver.cacheDensityProfileCoefficients(TestFixture::cacheDensityProfileCoefficients); + solver.cacheDomainGeometry(TestFixture::cacheDomainGeometry); + // --- Multigrid settings --- // + solver.extrapolation(TestFixture::extrapolation); + solver.maxLevels(TestFixture::maxLevels); + solver.preSmoothingSteps(preSmoothingSteps); + solver.postSmoothingSteps(postSmoothingSteps); + solver.multigridCycle(TestFixture::multigridCycle); + solver.FMG(TestFixture::FMG); + solver.FMG_iterations(TestFixture::FMG_iterations); + solver.FMG_cycle(TestFixture::FMG_cycle); + solver.PCG(TestFixture::PCG); + solver.PCG_FMG(TestFixture::PCG_FMG); + solver.PCG_FMG_iterations(TestFixture::PCG_FMG_iterations); + solver.PCG_FMG_cycle(TestFixture::PCG_FMG_cycle); + solver.PCG_MG_iterations(TestFixture::PCG_MG_iterations); + solver.PCG_MG_cycle(TestFixture::PCG_MG_cycle); + // --- Iterative solver controls --- // + solver.maxIterations(TestFixture::maxIterations); + solver.residualNormType(TestFixture::residualNormType); + solver.absoluteTolerance(TestFixture::absoluteTolerance); + solver.relativeTolerance(TestFixture::relativeTolerance); + + // --- Finalize solver setup --- // + solver.setup(); // (allocates internal data, prepares operators, etc.) + + // --- Provide optional exact solution --- // + solver.setSolution(&exact_solution); + // --- Solve Phase --- // + solver.solve(boundary_conditions, source_term); + + // --- Retrieve solution and associated grid --- // + Vector solution = solver.solution(); + const PolarGrid& solution_grid = solver.grid(); + + if (TestFixture::verbose > 0) { + solver.printTimings(); + } + + EXPECT_EQ(solver.verbose(), TestFixture::verbose); + EXPECT_EQ(solver.paraview(), paraview); + EXPECT_EQ(solver.maxOpenMPThreads(), TestFixture::maxOpenMPThreads); + EXPECT_EQ(solver.DirBC_Interior(), TestFixture::DirBC_Interior); + EXPECT_EQ(solver.stencilDistributionMethod(), TestFixture::stencilDistributionMethod); + EXPECT_EQ(solver.cacheDensityProfileCoefficients(), TestFixture::cacheDensityProfileCoefficients); + EXPECT_EQ(solver.cacheDomainGeometry(), TestFixture::cacheDomainGeometry); + EXPECT_EQ(solver.extrapolation(), TestFixture::extrapolation); + EXPECT_EQ(solver.maxLevels(), TestFixture::maxLevels); + EXPECT_EQ(solver.preSmoothingSteps(), preSmoothingSteps); + EXPECT_EQ(solver.postSmoothingSteps(), postSmoothingSteps); + EXPECT_EQ(solver.multigridCycle(), TestFixture::multigridCycle); + EXPECT_EQ(solver.FMG(), TestFixture::FMG); + EXPECT_EQ(solver.FMG_iterations(), TestFixture::FMG_iterations); + EXPECT_EQ(solver.FMG_cycle(), TestFixture::FMG_cycle); + EXPECT_EQ(solver.PCG(), TestFixture::PCG); + EXPECT_EQ(solver.PCG_FMG(), TestFixture::PCG_FMG); + EXPECT_EQ(solver.PCG_FMG_iterations(), TestFixture::PCG_FMG_iterations); + EXPECT_EQ(solver.PCG_FMG_cycle(), TestFixture::PCG_FMG_cycle); + EXPECT_EQ(solver.PCG_MG_iterations(), TestFixture::PCG_MG_iterations); + EXPECT_EQ(solver.PCG_MG_cycle(), TestFixture::PCG_MG_cycle); + EXPECT_EQ(solver.maxIterations(), TestFixture::maxIterations); + EXPECT_EQ(solver.residualNormType(), TestFixture::residualNormType); + if (solver.absoluteTolerance().has_value()) { + EXPECT_DOUBLE_EQ(solver.absoluteTolerance().value(), TestFixture::absoluteTolerance); + } + if (solver.relativeTolerance().has_value()) { + EXPECT_DOUBLE_EQ(solver.relativeTolerance().value(), TestFixture::relativeTolerance); + } + + // --- Verify results --- + int number_of_iterations = solver.numberOfIterations(); + std::optional exact_error_weighted_euclidean = solver.exactErrorWeightedEuclidean(); + std::optional exact_infinity_error = solver.exactErrorInfinity(); + double reductionFactor = solver.meanResidualReductionFactor(); + + ASSERT_TRUE(exact_error_weighted_euclidean.has_value()); + ASSERT_TRUE(exact_infinity_error.has_value()); + + EXPECT_EQ(number_of_iterations, TestFixture::expected_iterations); + EXPECT_LE(exact_error_weighted_euclidean.value(), TestFixture::expected_l2_error); + EXPECT_LE(exact_infinity_error.value(), TestFixture::expected_inf_error); + if (solver.absoluteTolerance().has_value() || solver.relativeTolerance().has_value()) { + EXPECT_LT(reductionFactor, TestFixture::expected_residual_reduction); + } + + // --- Verify timings (all available must be non-negative) --- + EXPECT_GE(solver.timeSetupTotal(), 0.0); + EXPECT_GE(solver.timeSetupCreateLevels(), 0.0); + EXPECT_GE(solver.timeSetupRHS(), 0.0); + EXPECT_GE(solver.timeSetupSmoother(), 0.0); + EXPECT_GE(solver.timeSetupDirectSolver(), 0.0); + + EXPECT_GE(solver.timeSolveTotal(), 0.0); + EXPECT_GE(solver.timeSolveInitialApproximation(), 0.0); + EXPECT_GE(solver.timeSolveMultigridIterations(), 0.0); + EXPECT_GE(solver.timeCheckConvergence(), 0.0); + EXPECT_GE(solver.timeCheckExactError(), 0.0); + EXPECT_GE(solver.timeConjugateGradient(), 0.0); + + EXPECT_GE(solver.timeAvgMGCTotal(), 0.0); + EXPECT_GE(solver.timeAvgMGCPreSmoothing(), 0.0); + EXPECT_GE(solver.timeAvgMGCPostSmoothing(), 0.0); + EXPECT_GE(solver.timeAvgMGCResidual(), 0.0); + EXPECT_GE(solver.timeAvgMGCDirectSolver(), 0.0); +} + +TYPED_TEST(PCGTestCase, GeneralTest) +{ + run_gmgpolar(); +} From 58f200d1e1d56b805fa65fc78db60cd12232c7b9 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Fri, 27 Feb 2026 23:47:51 +0100 Subject: [PATCH 04/10] Move applyMG after initial FMG guess --- src/GMGPolar/solver.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/GMGPolar/solver.cpp b/src/GMGPolar/solver.cpp index 39d168c8..008504eb 100644 --- a/src/GMGPolar/solver.cpp +++ b/src/GMGPolar/solver.cpp @@ -332,12 +332,12 @@ void IGMGPolar::solvePCG(double& initial_residual_norm, double& current_residual if (PCG_FMG_) { initRhsHierarchy(level.rhs()); fullMultigridApproximation(PCG_FMG_cycle_, PCG_FMG_iterations_); - applyMultigridIterations(level, PCG_MG_cycle_, PCG_MG_iterations_); } else { // z = I^{-1} * r (no preconditioning) Kokkos::deep_copy(level.solution(), level.rhs()); } + applyMultigridIterations(level, PCG_MG_cycle_, PCG_MG_iterations_); // p = z Kokkos::deep_copy(pcg_search_direction_, level.solution()); @@ -404,12 +404,12 @@ void IGMGPolar::solvePCG(double& initial_residual_norm, double& current_residual if (PCG_FMG_) { initRhsHierarchy(level.rhs()); fullMultigridApproximation(PCG_FMG_cycle_, PCG_FMG_iterations_); - applyMultigridIterations(level, PCG_MG_cycle_, PCG_MG_iterations_); } else { // z = I^{-1} * r (no preconditioning) Kokkos::deep_copy(level.solution(), level.rhs()); } + applyMultigridIterations(level, PCG_MG_cycle_, PCG_MG_iterations_); // r^T * z double r_z_new = dot_product(ConstVector(level.rhs()), ConstVector(level.solution())); From cf3f21456cb27f351d9002c81ab8c764958fb095 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Sat, 28 Feb 2026 00:29:11 +0100 Subject: [PATCH 05/10] Always initialize residual after initial guess. --- include/GMGPolar/solver.h | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/include/GMGPolar/solver.h b/include/GMGPolar/solver.h index 2b750ff2..d26e210b 100644 --- a/include/GMGPolar/solver.h +++ b/include/GMGPolar/solver.h @@ -104,10 +104,9 @@ void IGMGPolar::solve(const BoundaryConditions& boundary_conditions, const Sourc /* ---------------------------- */ auto start_check_convergence = std::chrono::high_resolution_clock::now(); - if (absolute_tolerance_.has_value() || relative_tolerance_.has_value()) { - updateResidualNorms(level, number_of_iterations_, initial_residual_norm, current_residual_norm, - current_relative_residual_norm); - } + // Initializes level.residual() and sets up the convergence criteria. + updateResidualNorms(level, number_of_iterations_, initial_residual_norm, current_residual_norm, + current_relative_residual_norm); auto end_check_convergence = std::chrono::high_resolution_clock::now(); t_check_convergence_ += std::chrono::duration(end_check_convergence - start_check_convergence).count(); From e8d7055d9f87faed18301ff9d04d501015bd45fa Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Sat, 28 Feb 2026 01:21:39 +0100 Subject: [PATCH 06/10] Simplify initRhsHierarchy() --- src/GMGPolar/solver.cpp | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/GMGPolar/solver.cpp b/src/GMGPolar/solver.cpp index 008504eb..aeb83390 100644 --- a/src/GMGPolar/solver.cpp +++ b/src/GMGPolar/solver.cpp @@ -429,14 +429,11 @@ void IGMGPolar::solvePCG(double& initial_residual_norm, double& current_residual void IGMGPolar::initRhsHierarchy(Vector rhs) { - Level& level = levels_[0]; - Kokkos::deep_copy(level.rhs(), rhs); - for (int level_depth = 0; level_depth < number_of_levels_; ++level_depth) { + Kokkos::deep_copy(levels_[0].rhs(), rhs); + for (int level_depth = 0; level_depth < number_of_levels_ - 1; ++level_depth) { Level& current_level = levels_[level_depth]; - if (level_depth + 1 < number_of_levels_) { - Level& next_level = levels_[level_depth + 1]; - restriction(level_depth, next_level.rhs(), current_level.rhs()); - } + Level& next_level = levels_[level_depth + 1]; + restriction(level_depth, next_level.rhs(), current_level.rhs()); } } From bcd67b1ad17a9acbcd66be22ba261a56d103ca43 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Sat, 28 Feb 2026 13:51:58 +0100 Subject: [PATCH 07/10] Increase PCG_MG_iterations from 1 to 2 --- scripts/tutorial/run.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/tutorial/run.sh b/scripts/tutorial/run.sh index 45271a30..fbd8c446 100755 --- a/scripts/tutorial/run.sh +++ b/scripts/tutorial/run.sh @@ -77,7 +77,7 @@ PCG_FMG=1 PCG_FMG_iterations=1 PCG_FMG_cycle=0 # V-Cycle(0), W-Cycle(1), F-Cycle(2) # Additional multigrid iterations after initial approximation to solve the linear system in each PCG iteration -PCG_MG_iterations=1 +PCG_MG_iterations=2 PCG_MG_cycle=0 # V-Cycle(0), W-Cycle(1), F-Cycle(2) # Extrapolation Method: From 434f1fad31519504a7055dd1e5436a38b5b06b19 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Sat, 28 Feb 2026 13:55:02 +0100 Subject: [PATCH 08/10] Fix formatting of next_level assignment in solver.cpp --- src/GMGPolar/solver.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/GMGPolar/solver.cpp b/src/GMGPolar/solver.cpp index aeb83390..2b721354 100644 --- a/src/GMGPolar/solver.cpp +++ b/src/GMGPolar/solver.cpp @@ -432,7 +432,7 @@ void IGMGPolar::initRhsHierarchy(Vector rhs) Kokkos::deep_copy(levels_[0].rhs(), rhs); for (int level_depth = 0; level_depth < number_of_levels_ - 1; ++level_depth) { Level& current_level = levels_[level_depth]; - Level& next_level = levels_[level_depth + 1]; + Level& next_level = levels_[level_depth + 1]; restriction(level_depth, next_level.rhs(), current_level.rhs()); } } @@ -475,4 +475,4 @@ void IGMGPolar::applyExtrapolatedMultigridIterations(Level& level, MultigridCycl break; } } -} \ No newline at end of file +} From fb8ae900f64ba2da928bd39e62e2c47eee39147e Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Sat, 28 Feb 2026 13:55:44 +0100 Subject: [PATCH 09/10] Remove unnecessary output in setup.cpp --- src/GMGPolar/setup.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/GMGPolar/setup.cpp b/src/GMGPolar/setup.cpp index 188cbfdf..a2f32129 100644 --- a/src/GMGPolar/setup.cpp +++ b/src/GMGPolar/setup.cpp @@ -241,6 +241,4 @@ void IGMGPolar::printSettings() const else { std::cout << "Preconditioned Conjugate Gradient: Disabled\n"; } - - std::cout << "------------------------------\n"; } From 64cdf38b2d01c5a2ea76bbbb08492a27348a0372 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Sat, 28 Feb 2026 18:03:26 +0100 Subject: [PATCH 10/10] Reduce PCG_MG_iterations from 2 to 1 --- scripts/tutorial/run.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/tutorial/run.sh b/scripts/tutorial/run.sh index fbd8c446..45271a30 100755 --- a/scripts/tutorial/run.sh +++ b/scripts/tutorial/run.sh @@ -77,7 +77,7 @@ PCG_FMG=1 PCG_FMG_iterations=1 PCG_FMG_cycle=0 # V-Cycle(0), W-Cycle(1), F-Cycle(2) # Additional multigrid iterations after initial approximation to solve the linear system in each PCG iteration -PCG_MG_iterations=2 +PCG_MG_iterations=1 PCG_MG_cycle=0 # V-Cycle(0), W-Cycle(1), F-Cycle(2) # Extrapolation Method: