@@ -86,77 +86,51 @@ void IGMGPolar::solve(const BoundaryConditions& boundary_conditions, const Sourc
8686
8787 printIterationHeader (exact_solution_);
8888
89- while (number_of_iterations_ < max_iterations_) {
90- /* ---------------------------------------------- */
91- /* Test solution against exact solution if given. */
92- /* ---------------------------------------------- */
93- LIKWID_STOP (" Solver"
94- auto start_check_exact_error = std::chrono::high_resolution_clock::now ();
95-
96- if (exact_solution_ != nullptr )
97- evaluateExactError (level, *exact_solution_);
98-
99- auto end_check_exact_error = std::chrono::high_resolution_clock::now ();
100- t_check_exact_error_ += std::chrono::duration<double >(end_check_exact_error - start_check_exact_error).count ();
101- LIKWID_START (" Solver"
102-
103- /* ---------------------------- */
104- /* Compute convergence criteria */
105- /* ---------------------------- */
106- auto start_check_convergence = std::chrono::high_resolution_clock::now ();
107-
108- if (absolute_tolerance_.has_value () || relative_tolerance_.has_value ()) {
109- updateResidualNorms (level, number_of_iterations_, initial_residual_norm, current_residual_norm,
110- current_relative_residual_norm);
111- }
89+ /* ---------------------------------------------- */
90+ /* Test solution against exact solution if given. */
91+ /* ---------------------------------------------- */
92+ LIKWID_STOP (" Solver"
93+ auto start_check_exact_error = std::chrono::high_resolution_clock::now ();
94+
95+ if (exact_solution_ != nullptr )
96+ evaluateExactError (level, *exact_solution_);
97+
98+ auto end_check_exact_error = std::chrono::high_resolution_clock::now ();
99+ t_check_exact_error_ += std::chrono::duration<double >(end_check_exact_error - start_check_exact_error).count ();
100+ LIKWID_START (" Solver"
101+
102+ /* ---------------------------- */
103+ /* Compute convergence criteria */
104+ /* ---------------------------- */
105+ auto start_check_convergence = std::chrono::high_resolution_clock::now ();
112106
113- auto end_check_convergence = std::chrono::high_resolution_clock::now ();
114- t_check_convergence_ += std::chrono::duration<double >(end_check_convergence - start_check_convergence).count ();
115-
116- printIterationInfo (number_of_iterations_, current_residual_norm, current_relative_residual_norm,
117- exact_solution_);
118-
119- if (converged (current_residual_norm, current_relative_residual_norm))
120- break ;
121-
122- /* ----------------------- */
123- /* Perform Multigrid Cycle */
124- /* ----------------------- */
125- auto start_solve_multigrid_iterations = std::chrono::high_resolution_clock::now ();
126-
127- switch (multigrid_cycle_) {
128- case MultigridCycleType::V_CYCLE:
129- if (extrapolation_ == ExtrapolationType::NONE) {
130- multigrid_V_Cycle (level.level_depth (), level.solution (), level.rhs (), level.residual ());
131- }
132- else {
133- extrapolated_multigrid_V_Cycle (level.level_depth (), level.solution (), level.rhs (), level.residual ());
134- }
135- break ;
136- case MultigridCycleType::W_CYCLE:
137- if (extrapolation_ == ExtrapolationType::NONE) {
138- multigrid_W_Cycle (level.level_depth (), level.solution (), level.rhs (), level.residual ());
139- }
140- else {
141- extrapolated_multigrid_W_Cycle (level.level_depth (), level.solution (), level.rhs (), level.residual ());
142- }
143- break ;
144- case MultigridCycleType::F_CYCLE:
145- if (extrapolation_ == ExtrapolationType::NONE) {
146- multigrid_F_Cycle (level.level_depth (), level.solution (), level.rhs (), level.residual ());
147- }
148- else {
149- extrapolated_multigrid_F_Cycle (level.level_depth (), level.solution (), level.rhs (), level.residual ());
150- }
151- break ;
152- default :
153- throw std::invalid_argument (" Unknown MultigridCycleType"
107+ // Initializes level.residual() and sets up the convergence criteria.
108+ updateResidualNorms (level, number_of_iterations_, initial_residual_norm, current_residual_norm,
109+ current_relative_residual_norm);
110+
111+ auto end_check_convergence = std::chrono::high_resolution_clock::now ();
112+ t_check_convergence_ += std::chrono::duration<double >(end_check_convergence - start_check_convergence).count ();
113+
114+ printIterationInfo (number_of_iterations_, current_residual_norm, current_relative_residual_norm, exact_solution_);
115+
116+ if (!converged (current_residual_norm, current_relative_residual_norm)) {
117+ if (!PCG_) {
118+ // Solve A*x = b directly using multigrid cycles (V/W/F-cycle)
119+ // until convergence or max_iterations_ is reached.
120+ solveMultigrid (initial_residual_norm, current_residual_norm, current_relative_residual_norm);
154121 }
155- number_of_iterations_++;
122+ else {
123+ // Solve A*x = b using Preconditioned Conjugate Gradient (PCG),
124+ // with multigrid cycles as the preconditioner (i.e., one multigrid
125+ // cycle approximates the action of A^{-1} at each PCG iteration).
126+ auto start_conjugate_gradient = std::chrono::high_resolution_clock::now ();
156127
157- auto end_solve_multigrid_iterations = std::chrono::high_resolution_clock::now ();
158- t_solve_multigrid_iterations_ +=
159- std::chrono::duration<double >(end_solve_multigrid_iterations - start_solve_multigrid_iterations).count ();
128+ solvePCG (initial_residual_norm, current_residual_norm, current_relative_residual_norm);
129+
130+ auto end_conjugate_gradient = std::chrono::high_resolution_clock::now ();
131+ t_conjugate_gradient_ +=
132+ std::chrono::duration<double >(end_conjugate_gradient - start_conjugate_gradient).count ();
133+ }
160134 }
161135
162136 /* ---------------------- */
0 commit comments