Skip to content

Support single-level multigrid (maxLevels = 1) when extrapolation is disabled #193

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
Mar 17, 2026
Merged

Support single-level multigrid (maxLevels = 1) when extrapolation is disabled #193

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
12 changes: 9 additions & 3 deletions include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigrid_F_Cycle(int level_depth,
Vector<double> solution,
Vector<double> rhs,
ConstVector<double> rhs,
Vector<double> residual)
{
assert(0 <= level_depth && level_depth < number_of_levels_ - 1);
Expand Down Expand Up @@ -94,8 +94,14 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
assign(next_level.residual(), 0.0);

/* Step 3: Solve for the error by recursively calling the multigrid cycle. */
multigrid_F_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
multigrid_V_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
multigrid_F_Cycle(level_depth + 1,
next_level.residual(), // error (solution)
next_level.error_correction(), // coarse residual (rhs)
next_level.solution()); // workspace
multigrid_V_Cycle(level_depth + 1,
next_level.residual(), // error (solution)
next_level.error_correction(), // coarse residual (rhs)
next_level.solution()); // workspace
}

/* Interpolate the correction */
Expand Down
7 changes: 5 additions & 2 deletions include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigrid_V_Cycle(int level_depth,
Vector<double> solution,
Vector<double> rhs,
ConstVector<double> rhs,
Vector<double> residual)
{
assert(0 <= level_depth && level_depth < number_of_levels_ - 1);
Expand Down Expand Up @@ -94,7 +94,10 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
assign(next_level.residual(), 0.0);

/* Step 3: Solve for the error by recursively calling the multigrid cycle. */
multigrid_V_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
multigrid_V_Cycle(level_depth + 1,
next_level.residual(), // error (solution)
next_level.error_correction(), // coarse residual (rhs)
next_level.solution()); // workspace
}

/* Interpolate the correction */
Expand Down
12 changes: 9 additions & 3 deletions include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigrid_W_Cycle(int level_depth,
Vector<double> solution,
Vector<double> rhs,
ConstVector<double> rhs,
Vector<double> residual)
{
assert(0 <= level_depth && level_depth < number_of_levels_ - 1);
Expand Down Expand Up @@ -94,8 +94,14 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
assign(next_level.residual(), 0.0);

/* Step 3: Solve for the error by recursively calling the multigrid cycle. */
multigrid_W_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
multigrid_W_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
multigrid_W_Cycle(level_depth + 1,
next_level.residual(), // error (solution)
next_level.error_correction(), // coarse residual (rhs)
next_level.solution()); // workspace
multigrid_W_Cycle(level_depth + 1,
next_level.residual(), // error (solution)
next_level.error_correction(), // coarse residual (rhs)
next_level.solution()); // workspace
}

/* Interpolate the correction */
Expand Down
135 changes: 75 additions & 60 deletions include/GMGPolar/MultigridMethods/multigrid_F_Cycle.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,63 +1,71 @@
#pragma once

