Add Preconditioned Conjugate Gradient Method

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_;

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                                           */
     /* ---------------------------------------------------------------------- */
@@ -178,6 +192,7 @@ class IGMGPolar
     double timeSolveMultigridIterations() const;
     double timeCheckConvergence() const;
     double timeCheckExactError() const;
+    double timeConjugateGradient() const;
 
     double timeAvgMGCTotal() const;
     double timeAvgMGCPreSmoothing() const;
@@ -215,6 +230,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 +256,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<double> pcg_solution_; // x (solution)
+    AllocatableVector<double> pcg_search_direction_; // p (search direction)
+
     /* -------------------- */
     /* Convergence criteria */
     int number_of_iterations_;
@@ -263,10 +301,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<double> 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<double> rhs);
+    void applyMultigridIterations(Level& level, MultigridCycleType cycle, int iterations);
+    void applyExtrapolatedMultigridIterations(Level& level, MultigridCycleType cycle, int iterations);
 
     /* ----------------- */
     /* Print information */
@@ -319,6 +363,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_;

include/GMGPolar/setup.h

@@ -56,14 +56,15 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::setup()
 
     auto finest_levelCache = std::make_unique<LevelCache>(*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<PolarGrid>(coarseningGrid(levels_[level_depth - 1].grid()));
         auto current_levelCache =
             std::make_unique<LevelCache>(*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<DomainGeometry, DensityProfileCoefficients>::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<double>("pcg_solution", grid_size);
+        }
+        if (std::ssize(pcg_search_direction_) != grid_size) {
+            pcg_search_direction_ = Vector<double>("pcg_search_direction", grid_size);
+        }
+    }
+
     // ------------------------------------------------ //
     // Define residual operator on all multigrid levels //
     // ------------------------------------------------ //
@@ -129,11 +143,13 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::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_);
         }
-        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_);
         }

include/GMGPolar/solver.h

@@ -86,77 +86,51 @@ 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<double>(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();
+
+    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<double>(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_check_convergence = std::chrono::high_resolution_clock::now();
-        t_check_convergence_ += std::chrono::duration<double>(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");
+    // 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<double>(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);
         }
-        number_of_iterations_++;
+        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).
+            auto start_conjugate_gradient = std::chrono::high_resolution_clock::now();
 
-        auto end_solve_multigrid_iterations = std::chrono::high_resolution_clock::now();
-        t_solve_multigrid_iterations_ +=
-            std::chrono::duration<double>(end_solve_multigrid_iterations - start_solve_multigrid_iterations).count();
+            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<double>(end_conjugate_gradient - start_conjugate_gradient).count();
+        }
     }
 
     /* ---------------------- */

include/Level/level.h

@@ -43,8 +43,8 @@ class Level
     // ----------- //
     // Constructor //
     explicit Level(const int level_depth, std::unique_ptr<const PolarGrid> grid,
-                   std::unique_ptr<const LevelCache> level_cache, const ExtrapolationType extrapolation,
-                   const bool FMG);
+                   std::unique_ptr<const LevelCache> 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<double> result, ConstVector<double> rhs, ConstVector<double> x) const;
+    void applySystemOperator(Vector<double> result, ConstVector<double> x) const;
 
     // ------------------- //
     // Solve coarse System //

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

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-10
+relativeTolerance=1e-10
 
 # 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 \