Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
julianlitz merged 1 commit into from
Feb 27, 2026
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion include/GMGPolar/igmgpolar.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
26 changes: 25 additions & 1 deletion include/GMGPolar/solver.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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_ =
Expand Down
110 changes: 39 additions & 71 deletions src/GMGPolar/solver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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;
}
}
}
Expand Down
14 changes: 7 additions & 7 deletions tests/GMGPolar/solve_tests.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down