Skip to content

Commit cefd330

Browse files
julianlitzjulianlitz
authored
Support single-level multigrid (maxLevels = 1) when extrapolation is disabled (#193)
1 parent d4c3359 commit cefd330

File tree

9 files changed

+349
-190
lines changed

9 files changed

+349
-190
lines changed

‎include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h‎

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
22
void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigrid_F_Cycle(int level_depth,
33
Vector<double> solution,
4-
Vector<double> rhs,
4+
ConstVector<double> rhs,
55
Vector<double> residual)
66
{
77
assert(0 <= level_depth && level_depth < number_of_levels_ - 1);
@@ -94,8 +94,14 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
9494
assign(next_level.residual(), 0.0);
9595

9696
/* Step 3: Solve for the error by recursively calling the multigrid cycle. */
97-
multigrid_F_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
98-
multigrid_V_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
97+
multigrid_F_Cycle(level_depth + 1,
98+
next_level.residual(), // error (solution)
99+
next_level.error_correction(), // coarse residual (rhs)
100+
next_level.solution()); // workspace
101+
multigrid_V_Cycle(level_depth + 1,
102+
next_level.residual(), // error (solution)
103+
next_level.error_correction(), // coarse residual (rhs)
104+
next_level.solution()); // workspace
99105
}
100106

101107
/* Interpolate the correction */

‎include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
22
void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigrid_V_Cycle(int level_depth,
33
Vector<double> solution,
4-
Vector<double> rhs,
4+
ConstVector<double> rhs,
55
Vector<double> residual)
66
{
77
assert(0 <= level_depth && level_depth < number_of_levels_ - 1);
@@ -94,7 +94,10 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
9494
assign(next_level.residual(), 0.0);
9595

9696
/* Step 3: Solve for the error by recursively calling the multigrid cycle. */
97-
multigrid_V_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
97+
multigrid_V_Cycle(level_depth + 1,
98+
next_level.residual(), // error (solution)
99+
next_level.error_correction(), // coarse residual (rhs)
100+
next_level.solution()); // workspace
98101
}
99102

100103
/* Interpolate the correction */

‎include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h‎

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
22
void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigrid_W_Cycle(int level_depth,
33
Vector<double> solution,
4-
Vector<double> rhs,
4+
ConstVector<double> rhs,
55
Vector<double> residual)
66
{
77
assert(0 <= level_depth && level_depth < number_of_levels_ - 1);
@@ -94,8 +94,14 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
9494
assign(next_level.residual(), 0.0);
9595

9696
/* Step 3: Solve for the error by recursively calling the multigrid cycle. */
97-
multigrid_W_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
98-
multigrid_W_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
97+
multigrid_W_Cycle(level_depth + 1,
98+
next_level.residual(), // error (solution)
99+
next_level.error_correction(), // coarse residual (rhs)
100+
next_level.solution()); // workspace
101+
multigrid_W_Cycle(level_depth + 1,
102+
next_level.residual(), // error (solution)
103+
next_level.error_correction(), // coarse residual (rhs)
104+
next_level.solution()); // workspace
99105
}
100106

101107
/* Interpolate the correction */

‎include/GMGPolar/MultigridMethods/multigrid_F_Cycle.h‎

Lines changed: 75 additions & 60 deletions
Original file line numberDiff line numberDiff line change
@@ -1,63 +1,71 @@
1+
#pragma once
2+
13
template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
24
void GMGPolar<DomainGeometry, DensityProfileCoefficients>::multigrid_F_Cycle(int level_depth, Vector<double> solution,
3-
Vector<double> rhs,
5+
ConstVector<double> rhs,
46
Vector<double> residual)
57
{
6-
assert(0 <= level_depth && level_depth < number_of_levels_ - 1);
8+
assert(0 <= level_depth && level_depth < number_of_levels_);
79

810
std::chrono::high_resolution_clock::time_point start_MGC;
911
if (level_depth == 0) {
1012
start_MGC = std::chrono::high_resolution_clock::now();
1113
}
1214

13-
Level<DomainGeometry>& level = levels_[level_depth];
14-
Level<DomainGeometry>& next_level = levels_[level_depth + 1];
15-
16-
auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
17-
18-
/* ------------ */
19-
/* Presmoothing */
20-
for (int i = 0; i < pre_smoothing_steps_; i++) {
21-
level.smoothing(solution, rhs, residual);
22-
}
23-
24-
auto end_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
25-
t_avg_MGC_preSmoothing_ += std::chrono::duration<double>(end_MGC_preSmoothing - start_MGC_preSmoothing).count();
26-
27-
/* ---------------------- */
28-
/* Coarse grid correction */
29-
/* ---------------------- */
30-
31-
auto start_MGC_residual = std::chrono::high_resolution_clock::now();
32-
33-
/* Compute the residual */
34-
level.computeResidual(residual, rhs, solution);
15+
/* ------------------------ */
16+
/* Solve A * solution = rhs */
17+
/* ------------------------ */
18+
if (level_depth == number_of_levels_ - 1) {
19+
/* ------------------------------ */
20+
/* Coarsest level: solve directly */
21+
/* ------------------------------ */
22+
Level<DomainGeometry>& level = levels_[level_depth];
3523

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

39-
/* -------------------------- */
40-
/* Solve A * error = residual */
41-
if (level_depth + 1 == number_of_levels_ - 1) {
42-
/* --------------------- */
43-
/* Using a direct solver */
44-
/* --------------------- */
45-
46-
/* Step 1: Restrict the residual */
47-
restriction(level_depth, next_level.residual(), residual);
48-
49-
/* Step 2: Solve for the error in place */
27+
/* Step 2: Solve for the solution in place */
5028
auto start_MGC_directSolver = std::chrono::high_resolution_clock::now();
5129

52-
next_level.directSolveInPlace(next_level.residual());
30+
level.directSolveInPlace(solution);
5331

5432
auto end_MGC_directSolver = std::chrono::high_resolution_clock::now();
5533
t_avg_MGC_directSolver_ += std::chrono::duration<double>(end_MGC_directSolver - start_MGC_directSolver).count();
5634
}
5735
else {
58-
/* ------------------------------------------ */
36+
/* ------------------------------------------------------- */
37+
/* Multigrid V-cycle on current level: */
38+
/* presmoothing, coarse-grid correction, and postsmoothing */
39+
/* ------------------------------------------------------- */
40+
Level<DomainGeometry>& level = levels_[level_depth];
41+
Level<DomainGeometry>& next_level = levels_[level_depth + 1];
42+
43+
auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
44+
45+
/* ------------ */
46+
/* Presmoothing */
47+
for (int i = 0; i < pre_smoothing_steps_; i++) {
48+
level.smoothing(solution, rhs, residual);
49+
}
50+
51+
auto end_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
52+
t_avg_MGC_preSmoothing_ += std::chrono::duration<double>(end_MGC_preSmoothing - start_MGC_preSmoothing).count();
53+
54+
/* ----------------------------- */
55+
/* Compute fine-grid residual */
56+
/* residual = rhs - A * solution */
57+
/* ----------------------------- */
58+
auto start_MGC_residual = std::chrono::high_resolution_clock::now();
59+
60+
level.computeResidual(residual, rhs, solution);
61+
62+
auto end_MGC_residual = std::chrono::high_resolution_clock::now();
63+
t_avg_MGC_residual_ += std::chrono::duration<double>(end_MGC_residual - start_MGC_residual).count();
64+
65+
/* -------------------------- */
66+
/* Solve A * error = residual */
67+
/* -------------------------- */
5968
/* By recursively calling the multigrid cycle */
60-
/* ------------------------------------------ */
6169

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

6876
/* Step 3: Solve for the error by recursively calling the multigrid cycle. */
69-
multigrid_F_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
70-
multigrid_V_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
71-
}
72-
73-
/* Interpolate the correction */
74-
prolongation(level_depth + 1, residual, next_level.residual());
75-
76-
/* Compute the corrected approximation: u = u + error */
77-
add(solution, ConstVector<double>(residual));
78-
79-
auto start_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
80-
81-
/* ------------- */
82-
/* Postsmoothing */
83-
for (int i = 0; i < post_smoothing_steps_; i++) {
84-
level.smoothing(solution, rhs, residual);
77+
multigrid_F_Cycle(level_depth + 1,
78+
next_level.residual(), // error (solution)
79+
next_level.error_correction(), // coarse residual (rhs)
80+
next_level.solution()); // workspace
81+
multigrid_V_Cycle(level_depth + 1,
82+
next_level.residual(), // error (solution)
83+
next_level.error_correction(), // coarse residual (rhs)
84+
next_level.solution()); // workspace
85+
86+
/* Interpolate the correction */
87+
prolongation(level_depth + 1, residual, next_level.residual());
88+
89+
/* Compute the corrected approximation: u = u + error */
90+
add(solution, ConstVector<double>(residual));
91+
92+
auto start_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
93+
94+
/* ------------- */
95+
/* Postsmoothing */
96+
for (int i = 0; i < post_smoothing_steps_; i++) {
97+
level.smoothing(solution, rhs, residual);
98+
}
99+
100+
auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
101+
t_avg_MGC_postSmoothing_ +=
102+
std::chrono::duration<double>(end_MGC_postSmoothing - start_MGC_postSmoothing).count();
85103
}
86104

87-
auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
88-
t_avg_MGC_postSmoothing_ += std::chrono::duration<double>(end_MGC_postSmoothing - start_MGC_postSmoothing).count();
89-
90105
if (level_depth == 0) {
91106
auto end_MGC = std::chrono::high_resolution_clock::now();
92107
t_avg_MGC_total_ += std::chrono::duration<double>(end_MGC - start_MGC).count();
93108
}
94-
}
109+
}

0 commit comments

Comments
 (0)