diff --git a/include/GMGPolar/igmgpolar.h b/include/GMGPolar/igmgpolar.h index 9e0d4bce..52c317ef 100644 --- a/include/GMGPolar/igmgpolar.h +++ b/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 residual) const; void evaluateExactError(Level& level, const ExactSolution& exact_solution); void updateResidualNorms(Level& level, int iteration, double& initial_residual_norm, double& current_residual_norm, diff --git a/include/GMGPolar/solver.h b/include/GMGPolar/solver.h index 161e9bc5..66bcaee9 100644 --- a/include/GMGPolar/solver.h +++ b/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_ = diff --git a/src/GMGPolar/solver.cpp b/src/GMGPolar/solver.cpp index ebc6d342..b56ae2d5 100644 --- a/src/GMGPolar/solver.cpp +++ b/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; } } } diff --git a/tests/GMGPolar/solve_tests.cpp b/tests/GMGPolar/solve_tests.cpp index f7be2f40..dbf55f6c 100644 --- a/tests/GMGPolar/solve_tests.cpp +++ b/tests/GMGPolar/solve_tests.cpp @@ -190,7 +190,7 @@ using gmgpolar_test_suite = testing::Types< std::integral_constant, // residualNormType std::integral_constant, // absoluteTolerance std::integral_constant, // relativeTolerance - std::integral_constant, // expected_iterations + std::integral_constant, // expected_iterations std::integral_constant, // expected_l2_error std::integral_constant, // expected_inf_error std::integral_constant // expected_residual_reduction @@ -227,7 +227,7 @@ using gmgpolar_test_suite = testing::Types< std::integral_constant, // residualNormType std::integral_constant, // absoluteTolerance std::integral_constant, // relativeTolerance - std::integral_constant, // expected_iterations + std::integral_constant, // expected_iterations std::integral_constant, // expected_l2_error std::integral_constant, // expected_inf_error std::integral_constant // expected_residual_reduction @@ -264,7 +264,7 @@ using gmgpolar_test_suite = testing::Types< std::integral_constant, // residualNormType std::integral_constant, // absoluteTolerance std::integral_constant, // relativeTolerance - std::integral_constant, // expected_iterations + std::integral_constant, // expected_iterations std::integral_constant, // expected_l2_error std::integral_constant, // expected_inf_error std::integral_constant // expected_residual_reduction @@ -299,7 +299,7 @@ using gmgpolar_test_suite = testing::Types< std::integral_constant, // residualNormType std::integral_constant, // absoluteTolerance std::integral_constant, // relativeTolerance - std::integral_constant, // expected_iterations + std::integral_constant, // expected_iterations std::integral_constant, // expected_l2_error std::integral_constant, // expected_inf_error std::integral_constant // expected_residual_reduction @@ -334,7 +334,7 @@ using gmgpolar_test_suite = testing::Types< std::integral_constant, // residualNormType std::integral_constant, // absoluteTolerance std::integral_constant, // relativeTolerance - std::integral_constant, // expected_iterations + std::integral_constant, // expected_iterations std::integral_constant, // expected_l2_error std::integral_constant, // expected_inf_error std::integral_constant // expected_residual_reduction @@ -472,7 +472,7 @@ using gmgpolar_test_suite = testing::Types< std::integral_constant, // residualNormType std::integral_constant, // absoluteTolerance std::integral_constant, // relativeTolerance - std::integral_constant, // expected_iterations + std::integral_constant, // expected_iterations std::integral_constant, // expected_l2_error std::integral_constant, // expected_inf_error std::integral_constant // expected_residual_reduction @@ -506,7 +506,7 @@ using gmgpolar_test_suite = testing::Types< std::integral_constant, // residualNormType std::integral_constant, // absoluteTolerance std::integral_constant, // relativeTolerance - std::integral_constant, // expected_iterations + std::integral_constant, // expected_iterations std::integral_constant, // expected_l2_error std::integral_constant, // expected_inf_error std::integral_constant // expected_residual_reduction