Adjust FMG approximation - Skip costly multigrid cycles on finest level (#185)

Author: julianlitz

include/GMGPolar/igmgpolar.h

@@ -262,7 +262,7 @@ class IGMGPolar
 
     /* --------------- */
     /* Solve Functions */
-    void initializeSolution();
+    void fullMultigridApproximation(MultigridCycleType FMG_cycle, int FMG_iterations);
     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,

include/GMGPolar/solver.h

@@ -41,7 +41,31 @@ void IGMGPolar::solve(const BoundaryConditions& boundary_conditions, const Sourc
     /* ---------------------------- */
     auto start_initial_approximation = std::chrono::high_resolution_clock::now();
 
-    initializeSolution();
+    if (!FMG_) {
+        // Assign zero initial guess if not using FMG
+        assign(levels_[0].solution(), 0.0);
+    }
+    else {
+        // Compute an initial approximation to the discrete system
+        //
+        //     A * levels_[0].solution() = levels_[0].rhs()
+        //
+        // using the Full Multigrid (FMG) algorithm.
+        //
+        // Prerequisite:
+        // The right-hand side must already be properly initialized on all levels,
+        // i.e., constructed on the finest level and transferred to coarser levels
+        // via injection and discretization. This ensures consistency of the coarse-
+        // grid operators and guarantees correctness of the multigrid hierarchy.
+        //
+        // The FMG algorithm performs an exact solve on the coarsest grid and then
+        // successively prolongates the solution to finer grids, applying multigrid
+        // cycles on each level to reduce the error. This produces a high-quality
+        // initial guess on the finest level, typically reducing the error to the
+        // order of the discretization error and significantly accelerating convergence
+        // of the subsequent solve phase.
+        fullMultigridApproximation(FMG_cycle_, FMG_iterations_);
+    }
 
     auto end_initial_approximation = std::chrono::high_resolution_clock::now();
     t_solve_initial_approximation_ =

src/GMGPolar/solver.cpp

@@ -1,80 +1,48 @@
 #include "../../include/GMGPolar/gmgpolar.h"
 
 // =============================================================================
-//   Solution Initialization
+//   Full Multigrid Approximation
 // =============================================================================
 
-void IGMGPolar::initializeSolution()
+void IGMGPolar::fullMultigridApproximation(MultigridCycleType FMG_cycle, int FMG_iterations)
 {
-    if (!FMG_) {
-        int start_level_depth = 0;
-        Level& level          = levels_[start_level_depth];
-        assign(level.solution(), 0.0); // Assign zero initial guess if not using FMG
-    }
-    else {
-        // Start from the coarsest level
-        int coarsest_depth    = number_of_levels_ - 1;
-        Level& coarsest_level = levels_[coarsest_depth];
-
-        // Solve directly on the coarsest level
-        Kokkos::deep_copy(coarsest_level.solution(), coarsest_level.rhs());
-        coarsest_level.directSolveInPlace(coarsest_level.solution()); // Direct solve on coarsest grid
-
-        // Prolongate the solution from the coarsest level up to the finest, while applying Multigrid Cycles on each level
-        for (int depth = coarsest_depth; depth > 0; --depth) {
-            Level& coarse_level = levels_[depth]; // Current coarse level
-            Level& fine_level   = levels_[depth - 1]; // Next finer level
-
-            // The bi-cubic FMG interpolation is of higher order
-            FMGInterpolation(coarse_level.level_depth(), fine_level.solution(), coarse_level.solution());
-
-            // Apply some FMG iterations
-            for (int i = 0; i < FMG_iterations_; i++) {
-                if (fine_level.level_depth() == 0 && (extrapolation_ != ExtrapolationType::NONE)) {
-                    switch (FMG_cycle_) {
-                    case MultigridCycleType::V_CYCLE:
-                        extrapolated_multigrid_V_Cycle(fine_level.level_depth(), fine_level.solution(),
-                                                       fine_level.rhs(), fine_level.residual());
-                        break;
-
-                    case MultigridCycleType::W_CYCLE:
-                        extrapolated_multigrid_W_Cycle(fine_level.level_depth(), fine_level.solution(),
-                                                       fine_level.rhs(), fine_level.residual());
-                        break;
-
-                    case MultigridCycleType::F_CYCLE:
-                        extrapolated_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;
-                    }
-                }
-                else {
-                    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;
-                    }
+    // Start from the coarsest level
+    int coarsest_depth    = number_of_levels_ - 1;
+    Level& coarsest_level = levels_[coarsest_depth];
+
+    // Solve directly on the coarsest level
+    Kokkos::deep_copy(coarsest_level.solution(), coarsest_level.rhs());
+    coarsest_level.directSolveInPlace(coarsest_level.solution()); // Direct solve on coarsest grid
+
+    // Prolongate the solution from the coarsest level up to the finest, while applying Multigrid Cycles on each level
+    for (int depth = coarsest_depth; depth > 0; --depth) {
+        Level& coarse_level = levels_[depth]; // Current coarse level
+        Level& fine_level   = levels_[depth - 1]; // Next finer level
+
+        // The bi-cubic FMG interpolation is of higher order
+        FMGInterpolation(coarse_level.level_depth(), fine_level.solution(), coarse_level.solution());
+
+        // 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;
                 }
             }
         }

tests/GMGPolar/solve_tests.cpp

@@ -190,7 +190,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<ResidualNormType, ResidualNormType::EUCLIDEAN>, // residualNormType
         std::integral_constant<double, 1e-8>, // absoluteTolerance
         std::integral_constant<double, 1e-8>, // relativeTolerance
-        std::integral_constant<int, 26>, // expected_iterations
+        std::integral_constant<int, 27>, // expected_iterations
         std::integral_constant<double, 2e-6>, // expected_l2_error
         std::integral_constant<double, 9e-6>, // expected_inf_error
         std::integral_constant<double, 0.7> // expected_residual_reduction
@@ -227,7 +227,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<ResidualNormType, ResidualNormType::EUCLIDEAN>, // residualNormType
         std::integral_constant<double, 1e-8>, // absoluteTolerance
         std::integral_constant<double, 1e-8>, // relativeTolerance
-        std::integral_constant<int, 26>, // expected_iterations
+        std::integral_constant<int, 27>, // expected_iterations
         std::integral_constant<double, 2e-6>, // expected_l2_error
         std::integral_constant<double, 9e-6>, // expected_inf_error
         std::integral_constant<double, 0.7> // expected_residual_reduction
@@ -264,7 +264,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<ResidualNormType, ResidualNormType::EUCLIDEAN>, // residualNormType
         std::integral_constant<double, 1e-8>, // absoluteTolerance
         std::integral_constant<double, 1e-8>, // relativeTolerance
-        std::integral_constant<int, 26>, // expected_iterations
+        std::integral_constant<int, 27>, // expected_iterations
         std::integral_constant<double, 2e-6>, // expected_l2_error
         std::integral_constant<double, 9e-6>, // expected_inf_error
         std::integral_constant<double, 0.7> // expected_residual_reduction
@@ -299,7 +299,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<ResidualNormType, ResidualNormType::INFINITY_NORM>, // residualNormType
         std::integral_constant<double, 1e-12>, // absoluteTolerance
         std::integral_constant<double, 1e-10>, // relativeTolerance
-        std::integral_constant<int, 12>, // expected_iterations
+        std::integral_constant<int, 13>, // expected_iterations
         std::integral_constant<double, 6e-6>, // expected_l2_error
         std::integral_constant<double, 2e-5>, // expected_inf_error
         std::integral_constant<double, 0.3> // expected_residual_reduction
@@ -334,7 +334,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<ResidualNormType, ResidualNormType::EUCLIDEAN>, // residualNormType
         std::integral_constant<double, 1e-8>, // absoluteTolerance
         std::integral_constant<double, 1e-6>, // relativeTolerance
-        std::integral_constant<int, 32>, // expected_iterations
+        std::integral_constant<int, 34>, // expected_iterations
         std::integral_constant<double, 5e-4>, // expected_l2_error
         std::integral_constant<double, 2e-3>, // expected_inf_error
         std::integral_constant<double, 0.7> // expected_residual_reduction
@@ -472,7 +472,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<ResidualNormType, ResidualNormType::WEIGHTED_EUCLIDEAN>, // residualNormType
         std::integral_constant<double, 1e-8>, // absoluteTolerance
         std::integral_constant<double, -1.0>, // relativeTolerance
-        std::integral_constant<int, 14>, // expected_iterations
+        std::integral_constant<int, 16>, // expected_iterations
         std::integral_constant<double, 9e-5>, // expected_l2_error
         std::integral_constant<double, 3e-4>, // expected_inf_error
         std::integral_constant<double, 0.6> // expected_residual_reduction
@@ -506,7 +506,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<ResidualNormType, ResidualNormType::EUCLIDEAN>, // residualNormType
         std::integral_constant<double, 1e-9>, // absoluteTolerance
         std::integral_constant<double, 1e-8>, // relativeTolerance
-        std::integral_constant<int, 7>, // expected_iterations
+        std::integral_constant<int, 8>, // expected_iterations
         std::integral_constant<double, 5e-6>, // expected_l2_error
         std::integral_constant<double, 2e-5>, // expected_inf_error
         std::integral_constant<double, 0.2> // expected_residual_reduction