diff --git a/include/GMGPolar/igmgpolar.h b/include/GMGPolar/igmgpolar.h index 9e0d4bce..adb8e2a9 100644 --- a/include/GMGPolar/igmgpolar.h +++ b/include/GMGPolar/igmgpolar.h @@ -277,14 +277,14 @@ class IGMGPolar /* ------------------- */ /* Multigrid Functions */ - void multigrid_V_Cycle(int level_depth, Vector solution, Vector rhs, Vector residual); - void multigrid_W_Cycle(int level_depth, Vector solution, Vector rhs, Vector residual); - void multigrid_F_Cycle(int level_depth, Vector solution, Vector rhs, Vector residual); - void extrapolated_multigrid_V_Cycle(int level_depth, Vector solution, Vector rhs, + void multigrid_V_Cycle(int level_depth, Vector solution, ConstVector rhs, Vector residual); + void multigrid_W_Cycle(int level_depth, Vector solution, ConstVector rhs, Vector residual); + void multigrid_F_Cycle(int level_depth, Vector solution, ConstVector rhs, Vector residual); + void extrapolated_multigrid_V_Cycle(int level_depth, Vector solution, ConstVector rhs, Vector residual); - void extrapolated_multigrid_W_Cycle(int level_depth, Vector solution, Vector rhs, + void extrapolated_multigrid_W_Cycle(int level_depth, Vector solution, ConstVector rhs, Vector residual); - void extrapolated_multigrid_F_Cycle(int level_depth, Vector solution, Vector rhs, + void extrapolated_multigrid_F_Cycle(int level_depth, Vector solution, ConstVector rhs, Vector residual); /* ----------------------- */ diff --git a/include/GMGPolar/setup.h b/include/GMGPolar/setup.h index 8eb78370..f3ca6e55 100644 --- a/include/GMGPolar/setup.h +++ b/include/GMGPolar/setup.h @@ -42,24 +42,22 @@ void GMGPolar::setup() } if (!is_uniform_refinement) { - throw std::runtime_error( - "Extrapolation Error: Finest PolarGrid does not originate from a single uniform refinement."); + std::cerr << "Warning: Finest PolarGrid does not originate from a single uniform refinement.\n"; } } // ---------------------------------------------------------- // // Building PolarGrid and LevelCache for all multigrid levels // // ---------------------------------------------------------- // - number_of_levels_ = chooseNumberOfLevels(*finest_grid); /* Implementation below */ + number_of_levels_ = chooseNumberOfLevels(*finest_grid); levels_.clear(); levels_.reserve(number_of_levels_); - int level_depth = 0; auto finest_levelCache = std::make_unique(*finest_grid, density_profile_coefficients_, domain_geometry_, cache_density_profile_coefficients_, cache_domain_geometry_); - levels_.emplace_back(level_depth, std::move(finest_grid), std::move(finest_levelCache), extrapolation_, FMG_); + levels_.emplace_back(0, std::move(finest_grid), std::move(finest_levelCache), extrapolation_, FMG_); - for (level_depth = 1; level_depth < number_of_levels_; level_depth++) { + for (int level_depth = 1; level_depth < number_of_levels_; level_depth++) { auto current_grid = std::make_unique(coarseningGrid(levels_[level_depth - 1].grid())); auto current_levelCache = std::make_unique(*current_grid, density_profile_coefficients_, domain_geometry_, @@ -81,84 +79,70 @@ void GMGPolar::setup() if (verbose_ > 0) printSettings(); - // ------------------------------------------------------- - // Initializing various operators based on the level index + // ------------------------------------------------ // + // Define residual operator on all multigrid levels // + // ------------------------------------------------ // for (int level_depth = 0; level_depth < number_of_levels_; level_depth++) { - // ---------------------- // - // Level 0 (finest Level) // - // ---------------------- // - if (level_depth == 0) { - auto start_setup_smoother = std::chrono::high_resolution_clock::now(); - switch (extrapolation_) { - case ExtrapolationType::NONE: - full_grid_smoothing_ = true; - levels_[level_depth].initializeSmoothing(domain_geometry_, density_profile_coefficients_, - DirBC_Interior_, max_omp_threads_, - stencil_distribution_method_); - break; - case ExtrapolationType::IMPLICIT_EXTRAPOLATION: - full_grid_smoothing_ = false; - levels_[level_depth].initializeExtrapolatedSmoothing(domain_geometry_, density_profile_coefficients_, - DirBC_Interior_, max_omp_threads_, - stencil_distribution_method_); - break; - case ExtrapolationType::IMPLICIT_FULL_GRID_SMOOTHING: - full_grid_smoothing_ = true; - levels_[level_depth].initializeSmoothing(domain_geometry_, density_profile_coefficients_, - DirBC_Interior_, max_omp_threads_, - stencil_distribution_method_); - break; - case ExtrapolationType::COMBINED: - full_grid_smoothing_ = true; - levels_[level_depth].initializeSmoothing(domain_geometry_, density_profile_coefficients_, - DirBC_Interior_, max_omp_threads_, - stencil_distribution_method_); - levels_[level_depth].initializeExtrapolatedSmoothing(domain_geometry_, density_profile_coefficients_, - DirBC_Interior_, max_omp_threads_, - stencil_distribution_method_); - break; - default: - full_grid_smoothing_ = false; - levels_[level_depth].initializeSmoothing(domain_geometry_, density_profile_coefficients_, - DirBC_Interior_, max_omp_threads_, - stencil_distribution_method_); - levels_[level_depth].initializeExtrapolatedSmoothing(domain_geometry_, density_profile_coefficients_, - DirBC_Interior_, max_omp_threads_, - stencil_distribution_method_); - break; - } - auto end_setup_smoother = std::chrono::high_resolution_clock::now(); - t_setup_smoother_ += std::chrono::duration(end_setup_smoother - start_setup_smoother).count(); - levels_[level_depth].initializeResidual(domain_geometry_, density_profile_coefficients_, DirBC_Interior_, - max_omp_threads_, stencil_distribution_method_); + levels_[level_depth].initializeResidual(domain_geometry_, density_profile_coefficients_, DirBC_Interior_, + max_omp_threads_, stencil_distribution_method_); + } + + // ----------------------------------------- // + // Build direct solver on the coarsest level // + // ----------------------------------------- // + auto start_setup_directSolver = std::chrono::high_resolution_clock::now(); + levels_[number_of_levels_ - 1].initializeDirectSolver(domain_geometry_, density_profile_coefficients_, + DirBC_Interior_, max_omp_threads_, + stencil_distribution_method_); + auto end_setup_directSolver = std::chrono::high_resolution_clock::now(); + t_setup_directSolver_ += std::chrono::duration(end_setup_directSolver - start_setup_directSolver).count(); + + // ---------------------------------------------------------- // + // Build the full-grid smoother and the extrapolated smoother // + // ---------------------------------------------------------- // + auto start_setup_smoother = std::chrono::high_resolution_clock::now(); + + bool do_full_grid_smoothing = false; + bool do_extrapolated_smoothing = false; + + switch (extrapolation_) { + + case ExtrapolationType::NONE: + do_full_grid_smoothing = true; + break; + + case ExtrapolationType::IMPLICIT_EXTRAPOLATION: + do_extrapolated_smoothing = true; + break; + + case ExtrapolationType::IMPLICIT_FULL_GRID_SMOOTHING: + do_full_grid_smoothing = true; + break; + + case ExtrapolationType::COMBINED: + do_full_grid_smoothing = true; + do_extrapolated_smoothing = true; + break; + } + + full_grid_smoothing_ = do_full_grid_smoothing; + + if (number_of_levels_ > 1) { + if (do_full_grid_smoothing) { + levels_[0].initializeSmoothing(domain_geometry_, density_profile_coefficients_, DirBC_Interior_, + max_omp_threads_, stencil_distribution_method_); } - // -------------------------- // - // Level n-1 (coarsest Level) // - // -------------------------- // - else if (level_depth == number_of_levels_ - 1) { - auto start_setup_directSolver = std::chrono::high_resolution_clock::now(); - levels_[level_depth].initializeDirectSolver(domain_geometry_, density_profile_coefficients_, - DirBC_Interior_, max_omp_threads_, - stencil_distribution_method_); - auto end_setup_directSolver = std::chrono::high_resolution_clock::now(); - t_setup_directSolver_ += - std::chrono::duration(end_setup_directSolver - start_setup_directSolver).count(); - levels_[level_depth].initializeResidual(domain_geometry_, density_profile_coefficients_, DirBC_Interior_, - max_omp_threads_, stencil_distribution_method_); + if (do_extrapolated_smoothing) { + levels_[0].initializeExtrapolatedSmoothing(domain_geometry_, density_profile_coefficients_, DirBC_Interior_, + max_omp_threads_, stencil_distribution_method_); } - // ------------------- // - // Intermediate levels // - // ------------------- // - else { - auto start_setup_smoother = std::chrono::high_resolution_clock::now(); + for (int level_depth = 1; level_depth < number_of_levels_ - 1; level_depth++) { levels_[level_depth].initializeSmoothing(domain_geometry_, density_profile_coefficients_, DirBC_Interior_, max_omp_threads_, stencil_distribution_method_); - auto end_setup_smoother = std::chrono::high_resolution_clock::now(); - t_setup_smoother_ += std::chrono::duration(end_setup_smoother - start_setup_smoother).count(); - levels_[level_depth].initializeResidual(domain_geometry_, density_profile_coefficients_, DirBC_Interior_, - max_omp_threads_, stencil_distribution_method_); } } + auto end_setup_smoother = std::chrono::high_resolution_clock::now(); + t_setup_smoother_ += std::chrono::duration(end_setup_smoother - start_setup_smoother).count(); auto end_setup = std::chrono::high_resolution_clock::now(); t_setup_total_ = std::chrono::duration(end_setup - start_setup).count(); diff --git a/src/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.cpp b/src/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.cpp index fa0a138b..e901c7aa 100644 --- a/src/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.cpp @@ -1,6 +1,6 @@ #include "../../../include/GMGPolar/gmgpolar.h" -void IGMGPolar::extrapolated_multigrid_F_Cycle(int level_depth, Vector solution, Vector rhs, +void IGMGPolar::extrapolated_multigrid_F_Cycle(int level_depth, Vector solution, ConstVector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); @@ -93,8 +93,14 @@ void IGMGPolar::extrapolated_multigrid_F_Cycle(int level_depth, Vector s 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 */ diff --git a/src/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.cpp b/src/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.cpp index 725d8085..1d7e8fbb 100644 --- a/src/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.cpp @@ -1,6 +1,6 @@ #include "../../../include/GMGPolar/gmgpolar.h" -void IGMGPolar::extrapolated_multigrid_V_Cycle(int level_depth, Vector solution, Vector rhs, +void IGMGPolar::extrapolated_multigrid_V_Cycle(int level_depth, Vector solution, ConstVector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); @@ -93,7 +93,10 @@ void IGMGPolar::extrapolated_multigrid_V_Cycle(int level_depth, Vector s 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 */ diff --git a/src/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.cpp b/src/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.cpp index cda07568..905f2833 100644 --- a/src/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.cpp @@ -1,6 +1,6 @@ #include "../../../include/GMGPolar/gmgpolar.h" -void IGMGPolar::extrapolated_multigrid_W_Cycle(int level_depth, Vector solution, Vector rhs, +void IGMGPolar::extrapolated_multigrid_W_Cycle(int level_depth, Vector solution, ConstVector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); @@ -93,8 +93,14 @@ void IGMGPolar::extrapolated_multigrid_W_Cycle(int level_depth, Vector s 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 */ diff --git a/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp b/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp index 6c497d6c..9d5bb633 100644 --- a/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp @@ -1,62 +1,69 @@ #include "../../../include/GMGPolar/gmgpolar.h" -void IGMGPolar::multigrid_F_Cycle(int level_depth, Vector solution, Vector rhs, Vector residual) +void IGMGPolar::multigrid_F_Cycle(int level_depth, Vector solution, ConstVector rhs, + Vector 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& level = levels_[level_depth]; - Level& next_level = levels_[level_depth + 1]; + /* ------------------------ */ + /* Solve A * solution = rhs */ + /* ------------------------ */ + if (level_depth == number_of_levels_ - 1) { + /* ------------------------------ */ + /* Coarsest level: solve directly */ + /* ------------------------------ */ + Level& level = levels_[level_depth]; - auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now(); + /* Step 1: Copy rhs in solution */ + Kokkos::deep_copy(solution, rhs); - /* ------------ */ - /* 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(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); - - auto end_MGC_residual = std::chrono::high_resolution_clock::now(); - t_avg_MGC_residual_ += std::chrono::duration(end_MGC_residual - start_MGC_residual).count(); - - /* -------------------------- */ - /* 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(end_MGC_directSolver - start_MGC_directSolver).count(); } else { - /* ------------------------------------------ */ + /* ------------------------------------------------------- */ + /* Multigrid V-cycle on current level: */ + /* presmoothing, coarse-grid correction, and postsmoothing */ + /* ------------------------------------------------------- */ + Level& level = levels_[level_depth]; + Level& 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(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(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); @@ -65,27 +72,34 @@ void IGMGPolar::multigrid_F_Cycle(int level_depth, Vector solution, Vect 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 */ + prolongation(level_depth + 1, residual, next_level.residual()); + + /* Compute the corrected approximation: u = u + error */ + add(solution, ConstVector(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(end_MGC_postSmoothing - start_MGC_postSmoothing).count(); } - /* Interpolate the correction */ - prolongation(level_depth + 1, residual, next_level.residual()); - - /* Compute the corrected approximation: u = u + error */ - add(solution, ConstVector(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(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(end_MGC - start_MGC).count(); diff --git a/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp b/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp index 7bb39c06..64790649 100644 --- a/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp @@ -1,62 +1,69 @@ #include "../../../include/GMGPolar/gmgpolar.h" -void IGMGPolar::multigrid_V_Cycle(int level_depth, Vector solution, Vector rhs, Vector residual) +void IGMGPolar::multigrid_V_Cycle(int level_depth, Vector solution, ConstVector rhs, + Vector 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& level = levels_[level_depth]; - Level& next_level = levels_[level_depth + 1]; + /* ------------------------ */ + /* Solve A * solution = rhs */ + /* ------------------------ */ + if (level_depth == number_of_levels_ - 1) { + /* ------------------------------ */ + /* Coarsest level: solve directly */ + /* ------------------------------ */ + Level& level = levels_[level_depth]; - auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now(); + /* Step 1: Copy rhs in solution */ + Kokkos::deep_copy(solution, rhs); - /* ------------ */ - /* 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(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); - - auto end_MGC_residual = std::chrono::high_resolution_clock::now(); - t_avg_MGC_residual_ += std::chrono::duration(end_MGC_residual - start_MGC_residual).count(); - - /* -------------------------- */ - /* 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(end_MGC_directSolver - start_MGC_directSolver).count(); } else { - /* ------------------------------------------ */ + /* ------------------------------------------------------- */ + /* Multigrid V-cycle on current level: */ + /* presmoothing, coarse-grid correction, and postsmoothing */ + /* ------------------------------------------------------- */ + Level& level = levels_[level_depth]; + Level& 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(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(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); @@ -65,25 +72,29 @@ void IGMGPolar::multigrid_V_Cycle(int level_depth, Vector solution, Vect 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 */ - prolongation(level_depth + 1, residual, next_level.residual()); + /* Interpolate the correction */ + prolongation(level_depth + 1, residual, next_level.residual()); - /* Compute the corrected approximation: u = u + error */ - add(solution, ConstVector(residual)); + /* Compute the corrected approximation: u = u + error */ + add(solution, ConstVector(residual)); - auto start_MGC_postSmoothing = std::chrono::high_resolution_clock::now(); + 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); - } + /* ------------- */ + /* 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(end_MGC_postSmoothing - start_MGC_postSmoothing).count(); + auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now(); + t_avg_MGC_postSmoothing_ += + std::chrono::duration(end_MGC_postSmoothing - start_MGC_postSmoothing).count(); + } if (level_depth == 0) { auto end_MGC = std::chrono::high_resolution_clock::now(); diff --git a/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp b/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp index 65fd22af..d6ca6d81 100644 --- a/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp @@ -1,62 +1,69 @@ #include "../../../include/GMGPolar/gmgpolar.h" -void IGMGPolar::multigrid_W_Cycle(int level_depth, Vector solution, Vector rhs, Vector residual) +void IGMGPolar::multigrid_W_Cycle(int level_depth, Vector solution, ConstVector rhs, + Vector 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& level = levels_[level_depth]; - Level& next_level = levels_[level_depth + 1]; + /* ------------------------ */ + /* Solve A * solution = rhs */ + /* ------------------------ */ + if (level_depth == number_of_levels_ - 1) { + /* ------------------------------ */ + /* Coarsest level: solve directly */ + /* ------------------------------ */ + Level& level = levels_[level_depth]; - auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now(); + /* Step 1: Copy rhs in solution */ + Kokkos::deep_copy(solution, rhs); - /* ------------ */ - /* 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(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); - - auto end_MGC_residual = std::chrono::high_resolution_clock::now(); - t_avg_MGC_residual_ += std::chrono::duration(end_MGC_residual - start_MGC_residual).count(); - - /* -------------------------- */ - /* 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(end_MGC_directSolver - start_MGC_directSolver).count(); } else { - /* ------------------------------------------ */ + /* ------------------------------------------------------- */ + /* Multigrid V-cycle on current level: */ + /* presmoothing, coarse-grid correction, and postsmoothing */ + /* ------------------------------------------------------- */ + Level& level = levels_[level_depth]; + Level& 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(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(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); @@ -65,27 +72,34 @@ void IGMGPolar::multigrid_W_Cycle(int level_depth, Vector solution, Vect 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 */ + prolongation(level_depth + 1, residual, next_level.residual()); + + /* Compute the corrected approximation: u = u + error */ + add(solution, ConstVector(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(end_MGC_postSmoothing - start_MGC_postSmoothing).count(); } - /* Interpolate the correction */ - prolongation(level_depth + 1, residual, next_level.residual()); - - /* Compute the corrected approximation: u = u + error */ - add(solution, ConstVector(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(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(end_MGC - start_MGC).count(); diff --git a/src/GMGPolar/setup.cpp b/src/GMGPolar/setup.cpp index ff896ba1..60ec4f0b 100644 --- a/src/GMGPolar/setup.cpp +++ b/src/GMGPolar/setup.cpp @@ -6,7 +6,7 @@ int IGMGPolar::chooseNumberOfLevels(const PolarGrid& finestGrid) const int minAngularDivisions = 4; // Minimum level for Multigrid - const int multigridMinLevel = 2; + const int multigridMinLevel = (extrapolation_ == ExtrapolationType::NONE) ? 1 : 2; // Calculate radial maximum level int radialNodes = finestGrid.nr();