template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
void GMGPolar<DomainGeometry, DensityProfileCoefficients>::multigrid_F_Cycle(int level_depth, Vector<double> solution,
Vector<double> rhs,
ConstVector<double> rhs,
Vector<double> residual)
{
assert(0 <= level_depth && level_depth < number_of_levels_ - 1);
assert(0 <= level_depth && level_depth < number_of_levels_);

std::chrono::high_resolution_clock::time_point start_MGC;
if (level_depth == 0) {
start_MGC = std::chrono::high_resolution_clock::now();
}

Level<DomainGeometry>& level = levels_[level_depth];
Level<DomainGeometry>& next_level = levels_[level_depth + 1];

auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now();

/* ------------ */
/* Presmoothing */
for (int i = 0; i < pre_smoothing_steps_; i++) {
level.smoothing(solution, rhs, residual);
}

auto end_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
t_avg_MGC_preSmoothing_ += std::chrono::duration<double>(end_MGC_preSmoothing - start_MGC_preSmoothing).count();

/* ---------------------- */
/* Coarse grid correction */
/* ---------------------- */

auto start_MGC_residual = std::chrono::high_resolution_clock::now();

/* Compute the residual */
level.computeResidual(residual, rhs, solution);
/* ------------------------ */
/* Solve A * solution = rhs */
/* ------------------------ */
if (level_depth == number_of_levels_ - 1) {
/* ------------------------------ */
/* Coarsest level: solve directly */
/* ------------------------------ */
Level<DomainGeometry>& level = levels_[level_depth];

auto end_MGC_residual = std::chrono::high_resolution_clock::now();
t_avg_MGC_residual_ += std::chrono::duration<double>(end_MGC_residual - start_MGC_residual).count();
/* Step 1: Copy rhs in solution */
Kokkos::deep_copy(solution, rhs);

/* -------------------------- */
/* Solve A * error = residual */
if (level_depth + 1 == number_of_levels_ - 1) {
/* --------------------- */
/* Using a direct solver */
/* --------------------- */

/* Step 1: Restrict the residual */
restriction(level_depth, next_level.residual(), residual);

/* Step 2: Solve for the error in place */
/* Step 2: Solve for the solution in place */
auto start_MGC_directSolver = std::chrono::high_resolution_clock::now();

next_level.directSolveInPlace(next_level.residual());
level.directSolveInPlace(solution);

auto end_MGC_directSolver = std::chrono::high_resolution_clock::now();
t_avg_MGC_directSolver_ += std::chrono::duration<double>(end_MGC_directSolver - start_MGC_directSolver).count();
}
else {
/* ------------------------------------------ */
/* ------------------------------------------------------- */
/* Multigrid V-cycle on current level: */
/* presmoothing, coarse-grid correction, and postsmoothing */
/* ------------------------------------------------------- */
Level<DomainGeometry>& level = levels_[level_depth];
Level<DomainGeometry>& next_level = levels_[level_depth + 1];

auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now();

/* ------------ */
/* Presmoothing */
for (int i = 0; i < pre_smoothing_steps_; i++) {
level.smoothing(solution, rhs, residual);
}

auto end_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
t_avg_MGC_preSmoothing_ += std::chrono::duration<double>(end_MGC_preSmoothing - start_MGC_preSmoothing).count();

/* ----------------------------- */
/* Compute fine-grid residual */
/* residual = rhs - A * solution */
/* ----------------------------- */
auto start_MGC_residual = std::chrono::high_resolution_clock::now();

level.computeResidual(residual, rhs, solution);

auto end_MGC_residual = std::chrono::high_resolution_clock::now();
t_avg_MGC_residual_ += std::chrono::duration<double>(end_MGC_residual - start_MGC_residual).count();

/* -------------------------- */
/* Solve A * error = residual */
/* -------------------------- */
/* By recursively calling the multigrid cycle */
/* ------------------------------------------ */

/* Step 1: Restrict the residual. */
restriction(level_depth, next_level.error_correction(), residual);
Expand All @@ -66,29 +74,36 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::multigrid_F_Cycle(int
assign(next_level.residual(), 0.0);

/* Step 3: Solve for the error by recursively calling the multigrid cycle. */
multigrid_F_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
multigrid_V_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
}

/* Interpolate the correction */
prolongation(level_depth + 1, residual, next_level.residual());

/* Compute the corrected approximation: u = u + error */
add(solution, ConstVector<double>(residual));

auto start_MGC_postSmoothing = std::chrono::high_resolution_clock::now();

/* ------------- */
/* Postsmoothing */
for (int i = 0; i < post_smoothing_steps_; i++) {
level.smoothing(solution, rhs, residual);
multigrid_F_Cycle(level_depth + 1,
next_level.residual(), // error (solution)
next_level.error_correction(), // coarse residual (rhs)
next_level.solution()); // workspace
multigrid_V_Cycle(level_depth + 1,
next_level.residual(), // error (solution)
next_level.error_correction(), // coarse residual (rhs)
next_level.solution()); // workspace

/* Interpolate the correction */
prolongation(level_depth + 1, residual, next_level.residual());

/* Compute the corrected approximation: u = u + error */
add(solution, ConstVector<double>(residual));

auto start_MGC_postSmoothing = std::chrono::high_resolution_clock::now();

/* ------------- */
/* Postsmoothing */
for (int i = 0; i < post_smoothing_steps_; i++) {
level.smoothing(solution, rhs, residual);
}

auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
t_avg_MGC_postSmoothing_ +=
std::chrono::duration<double>(end_MGC_postSmoothing - start_MGC_postSmoothing).count();
}

auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
t_avg_MGC_postSmoothing_ += std::chrono::duration<double>(end_MGC_postSmoothing - start_MGC_postSmoothing).count();

if (level_depth == 0) {
auto end_MGC = std::chrono::high_resolution_clock::now();
t_avg_MGC_total_ += std::chrono::duration<double>(end_MGC - start_MGC).count();
}
}
}
Loading
Loading