From 0defbc99678fb0900d225ab374bd83eefe024c9d Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Mon, 15 Dec 2025 14:04:18 +0100 Subject: [PATCH 01/13] Add a concept describing a DomainGeometry --- include/InputFunctions/domainGeometry.h | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/include/InputFunctions/domainGeometry.h b/include/InputFunctions/domainGeometry.h index 53899726..81069827 100644 --- a/include/InputFunctions/domainGeometry.h +++ b/include/InputFunctions/domainGeometry.h @@ -29,3 +29,13 @@ class DomainGeometry virtual double dFx_dt(double r, double theta) const = 0; virtual double dFy_dt(double r, double theta) const = 0; }; + +template +concept DomainGeometryConcept = requires(const T geom, double r, double theta) { + { geom.Fx(r, theta) } -> std::convertible_to; + { geom.Fy(r, theta) } -> std::convertible_to; + { geom.dFx_dr(r, theta) } -> std::convertible_to; + { geom.dFy_dr(r, theta) } -> std::convertible_to; + { geom.dFx_dt(r, theta) } -> std::convertible_to; + { geom.dFy_dt(r, theta) } -> std::convertible_to; +}; From 45e599a9159165939d14824de4b596a1fa4f9469 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Mon, 15 Dec 2025 14:04:30 +0100 Subject: [PATCH 02/13] Template GMGPolar on DomainGeometry --- include/GMGPolar/gmgpolar.h | 1 + 1 file changed, 1 insertion(+) diff --git a/include/GMGPolar/gmgpolar.h b/include/GMGPolar/gmgpolar.h index 0bc932a5..6c8956ad 100644 --- a/include/GMGPolar/gmgpolar.h +++ b/include/GMGPolar/gmgpolar.h @@ -24,6 +24,7 @@ class LevelCache; #include "../common/global_definitions.h" #include "test_cases.h" +template class GMGPolar { public: From 880415b82bf4e7b2467251337292ea4e058c68d1 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Mon, 15 Dec 2025 16:28:30 +0100 Subject: [PATCH 03/13] Add templating on class methods --- include/GMGPolar/gmgpolar.h | 13 ++ src/CMakeLists.txt | 22 -- ...licitly_extrapolated_multigrid_F_Cycle.cpp | 6 +- ...licitly_extrapolated_multigrid_V_Cycle.cpp | 6 +- ...licitly_extrapolated_multigrid_W_Cycle.cpp | 6 +- .../MultigridMethods/multigrid_F_Cycle.cpp | 6 +- .../MultigridMethods/multigrid_V_Cycle.cpp | 6 +- .../MultigridMethods/multigrid_W_Cycle.cpp | 6 +- src/GMGPolar/build_rhs_f.cpp | 7 +- src/GMGPolar/gmgpolar.cpp | 208 ++++++++++++------ src/GMGPolar/level_interpolation.cpp | 21 +- src/GMGPolar/setup.cpp | 10 +- src/GMGPolar/solver.cpp | 32 ++- src/GMGPolar/writeToVTK.cpp | 7 +- 14 files changed, 216 insertions(+), 140 deletions(-) diff --git a/include/GMGPolar/gmgpolar.h b/include/GMGPolar/gmgpolar.h index 6c8956ad..cacbd137 100644 --- a/include/GMGPolar/gmgpolar.h +++ b/include/GMGPolar/gmgpolar.h @@ -337,3 +337,16 @@ class GMGPolar double t_avg_MGC_residual_; double t_avg_MGC_directSolver_; }; + +#include "build_rhs_f.h" +#include "gmgpolar_methods.h" +#include "level_interpolation.h" +#include "setup.h" +#include "solver.h" +#include "writeToVTK.h" +#include "MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.h" +#include "MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.h" +#include "MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.h" +#include "MultigridMethods/multigrid_F_Cycle.h" +#include "MultigridMethods/multigrid_V_Cycle.h" +#include "MultigridMethods/multigrid_W_Cycle.h" diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1b7ea608..dab54db3 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -12,26 +12,6 @@ set(POLAR_GRID_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/PolarGrid/polargrid.cpp ) -# file(GLOB_RECURSE GMG_POLAR_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/*.cpp) -set(GMG_POLAR_SOURCES - ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/build_rhs_f.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/gmgpolar.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/level_interpolation.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/setup.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/solver.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/writeToVTK.cpp -) - -# file(GLOB_RECURSE MULTIGRID_METHODS_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/MultigridMethods/*.cpp) -set(MULTIGRID_METHODS_SOURCES - ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp -) - #file(GLOB_RECURSE LEVEL_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/Level/*.cpp) set(LEVEL_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/Level/level.cpp @@ -162,8 +142,6 @@ set(CONFIG_PARSER_SOURCES # Create the main library add_library(GMGPolarLib STATIC ${POLAR_GRID_SOURCES} - ${GMG_POLAR_SOURCES} - ${MULTIGRID_METHODS_SOURCES} ${LEVEL_SOURCES} ${STENCIL_SOURCES} ${INTERPOLATION_SOURCES} diff --git a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp index 63b00848..a9b86cf5 100644 --- a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp @@ -1,6 +1,6 @@ -#include "../../../include/GMGPolar/gmgpolar.h" -void GMGPolar::implicitlyExtrapolatedMultigrid_F_Cycle(const int level_depth, Vector solution, +template +void GMGPolar::implicitlyExtrapolatedMultigrid_F_Cycle(const int level_depth, Vector solution, Vector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); @@ -123,4 +123,4 @@ void GMGPolar::implicitlyExtrapolatedMultigrid_F_Cycle(const int level_depth, Ve auto end_MGC = std::chrono::high_resolution_clock::now(); t_avg_MGC_total_ += std::chrono::duration(end_MGC - start_MGC).count(); } -} \ No newline at end of file +} diff --git a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp index 8ef84850..a19eb16c 100644 --- a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp @@ -1,6 +1,6 @@ -#include "../../../include/GMGPolar/gmgpolar.h" -void GMGPolar::implicitlyExtrapolatedMultigrid_V_Cycle(const int level_depth, Vector solution, +template +void GMGPolar::implicitlyExtrapolatedMultigrid_V_Cycle(const int level_depth, Vector solution, Vector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); @@ -122,4 +122,4 @@ void GMGPolar::implicitlyExtrapolatedMultigrid_V_Cycle(const int level_depth, Ve auto end_MGC = std::chrono::high_resolution_clock::now(); t_avg_MGC_total_ += std::chrono::duration(end_MGC - start_MGC).count(); } -} \ No newline at end of file +} diff --git a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp index cffcccc7..d33c1f8e 100644 --- a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp @@ -1,6 +1,6 @@ -#include "../../../include/GMGPolar/gmgpolar.h" -void GMGPolar::implicitlyExtrapolatedMultigrid_W_Cycle(const int level_depth, Vector solution, +template +void GMGPolar::implicitlyExtrapolatedMultigrid_W_Cycle(const int level_depth, Vector solution, Vector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); @@ -123,4 +123,4 @@ void GMGPolar::implicitlyExtrapolatedMultigrid_W_Cycle(const int level_depth, Ve auto end_MGC = std::chrono::high_resolution_clock::now(); t_avg_MGC_total_ += std::chrono::duration(end_MGC - start_MGC).count(); } -} \ No newline at end of file +} diff --git a/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp b/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp index aa9ea2d5..df3444a9 100644 --- a/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp @@ -1,6 +1,6 @@ -#include "../../../include/GMGPolar/gmgpolar.h" -void GMGPolar::multigrid_F_Cycle(const int level_depth, Vector solution, Vector rhs, +template +void GMGPolar::multigrid_F_Cycle(const int level_depth, Vector solution, Vector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); @@ -91,4 +91,4 @@ void GMGPolar::multigrid_F_Cycle(const int level_depth, Vector solution, auto end_MGC = std::chrono::high_resolution_clock::now(); t_avg_MGC_total_ += std::chrono::duration(end_MGC - start_MGC).count(); } -} \ No newline at end of file +} diff --git a/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp b/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp index 3dbec4fd..1159c753 100644 --- a/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp @@ -1,6 +1,6 @@ -#include "../../../include/GMGPolar/gmgpolar.h" -void GMGPolar::multigrid_V_Cycle(const int level_depth, Vector solution, Vector rhs, +template +void GMGPolar::multigrid_V_Cycle(const int level_depth, Vector solution, Vector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); @@ -90,4 +90,4 @@ void GMGPolar::multigrid_V_Cycle(const int level_depth, Vector solution, auto end_MGC = std::chrono::high_resolution_clock::now(); t_avg_MGC_total_ += std::chrono::duration(end_MGC - start_MGC).count(); } -} \ No newline at end of file +} diff --git a/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp b/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp index f9165b29..e19628a0 100644 --- a/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp @@ -1,6 +1,6 @@ -#include "../../../include/GMGPolar/gmgpolar.h" -void GMGPolar::multigrid_W_Cycle(const int level_depth, Vector solution, Vector rhs, +template +void GMGPolar::multigrid_W_Cycle(const int level_depth, Vector solution, Vector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); @@ -91,4 +91,4 @@ void GMGPolar::multigrid_W_Cycle(const int level_depth, Vector solution, auto end_MGC = std::chrono::high_resolution_clock::now(); t_avg_MGC_total_ += std::chrono::duration(end_MGC - start_MGC).count(); } -} \ No newline at end of file +} diff --git a/src/GMGPolar/build_rhs_f.cpp b/src/GMGPolar/build_rhs_f.cpp index c97b640e..a3c1287e 100644 --- a/src/GMGPolar/build_rhs_f.cpp +++ b/src/GMGPolar/build_rhs_f.cpp @@ -1,6 +1,6 @@ -#include "../../include/GMGPolar/gmgpolar.h" -void GMGPolar::build_rhs_f(const Level& level, Vector rhs_f, const BoundaryConditions& boundary_conditions, +template +void GMGPolar::build_rhs_f(const Level& level, Vector rhs_f, const BoundaryConditions& boundary_conditions, const SourceTerm& source_term) { const PolarGrid& grid = level.grid(); @@ -52,7 +52,8 @@ void GMGPolar::build_rhs_f(const Level& level, Vector rhs_f, const Bound } } -void GMGPolar::discretize_rhs_f(const Level& level, Vector rhs_f) +template +void GMGPolar::discretize_rhs_f(const Level& level, Vector rhs_f) { const PolarGrid& grid = level.grid(); assert(rhs_f.size() == static_cast(grid.numberOfNodes())); diff --git a/src/GMGPolar/gmgpolar.cpp b/src/GMGPolar/gmgpolar.cpp index 78ab4b8b..8f8624a7 100644 --- a/src/GMGPolar/gmgpolar.cpp +++ b/src/GMGPolar/gmgpolar.cpp @@ -1,9 +1,9 @@ -#include "../../include/GMGPolar/gmgpolar.h" /* ---------------------------------------------------------------------- */ /* Constructor & Initialization */ /* ---------------------------------------------------------------------- */ -GMGPolar::GMGPolar(const PolarGrid& grid, const DomainGeometry& domain_geometry, +template +GMGPolar::GMGPolar(const PolarGrid& grid, const DomainGeometry& domain_geometry, const DensityProfileCoefficients& density_profile_coefficients) : grid_(grid) , domain_geometry_(domain_geometry) @@ -47,7 +47,8 @@ GMGPolar::GMGPolar(const PolarGrid& grid, const DomainGeometry& domain_geometry, LIKWID_REGISTER("Solve"); } -void GMGPolar::setSolution(const ExactSolution* exact_solution) +template +void GMGPolar::setSolution(const ExactSolution* exact_solution) { exact_solution_ = exact_solution; } @@ -55,20 +56,24 @@ void GMGPolar::setSolution(const ExactSolution* exact_solution) /* ---------------------------------------------------------------------- */ /* General output & visualization */ /* ---------------------------------------------------------------------- */ -int GMGPolar::verbose() const +template +int GMGPolar::verbose() const { return verbose_; } -void GMGPolar::verbose(int verbose) +template +void GMGPolar::verbose(int verbose) { verbose_ = verbose; } -bool GMGPolar::paraview() const +template +bool GMGPolar::paraview() const { return paraview_; } -void GMGPolar::paraview(bool paraview) +template +void GMGPolar::paraview(bool paraview) { paraview_ = paraview; } @@ -76,20 +81,24 @@ void GMGPolar::paraview(bool paraview) /* ---------------------------------------------------------------------- */ /* Parallelization & threading */ /* ---------------------------------------------------------------------- */ -int GMGPolar::maxOpenMPThreads() const +template +int GMGPolar::maxOpenMPThreads() const { return max_omp_threads_; } -void GMGPolar::maxOpenMPThreads(int max_omp_threads) +template +void GMGPolar::maxOpenMPThreads(int max_omp_threads) { max_omp_threads_ = max_omp_threads; } -double GMGPolar::threadReductionFactor() const +template +double GMGPolar::threadReductionFactor() const { return thread_reduction_factor_; } -void GMGPolar::threadReductionFactor(double thread_reduction_factor) +template +void GMGPolar::threadReductionFactor(double thread_reduction_factor) { thread_reduction_factor_ = thread_reduction_factor; } @@ -97,38 +106,46 @@ void GMGPolar::threadReductionFactor(double thread_reduction_factor) /* ---------------------------------------------------------------------- */ /* Numerical method options */ /* ---------------------------------------------------------------------- */ -bool GMGPolar::DirBC_Interior() const +template +bool GMGPolar::DirBC_Interior() const { return DirBC_Interior_; } -void GMGPolar::DirBC_Interior(bool DirBC_Interior) +template +void GMGPolar::DirBC_Interior(bool DirBC_Interior) { DirBC_Interior_ = DirBC_Interior; } -StencilDistributionMethod GMGPolar::stencilDistributionMethod() const +template +StencilDistributionMethod GMGPolar::stencilDistributionMethod() const { return stencil_distribution_method_; } -void GMGPolar::stencilDistributionMethod(StencilDistributionMethod stencil_distribution_method) +template +void GMGPolar::stencilDistributionMethod(StencilDistributionMethod stencil_distribution_method) { stencil_distribution_method_ = stencil_distribution_method; } -bool GMGPolar::cacheDensityProfileCoefficients() const +template +bool GMGPolar::cacheDensityProfileCoefficients() const { return cache_density_profile_coefficients_; } -void GMGPolar::cacheDensityProfileCoefficients(bool cache_density_profile_coefficients) +template +void GMGPolar::cacheDensityProfileCoefficients(bool cache_density_profile_coefficients) { cache_density_profile_coefficients_ = cache_density_profile_coefficients; } -bool GMGPolar::cacheDomainGeometry() const +template +bool GMGPolar::cacheDomainGeometry() const { return cache_domain_geometry_; } -void GMGPolar::cacheDomainGeometry(bool cache_domain_geometry) +template +void GMGPolar::cacheDomainGeometry(bool cache_domain_geometry) { cache_domain_geometry_ = cache_domain_geometry; } @@ -136,74 +153,90 @@ void GMGPolar::cacheDomainGeometry(bool cache_domain_geometry) /* ---------------------------------------------------------------------- */ /* Multigrid controls */ /* ---------------------------------------------------------------------- */ -ExtrapolationType GMGPolar::extrapolation() const +template +ExtrapolationType GMGPolar::extrapolation() const { return extrapolation_; } -void GMGPolar::extrapolation(ExtrapolationType extrapolation) +template +void GMGPolar::extrapolation(ExtrapolationType extrapolation) { extrapolation_ = extrapolation; } -int GMGPolar::maxLevels() const +template +int GMGPolar::maxLevels() const { return max_levels_; } -void GMGPolar::maxLevels(int max_levels) +template +void GMGPolar::maxLevels(int max_levels) { max_levels_ = max_levels; } -MultigridCycleType GMGPolar::multigridCycle() const +template +MultigridCycleType GMGPolar::multigridCycle() const { return multigrid_cycle_; } -void GMGPolar::multigridCycle(MultigridCycleType multigrid_cycle) +template +void GMGPolar::multigridCycle(MultigridCycleType multigrid_cycle) { multigrid_cycle_ = multigrid_cycle; } -int GMGPolar::preSmoothingSteps() const +template +int GMGPolar::preSmoothingSteps() const { return pre_smoothing_steps_; } -void GMGPolar::preSmoothingSteps(int pre_smoothing_steps) +template +void GMGPolar::preSmoothingSteps(int pre_smoothing_steps) { pre_smoothing_steps_ = pre_smoothing_steps; } -int GMGPolar::postSmoothingSteps() const +template +int GMGPolar::postSmoothingSteps() const { return post_smoothing_steps_; } -void GMGPolar::postSmoothingSteps(int post_smoothing_steps) +template +void GMGPolar::postSmoothingSteps(int post_smoothing_steps) { post_smoothing_steps_ = post_smoothing_steps; } -bool GMGPolar::FMG() const +template +bool GMGPolar::FMG() const { return FMG_; } -void GMGPolar::FMG(bool FMG) +template +void GMGPolar::FMG(bool FMG) { FMG_ = FMG; } -int GMGPolar::FMG_iterations() const +template +int GMGPolar::FMG_iterations() const { return FMG_iterations_; } -void GMGPolar::FMG_iterations(int FMG_iterations) +template +void GMGPolar::FMG_iterations(int FMG_iterations) { FMG_iterations_ = FMG_iterations; } -MultigridCycleType GMGPolar::FMG_cycle() const +template +MultigridCycleType GMGPolar::FMG_cycle() const { return FMG_cycle_; } -void GMGPolar::FMG_cycle(MultigridCycleType FMG_cycle) +template +void GMGPolar::FMG_cycle(MultigridCycleType FMG_cycle) { FMG_cycle_ = FMG_cycle; } @@ -212,38 +245,46 @@ void GMGPolar::FMG_cycle(MultigridCycleType FMG_cycle) /* Iterative solver termination */ /* ---------------------------------------------------------------------- */ -int GMGPolar::maxIterations() const +template +int GMGPolar::maxIterations() const { return max_iterations_; } -void GMGPolar::maxIterations(int maxIterations) +template +void GMGPolar::maxIterations(int maxIterations) { max_iterations_ = maxIterations; } -ResidualNormType GMGPolar::residualNormType() const +template +ResidualNormType GMGPolar::residualNormType() const { return residual_norm_type_; } -void GMGPolar::residualNormType(ResidualNormType residualNormType) +template +void GMGPolar::residualNormType(ResidualNormType residualNormType) { residual_norm_type_ = residualNormType; } -std::optional GMGPolar::absoluteTolerance() const +template +std::optional GMGPolar::absoluteTolerance() const { return absolute_tolerance_; } -void GMGPolar::absoluteTolerance(std::optional tol) +template +void GMGPolar::absoluteTolerance(std::optional tol) { absolute_tolerance_ = tol.has_value() && tol.value() >= 0.0 ? tol : std::nullopt; } -std::optional GMGPolar::relativeTolerance() const +template +std::optional GMGPolar::relativeTolerance() const { return relative_tolerance_; } -void GMGPolar::relativeTolerance(std::optional tol) +template +void GMGPolar::relativeTolerance(std::optional tol) { relative_tolerance_ = tol.has_value() && tol.value() >= 0.0 ? tol : std::nullopt; } @@ -251,18 +292,21 @@ void GMGPolar::relativeTolerance(std::optional tol) /* ---------------------------------------------------------------------- */ /* Solution & Grid Access */ /* ---------------------------------------------------------------------- */ -Vector GMGPolar::solution() +template +Vector GMGPolar::solution() { int level_depth = 0; return levels_[level_depth].solution(); } -ConstVector GMGPolar::solution() const +template +ConstVector GMGPolar::solution() const { int level_depth = 0; return levels_[level_depth].solution(); } -const PolarGrid& GMGPolar::grid() const +template +const PolarGrid& GMGPolar::grid() const { return grid_; } @@ -270,23 +314,28 @@ const PolarGrid& GMGPolar::grid() const /* ---------------------------------------------------------------------- */ /* Setup timings */ /* ---------------------------------------------------------------------- */ -double GMGPolar::timeSetupTotal() const +template +double GMGPolar::timeSetupTotal() const { return t_setup_total_; } -double GMGPolar::timeSetupCreateLevels() const +template +double GMGPolar::timeSetupCreateLevels() const { return t_setup_createLevels_; } -double GMGPolar::timeSetupRHS() const +template +double GMGPolar::timeSetupRHS() const { return t_setup_rhs_; } -double GMGPolar::timeSetupSmoother() const +template +double GMGPolar::timeSetupSmoother() const { return t_setup_smoother_; } -double GMGPolar::timeSetupDirectSolver() const +template +double GMGPolar::timeSetupDirectSolver() const { return t_setup_directSolver_; } @@ -294,23 +343,28 @@ double GMGPolar::timeSetupDirectSolver() const /* ---------------------------------------------------------------------- */ /* Solve timings */ /* ---------------------------------------------------------------------- */ -double GMGPolar::timeSolveTotal() const +template +double GMGPolar::timeSolveTotal() const { return t_solve_total_; } -double GMGPolar::timeSolveInitialApproximation() const +template +double GMGPolar::timeSolveInitialApproximation() const { return t_solve_initial_approximation_; } -double GMGPolar::timeSolveMultigridIterations() const +template +double GMGPolar::timeSolveMultigridIterations() const { return t_solve_multigrid_iterations_; } -double GMGPolar::timeCheckConvergence() const +template +double GMGPolar::timeCheckConvergence() const { return t_check_convergence_; } -double GMGPolar::timeCheckExactError() const +template +double GMGPolar::timeCheckExactError() const { return t_check_exact_error_; } @@ -318,23 +372,28 @@ double GMGPolar::timeCheckExactError() const /* ---------------------------------------------------------------------- */ /* Average Multigrid Cycle timings */ /* ---------------------------------------------------------------------- */ -double GMGPolar::timeAvgMGCTotal() const +template +double GMGPolar::timeAvgMGCTotal() const { return t_avg_MGC_total_; } -double GMGPolar::timeAvgMGCPreSmoothing() const +template +double GMGPolar::timeAvgMGCPreSmoothing() const { return t_avg_MGC_preSmoothing_; } -double GMGPolar::timeAvgMGCPostSmoothing() const +template +double GMGPolar::timeAvgMGCPostSmoothing() const { return t_avg_MGC_postSmoothing_; } -double GMGPolar::timeAvgMGCResidual() const +template +double GMGPolar::timeAvgMGCResidual() const { return t_avg_MGC_residual_; } -double GMGPolar::timeAvgMGCDirectSolver() const +template +double GMGPolar::timeAvgMGCDirectSolver() const { return t_avg_MGC_directSolver_; } @@ -343,14 +402,16 @@ double GMGPolar::timeAvgMGCDirectSolver() const /* Reset timings */ /* ---------------------------------------------------------------------- */ -void GMGPolar::resetAllTimings() +template +void GMGPolar::resetAllTimings() { resetSetupPhaseTimings(); resetSolvePhaseTimings(); resetAvgMultigridCycleTimings(); } -void GMGPolar::resetSetupPhaseTimings() +template +void GMGPolar::resetSetupPhaseTimings() { t_setup_total_ = 0.0; t_setup_createLevels_ = 0.0; @@ -359,7 +420,8 @@ void GMGPolar::resetSetupPhaseTimings() t_setup_directSolver_ = 0.0; } -void GMGPolar::resetSolvePhaseTimings() +template +void GMGPolar::resetSolvePhaseTimings() { t_solve_total_ = 0.0; t_solve_initial_approximation_ = 0.0; @@ -368,7 +430,8 @@ void GMGPolar::resetSolvePhaseTimings() t_check_exact_error_ = 0.0; } -void GMGPolar::resetAvgMultigridCycleTimings() +template +void GMGPolar::resetAvgMultigridCycleTimings() { t_avg_MGC_total_ = 0.0; t_avg_MGC_preSmoothing_ = 0.0; @@ -381,7 +444,8 @@ void GMGPolar::resetAvgMultigridCycleTimings() /* Diagnostics & statistics */ /* ---------------------------------------------------------------------- */ // Print timing breakdown for setup, smoothing, coarse solve, etc. -void GMGPolar::printTimings() const +template +void GMGPolar::printTimings() const { // t_setup_rhs_ is neither included in t_setup_total_ and t_solve_total_. std::cout << "\n------------------" << std::endl; @@ -410,26 +474,30 @@ void GMGPolar::printTimings() const } // Number of iterations taken by last solve. -int GMGPolar::numberOfIterations() const +template +int GMGPolar::numberOfIterations() const { return number_of_iterations_; } // Mean residual reduction factor per iteration. -double GMGPolar::meanResidualReductionFactor() const +template +double GMGPolar::meanResidualReductionFactor() const { return mean_residual_reduction_factor_; } // Error norms (only available if exact solution was set). -std::optional GMGPolar::exactErrorWeightedEuclidean() const +template +std::optional GMGPolar::exactErrorWeightedEuclidean() const { if (exact_solution_) { return exact_errors_.back().first; } return std::nullopt; } -std::optional GMGPolar::exactErrorInfinity() const +template +std::optional GMGPolar::exactErrorInfinity() const { if (exact_solution_) { return exact_errors_.back().second; diff --git a/src/GMGPolar/level_interpolation.cpp b/src/GMGPolar/level_interpolation.cpp index de6aa8b8..86a3fa2b 100644 --- a/src/GMGPolar/level_interpolation.cpp +++ b/src/GMGPolar/level_interpolation.cpp @@ -1,6 +1,6 @@ -#include "../../include/GMGPolar/gmgpolar.h" -void GMGPolar::prolongation(const int current_level, Vector result, ConstVector x) const +template +void GMGPolar::prolongation(const int current_level, Vector result, ConstVector x) const { assert(current_level < number_of_levels_ && 1 <= current_level); if (!interpolation_) @@ -9,7 +9,8 @@ void GMGPolar::prolongation(const int current_level, Vector result, Cons interpolation_->applyProlongation(levels_[current_level], levels_[current_level - 1], result, x); } -void GMGPolar::restriction(const int current_level, Vector result, ConstVector x) const +template +void GMGPolar::restriction(const int current_level, Vector result, ConstVector x) const { assert(current_level < number_of_levels_ - 1 && 0 <= current_level); if (!interpolation_) @@ -18,7 +19,8 @@ void GMGPolar::restriction(const int current_level, Vector result, Const interpolation_->applyRestriction(levels_[current_level], levels_[current_level + 1], result, x); } -void GMGPolar::injection(const int current_level, Vector result, ConstVector x) const +template +void GMGPolar::injection(const int current_level, Vector result, ConstVector x) const { assert(current_level < number_of_levels_ - 1 && 0 <= current_level); if (!interpolation_) @@ -27,7 +29,8 @@ void GMGPolar::injection(const int current_level, Vector result, ConstVe interpolation_->applyInjection(levels_[current_level], levels_[current_level + 1], result, x); } -void GMGPolar::extrapolatedProlongation(const int current_level, Vector result, ConstVector x) const +template +void GMGPolar::extrapolatedProlongation(const int current_level, Vector result, ConstVector x) const { assert(current_level < number_of_levels_ && 1 <= current_level); if (!interpolation_) @@ -36,7 +39,8 @@ void GMGPolar::extrapolatedProlongation(const int current_level, Vector interpolation_->applyExtrapolatedProlongation(levels_[current_level], levels_[current_level - 1], result, x); } -void GMGPolar::extrapolatedRestriction(const int current_level, Vector result, ConstVector x) const +template +void GMGPolar::extrapolatedRestriction(const int current_level, Vector result, ConstVector x) const { assert(current_level < number_of_levels_ - 1 && 0 <= current_level); if (!interpolation_) @@ -45,11 +49,12 @@ void GMGPolar::extrapolatedRestriction(const int current_level, Vector r interpolation_->applyExtrapolatedRestriction(levels_[current_level], levels_[current_level + 1], result, x); } -void GMGPolar::FMGInterpolation(const int current_level, Vector result, ConstVector x) const +template +void GMGPolar::FMGInterpolation(const int current_level, Vector result, ConstVector x) const { assert(current_level < number_of_levels_ && 1 <= current_level); if (!interpolation_) throw std::runtime_error("Interpolation not initialized."); interpolation_->applyFMGInterpolation(levels_[current_level], levels_[current_level - 1], result, x); -} \ No newline at end of file +} diff --git a/src/GMGPolar/setup.cpp b/src/GMGPolar/setup.cpp index 9aeb651b..a1435a4a 100644 --- a/src/GMGPolar/setup.cpp +++ b/src/GMGPolar/setup.cpp @@ -1,6 +1,6 @@ -#include "../../include/GMGPolar/gmgpolar.h" -void GMGPolar::setup() +template +void GMGPolar::setup() { LIKWID_START("Setup"); auto start_setup = std::chrono::high_resolution_clock::now(); @@ -172,7 +172,8 @@ void GMGPolar::setup() LIKWID_STOP("Setup"); } -int GMGPolar::chooseNumberOfLevels(const PolarGrid& finestGrid) +template +int GMGPolar::chooseNumberOfLevels(const PolarGrid& finestGrid) { const int minRadialNodes = 5; const int minAngularDivisions = 4; @@ -215,7 +216,8 @@ int GMGPolar::chooseNumberOfLevels(const PolarGrid& finestGrid) return levels; } -void GMGPolar::printSettings() const +template +void GMGPolar::printSettings() const { std::cout << "------------------------------\n"; diff --git a/src/GMGPolar/solver.cpp b/src/GMGPolar/solver.cpp index 2bcf800a..4bf83147 100644 --- a/src/GMGPolar/solver.cpp +++ b/src/GMGPolar/solver.cpp @@ -1,12 +1,11 @@ -#include "../../include/GMGPolar/gmgpolar.h" - #include // ============================================================================= // Main Solver Routine // ============================================================================= -void GMGPolar::solve(const BoundaryConditions& boundary_conditions, const SourceTerm& source_term) +template +void GMGPolar::solve(const BoundaryConditions& boundary_conditions, const SourceTerm& source_term) { auto start_setup_rhs = std::chrono::high_resolution_clock::now(); @@ -185,7 +184,8 @@ void GMGPolar::solve(const BoundaryConditions& boundary_conditions, const Source // Solution Initialization // ============================================================================= -void GMGPolar::initializeSolution() +template +void GMGPolar::initializeSolution() { if (!FMG_) { int start_level_depth = 0; @@ -266,7 +266,8 @@ void GMGPolar::initializeSolution() // Residual Handling Functions // ============================================================================= -double GMGPolar::residualNorm(const ResidualNormType& norm_type, const Level& level, ConstVector residual) const +template +double GMGPolar::residualNorm(const ResidualNormType& norm_type, const Level& level, ConstVector residual) const { switch (norm_type) { case ResidualNormType::EUCLIDEAN: @@ -280,7 +281,8 @@ double GMGPolar::residualNorm(const ResidualNormType& norm_type, const Level& le } } -void GMGPolar::updateResidualNorms(Level& level, int iteration, double& initial_residual_norm, +template +void GMGPolar::updateResidualNorms(Level& level, int iteration, double& initial_residual_norm, double& current_residual_norm, double& current_relative_residual_norm) { level.computeResidual(level.residual(), level.rhs(), level.solution()); @@ -313,7 +315,8 @@ void GMGPolar::updateResidualNorms(Level& level, int iteration, double& initial_ } } -void GMGPolar::extrapolatedResidual(const int current_level, Vector residual, +template +void GMGPolar::extrapolatedResidual(const int current_level, Vector residual, ConstVector residual_next_level) { const PolarGrid& fineGrid = levels_[current_level].grid(); @@ -368,7 +371,8 @@ void GMGPolar::extrapolatedResidual(const int current_level, Vector resi // Convergence and Error Analysis Functions // ============================================================================= -bool GMGPolar::converged(double residual_norm, double relative_residual_norm) +template +bool GMGPolar::converged(double residual_norm, double relative_residual_norm) { if (relative_tolerance_.has_value()) { if (!(relative_residual_norm > relative_tolerance_.value())) { @@ -383,14 +387,16 @@ bool GMGPolar::converged(double residual_norm, double relative_residual_norm) return false; } -void GMGPolar::evaluateExactError(Level& level, const ExactSolution& exact_solution) +template +void GMGPolar::evaluateExactError(Level& level, const ExactSolution& exact_solution) { // Compute the weighted L2 norm and infinity norm of the error between the numerical and exact solution. // The results are stored as a pair: (weighted L2 error, infinity error). exact_errors_.push_back(computeExactError(level, level.solution(), level.residual(), exact_solution)); } -std::pair GMGPolar::computeExactError(Level& level, ConstVector solution, Vector error, +template +std::pair GMGPolar::computeExactError(Level& level, ConstVector solution, Vector error, const ExactSolution& exact_solution) { const PolarGrid& grid = level.grid(); @@ -431,7 +437,8 @@ std::pair GMGPolar::computeExactError(Level& level, ConstVector< // Output and Logging Functions // ============================================================================= -void GMGPolar::printIterationHeader(const ExactSolution* exact_solution) +template +void GMGPolar::printIterationHeader(const ExactSolution* exact_solution) { if (verbose_ <= 0) return; @@ -457,7 +464,8 @@ void GMGPolar::printIterationHeader(const ExactSolution* exact_solution) std::cout << std::right << std::setfill(' '); } -void GMGPolar::printIterationInfo(int iteration, double current_residual_norm, double current_relative_residual_norm, +template +void GMGPolar::printIterationInfo(int iteration, double current_residual_norm, double current_relative_residual_norm, const ExactSolution* exact_solution) { if (verbose_ <= 0) diff --git a/src/GMGPolar/writeToVTK.cpp b/src/GMGPolar/writeToVTK.cpp index 48d1fa32..41280a0e 100644 --- a/src/GMGPolar/writeToVTK.cpp +++ b/src/GMGPolar/writeToVTK.cpp @@ -1,6 +1,6 @@ -#include "../../include/GMGPolar/gmgpolar.h" -void GMGPolar::writeToVTK(const std::filesystem::path& file_path, const PolarGrid& grid) +template +void GMGPolar::writeToVTK(const std::filesystem::path& file_path, const PolarGrid& grid) { const auto filename = file_path.stem().string() + ".vtu"; @@ -58,7 +58,8 @@ void GMGPolar::writeToVTK(const std::filesystem::path& file_path, const PolarGri << "\n"; } -void GMGPolar::writeToVTK(const std::filesystem::path& file_path, const Level& level, ConstVector grid_function) +template +void GMGPolar::writeToVTK(const std::filesystem::path& file_path, const Level& level, ConstVector grid_function) { const PolarGrid& grid = level.grid(); const LevelCache& level_cache = level.levelCache(); From fd3152c43c89482aad97bd10ed56bd8d5ff7fa5e Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 17 Dec 2025 13:40:46 +0100 Subject: [PATCH 04/13] Modify concept so abstract class does not match --- include/InputFunctions/domainGeometry.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/InputFunctions/domainGeometry.h b/include/InputFunctions/domainGeometry.h index 81069827..7a7f83fb 100644 --- a/include/InputFunctions/domainGeometry.h +++ b/include/InputFunctions/domainGeometry.h @@ -31,7 +31,7 @@ class DomainGeometry }; template -concept DomainGeometryConcept = requires(const T geom, double r, double theta) { +concept DomainGeometryConcept = !std::same_as && requires(const T geom, double r, double theta) { { geom.Fx(r, theta) } -> std::convertible_to; { geom.Fy(r, theta) } -> std::convertible_to; { geom.dFx_dr(r, theta) } -> std::convertible_to; From 7964abd8dc1ce5709aa38c6c64d936aafc5e1feb Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 17 Dec 2025 13:41:27 +0100 Subject: [PATCH 05/13] Clang format --- .../implicitly_extrapolated_multigrid_F_Cycle.cpp | 4 ++-- .../implicitly_extrapolated_multigrid_V_Cycle.cpp | 4 ++-- .../implicitly_extrapolated_multigrid_W_Cycle.cpp | 4 ++-- src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp | 4 ++-- src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp | 4 ++-- src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp | 4 ++-- 6 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp index a9b86cf5..72ad172d 100644 --- a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp @@ -1,7 +1,7 @@ -template +template void GMGPolar::implicitlyExtrapolatedMultigrid_F_Cycle(const int level_depth, Vector solution, - Vector rhs, Vector residual) + Vector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); diff --git a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp index a19eb16c..dadfe565 100644 --- a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp @@ -1,7 +1,7 @@ -template +template void GMGPolar::implicitlyExtrapolatedMultigrid_V_Cycle(const int level_depth, Vector solution, - Vector rhs, Vector residual) + Vector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); diff --git a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp index d33c1f8e..8fe51fad 100644 --- a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp @@ -1,7 +1,7 @@ -template +template void GMGPolar::implicitlyExtrapolatedMultigrid_W_Cycle(const int level_depth, Vector solution, - Vector rhs, Vector residual) + Vector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); diff --git a/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp b/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp index df3444a9..a0bc1396 100644 --- a/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp @@ -1,7 +1,7 @@ -template +template void GMGPolar::multigrid_F_Cycle(const int level_depth, Vector solution, Vector rhs, - Vector residual) + Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); diff --git a/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp b/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp index 1159c753..a47a197a 100644 --- a/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp @@ -1,7 +1,7 @@ -template +template void GMGPolar::multigrid_V_Cycle(const int level_depth, Vector solution, Vector rhs, - Vector residual) + Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); diff --git a/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp b/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp index e19628a0..52b95069 100644 --- a/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp @@ -1,7 +1,7 @@ -template +template void GMGPolar::multigrid_W_Cycle(const int level_depth, Vector solution, Vector rhs, - Vector residual) + Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); From 9b230dac6a3c7c2ecf1820a387a93866867c352a Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 17 Dec 2025 14:16:46 +0100 Subject: [PATCH 06/13] Use inheritance for GMGPolar to allow main declaration with ConfigParser --- include/GMGPolar/build_rhs_f.h | 141 ++++++++ include/GMGPolar/gmgpolar.h | 304 +--------------- include/GMGPolar/igmgpolar.h | 336 ++++++++++++++++++ include/GMGPolar/setup.h | 173 +++++++++ .../GMGPolar/writeToVTK.h | 1 + src/CMakeLists.txt | 21 ++ ...licitly_extrapolated_multigrid_F_Cycle.cpp | 8 +- ...licitly_extrapolated_multigrid_V_Cycle.cpp | 8 +- ...licitly_extrapolated_multigrid_W_Cycle.cpp | 8 +- .../MultigridMethods/multigrid_F_Cycle.cpp | 8 +- .../MultigridMethods/multigrid_V_Cycle.cpp | 8 +- .../MultigridMethods/multigrid_W_Cycle.cpp | 8 +- src/GMGPolar/build_rhs_f.cpp | 144 +------- src/GMGPolar/gmgpolar.cpp | 209 ++++------- src/GMGPolar/level_interpolation.cpp | 21 +- src/GMGPolar/setup.cpp | 180 +--------- src/GMGPolar/solver.cpp | 32 +- 17 files changed, 804 insertions(+), 806 deletions(-) create mode 100644 include/GMGPolar/build_rhs_f.h create mode 100644 include/GMGPolar/igmgpolar.h create mode 100644 include/GMGPolar/setup.h rename src/GMGPolar/writeToVTK.cpp => include/GMGPolar/writeToVTK.h (99%) diff --git a/include/GMGPolar/build_rhs_f.h b/include/GMGPolar/build_rhs_f.h new file mode 100644 index 00000000..81342f7a --- /dev/null +++ b/include/GMGPolar/build_rhs_f.h @@ -0,0 +1,141 @@ +#include "../../include/GMGPolar/gmgpolar.h" + +template +void GMGPolar::discretize_rhs_f(const Level& level, Vector rhs_f) +{ + const PolarGrid& grid = level.grid(); + assert(rhs_f.size() == static_cast(grid.numberOfNodes())); + + if (level.levelCache().cacheDomainGeometry()) { + /* DomainGeometry is cached */ + const auto& detDF_cache = level.levelCache().detDF(); +#pragma omp parallel + { +// ---------------------------------------------- // +// Discretize rhs values (circular index section) // +// ---------------------------------------------- // +#pragma omp for nowait + for (int i_r = 0; i_r < grid.numberSmootherCircles(); i_r++) { + double r = grid.radius(i_r); + for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { + double theta = grid.theta(i_theta); + if ((0 < i_r && i_r < grid.nr() - 1) || (i_r == 0 && !DirBC_Interior_)) { + double h1 = (i_r == 0) ? 2.0 * grid.radius(0) : grid.radialSpacing(i_r - 1); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + const double detDF = detDF_cache[grid.index(i_r, i_theta)]; + rhs_f[grid.index(i_r, i_theta)] *= 0.25 * (h1 + h2) * (k1 + k2) * fabs(detDF); + } + else if (i_r == 0 && DirBC_Interior_) { + rhs_f[grid.index(i_r, i_theta)] *= 1.0; + } + else if (i_r == grid.nr() - 1) { + rhs_f[grid.index(i_r, i_theta)] *= 1.0; + } + } + } + +// -------------------------------------------- // +// Discretize rhs values (radial index section) // +// -------------------------------------------- // +#pragma omp for nowait + for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { + double theta = grid.theta(i_theta); + for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { + double r = grid.radius(i_r); + if ((0 < i_r && i_r < grid.nr() - 1) || (i_r == 0 && !DirBC_Interior_)) { + double h1 = (i_r == 0) ? 2.0 * grid.radius(0) : grid.radialSpacing(i_r - 1); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + const double detDF = detDF_cache[grid.index(i_r, i_theta)]; + rhs_f[grid.index(i_r, i_theta)] *= 0.25 * (h1 + h2) * (k1 + k2) * fabs(detDF); + } + else if (i_r == 0 && DirBC_Interior_) { + rhs_f[grid.index(i_r, i_theta)] *= 1.0; + } + else if (i_r == grid.nr() - 1) { + rhs_f[grid.index(i_r, i_theta)] *= 1.0; + } + } + } + } + } + else { + /* DomainGeometry is not cached */ + +#pragma omp parallel + { +// ---------------------------------------------- // +// Discretize rhs values (circular index section) // +// ---------------------------------------------- // +#pragma omp for nowait + for (int i_r = 0; i_r < grid.numberSmootherCircles(); i_r++) { + double r = grid.radius(i_r); + for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { + double theta = grid.theta(i_theta); + + if ((0 < i_r && i_r < grid.nr() - 1) || (i_r == 0 && !DirBC_Interior_)) { + double h1 = (i_r == 0) ? 2.0 * grid.radius(0) : grid.radialSpacing(i_r - 1); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + /* Calculate the elements of the Jacobian matrix for the transformation mapping */ + /* The Jacobian matrix is: */ + /* [Jrr, Jrt] */ + /* [Jtr, Jtt] */ + double Jrr = domain_geometry_.dFx_dr(r, theta); + double Jtr = domain_geometry_.dFy_dr(r, theta); + double Jrt = domain_geometry_.dFx_dt(r, theta); + double Jtt = domain_geometry_.dFy_dt(r, theta); + /* Compute the determinant of the Jacobian matrix */ + double detDF = Jrr * Jtt - Jrt * Jtr; + rhs_f[grid.index(i_r, i_theta)] *= 0.25 * (h1 + h2) * (k1 + k2) * fabs(detDF); + } + else if (i_r == 0 && DirBC_Interior_) { + rhs_f[grid.index(i_r, i_theta)] *= 1.0; + } + else if (i_r == grid.nr() - 1) { + rhs_f[grid.index(i_r, i_theta)] *= 1.0; + } + } + } + +// -------------------------------------------- // +// Discretize rhs values (radial index section) // +// -------------------------------------------- // +#pragma omp for nowait + for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { + double theta = grid.theta(i_theta); + + for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { + double r = grid.radius(i_r); + if ((0 < i_r && i_r < grid.nr() - 1) || (i_r == 0 && !DirBC_Interior_)) { + double h1 = (i_r == 0) ? 2.0 * grid.radius(0) : grid.radialSpacing(i_r - 1); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + /* Calculate the elements of the Jacobian matrix for the transformation mapping */ + /* The Jacobian matrix is: */ + /* [Jrr, Jrt] */ + /* [Jtr, Jtt] */ + double Jrr = domain_geometry_.dFx_dr(r, theta); + double Jtr = domain_geometry_.dFy_dr(r, theta); + double Jrt = domain_geometry_.dFx_dt(r, theta); + double Jtt = domain_geometry_.dFy_dt(r, theta); + /* Compute the determinant of the Jacobian matrix */ + double detDF = Jrr * Jtt - Jrt * Jtr; + rhs_f[grid.index(i_r, i_theta)] *= 0.25 * (h1 + h2) * (k1 + k2) * fabs(detDF); + } + else if (i_r == 0 && DirBC_Interior_) { + rhs_f[grid.index(i_r, i_theta)] *= 1.0; + } + else if (i_r == grid.nr() - 1) { + rhs_f[grid.index(i_r, i_theta)] *= 1.0; + } + } + } + } + } +} diff --git a/include/GMGPolar/gmgpolar.h b/include/GMGPolar/gmgpolar.h index cacbd137..81ec2e38 100644 --- a/include/GMGPolar/gmgpolar.h +++ b/include/GMGPolar/gmgpolar.h @@ -23,9 +23,10 @@ class LevelCache; #include "../PolarGrid/polargrid.h" #include "../common/global_definitions.h" #include "test_cases.h" +#include "igmgpolar.h" -template -class GMGPolar +template +class GMGPolar : public IGMGPolar { public: /* ---------------------------------------------------------------------- */ @@ -39,314 +40,35 @@ class GMGPolar // - domain_geometry: Mapping from the reference domain to the physical domain \Omega. // - density_profile_coefficients: Coefficients \alpha and \beta defining the PDE. GMGPolar(const PolarGrid& grid, const DomainGeometry& domain_geometry, - const DensityProfileCoefficients& density_profile_coefficients); - - /* ---------------------------------------------------------------------- */ - /* General output & visualization */ - /* ---------------------------------------------------------------------- */ - // Verbose level (0 = silent, >0 = increasingly detailed logs). - int verbose() const; - void verbose(int verbose); - - // Enable/disable ParaView output for visualization. - bool paraview() const; - void paraview(bool paraview); - - /* ---------------------------------------------------------------------- */ - /* Parallelization & threading */ - /* ---------------------------------------------------------------------- */ - // Maximum number of OpenMP threads to use. - int maxOpenMPThreads() const; - void maxOpenMPThreads(int max_omp_threads); - - // Thread reduction factor on coarser grids (e.g., 0.5 halves threads each level). - double threadReductionFactor() const; - void threadReductionFactor(double thread_reduction_factor); - - /* ---------------------------------------------------------------------- */ - /* Numerical method options */ - /* ---------------------------------------------------------------------- */ - // Treatment of the interior boundary at the origin: - // true -> use Dirichlet boundary condition - // false -> use discretization across the origin - bool DirBC_Interior() const; - void DirBC_Interior(bool DirBC_Interior); - - // Strategy for distributing the stencil (Take, Give). - StencilDistributionMethod stencilDistributionMethod() const; - void stencilDistributionMethod(StencilDistributionMethod stencil_distribution_method); - - // Cache density profile coefficients (alpha, beta). - bool cacheDensityProfileCoefficients() const; - void cacheDensityProfileCoefficients(bool cache_density_profile_coefficients); - - // Cache domain geometry data (arr, att, art, detDF). - bool cacheDomainGeometry() const; - void cacheDomainGeometry(bool cache_domain_geometry); - - /* ---------------------------------------------------------------------- */ - /* Multigrid controls */ - /* ---------------------------------------------------------------------- */ - // Implicit extrapolation to increase the order of convergence - ExtrapolationType extrapolation() const; - void extrapolation(ExtrapolationType extrapolation); - - // Maximum multigrid levels (-1 = use deepest possible). - int maxLevels() const; - void maxLevels(int max_levels); - - // Multigrid cycle type (V-cycle, W-cycle, F-cycle). - MultigridCycleType multigridCycle() const; - void multigridCycle(MultigridCycleType multigrid_cycle); - - // Pre-/post-smoothing steps per level. - int preSmoothingSteps() const; - void preSmoothingSteps(int pre_smoothing_steps); - int postSmoothingSteps() const; - void postSmoothingSteps(int post_smoothing_steps); - - // Full Multigrid (FMG) control. - bool FMG() const; - void FMG(bool FMG); - int FMG_iterations() const; - void FMG_iterations(int FMG_iterations); - MultigridCycleType FMG_cycle() const; - void FMG_cycle(MultigridCycleType FMG_cycle); - - /* ---------------------------------------------------------------------- */ - /* Iterative solver termination */ - /* ---------------------------------------------------------------------- */ - // Maximum number of iterations for the solver. - int maxIterations() const; - void maxIterations(int max_iterations); - - // Type of residual norm used to check convergence. - ResidualNormType residualNormType() const; - void residualNormType(ResidualNormType residual_norm_type); - - // Absolute residual tolerance for convergence. - std::optional absoluteTolerance() const; - void absoluteTolerance(std::optional tol); - - // Relative residual tolerance (relative to initial residual). - std::optional relativeTolerance() const; - void relativeTolerance(std::optional tol); + const DensityProfileCoefficients& density_profile_coefficients) + : IGMGPolar(grid, density_profile_coefficients) + , domain_geometry_(domain_geometry) + { + } /* ---------------------------------------------------------------------- */ /* Setup & Solve */ /* ---------------------------------------------------------------------- */ // Finalize solver setup (allocate data, build operators, etc.). - void setup(); - - // If an exact solution is provided, the solver will compute the exact error at each iteration. - void setSolution(const ExactSolution* exact_solution); - - // Solve system with given boundary conditions and source term. - // Multiple solves with different inputs are supported. - void solve(const BoundaryConditions& boundary_conditions, const SourceTerm& source_term); - - /* ---------------------------------------------------------------------- */ - /* Solution & Grid Access */ - /* ---------------------------------------------------------------------- */ - // Return a reference to the computed solution vector. - Vector solution(); - ConstVector solution() const; - - // Return the underlying cartesian mesh used for discretization. - const PolarGrid& grid() const; - - /* ---------------------------------------------------------------------- */ - /* Diagnostics & statistics */ - /* ---------------------------------------------------------------------- */ - // Print timing breakdown for setup, smoothing, coarse solve, etc. - void printTimings() const; - - // Number of iterations taken by last solve. - int numberOfIterations() const; - - // Mean residual reduction factor per iteration. - double meanResidualReductionFactor() const; - - // Error norms (only available if exact solution was set). - std::optional exactErrorWeightedEuclidean() const; - std::optional exactErrorInfinity() const; - - /* ---------------------------------------------------------------------- */ - /* Timing Statistics */ - /* ---------------------------------------------------------------------- */ - double timeSetupTotal() const; - double timeSetupCreateLevels() const; - double timeSetupRHS() const; - double timeSetupSmoother() const; - double timeSetupDirectSolver() const; - - double timeSolveTotal() const; - double timeSolveInitialApproximation() const; - double timeSolveMultigridIterations() const; - double timeCheckConvergence() const; - double timeCheckExactError() const; - - double timeAvgMGCTotal() const; - double timeAvgMGCPreSmoothing() const; - double timeAvgMGCPostSmoothing() const; - double timeAvgMGCResidual() const; - double timeAvgMGCDirectSolver() const; + void setup() final; private: /* ------------------------------------ */ /* Grid Configuration & Input Functions */ /* ------------------------------------ */ - const PolarGrid& grid_; const DomainGeometry& domain_geometry_; - const DensityProfileCoefficients& density_profile_coefficients_; - const ExactSolution* exact_solution_; // Optional exact solution for validation - - /* ------------------ */ - /* Control Parameters */ - /* ------------------ */ - // General solver output and visualization settings - int verbose_; - bool paraview_; - // Parallelization and threading settings - int max_omp_threads_; - double thread_reduction_factor_; - // Numerical method setup - bool DirBC_Interior_; - StencilDistributionMethod stencil_distribution_method_; - bool cache_density_profile_coefficients_; - bool cache_domain_geometry_; - // Multigrid settings - ExtrapolationType extrapolation_; - int max_levels_; - int pre_smoothing_steps_; - int post_smoothing_steps_; - MultigridCycleType multigrid_cycle_; - // FMG settings - bool FMG_; - int FMG_iterations_; - MultigridCycleType FMG_cycle_; - // Convergence settings - int max_iterations_; - ResidualNormType residual_norm_type_; - std::optional absolute_tolerance_; - std::optional relative_tolerance_; - - /* ---------------- */ - /* Multigrid levels */ - int number_of_levels_; - std::vector levels_; - std::vector threads_per_level_; - - /* ---------------------- */ - /* Interpolation operator */ - std::unique_ptr interpolation_; - - /* ------------------------------------------------------------------------- */ - /* Chooses if full grid smoothing is active on level 0 for extrapolation > 0 */ - bool full_grid_smoothing_ = false; - - /* -------------------- */ - /* Convergence criteria */ - int number_of_iterations_; - std::vector residual_norms_; - double mean_residual_reduction_factor_; - bool converged(double current_residual_norm, double first_residual_norm); - - /* ---------------------------------------------------- */ - /* Compute exact error if an exact solution is provided */ - // The results are stored as a pair: (weighted L2 error, infinity error). - std::vector> exact_errors_; - std::pair computeExactError(Level& level, ConstVector solution, Vector error, - const ExactSolution& exact_solution); - - /* ------------------------------------------------------------------------- */ - /* Compute the extrapolated residual: res_ex = 4/3 res_fine - 1/3 res_coarse */ - void extrapolatedResidual(const int current_level, Vector residual, - ConstVector residual_next_level); /* --------------- */ /* Setup Functions */ - int chooseNumberOfLevels(const PolarGrid& finest_grid); - void build_rhs_f(const Level& level, Vector rhs_f, const BoundaryConditions& boundary_conditions, - const SourceTerm& source_term); - void discretize_rhs_f(const Level& level, Vector rhs_f); - - /* --------------- */ - /* Solve Functions */ - void initializeSolution(); - double residualNorm(const ResidualNormType& norm_type, const Level& level, ConstVector residual) const; - void evaluateExactError(Level& level, const ExactSolution& exact_solution); - void updateResidualNorms(Level& level, int iteration, double& initial_residual_norm, double& current_residual_norm, - double& current_relative_residual_norm); - - /* ----------------- */ - /* Print information */ - void printSettings() const; - void printIterationHeader(const ExactSolution* exact_solution); - void printIterationInfo(int iteration, double current_residual_norm, double current_relative_residual_norm, - const ExactSolution* exact_solution); - - /* ------------------- */ - /* Multigrid Functions */ - void multigrid_V_Cycle(const int level_depth, Vector solution, Vector rhs, Vector residual); - void multigrid_W_Cycle(const int level_depth, Vector solution, Vector rhs, Vector residual); - void multigrid_F_Cycle(const int level_depth, Vector solution, Vector rhs, Vector residual); - void implicitlyExtrapolatedMultigrid_V_Cycle(const int level_depth, Vector solution, Vector rhs, - Vector residual); - void implicitlyExtrapolatedMultigrid_W_Cycle(const int level_depth, Vector solution, Vector rhs, - Vector residual); - void implicitlyExtrapolatedMultigrid_F_Cycle(const int level_depth, Vector solution, Vector rhs, - Vector residual); - - /* ----------------------- */ - /* Interpolation functions */ - void prolongation(const int current_level, Vector result, ConstVector x) const; - void restriction(const int current_level, Vector result, ConstVector x) const; - void injection(const int current_level, Vector result, ConstVector x) const; - void extrapolatedProlongation(const int current_level, Vector result, ConstVector x) const; - void extrapolatedRestriction(const int current_level, Vector result, ConstVector x) const; - void FMGInterpolation(const int current_level, Vector result, ConstVector x) const; + void discretize_rhs_f(const Level& level, Vector rhs_f) final; /* ------------- */ /* Visualization */ - void writeToVTK(const std::filesystem::path& file_path, const PolarGrid& grid); - void writeToVTK(const std::filesystem::path& file_path, const Level& level, ConstVector grid_function); - - /* ------------------------------ */ - /* Timing statistics for GMGPolar */ - void resetAllTimings(); - - void resetSetupPhaseTimings(); - double t_setup_total_; - double t_setup_createLevels_; - double t_setup_rhs_; - double t_setup_smoother_; - double t_setup_directSolver_; - - void resetSolvePhaseTimings(); - double t_solve_total_; - double t_solve_initial_approximation_; - double t_solve_multigrid_iterations_; - double t_check_convergence_; - double t_check_exact_error_; - - void resetAvgMultigridCycleTimings(); - double t_avg_MGC_total_; - double t_avg_MGC_preSmoothing_; - double t_avg_MGC_postSmoothing_; - double t_avg_MGC_residual_; - double t_avg_MGC_directSolver_; + void writeToVTK(const std::filesystem::path& file_path, const PolarGrid& grid) final; + void writeToVTK(const std::filesystem::path& file_path, const Level& level, + ConstVector grid_function) final; }; #include "build_rhs_f.h" -#include "gmgpolar_methods.h" -#include "level_interpolation.h" #include "setup.h" -#include "solver.h" #include "writeToVTK.h" -#include "MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.h" -#include "MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.h" -#include "MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.h" -#include "MultigridMethods/multigrid_F_Cycle.h" -#include "MultigridMethods/multigrid_V_Cycle.h" -#include "MultigridMethods/multigrid_W_Cycle.h" diff --git a/include/GMGPolar/igmgpolar.h b/include/GMGPolar/igmgpolar.h new file mode 100644 index 00000000..60d8585d --- /dev/null +++ b/include/GMGPolar/igmgpolar.h @@ -0,0 +1,336 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +class Level; +class LevelCache; + +#include "../InputFunctions/boundaryConditions.h" +#include "../InputFunctions/densityProfileCoefficients.h" +#include "../InputFunctions/domainGeometry.h" +#include "../InputFunctions/exactSolution.h" +#include "../InputFunctions/sourceTerm.h" +#include "../Interpolation/interpolation.h" +#include "../Level/level.h" +#include "../LinearAlgebra/vector.h" +#include "../LinearAlgebra/vector_operations.h" +#include "../PolarGrid/polargrid.h" +#include "../common/global_definitions.h" +#include "test_cases.h" + +class IGMGPolar +{ +public: + /* ---------------------------------------------------------------------- */ + /* Constructor & Initialization */ + /* ---------------------------------------------------------------------- */ + // Construct a polar PDE multigrid solver for the Poisson-like equation: + // - \nabla \cdot (\alpha \nabla u) + \beta u = rhs_f in \Omega, + // with Dirichlet boundary condition u = u_D on \partial \Omega. + // Parameters: + // - grid: Cartesian mesh discretizing the computational domain. + // - density_profile_coefficients: Coefficients \alpha and \beta defining the PDE. + IGMGPolar(const PolarGrid& grid, + const DensityProfileCoefficients& density_profile_coefficients); + + /* ---------------------------------------------------------------------- */ + /* General output & visualization */ + /* ---------------------------------------------------------------------- */ + // Verbose level (0 = silent, >0 = increasingly detailed logs). + int verbose() const; + void verbose(int verbose); + + // Enable/disable ParaView output for visualization. + bool paraview() const; + void paraview(bool paraview); + + /* ---------------------------------------------------------------------- */ + /* Parallelization & threading */ + /* ---------------------------------------------------------------------- */ + // Maximum number of OpenMP threads to use. + int maxOpenMPThreads() const; + void maxOpenMPThreads(int max_omp_threads); + + // Thread reduction factor on coarser grids (e.g., 0.5 halves threads each level). + double threadReductionFactor() const; + void threadReductionFactor(double thread_reduction_factor); + + /* ---------------------------------------------------------------------- */ + /* Numerical method options */ + /* ---------------------------------------------------------------------- */ + // Treatment of the interior boundary at the origin: + // true -> use Dirichlet boundary condition + // false -> use discretization across the origin + bool DirBC_Interior() const; + void DirBC_Interior(bool DirBC_Interior); + + // Strategy for distributing the stencil (Take, Give). + StencilDistributionMethod stencilDistributionMethod() const; + void stencilDistributionMethod(StencilDistributionMethod stencil_distribution_method); + + // Cache density profile coefficients (alpha, beta). + bool cacheDensityProfileCoefficients() const; + void cacheDensityProfileCoefficients(bool cache_density_profile_coefficients); + + // Cache domain geometry data (arr, att, art, detDF). + bool cacheDomainGeometry() const; + void cacheDomainGeometry(bool cache_domain_geometry); + + /* ---------------------------------------------------------------------- */ + /* Multigrid controls */ + /* ---------------------------------------------------------------------- */ + // Implicit extrapolation to increase the order of convergence + ExtrapolationType extrapolation() const; + void extrapolation(ExtrapolationType extrapolation); + + // Maximum multigrid levels (-1 = use deepest possible). + int maxLevels() const; + void maxLevels(int max_levels); + + // Multigrid cycle type (V-cycle, W-cycle, F-cycle). + MultigridCycleType multigridCycle() const; + void multigridCycle(MultigridCycleType multigrid_cycle); + + // Pre-/post-smoothing steps per level. + int preSmoothingSteps() const; + void preSmoothingSteps(int pre_smoothing_steps); + int postSmoothingSteps() const; + void postSmoothingSteps(int post_smoothing_steps); + + // Full Multigrid (FMG) control. + bool FMG() const; + void FMG(bool FMG); + int FMG_iterations() const; + void FMG_iterations(int FMG_iterations); + MultigridCycleType FMG_cycle() const; + void FMG_cycle(MultigridCycleType FMG_cycle); + + /* ---------------------------------------------------------------------- */ + /* Iterative solver termination */ + /* ---------------------------------------------------------------------- */ + // Maximum number of iterations for the solver. + int maxIterations() const; + void maxIterations(int max_iterations); + + // Type of residual norm used to check convergence. + ResidualNormType residualNormType() const; + void residualNormType(ResidualNormType residual_norm_type); + + // Absolute residual tolerance for convergence. + std::optional absoluteTolerance() const; + void absoluteTolerance(std::optional tol); + + // Relative residual tolerance (relative to initial residual). + std::optional relativeTolerance() const; + void relativeTolerance(std::optional tol); + + /* ---------------------------------------------------------------------- */ + /* Setup & Solve */ + /* ---------------------------------------------------------------------- */ + // Finalize solver setup (allocate data, build operators, etc.). + virtual void setup() = 0; + + // If an exact solution is provided, the solver will compute the exact error at each iteration. + void setSolution(const ExactSolution* exact_solution); + + // Solve system with given boundary conditions and source term. + // Multiple solves with different inputs are supported. + void solve(const BoundaryConditions& boundary_conditions, const SourceTerm& source_term); + + /* ---------------------------------------------------------------------- */ + /* Solution & Grid Access */ + /* ---------------------------------------------------------------------- */ + // Return a reference to the computed solution vector. + Vector solution(); + ConstVector solution() const; + + // Return the underlying cartesian mesh used for discretization. + const PolarGrid& grid() const; + + /* ---------------------------------------------------------------------- */ + /* Diagnostics & statistics */ + /* ---------------------------------------------------------------------- */ + // Print timing breakdown for setup, smoothing, coarse solve, etc. + void printTimings() const; + + // Number of iterations taken by last solve. + int numberOfIterations() const; + + // Mean residual reduction factor per iteration. + double meanResidualReductionFactor() const; + + // Error norms (only available if exact solution was set). + std::optional exactErrorWeightedEuclidean() const; + std::optional exactErrorInfinity() const; + + /* ---------------------------------------------------------------------- */ + /* Timing Statistics */ + /* ---------------------------------------------------------------------- */ + double timeSetupTotal() const; + double timeSetupCreateLevels() const; + double timeSetupRHS() const; + double timeSetupSmoother() const; + double timeSetupDirectSolver() const; + + double timeSolveTotal() const; + double timeSolveInitialApproximation() const; + double timeSolveMultigridIterations() const; + double timeCheckConvergence() const; + double timeCheckExactError() const; + + double timeAvgMGCTotal() const; + double timeAvgMGCPreSmoothing() const; + double timeAvgMGCPostSmoothing() const; + double timeAvgMGCResidual() const; + double timeAvgMGCDirectSolver() const; + +protected: + /* ------------------------------------ */ + /* Grid Configuration & Input Functions */ + /* ------------------------------------ */ + const PolarGrid& grid_; + const DensityProfileCoefficients& density_profile_coefficients_; + const ExactSolution* exact_solution_; // Optional exact solution for validation + + /* ------------------ */ + /* Control Parameters */ + /* ------------------ */ + // General solver output and visualization settings + int verbose_; + bool paraview_; + // Parallelization and threading settings + int max_omp_threads_; + double thread_reduction_factor_; + // Numerical method setup + bool DirBC_Interior_; + StencilDistributionMethod stencil_distribution_method_; + bool cache_density_profile_coefficients_; + bool cache_domain_geometry_; + // Multigrid settings + ExtrapolationType extrapolation_; + int max_levels_; + int pre_smoothing_steps_; + int post_smoothing_steps_; + MultigridCycleType multigrid_cycle_; + // FMG settings + bool FMG_; + int FMG_iterations_; + MultigridCycleType FMG_cycle_; + // Convergence settings + int max_iterations_; + ResidualNormType residual_norm_type_; + std::optional absolute_tolerance_; + std::optional relative_tolerance_; + + /* ---------------- */ + /* Multigrid levels */ + int number_of_levels_; + std::vector levels_; + std::vector threads_per_level_; + + /* ---------------------- */ + /* Interpolation operator */ + std::unique_ptr interpolation_; + + /* ------------------------------------------------------------------------- */ + /* Chooses if full grid smoothing is active on level 0 for extrapolation > 0 */ + bool full_grid_smoothing_ = false; + + /* -------------------- */ + /* Convergence criteria */ + int number_of_iterations_; + std::vector residual_norms_; + double mean_residual_reduction_factor_; + bool converged(double current_residual_norm, double first_residual_norm); + + /* ---------------------------------------------------- */ + /* Compute exact error if an exact solution is provided */ + // The results are stored as a pair: (weighted L2 error, infinity error). + std::vector> exact_errors_; + std::pair computeExactError(Level& level, ConstVector solution, Vector error, + const ExactSolution& exact_solution); + + /* ------------------------------------------------------------------------- */ + /* Compute the extrapolated residual: res_ex = 4/3 res_fine - 1/3 res_coarse */ + void extrapolatedResidual(const int current_level, Vector residual, + ConstVector residual_next_level); + + /* --------------- */ + /* Setup Functions */ + int chooseNumberOfLevels(const PolarGrid& finest_grid); + void build_rhs_f(const Level& level, Vector rhs_f, const BoundaryConditions& boundary_conditions, + const SourceTerm& source_term); + virtual void discretize_rhs_f(const Level& level, Vector rhs_f) = 0; + + /* --------------- */ + /* Solve Functions */ + void initializeSolution(); + double residualNorm(const ResidualNormType& norm_type, const Level& level, ConstVector residual) const; + void evaluateExactError(Level& level, const ExactSolution& exact_solution); + void updateResidualNorms(Level& level, int iteration, double& initial_residual_norm, double& current_residual_norm, + double& current_relative_residual_norm); + + /* ----------------- */ + /* Print information */ + void printSettings() const; + void printIterationHeader(const ExactSolution* exact_solution); + void printIterationInfo(int iteration, double current_residual_norm, double current_relative_residual_norm, + const ExactSolution* exact_solution); + + /* ------------------- */ + /* Multigrid Functions */ + void multigrid_V_Cycle(const int level_depth, Vector solution, Vector rhs, Vector residual); + void multigrid_W_Cycle(const int level_depth, Vector solution, Vector rhs, Vector residual); + void multigrid_F_Cycle(const int level_depth, Vector solution, Vector rhs, Vector residual); + void implicitlyExtrapolatedMultigrid_V_Cycle(const int level_depth, Vector solution, Vector rhs, + Vector residual); + void implicitlyExtrapolatedMultigrid_W_Cycle(const int level_depth, Vector solution, Vector rhs, + Vector residual); + void implicitlyExtrapolatedMultigrid_F_Cycle(const int level_depth, Vector solution, Vector rhs, + Vector residual); + + /* ----------------------- */ + /* Interpolation functions */ + void prolongation(const int current_level, Vector result, ConstVector x) const; + void restriction(const int current_level, Vector result, ConstVector x) const; + void injection(const int current_level, Vector result, ConstVector x) const; + void extrapolatedProlongation(const int current_level, Vector result, ConstVector x) const; + void extrapolatedRestriction(const int current_level, Vector result, ConstVector x) const; + void FMGInterpolation(const int current_level, Vector result, ConstVector x) const; + + /* ------------- */ + /* Visualization */ + virtual void writeToVTK(const std::filesystem::path& file_path, const PolarGrid& grid) = 0; + virtual void writeToVTK(const std::filesystem::path& file_path, const Level& level, ConstVector grid_function) = 0; + + /* ------------------------------ */ + /* Timing statistics for GMGPolar */ + void resetAllTimings(); + + void resetSetupPhaseTimings(); + double t_setup_total_; + double t_setup_createLevels_; + double t_setup_rhs_; + double t_setup_smoother_; + double t_setup_directSolver_; + + void resetSolvePhaseTimings(); + double t_solve_total_; + double t_solve_initial_approximation_; + double t_solve_multigrid_iterations_; + double t_check_convergence_; + double t_check_exact_error_; + + void resetAvgMultigridCycleTimings(); + double t_avg_MGC_total_; + double t_avg_MGC_preSmoothing_; + double t_avg_MGC_postSmoothing_; + double t_avg_MGC_residual_; + double t_avg_MGC_directSolver_; +}; diff --git a/include/GMGPolar/setup.h b/include/GMGPolar/setup.h new file mode 100644 index 00000000..e17a11cf --- /dev/null +++ b/include/GMGPolar/setup.h @@ -0,0 +1,173 @@ + +template +void GMGPolar::setup() +{ + LIKWID_START("Setup"); + auto start_setup = std::chrono::high_resolution_clock::now(); + + resetSetupPhaseTimings(); + + auto start_setup_createLevels = std::chrono::high_resolution_clock::now(); + + if (stencil_distribution_method_ == StencilDistributionMethod::CPU_TAKE) { + if (!cache_density_profile_coefficients_ || !cache_domain_geometry_) { + throw std::runtime_error("Error: Caching must be enabled for both density profile coefficients and domain " + "geometry in 'Take' implementation strategy."); + } + } + + // -------------------------------- // + // Create the finest mesh (level 0) // + // -------------------------------- // + auto finest_grid = std::make_unique(grid_); + if (paraview_) + writeToVTK("output_finest_grid", *finest_grid); + + if (extrapolation_ != ExtrapolationType::NONE) { + // Check if grid comes from a single uniform refinement + bool is_uniform_refinement = true; + + for (int i_r = 1; i_r < finest_grid->nr() - 1; i_r += 2) { + double mid = 0.5 * (finest_grid->radius(i_r - 1) + finest_grid->radius(i_r + 1)); + if (std::abs(mid - finest_grid->radius(i_r)) > 1e-12) { + is_uniform_refinement = false; + break; + } + } + for (int i_theta = 1; i_theta < finest_grid->ntheta(); i_theta += 2) { + double mid = 0.5 * (finest_grid->theta(i_theta - 1) + finest_grid->theta(i_theta + 1)); + if (std::abs(mid - finest_grid->theta(i_theta)) > 1e-12) { + is_uniform_refinement = false; + break; + } + } + + if (!is_uniform_refinement) { + throw std::runtime_error( + "Extrapolation Error: Finest PolarGrid does not originate from a single uniform refinement."); + } + } + + // ---------------------------------------------------------- // + // Building PolarGrid and LevelCache for all multigrid levels // + // ---------------------------------------------------------- // + number_of_levels_ = chooseNumberOfLevels(*finest_grid); /* Implementation below */ + 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_); + + for (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(levels_[level_depth - 1], *current_grid); + levels_.emplace_back(level_depth, std::move(current_grid), std::move(current_levelCache), extrapolation_, FMG_); + } + + auto end_setup_createLevels = std::chrono::high_resolution_clock::now(); + t_setup_createLevels_ = std::chrono::duration(end_setup_createLevels - start_setup_createLevels).count(); + + if (paraview_) + writeToVTK("output_coarsest_grid", levels_.back().grid()); + + // ----------------------------------------------------------- // + // Initializing the optimal number of threads for OpenMP tasks // + // ----------------------------------------------------------- // + threads_per_level_.resize(number_of_levels_, max_omp_threads_); + for (int level_depth = 0; level_depth < number_of_levels_; level_depth++) { + threads_per_level_[level_depth] = std::max( + 1, + std::min(max_omp_threads_, + static_cast(std::floor(max_omp_threads_ * std::pow(thread_reduction_factor_, level_depth))))); + } + + interpolation_ = std::make_unique(threads_per_level_, DirBC_Interior_); + + if (verbose_ > 0) + printSettings(); + + // ------------------------------------------------------- + // Initializing various operators based on the level index + 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_, threads_per_level_[level_depth], + stencil_distribution_method_); + break; + case ExtrapolationType::IMPLICIT_EXTRAPOLATION: + full_grid_smoothing_ = false; + levels_[level_depth].initializeExtrapolatedSmoothing(domain_geometry_, density_profile_coefficients_, + DirBC_Interior_, threads_per_level_[level_depth], + 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_, threads_per_level_[level_depth], + stencil_distribution_method_); + break; + case ExtrapolationType::COMBINED: + full_grid_smoothing_ = true; + levels_[level_depth].initializeSmoothing(domain_geometry_, density_profile_coefficients_, + DirBC_Interior_, threads_per_level_[level_depth], + stencil_distribution_method_); + levels_[level_depth].initializeExtrapolatedSmoothing(domain_geometry_, density_profile_coefficients_, + DirBC_Interior_, threads_per_level_[level_depth], + stencil_distribution_method_); + break; + default: + full_grid_smoothing_ = false; + levels_[level_depth].initializeSmoothing(domain_geometry_, density_profile_coefficients_, + DirBC_Interior_, threads_per_level_[level_depth], + stencil_distribution_method_); + levels_[level_depth].initializeExtrapolatedSmoothing(domain_geometry_, density_profile_coefficients_, + DirBC_Interior_, threads_per_level_[level_depth], + 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_, + threads_per_level_[level_depth], 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_, threads_per_level_[level_depth], + 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_, + threads_per_level_[level_depth], stencil_distribution_method_); + } + // ------------------- // + // Intermediate levels // + // ------------------- // + else { + auto start_setup_smoother = std::chrono::high_resolution_clock::now(); + levels_[level_depth].initializeSmoothing(domain_geometry_, density_profile_coefficients_, DirBC_Interior_, + threads_per_level_[level_depth], 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_, + threads_per_level_[level_depth], stencil_distribution_method_); + } + } + + auto end_setup = std::chrono::high_resolution_clock::now(); + t_setup_total_ = std::chrono::duration(end_setup - start_setup).count(); + LIKWID_STOP("Setup"); +} diff --git a/src/GMGPolar/writeToVTK.cpp b/include/GMGPolar/writeToVTK.h similarity index 99% rename from src/GMGPolar/writeToVTK.cpp rename to include/GMGPolar/writeToVTK.h index 41280a0e..88a7fc9f 100644 --- a/src/GMGPolar/writeToVTK.cpp +++ b/include/GMGPolar/writeToVTK.h @@ -1,3 +1,4 @@ +#include "../../include/GMGPolar/gmgpolar.h" template void GMGPolar::writeToVTK(const std::filesystem::path& file_path, const PolarGrid& grid) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index dab54db3..1906732b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -12,6 +12,25 @@ set(POLAR_GRID_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/PolarGrid/polargrid.cpp ) +# file(GLOB_RECURSE GMG_POLAR_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/*.cpp) +set(GMG_POLAR_SOURCES + ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/build_rhs_f.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/gmgpolar.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/level_interpolation.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/setup.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/solver.cpp +) + +# file(GLOB_RECURSE MULTIGRID_METHODS_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/MultigridMethods/*.cpp) +set(MULTIGRID_METHODS_SOURCES + ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp +) + #file(GLOB_RECURSE LEVEL_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/Level/*.cpp) set(LEVEL_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/Level/level.cpp @@ -142,6 +161,8 @@ set(CONFIG_PARSER_SOURCES # Create the main library add_library(GMGPolarLib STATIC ${POLAR_GRID_SOURCES} + ${GMG_POLAR_SOURCES} + ${MULTIGRID_METHODS_SOURCES} ${LEVEL_SOURCES} ${STENCIL_SOURCES} ${INTERPOLATION_SOURCES} diff --git a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp index 72ad172d..7a91d3c8 100644 --- a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp @@ -1,7 +1,7 @@ +#include "../../../include/GMGPolar/gmgpolar.h" -template -void GMGPolar::implicitlyExtrapolatedMultigrid_F_Cycle(const int level_depth, Vector solution, - Vector rhs, Vector residual) +void IGMGPolar::implicitlyExtrapolatedMultigrid_F_Cycle(const int level_depth, Vector solution, + Vector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); @@ -123,4 +123,4 @@ void GMGPolar::implicitlyExtrapolatedMultigrid_F_Cycle(const int auto end_MGC = std::chrono::high_resolution_clock::now(); t_avg_MGC_total_ += std::chrono::duration(end_MGC - start_MGC).count(); } -} +} \ No newline at end of file diff --git a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp index dadfe565..00e5443b 100644 --- a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp @@ -1,7 +1,7 @@ +#include "../../../include/GMGPolar/gmgpolar.h" -template -void GMGPolar::implicitlyExtrapolatedMultigrid_V_Cycle(const int level_depth, Vector solution, - Vector rhs, Vector residual) +void IGMGPolar::implicitlyExtrapolatedMultigrid_V_Cycle(const int level_depth, Vector solution, + Vector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); @@ -122,4 +122,4 @@ void GMGPolar::implicitlyExtrapolatedMultigrid_V_Cycle(const int auto end_MGC = std::chrono::high_resolution_clock::now(); t_avg_MGC_total_ += std::chrono::duration(end_MGC - start_MGC).count(); } -} +} \ No newline at end of file diff --git a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp index 8fe51fad..5863f356 100644 --- a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp @@ -1,7 +1,7 @@ +#include "../../../include/GMGPolar/gmgpolar.h" -template -void GMGPolar::implicitlyExtrapolatedMultigrid_W_Cycle(const int level_depth, Vector solution, - Vector rhs, Vector residual) +void IGMGPolar::implicitlyExtrapolatedMultigrid_W_Cycle(const int level_depth, Vector solution, + Vector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); @@ -123,4 +123,4 @@ void GMGPolar::implicitlyExtrapolatedMultigrid_W_Cycle(const int auto end_MGC = std::chrono::high_resolution_clock::now(); t_avg_MGC_total_ += std::chrono::duration(end_MGC - start_MGC).count(); } -} +} \ No newline at end of file diff --git a/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp b/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp index a0bc1396..60a65d22 100644 --- a/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp @@ -1,7 +1,7 @@ +#include "../../../include/GMGPolar/gmgpolar.h" -template -void GMGPolar::multigrid_F_Cycle(const int level_depth, Vector solution, Vector rhs, - Vector residual) +void IGMGPolar::multigrid_F_Cycle(const int level_depth, Vector solution, Vector rhs, + Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); @@ -91,4 +91,4 @@ void GMGPolar::multigrid_F_Cycle(const int level_depth, Vector(end_MGC - start_MGC).count(); } -} +} \ No newline at end of file diff --git a/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp b/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp index a47a197a..1e6ce437 100644 --- a/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp @@ -1,7 +1,7 @@ +#include "../../../include/GMGPolar/gmgpolar.h" -template -void GMGPolar::multigrid_V_Cycle(const int level_depth, Vector solution, Vector rhs, - Vector residual) +void IGMGPolar::multigrid_V_Cycle(const int level_depth, Vector solution, Vector rhs, + Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); @@ -90,4 +90,4 @@ void GMGPolar::multigrid_V_Cycle(const int level_depth, Vector(end_MGC - start_MGC).count(); } -} +} \ No newline at end of file diff --git a/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp b/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp index 52b95069..a6cfc599 100644 --- a/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp @@ -1,7 +1,7 @@ +#include "../../../include/GMGPolar/gmgpolar.h" -template -void GMGPolar::multigrid_W_Cycle(const int level_depth, Vector solution, Vector rhs, - Vector residual) +void IGMGPolar::multigrid_W_Cycle(const int level_depth, Vector solution, Vector rhs, + Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); @@ -91,4 +91,4 @@ void GMGPolar::multigrid_W_Cycle(const int level_depth, Vector(end_MGC - start_MGC).count(); } -} +} \ No newline at end of file diff --git a/src/GMGPolar/build_rhs_f.cpp b/src/GMGPolar/build_rhs_f.cpp index a3c1287e..ba3c5b85 100644 --- a/src/GMGPolar/build_rhs_f.cpp +++ b/src/GMGPolar/build_rhs_f.cpp @@ -1,6 +1,6 @@ +#include "../../include/GMGPolar/gmgpolar.h" -template -void GMGPolar::build_rhs_f(const Level& level, Vector rhs_f, const BoundaryConditions& boundary_conditions, +void IGMGPolar::build_rhs_f(const Level& level, Vector rhs_f, const BoundaryConditions& boundary_conditions, const SourceTerm& source_term) { const PolarGrid& grid = level.grid(); @@ -51,143 +51,3 @@ void GMGPolar::build_rhs_f(const Level& level, Vector rh } } } - -template -void GMGPolar::discretize_rhs_f(const Level& level, Vector rhs_f) -{ - const PolarGrid& grid = level.grid(); - assert(rhs_f.size() == static_cast(grid.numberOfNodes())); - - if (level.levelCache().cacheDomainGeometry()) { - /* DomainGeometry is cached */ - const auto& detDF_cache = level.levelCache().detDF(); -#pragma omp parallel - { -// ---------------------------------------------- // -// Discretize rhs values (circular index section) // -// ---------------------------------------------- // -#pragma omp for nowait - for (int i_r = 0; i_r < grid.numberSmootherCircles(); i_r++) { - double r = grid.radius(i_r); - for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { - double theta = grid.theta(i_theta); - if ((0 < i_r && i_r < grid.nr() - 1) || (i_r == 0 && !DirBC_Interior_)) { - double h1 = (i_r == 0) ? 2.0 * grid.radius(0) : grid.radialSpacing(i_r - 1); - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); - const double detDF = detDF_cache[grid.index(i_r, i_theta)]; - rhs_f[grid.index(i_r, i_theta)] *= 0.25 * (h1 + h2) * (k1 + k2) * fabs(detDF); - } - else if (i_r == 0 && DirBC_Interior_) { - rhs_f[grid.index(i_r, i_theta)] *= 1.0; - } - else if (i_r == grid.nr() - 1) { - rhs_f[grid.index(i_r, i_theta)] *= 1.0; - } - } - } - -// -------------------------------------------- // -// Discretize rhs values (radial index section) // -// -------------------------------------------- // -#pragma omp for nowait - for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { - double theta = grid.theta(i_theta); - for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { - double r = grid.radius(i_r); - if ((0 < i_r && i_r < grid.nr() - 1) || (i_r == 0 && !DirBC_Interior_)) { - double h1 = (i_r == 0) ? 2.0 * grid.radius(0) : grid.radialSpacing(i_r - 1); - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); - const double detDF = detDF_cache[grid.index(i_r, i_theta)]; - rhs_f[grid.index(i_r, i_theta)] *= 0.25 * (h1 + h2) * (k1 + k2) * fabs(detDF); - } - else if (i_r == 0 && DirBC_Interior_) { - rhs_f[grid.index(i_r, i_theta)] *= 1.0; - } - else if (i_r == grid.nr() - 1) { - rhs_f[grid.index(i_r, i_theta)] *= 1.0; - } - } - } - } - } - else { - /* DomainGeometry is not cached */ - -#pragma omp parallel - { -// ---------------------------------------------- // -// Discretize rhs values (circular index section) // -// ---------------------------------------------- // -#pragma omp for nowait - for (int i_r = 0; i_r < grid.numberSmootherCircles(); i_r++) { - double r = grid.radius(i_r); - for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { - double theta = grid.theta(i_theta); - - if ((0 < i_r && i_r < grid.nr() - 1) || (i_r == 0 && !DirBC_Interior_)) { - double h1 = (i_r == 0) ? 2.0 * grid.radius(0) : grid.radialSpacing(i_r - 1); - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); - /* Calculate the elements of the Jacobian matrix for the transformation mapping */ - /* The Jacobian matrix is: */ - /* [Jrr, Jrt] */ - /* [Jtr, Jtt] */ - double Jrr = domain_geometry_.dFx_dr(r, theta); - double Jtr = domain_geometry_.dFy_dr(r, theta); - double Jrt = domain_geometry_.dFx_dt(r, theta); - double Jtt = domain_geometry_.dFy_dt(r, theta); - /* Compute the determinant of the Jacobian matrix */ - double detDF = Jrr * Jtt - Jrt * Jtr; - rhs_f[grid.index(i_r, i_theta)] *= 0.25 * (h1 + h2) * (k1 + k2) * fabs(detDF); - } - else if (i_r == 0 && DirBC_Interior_) { - rhs_f[grid.index(i_r, i_theta)] *= 1.0; - } - else if (i_r == grid.nr() - 1) { - rhs_f[grid.index(i_r, i_theta)] *= 1.0; - } - } - } - -// -------------------------------------------- // -// Discretize rhs values (radial index section) // -// -------------------------------------------- // -#pragma omp for nowait - for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { - double theta = grid.theta(i_theta); - - for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { - double r = grid.radius(i_r); - if ((0 < i_r && i_r < grid.nr() - 1) || (i_r == 0 && !DirBC_Interior_)) { - double h1 = (i_r == 0) ? 2.0 * grid.radius(0) : grid.radialSpacing(i_r - 1); - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); - /* Calculate the elements of the Jacobian matrix for the transformation mapping */ - /* The Jacobian matrix is: */ - /* [Jrr, Jrt] */ - /* [Jtr, Jtt] */ - double Jrr = domain_geometry_.dFx_dr(r, theta); - double Jtr = domain_geometry_.dFy_dr(r, theta); - double Jrt = domain_geometry_.dFx_dt(r, theta); - double Jtt = domain_geometry_.dFy_dt(r, theta); - /* Compute the determinant of the Jacobian matrix */ - double detDF = Jrr * Jtt - Jrt * Jtr; - rhs_f[grid.index(i_r, i_theta)] *= 0.25 * (h1 + h2) * (k1 + k2) * fabs(detDF); - } - else if (i_r == 0 && DirBC_Interior_) { - rhs_f[grid.index(i_r, i_theta)] *= 1.0; - } - else if (i_r == grid.nr() - 1) { - rhs_f[grid.index(i_r, i_theta)] *= 1.0; - } - } - } - } - } -} diff --git a/src/GMGPolar/gmgpolar.cpp b/src/GMGPolar/gmgpolar.cpp index 8f8624a7..c74028ad 100644 --- a/src/GMGPolar/gmgpolar.cpp +++ b/src/GMGPolar/gmgpolar.cpp @@ -1,12 +1,11 @@ +#include "../../include/GMGPolar/gmgpolar.h" /* ---------------------------------------------------------------------- */ /* Constructor & Initialization */ /* ---------------------------------------------------------------------- */ -template -GMGPolar::GMGPolar(const PolarGrid& grid, const DomainGeometry& domain_geometry, +IGMGPolar::IGMGPolar(const PolarGrid& grid, const DensityProfileCoefficients& density_profile_coefficients) : grid_(grid) - , domain_geometry_(domain_geometry) , density_profile_coefficients_(density_profile_coefficients) , exact_solution_(nullptr) // General solver output and visualization settings @@ -47,8 +46,7 @@ GMGPolar::GMGPolar(const PolarGrid& grid, const DomainGeometry& LIKWID_REGISTER("Solve"); } -template -void GMGPolar::setSolution(const ExactSolution* exact_solution) +void IGMGPolar::setSolution(const ExactSolution* exact_solution) { exact_solution_ = exact_solution; } @@ -56,24 +54,20 @@ void GMGPolar::setSolution(const ExactSolution* exact_solution) /* ---------------------------------------------------------------------- */ /* General output & visualization */ /* ---------------------------------------------------------------------- */ -template -int GMGPolar::verbose() const +int IGMGPolar::verbose() const { return verbose_; } -template -void GMGPolar::verbose(int verbose) +void IGMGPolar::verbose(int verbose) { verbose_ = verbose; } -template -bool GMGPolar::paraview() const +bool IGMGPolar::paraview() const { return paraview_; } -template -void GMGPolar::paraview(bool paraview) +void IGMGPolar::paraview(bool paraview) { paraview_ = paraview; } @@ -81,24 +75,20 @@ void GMGPolar::paraview(bool paraview) /* ---------------------------------------------------------------------- */ /* Parallelization & threading */ /* ---------------------------------------------------------------------- */ -template -int GMGPolar::maxOpenMPThreads() const +int IGMGPolar::maxOpenMPThreads() const { return max_omp_threads_; } -template -void GMGPolar::maxOpenMPThreads(int max_omp_threads) +void IGMGPolar::maxOpenMPThreads(int max_omp_threads) { max_omp_threads_ = max_omp_threads; } -template -double GMGPolar::threadReductionFactor() const +double IGMGPolar::threadReductionFactor() const { return thread_reduction_factor_; } -template -void GMGPolar::threadReductionFactor(double thread_reduction_factor) +void IGMGPolar::threadReductionFactor(double thread_reduction_factor) { thread_reduction_factor_ = thread_reduction_factor; } @@ -106,46 +96,38 @@ void GMGPolar::threadReductionFactor(double thread_reduction_fac /* ---------------------------------------------------------------------- */ /* Numerical method options */ /* ---------------------------------------------------------------------- */ -template -bool GMGPolar::DirBC_Interior() const +bool IGMGPolar::DirBC_Interior() const { return DirBC_Interior_; } -template -void GMGPolar::DirBC_Interior(bool DirBC_Interior) +void IGMGPolar::DirBC_Interior(bool DirBC_Interior) { DirBC_Interior_ = DirBC_Interior; } -template -StencilDistributionMethod GMGPolar::stencilDistributionMethod() const +StencilDistributionMethod IGMGPolar::stencilDistributionMethod() const { return stencil_distribution_method_; } -template -void GMGPolar::stencilDistributionMethod(StencilDistributionMethod stencil_distribution_method) +void IGMGPolar::stencilDistributionMethod(StencilDistributionMethod stencil_distribution_method) { stencil_distribution_method_ = stencil_distribution_method; } -template -bool GMGPolar::cacheDensityProfileCoefficients() const +bool IGMGPolar::cacheDensityProfileCoefficients() const { return cache_density_profile_coefficients_; } -template -void GMGPolar::cacheDensityProfileCoefficients(bool cache_density_profile_coefficients) +void IGMGPolar::cacheDensityProfileCoefficients(bool cache_density_profile_coefficients) { cache_density_profile_coefficients_ = cache_density_profile_coefficients; } -template -bool GMGPolar::cacheDomainGeometry() const +bool IGMGPolar::cacheDomainGeometry() const { return cache_domain_geometry_; } -template -void GMGPolar::cacheDomainGeometry(bool cache_domain_geometry) +void IGMGPolar::cacheDomainGeometry(bool cache_domain_geometry) { cache_domain_geometry_ = cache_domain_geometry; } @@ -153,90 +135,74 @@ void GMGPolar::cacheDomainGeometry(bool cache_domain_geometry) /* ---------------------------------------------------------------------- */ /* Multigrid controls */ /* ---------------------------------------------------------------------- */ -template -ExtrapolationType GMGPolar::extrapolation() const +ExtrapolationType IGMGPolar::extrapolation() const { return extrapolation_; } -template -void GMGPolar::extrapolation(ExtrapolationType extrapolation) +void IGMGPolar::extrapolation(ExtrapolationType extrapolation) { extrapolation_ = extrapolation; } -template -int GMGPolar::maxLevels() const +int IGMGPolar::maxLevels() const { return max_levels_; } -template -void GMGPolar::maxLevels(int max_levels) +void IGMGPolar::maxLevels(int max_levels) { max_levels_ = max_levels; } -template -MultigridCycleType GMGPolar::multigridCycle() const +MultigridCycleType IGMGPolar::multigridCycle() const { return multigrid_cycle_; } -template -void GMGPolar::multigridCycle(MultigridCycleType multigrid_cycle) +void IGMGPolar::multigridCycle(MultigridCycleType multigrid_cycle) { multigrid_cycle_ = multigrid_cycle; } -template -int GMGPolar::preSmoothingSteps() const +int IGMGPolar::preSmoothingSteps() const { return pre_smoothing_steps_; } -template -void GMGPolar::preSmoothingSteps(int pre_smoothing_steps) +void IGMGPolar::preSmoothingSteps(int pre_smoothing_steps) { pre_smoothing_steps_ = pre_smoothing_steps; } -template -int GMGPolar::postSmoothingSteps() const +int IGMGPolar::postSmoothingSteps() const { return post_smoothing_steps_; } -template -void GMGPolar::postSmoothingSteps(int post_smoothing_steps) +void IGMGPolar::postSmoothingSteps(int post_smoothing_steps) { post_smoothing_steps_ = post_smoothing_steps; } -template -bool GMGPolar::FMG() const +bool IGMGPolar::FMG() const { return FMG_; } -template -void GMGPolar::FMG(bool FMG) +void IGMGPolar::FMG(bool FMG) { FMG_ = FMG; } -template -int GMGPolar::FMG_iterations() const +int IGMGPolar::FMG_iterations() const { return FMG_iterations_; } -template -void GMGPolar::FMG_iterations(int FMG_iterations) +void IGMGPolar::FMG_iterations(int FMG_iterations) { FMG_iterations_ = FMG_iterations; } -template -MultigridCycleType GMGPolar::FMG_cycle() const +MultigridCycleType IGMGPolar::FMG_cycle() const { return FMG_cycle_; } -template -void GMGPolar::FMG_cycle(MultigridCycleType FMG_cycle) +void IGMGPolar::FMG_cycle(MultigridCycleType FMG_cycle) { FMG_cycle_ = FMG_cycle; } @@ -245,46 +211,38 @@ void GMGPolar::FMG_cycle(MultigridCycleType FMG_cycle) /* Iterative solver termination */ /* ---------------------------------------------------------------------- */ -template -int GMGPolar::maxIterations() const +int IGMGPolar::maxIterations() const { return max_iterations_; } -template -void GMGPolar::maxIterations(int maxIterations) +void IGMGPolar::maxIterations(int maxIterations) { max_iterations_ = maxIterations; } -template -ResidualNormType GMGPolar::residualNormType() const +ResidualNormType IGMGPolar::residualNormType() const { return residual_norm_type_; } -template -void GMGPolar::residualNormType(ResidualNormType residualNormType) +void IGMGPolar::residualNormType(ResidualNormType residualNormType) { residual_norm_type_ = residualNormType; } -template -std::optional GMGPolar::absoluteTolerance() const +std::optional IGMGPolar::absoluteTolerance() const { return absolute_tolerance_; } -template -void GMGPolar::absoluteTolerance(std::optional tol) +void IGMGPolar::absoluteTolerance(std::optional tol) { absolute_tolerance_ = tol.has_value() && tol.value() >= 0.0 ? tol : std::nullopt; } -template -std::optional GMGPolar::relativeTolerance() const +std::optional IGMGPolar::relativeTolerance() const { return relative_tolerance_; } -template -void GMGPolar::relativeTolerance(std::optional tol) +void IGMGPolar::relativeTolerance(std::optional tol) { relative_tolerance_ = tol.has_value() && tol.value() >= 0.0 ? tol : std::nullopt; } @@ -292,21 +250,18 @@ void GMGPolar::relativeTolerance(std::optional tol) /* ---------------------------------------------------------------------- */ /* Solution & Grid Access */ /* ---------------------------------------------------------------------- */ -template -Vector GMGPolar::solution() +Vector IGMGPolar::solution() { int level_depth = 0; return levels_[level_depth].solution(); } -template -ConstVector GMGPolar::solution() const +ConstVector IGMGPolar::solution() const { int level_depth = 0; return levels_[level_depth].solution(); } -template -const PolarGrid& GMGPolar::grid() const +const PolarGrid& IGMGPolar::grid() const { return grid_; } @@ -314,28 +269,23 @@ const PolarGrid& GMGPolar::grid() const /* ---------------------------------------------------------------------- */ /* Setup timings */ /* ---------------------------------------------------------------------- */ -template -double GMGPolar::timeSetupTotal() const +double IGMGPolar::timeSetupTotal() const { return t_setup_total_; } -template -double GMGPolar::timeSetupCreateLevels() const +double IGMGPolar::timeSetupCreateLevels() const { return t_setup_createLevels_; } -template -double GMGPolar::timeSetupRHS() const +double IGMGPolar::timeSetupRHS() const { return t_setup_rhs_; } -template -double GMGPolar::timeSetupSmoother() const +double IGMGPolar::timeSetupSmoother() const { return t_setup_smoother_; } -template -double GMGPolar::timeSetupDirectSolver() const +double IGMGPolar::timeSetupDirectSolver() const { return t_setup_directSolver_; } @@ -343,28 +293,23 @@ double GMGPolar::timeSetupDirectSolver() const /* ---------------------------------------------------------------------- */ /* Solve timings */ /* ---------------------------------------------------------------------- */ -template -double GMGPolar::timeSolveTotal() const +double IGMGPolar::timeSolveTotal() const { return t_solve_total_; } -template -double GMGPolar::timeSolveInitialApproximation() const +double IGMGPolar::timeSolveInitialApproximation() const { return t_solve_initial_approximation_; } -template -double GMGPolar::timeSolveMultigridIterations() const +double IGMGPolar::timeSolveMultigridIterations() const { return t_solve_multigrid_iterations_; } -template -double GMGPolar::timeCheckConvergence() const +double IGMGPolar::timeCheckConvergence() const { return t_check_convergence_; } -template -double GMGPolar::timeCheckExactError() const +double IGMGPolar::timeCheckExactError() const { return t_check_exact_error_; } @@ -372,28 +317,23 @@ double GMGPolar::timeCheckExactError() const /* ---------------------------------------------------------------------- */ /* Average Multigrid Cycle timings */ /* ---------------------------------------------------------------------- */ -template -double GMGPolar::timeAvgMGCTotal() const +double IGMGPolar::timeAvgMGCTotal() const { return t_avg_MGC_total_; } -template -double GMGPolar::timeAvgMGCPreSmoothing() const +double IGMGPolar::timeAvgMGCPreSmoothing() const { return t_avg_MGC_preSmoothing_; } -template -double GMGPolar::timeAvgMGCPostSmoothing() const +double IGMGPolar::timeAvgMGCPostSmoothing() const { return t_avg_MGC_postSmoothing_; } -template -double GMGPolar::timeAvgMGCResidual() const +double IGMGPolar::timeAvgMGCResidual() const { return t_avg_MGC_residual_; } -template -double GMGPolar::timeAvgMGCDirectSolver() const +double IGMGPolar::timeAvgMGCDirectSolver() const { return t_avg_MGC_directSolver_; } @@ -402,16 +342,14 @@ double GMGPolar::timeAvgMGCDirectSolver() const /* Reset timings */ /* ---------------------------------------------------------------------- */ -template -void GMGPolar::resetAllTimings() +void IGMGPolar::resetAllTimings() { resetSetupPhaseTimings(); resetSolvePhaseTimings(); resetAvgMultigridCycleTimings(); } -template -void GMGPolar::resetSetupPhaseTimings() +void IGMGPolar::resetSetupPhaseTimings() { t_setup_total_ = 0.0; t_setup_createLevels_ = 0.0; @@ -420,8 +358,7 @@ void GMGPolar::resetSetupPhaseTimings() t_setup_directSolver_ = 0.0; } -template -void GMGPolar::resetSolvePhaseTimings() +void IGMGPolar::resetSolvePhaseTimings() { t_solve_total_ = 0.0; t_solve_initial_approximation_ = 0.0; @@ -430,8 +367,7 @@ void GMGPolar::resetSolvePhaseTimings() t_check_exact_error_ = 0.0; } -template -void GMGPolar::resetAvgMultigridCycleTimings() +void IGMGPolar::resetAvgMultigridCycleTimings() { t_avg_MGC_total_ = 0.0; t_avg_MGC_preSmoothing_ = 0.0; @@ -444,8 +380,7 @@ void GMGPolar::resetAvgMultigridCycleTimings() /* Diagnostics & statistics */ /* ---------------------------------------------------------------------- */ // Print timing breakdown for setup, smoothing, coarse solve, etc. -template -void GMGPolar::printTimings() const +void IGMGPolar::printTimings() const { // t_setup_rhs_ is neither included in t_setup_total_ and t_solve_total_. std::cout << "\n------------------" << std::endl; @@ -474,30 +409,26 @@ void GMGPolar::printTimings() const } // Number of iterations taken by last solve. -template -int GMGPolar::numberOfIterations() const +int IGMGPolar::numberOfIterations() const { return number_of_iterations_; } // Mean residual reduction factor per iteration. -template -double GMGPolar::meanResidualReductionFactor() const +double IGMGPolar::meanResidualReductionFactor() const { return mean_residual_reduction_factor_; } // Error norms (only available if exact solution was set). -template -std::optional GMGPolar::exactErrorWeightedEuclidean() const +std::optional IGMGPolar::exactErrorWeightedEuclidean() const { if (exact_solution_) { return exact_errors_.back().first; } return std::nullopt; } -template -std::optional GMGPolar::exactErrorInfinity() const +std::optional IGMGPolar::exactErrorInfinity() const { if (exact_solution_) { return exact_errors_.back().second; diff --git a/src/GMGPolar/level_interpolation.cpp b/src/GMGPolar/level_interpolation.cpp index 86a3fa2b..57d4bfab 100644 --- a/src/GMGPolar/level_interpolation.cpp +++ b/src/GMGPolar/level_interpolation.cpp @@ -1,6 +1,6 @@ +#include "../../include/GMGPolar/gmgpolar.h" -template -void GMGPolar::prolongation(const int current_level, Vector result, ConstVector x) const +void IGMGPolar::prolongation(const int current_level, Vector result, ConstVector x) const { assert(current_level < number_of_levels_ && 1 <= current_level); if (!interpolation_) @@ -9,8 +9,7 @@ void GMGPolar::prolongation(const int current_level, VectorapplyProlongation(levels_[current_level], levels_[current_level - 1], result, x); } -template -void GMGPolar::restriction(const int current_level, Vector result, ConstVector x) const +void IGMGPolar::restriction(const int current_level, Vector result, ConstVector x) const { assert(current_level < number_of_levels_ - 1 && 0 <= current_level); if (!interpolation_) @@ -19,8 +18,7 @@ void GMGPolar::restriction(const int current_level, VectorapplyRestriction(levels_[current_level], levels_[current_level + 1], result, x); } -template -void GMGPolar::injection(const int current_level, Vector result, ConstVector x) const +void IGMGPolar::injection(const int current_level, Vector result, ConstVector x) const { assert(current_level < number_of_levels_ - 1 && 0 <= current_level); if (!interpolation_) @@ -29,8 +27,7 @@ void GMGPolar::injection(const int current_level, Vector interpolation_->applyInjection(levels_[current_level], levels_[current_level + 1], result, x); } -template -void GMGPolar::extrapolatedProlongation(const int current_level, Vector result, ConstVector x) const +void IGMGPolar::extrapolatedProlongation(const int current_level, Vector result, ConstVector x) const { assert(current_level < number_of_levels_ && 1 <= current_level); if (!interpolation_) @@ -39,8 +36,7 @@ void GMGPolar::extrapolatedProlongation(const int current_level, interpolation_->applyExtrapolatedProlongation(levels_[current_level], levels_[current_level - 1], result, x); } -template -void GMGPolar::extrapolatedRestriction(const int current_level, Vector result, ConstVector x) const +void IGMGPolar::extrapolatedRestriction(const int current_level, Vector result, ConstVector x) const { assert(current_level < number_of_levels_ - 1 && 0 <= current_level); if (!interpolation_) @@ -49,12 +45,11 @@ void GMGPolar::extrapolatedRestriction(const int current_level, interpolation_->applyExtrapolatedRestriction(levels_[current_level], levels_[current_level + 1], result, x); } -template -void GMGPolar::FMGInterpolation(const int current_level, Vector result, ConstVector x) const +void IGMGPolar::FMGInterpolation(const int current_level, Vector result, ConstVector x) const { assert(current_level < number_of_levels_ && 1 <= current_level); if (!interpolation_) throw std::runtime_error("Interpolation not initialized."); interpolation_->applyFMGInterpolation(levels_[current_level], levels_[current_level - 1], result, x); -} +} \ No newline at end of file diff --git a/src/GMGPolar/setup.cpp b/src/GMGPolar/setup.cpp index a1435a4a..dc14afe9 100644 --- a/src/GMGPolar/setup.cpp +++ b/src/GMGPolar/setup.cpp @@ -1,179 +1,6 @@ +#include "../../include/GMGPolar/gmgpolar.h" -template -void GMGPolar::setup() -{ - LIKWID_START("Setup"); - auto start_setup = std::chrono::high_resolution_clock::now(); - - resetSetupPhaseTimings(); - - auto start_setup_createLevels = std::chrono::high_resolution_clock::now(); - - if (stencil_distribution_method_ == StencilDistributionMethod::CPU_TAKE) { - if (!cache_density_profile_coefficients_ || !cache_domain_geometry_) { - throw std::runtime_error("Error: Caching must be enabled for both density profile coefficients and domain " - "geometry in 'Take' implementation strategy."); - } - } - - // -------------------------------- // - // Create the finest mesh (level 0) // - // -------------------------------- // - auto finest_grid = std::make_unique(grid_); - if (paraview_) - writeToVTK("output_finest_grid", *finest_grid); - - if (extrapolation_ != ExtrapolationType::NONE) { - // Check if grid comes from a single uniform refinement - bool is_uniform_refinement = true; - - for (int i_r = 1; i_r < finest_grid->nr() - 1; i_r += 2) { - double mid = 0.5 * (finest_grid->radius(i_r - 1) + finest_grid->radius(i_r + 1)); - if (std::abs(mid - finest_grid->radius(i_r)) > 1e-12) { - is_uniform_refinement = false; - break; - } - } - for (int i_theta = 1; i_theta < finest_grid->ntheta(); i_theta += 2) { - double mid = 0.5 * (finest_grid->theta(i_theta - 1) + finest_grid->theta(i_theta + 1)); - if (std::abs(mid - finest_grid->theta(i_theta)) > 1e-12) { - is_uniform_refinement = false; - break; - } - } - - if (!is_uniform_refinement) { - throw std::runtime_error( - "Extrapolation Error: Finest PolarGrid does not originate from a single uniform refinement."); - } - } - - // ---------------------------------------------------------- // - // Building PolarGrid and LevelCache for all multigrid levels // - // ---------------------------------------------------------- // - number_of_levels_ = chooseNumberOfLevels(*finest_grid); /* Implementation below */ - 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_); - - for (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(levels_[level_depth - 1], *current_grid); - levels_.emplace_back(level_depth, std::move(current_grid), std::move(current_levelCache), extrapolation_, FMG_); - } - - auto end_setup_createLevels = std::chrono::high_resolution_clock::now(); - t_setup_createLevels_ = std::chrono::duration(end_setup_createLevels - start_setup_createLevels).count(); - - if (paraview_) - writeToVTK("output_coarsest_grid", levels_.back().grid()); - - // ----------------------------------------------------------- // - // Initializing the optimal number of threads for OpenMP tasks // - // ----------------------------------------------------------- // - threads_per_level_.resize(number_of_levels_, max_omp_threads_); - for (int level_depth = 0; level_depth < number_of_levels_; level_depth++) { - threads_per_level_[level_depth] = std::max( - 1, - std::min(max_omp_threads_, - static_cast(std::floor(max_omp_threads_ * std::pow(thread_reduction_factor_, level_depth))))); - } - - interpolation_ = std::make_unique(threads_per_level_, DirBC_Interior_); - - if (verbose_ > 0) - printSettings(); - - // ------------------------------------------------------- - // Initializing various operators based on the level index - 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_, threads_per_level_[level_depth], - stencil_distribution_method_); - break; - case ExtrapolationType::IMPLICIT_EXTRAPOLATION: - full_grid_smoothing_ = false; - levels_[level_depth].initializeExtrapolatedSmoothing(domain_geometry_, density_profile_coefficients_, - DirBC_Interior_, threads_per_level_[level_depth], - 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_, threads_per_level_[level_depth], - stencil_distribution_method_); - break; - case ExtrapolationType::COMBINED: - full_grid_smoothing_ = true; - levels_[level_depth].initializeSmoothing(domain_geometry_, density_profile_coefficients_, - DirBC_Interior_, threads_per_level_[level_depth], - stencil_distribution_method_); - levels_[level_depth].initializeExtrapolatedSmoothing(domain_geometry_, density_profile_coefficients_, - DirBC_Interior_, threads_per_level_[level_depth], - stencil_distribution_method_); - break; - default: - full_grid_smoothing_ = false; - levels_[level_depth].initializeSmoothing(domain_geometry_, density_profile_coefficients_, - DirBC_Interior_, threads_per_level_[level_depth], - stencil_distribution_method_); - levels_[level_depth].initializeExtrapolatedSmoothing(domain_geometry_, density_profile_coefficients_, - DirBC_Interior_, threads_per_level_[level_depth], - 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_, - threads_per_level_[level_depth], 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_, threads_per_level_[level_depth], - 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_, - threads_per_level_[level_depth], stencil_distribution_method_); - } - // ------------------- // - // Intermediate levels // - // ------------------- // - else { - auto start_setup_smoother = std::chrono::high_resolution_clock::now(); - levels_[level_depth].initializeSmoothing(domain_geometry_, density_profile_coefficients_, DirBC_Interior_, - threads_per_level_[level_depth], 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_, - threads_per_level_[level_depth], stencil_distribution_method_); - } - } - - auto end_setup = std::chrono::high_resolution_clock::now(); - t_setup_total_ = std::chrono::duration(end_setup - start_setup).count(); - LIKWID_STOP("Setup"); -} - -template -int GMGPolar::chooseNumberOfLevels(const PolarGrid& finestGrid) +int IGMGPolar::chooseNumberOfLevels(const PolarGrid& finestGrid) { const int minRadialNodes = 5; const int minAngularDivisions = 4; @@ -216,8 +43,7 @@ int GMGPolar::chooseNumberOfLevels(const PolarGrid& finestGrid) return levels; } -template -void GMGPolar::printSettings() const +void IGMGPolar::printSettings() const { std::cout << "------------------------------\n"; diff --git a/src/GMGPolar/solver.cpp b/src/GMGPolar/solver.cpp index 4bf83147..9de788fa 100644 --- a/src/GMGPolar/solver.cpp +++ b/src/GMGPolar/solver.cpp @@ -1,11 +1,12 @@ +#include "../../include/GMGPolar/gmgpolar.h" + #include // ============================================================================= // Main Solver Routine // ============================================================================= -template -void GMGPolar::solve(const BoundaryConditions& boundary_conditions, const SourceTerm& source_term) +void IGMGPolar::solve(const BoundaryConditions& boundary_conditions, const SourceTerm& source_term) { auto start_setup_rhs = std::chrono::high_resolution_clock::now(); @@ -184,8 +185,7 @@ void GMGPolar::solve(const BoundaryConditions& boundary_conditio // Solution Initialization // ============================================================================= -template -void GMGPolar::initializeSolution() +void IGMGPolar::initializeSolution() { if (!FMG_) { int start_level_depth = 0; @@ -266,8 +266,7 @@ void GMGPolar::initializeSolution() // Residual Handling Functions // ============================================================================= -template -double GMGPolar::residualNorm(const ResidualNormType& norm_type, const Level& level, ConstVector residual) const +double IGMGPolar::residualNorm(const ResidualNormType& norm_type, const Level& level, ConstVector residual) const { switch (norm_type) { case ResidualNormType::EUCLIDEAN: @@ -281,8 +280,7 @@ double GMGPolar::residualNorm(const ResidualNormType& norm_type, } } -template -void GMGPolar::updateResidualNorms(Level& level, int iteration, double& initial_residual_norm, +void IGMGPolar::updateResidualNorms(Level& level, int iteration, double& initial_residual_norm, double& current_residual_norm, double& current_relative_residual_norm) { level.computeResidual(level.residual(), level.rhs(), level.solution()); @@ -315,8 +313,7 @@ void GMGPolar::updateResidualNorms(Level& level, int iteration, } } -template -void GMGPolar::extrapolatedResidual(const int current_level, Vector residual, +void IGMGPolar::extrapolatedResidual(const int current_level, Vector residual, ConstVector residual_next_level) { const PolarGrid& fineGrid = levels_[current_level].grid(); @@ -371,8 +368,7 @@ void GMGPolar::extrapolatedResidual(const int current_level, Vec // Convergence and Error Analysis Functions // ============================================================================= -template -bool GMGPolar::converged(double residual_norm, double relative_residual_norm) +bool IGMGPolar::converged(double residual_norm, double relative_residual_norm) { if (relative_tolerance_.has_value()) { if (!(relative_residual_norm > relative_tolerance_.value())) { @@ -387,16 +383,14 @@ bool GMGPolar::converged(double residual_norm, double relative_r return false; } -template -void GMGPolar::evaluateExactError(Level& level, const ExactSolution& exact_solution) +void IGMGPolar::evaluateExactError(Level& level, const ExactSolution& exact_solution) { // Compute the weighted L2 norm and infinity norm of the error between the numerical and exact solution. // The results are stored as a pair: (weighted L2 error, infinity error). exact_errors_.push_back(computeExactError(level, level.solution(), level.residual(), exact_solution)); } -template -std::pair GMGPolar::computeExactError(Level& level, ConstVector solution, Vector error, +std::pair IGMGPolar::computeExactError(Level& level, ConstVector solution, Vector error, const ExactSolution& exact_solution) { const PolarGrid& grid = level.grid(); @@ -437,8 +431,7 @@ std::pair GMGPolar::computeExactError(Level& lev // Output and Logging Functions // ============================================================================= -template -void GMGPolar::printIterationHeader(const ExactSolution* exact_solution) +void IGMGPolar::printIterationHeader(const ExactSolution* exact_solution) { if (verbose_ <= 0) return; @@ -464,8 +457,7 @@ void GMGPolar::printIterationHeader(const ExactSolution* exact_s std::cout << std::right << std::setfill(' '); } -template -void GMGPolar::printIterationInfo(int iteration, double current_residual_norm, double current_relative_residual_norm, +void IGMGPolar::printIterationInfo(int iteration, double current_residual_norm, double current_relative_residual_norm, const ExactSolution* exact_solution) { if (verbose_ <= 0) From ae1aee5505099d43e74acd218d6dc69cbf92aa0c Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 17 Dec 2025 15:10:07 +0100 Subject: [PATCH 07/13] Add a method to get a std::unique_ptr from the ConfigParser --- include/ConfigParser/config_parser.h | 7 +++- src/ConfigParser/config_parser.cpp | 2 +- src/ConfigParser/select_test_case.cpp | 24 +++++++++--- src/main.cpp | 54 +++++++++++++-------------- 4 files changed, 52 insertions(+), 35 deletions(-) diff --git a/include/ConfigParser/config_parser.h b/include/ConfigParser/config_parser.h index 080d4490..c8d7f01f 100644 --- a/include/ConfigParser/config_parser.h +++ b/include/ConfigParser/config_parser.h @@ -9,6 +9,8 @@ #include "../../include/common/global_definitions.h" #include "../../include/PolarGrid/polargrid.h" #include "../../include/GMGPolar/test_cases.h" +#include "../../include/GMGPolar/igmgpolar.h" +#include "test_selection.h" class ConfigParser { @@ -22,11 +24,12 @@ class ConfigParser bool parse(int argc, char* argv[]); // Test Case - const DomainGeometry& domainGeometry() const; + const DomainGeometryVariant& domainGeometry() const; const DensityProfileCoefficients& densityProfileCoefficients() const; const BoundaryConditions& boundaryConditions() const; const SourceTerm& sourceTerm() const; const ExactSolution& exactSolution() const; + std::unique_ptr solver() const; // Control Parameters int verbose() const; @@ -58,7 +61,7 @@ class ConfigParser // Parse command-line arguments to extract problem configuration cmdline::parser parser_; // Input Functions - std::unique_ptr domain_geometry_; + std::unique_ptr domain_geometry_; std::unique_ptr density_profile_coefficients_; std::unique_ptr boundary_conditions_; std::unique_ptr source_term_; diff --git a/src/ConfigParser/config_parser.cpp b/src/ConfigParser/config_parser.cpp index 73d3e179..dd0227f9 100644 --- a/src/ConfigParser/config_parser.cpp +++ b/src/ConfigParser/config_parser.cpp @@ -373,7 +373,7 @@ std::optional ConfigParser::relativeTolerance() const return relative_tolerance_; } -const DomainGeometry& ConfigParser::domainGeometry() const +const DomainGeometryVariant& ConfigParser::domainGeometry() const { return *domain_geometry_.get(); } diff --git a/src/ConfigParser/select_test_case.cpp b/src/ConfigParser/select_test_case.cpp index 821afe02..08e62292 100644 --- a/src/ConfigParser/select_test_case.cpp +++ b/src/ConfigParser/select_test_case.cpp @@ -1,4 +1,18 @@ #include "../include/ConfigParser/config_parser.h" +#include "../../include/GMGPolar/gmgpolar.h" + +std::unique_ptr ConfigParser::solver() const +{ + const PolarGrid& grid = grid_; + const DensityProfileCoefficients& density_profile_coefficients = *density_profile_coefficients_; + return std::visit( + [&grid, &density_profile_coefficients](auto const& domain_geometry) { + using DomainGeomType = std::decay_t; + return static_cast>( + std::make_unique>(grid, domain_geometry, density_profile_coefficients)); + }, + *domain_geometry_); +} void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType problem_type, AlphaCoeff alpha_type, BetaCoeff beta_type, double Rmax, double kappa_eps, double delta_e, double alpha_jump) @@ -7,19 +21,19 @@ void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType proble /* Domain Geometry */ switch (geometry_type) { case GeometryType::CIRCULAR: - domain_geometry_ = std::make_unique(Rmax); + domain_geometry_ = std::make_unique(CircularGeometry(Rmax)); break; case GeometryType::SHAFRANOV: - domain_geometry_ = std::make_unique(Rmax, kappa_eps, delta_e); + domain_geometry_ = std::make_unique(ShafranovGeometry(Rmax, kappa_eps, delta_e)); break; case GeometryType::CZARNY: - domain_geometry_ = std::make_unique(Rmax, kappa_eps, delta_e); + domain_geometry_ = std::make_unique(CzarnyGeometry(Rmax, kappa_eps, delta_e)); break; case GeometryType::CULHAM: - domain_geometry_ = std::make_unique(Rmax); + domain_geometry_ = std::make_unique(CulhamGeometry(Rmax)); break; default: @@ -727,4 +741,4 @@ void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType proble default: throw std::runtime_error("Invalid problem.\n"); } -} \ No newline at end of file +} diff --git a/src/main.cpp b/src/main.cpp index b375a46d..4063cd1f 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -14,58 +14,58 @@ int main(int argc, char* argv[]) parser.parse(argc, argv); // Create GMGPolar solver - GMGPolar solver(parser.grid(), parser.domainGeometry(), parser.densityProfileCoefficients()); + std::unique_ptr solver(parser.solver()); // --- General solver output and visualization settings --- // - solver.verbose(parser.verbose()); // Enable/disable verbose output - solver.paraview(parser.paraview()); // Enable/disable ParaView output + solver->verbose(parser.verbose()); // Enable/disable verbose output + solver->paraview(parser.paraview()); // Enable/disable ParaView output // --- Parallelization and threading settings --- // - solver.maxOpenMPThreads(parser.maxOpenMPThreads()); // Maximum OpenMP threads to use - solver.threadReductionFactor(parser.threadReductionFactor()); // Reduce threads on coarser grids + solver->maxOpenMPThreads(parser.maxOpenMPThreads()); // Maximum OpenMP threads to use + solver->threadReductionFactor(parser.threadReductionFactor()); // Reduce threads on coarser grids omp_set_num_threads(parser.maxOpenMPThreads()); // Global OpenMP thread limit // --- Numerical method setup --- // - solver.DirBC_Interior(parser.DirBC_Interior()); // Interior boundary conditions: Dirichlet, Across-the-origin, - solver.stencilDistributionMethod(parser.stencilDistributionMethod()); // Stencil distribution strategy: Take, Give - solver.cacheDensityProfileCoefficients( + solver->DirBC_Interior(parser.DirBC_Interior()); // Interior boundary conditions: Dirichlet, Across-the-origin, + solver->stencilDistributionMethod(parser.stencilDistributionMethod()); // Stencil distribution strategy: Take, Give + solver->cacheDensityProfileCoefficients( parser.cacheDensityProfileCoefficients()); // Cache density profile coefficients: alpha, beta - solver.cacheDomainGeometry(parser.cacheDomainGeometry()); // Cache domain geometry data: arr, att, art, detDF + solver->cacheDomainGeometry(parser.cacheDomainGeometry()); // Cache domain geometry data: arr, att, art, detDF // --- Multigrid settings --- // - solver.extrapolation(parser.extrapolation()); // Enable/disable extrapolation - solver.maxLevels(parser.maxLevels()); // Max multigrid levels (-1 = use deepest possible) - solver.preSmoothingSteps(parser.preSmoothingSteps()); // Smoothing before coarse-grid correction - solver.postSmoothingSteps(parser.postSmoothingSteps()); // Smoothing after coarse-grid correction - solver.multigridCycle(parser.multigridCycle()); // Multigrid cycle type - solver.FMG(parser.FMG()); // Full Multigrid mode on/off - solver.FMG_iterations(parser.FMG_iterations()); // FMG iteration count - solver.FMG_cycle(parser.FMG_cycle()); // FMG cycle type + solver->extrapolation(parser.extrapolation()); // Enable/disable extrapolation + solver->maxLevels(parser.maxLevels()); // Max multigrid levels (-1 = use deepest possible) + solver->preSmoothingSteps(parser.preSmoothingSteps()); // Smoothing before coarse-grid correction + solver->postSmoothingSteps(parser.postSmoothingSteps()); // Smoothing after coarse-grid correction + solver->multigridCycle(parser.multigridCycle()); // Multigrid cycle type + solver->FMG(parser.FMG()); // Full Multigrid mode on/off + solver->FMG_iterations(parser.FMG_iterations()); // FMG iteration count + solver->FMG_cycle(parser.FMG_cycle()); // FMG cycle type // --- Iterative solver controls --- // - solver.maxIterations(parser.maxIterations()); // Max number of iterations - solver.residualNormType(parser.residualNormType()); // Residual norm type (L2, weighted-L2, L∞) - solver.absoluteTolerance(parser.absoluteTolerance()); // Absolute residual tolerance - solver.relativeTolerance(parser.relativeTolerance()); // Relative residual tolerance + solver->maxIterations(parser.maxIterations()); // Max number of iterations + solver->residualNormType(parser.residualNormType()); // Residual norm type (L2, weighted-L2, L∞) + solver->absoluteTolerance(parser.absoluteTolerance()); // Absolute residual tolerance + solver->relativeTolerance(parser.relativeTolerance()); // Relative residual tolerance // --- Finalize solver setup --- // - solver.setup(); // (allocates internal data, prepares operators, etc.) + solver->setup(); // (allocates internal data, prepares operators, etc.) // --- Provide optional exact solution --- // - solver.setSolution(&parser.exactSolution()); + solver->setSolution(&parser.exactSolution()); // --- Solve Phase --- // - solver.solve(parser.boundaryConditions(), parser.sourceTerm()); + solver->solve(parser.boundaryConditions(), parser.sourceTerm()); // --- Retrieve solution and associated grid --- // - Vector solution = solver.solution(); - const PolarGrid& grid = solver.grid(); + Vector solution = solver->solution(); + const PolarGrid& grid = solver->grid(); // Finalize LIKWID performance markers LIKWID_CLOSE(); // Print timing statistics for each solver phase - solver.printTimings(); + solver->printTimings(); return 0; } From bf8cc11833de2f3ff34ed6704666adb128be53f1 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 17 Dec 2025 15:13:06 +0100 Subject: [PATCH 08/13] Clang format --- include/GMGPolar/build_rhs_f.h | 2 +- include/GMGPolar/igmgpolar.h | 6 +++--- include/GMGPolar/setup.h | 2 +- include/GMGPolar/writeToVTK.h | 7 ++++--- .../implicitly_extrapolated_multigrid_F_Cycle.cpp | 2 +- .../implicitly_extrapolated_multigrid_V_Cycle.cpp | 2 +- .../implicitly_extrapolated_multigrid_W_Cycle.cpp | 2 +- src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp | 2 +- src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp | 2 +- src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp | 2 +- src/GMGPolar/build_rhs_f.cpp | 2 +- src/GMGPolar/gmgpolar.cpp | 3 +-- src/GMGPolar/solver.cpp | 11 ++++++----- 13 files changed, 23 insertions(+), 22 deletions(-) diff --git a/include/GMGPolar/build_rhs_f.h b/include/GMGPolar/build_rhs_f.h index 81342f7a..ff2d5f0a 100644 --- a/include/GMGPolar/build_rhs_f.h +++ b/include/GMGPolar/build_rhs_f.h @@ -1,6 +1,6 @@ #include "../../include/GMGPolar/gmgpolar.h" -template +template void GMGPolar::discretize_rhs_f(const Level& level, Vector rhs_f) { const PolarGrid& grid = level.grid(); diff --git a/include/GMGPolar/igmgpolar.h b/include/GMGPolar/igmgpolar.h index 60d8585d..7d72864e 100644 --- a/include/GMGPolar/igmgpolar.h +++ b/include/GMGPolar/igmgpolar.h @@ -36,8 +36,7 @@ class IGMGPolar // Parameters: // - grid: Cartesian mesh discretizing the computational domain. // - density_profile_coefficients: Coefficients \alpha and \beta defining the PDE. - IGMGPolar(const PolarGrid& grid, - const DensityProfileCoefficients& density_profile_coefficients); + IGMGPolar(const PolarGrid& grid, const DensityProfileCoefficients& density_profile_coefficients); /* ---------------------------------------------------------------------- */ /* General output & visualization */ @@ -307,7 +306,8 @@ class IGMGPolar /* ------------- */ /* Visualization */ virtual void writeToVTK(const std::filesystem::path& file_path, const PolarGrid& grid) = 0; - virtual void writeToVTK(const std::filesystem::path& file_path, const Level& level, ConstVector grid_function) = 0; + virtual void writeToVTK(const std::filesystem::path& file_path, const Level& level, + ConstVector grid_function) = 0; /* ------------------------------ */ /* Timing statistics for GMGPolar */ diff --git a/include/GMGPolar/setup.h b/include/GMGPolar/setup.h index e17a11cf..4ac84235 100644 --- a/include/GMGPolar/setup.h +++ b/include/GMGPolar/setup.h @@ -1,5 +1,5 @@ -template +template void GMGPolar::setup() { LIKWID_START("Setup"); diff --git a/include/GMGPolar/writeToVTK.h b/include/GMGPolar/writeToVTK.h index 88a7fc9f..a3cc2b59 100644 --- a/include/GMGPolar/writeToVTK.h +++ b/include/GMGPolar/writeToVTK.h @@ -1,6 +1,6 @@ #include "../../include/GMGPolar/gmgpolar.h" -template +template void GMGPolar::writeToVTK(const std::filesystem::path& file_path, const PolarGrid& grid) { const auto filename = file_path.stem().string() + ".vtu"; @@ -59,8 +59,9 @@ void GMGPolar::writeToVTK(const std::filesystem::path& file_path << "\n"; } -template -void GMGPolar::writeToVTK(const std::filesystem::path& file_path, const Level& level, ConstVector grid_function) +template +void GMGPolar::writeToVTK(const std::filesystem::path& file_path, const Level& level, + ConstVector grid_function) { const PolarGrid& grid = level.grid(); const LevelCache& level_cache = level.levelCache(); diff --git a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp index 7a91d3c8..e7b4acb6 100644 --- a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_F_Cycle.cpp @@ -1,7 +1,7 @@ #include "../../../include/GMGPolar/gmgpolar.h" void IGMGPolar::implicitlyExtrapolatedMultigrid_F_Cycle(const int level_depth, Vector solution, - Vector rhs, Vector residual) + Vector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); diff --git a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp index 00e5443b..66673392 100644 --- a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_V_Cycle.cpp @@ -1,7 +1,7 @@ #include "../../../include/GMGPolar/gmgpolar.h" void IGMGPolar::implicitlyExtrapolatedMultigrid_V_Cycle(const int level_depth, Vector solution, - Vector rhs, Vector residual) + Vector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); diff --git a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp index 5863f356..655efe5a 100644 --- a/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/implicitly_extrapolated_multigrid_W_Cycle.cpp @@ -1,7 +1,7 @@ #include "../../../include/GMGPolar/gmgpolar.h" void IGMGPolar::implicitlyExtrapolatedMultigrid_W_Cycle(const int level_depth, Vector solution, - Vector rhs, Vector residual) + Vector rhs, Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); diff --git a/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp b/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp index 60a65d22..6023aeaa 100644 --- a/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/multigrid_F_Cycle.cpp @@ -1,7 +1,7 @@ #include "../../../include/GMGPolar/gmgpolar.h" void IGMGPolar::multigrid_F_Cycle(const int level_depth, Vector solution, Vector rhs, - Vector residual) + Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); diff --git a/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp b/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp index 1e6ce437..8026c26e 100644 --- a/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/multigrid_V_Cycle.cpp @@ -1,7 +1,7 @@ #include "../../../include/GMGPolar/gmgpolar.h" void IGMGPolar::multigrid_V_Cycle(const int level_depth, Vector solution, Vector rhs, - Vector residual) + Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); diff --git a/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp b/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp index a6cfc599..39f727a0 100644 --- a/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp +++ b/src/GMGPolar/MultigridMethods/multigrid_W_Cycle.cpp @@ -1,7 +1,7 @@ #include "../../../include/GMGPolar/gmgpolar.h" void IGMGPolar::multigrid_W_Cycle(const int level_depth, Vector solution, Vector rhs, - Vector residual) + Vector residual) { assert(0 <= level_depth && level_depth < number_of_levels_ - 1); diff --git a/src/GMGPolar/build_rhs_f.cpp b/src/GMGPolar/build_rhs_f.cpp index ba3c5b85..9bd66462 100644 --- a/src/GMGPolar/build_rhs_f.cpp +++ b/src/GMGPolar/build_rhs_f.cpp @@ -1,7 +1,7 @@ #include "../../include/GMGPolar/gmgpolar.h" void IGMGPolar::build_rhs_f(const Level& level, Vector rhs_f, const BoundaryConditions& boundary_conditions, - const SourceTerm& source_term) + const SourceTerm& source_term) { const PolarGrid& grid = level.grid(); assert(rhs_f.size() == static_cast(grid.numberOfNodes())); diff --git a/src/GMGPolar/gmgpolar.cpp b/src/GMGPolar/gmgpolar.cpp index c74028ad..4c71c680 100644 --- a/src/GMGPolar/gmgpolar.cpp +++ b/src/GMGPolar/gmgpolar.cpp @@ -3,8 +3,7 @@ /* ---------------------------------------------------------------------- */ /* Constructor & Initialization */ /* ---------------------------------------------------------------------- */ -IGMGPolar::IGMGPolar(const PolarGrid& grid, - const DensityProfileCoefficients& density_profile_coefficients) +IGMGPolar::IGMGPolar(const PolarGrid& grid, const DensityProfileCoefficients& density_profile_coefficients) : grid_(grid) , density_profile_coefficients_(density_profile_coefficients) , exact_solution_(nullptr) diff --git a/src/GMGPolar/solver.cpp b/src/GMGPolar/solver.cpp index 9de788fa..e3767cb2 100644 --- a/src/GMGPolar/solver.cpp +++ b/src/GMGPolar/solver.cpp @@ -266,7 +266,8 @@ void IGMGPolar::initializeSolution() // Residual Handling Functions // ============================================================================= -double IGMGPolar::residualNorm(const ResidualNormType& norm_type, const Level& level, ConstVector residual) const +double IGMGPolar::residualNorm(const ResidualNormType& norm_type, const Level& level, + ConstVector residual) const { switch (norm_type) { case ResidualNormType::EUCLIDEAN: @@ -281,7 +282,7 @@ double IGMGPolar::residualNorm(const ResidualNormType& norm_type, const Level& l } void IGMGPolar::updateResidualNorms(Level& level, int iteration, double& initial_residual_norm, - double& current_residual_norm, double& current_relative_residual_norm) + double& current_residual_norm, double& current_relative_residual_norm) { level.computeResidual(level.residual(), level.rhs(), level.solution()); if (extrapolation_ != ExtrapolationType::NONE) { @@ -314,7 +315,7 @@ void IGMGPolar::updateResidualNorms(Level& level, int iteration, double& initial } void IGMGPolar::extrapolatedResidual(const int current_level, Vector residual, - ConstVector residual_next_level) + ConstVector residual_next_level) { const PolarGrid& fineGrid = levels_[current_level].grid(); const PolarGrid& coarseGrid = levels_[current_level + 1].grid(); @@ -391,7 +392,7 @@ void IGMGPolar::evaluateExactError(Level& level, const ExactSolution& exact_solu } std::pair IGMGPolar::computeExactError(Level& level, ConstVector solution, Vector error, - const ExactSolution& exact_solution) + const ExactSolution& exact_solution) { const PolarGrid& grid = level.grid(); const LevelCache& levelCache = level.levelCache(); @@ -458,7 +459,7 @@ void IGMGPolar::printIterationHeader(const ExactSolution* exact_solution) } void IGMGPolar::printIterationInfo(int iteration, double current_residual_norm, double current_relative_residual_norm, - const ExactSolution* exact_solution) + const ExactSolution* exact_solution) { if (verbose_ <= 0) return; From 880d87be72cd37224d95744bf0f7cab65c419fdf Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 17 Dec 2025 15:18:44 +0100 Subject: [PATCH 09/13] Use more common naming strategy --- include/GMGPolar/build_rhs_f.h | 2 +- include/GMGPolar/gmgpolar.h | 2 +- include/GMGPolar/setup.h | 2 +- include/GMGPolar/writeToVTK.h | 4 ++-- include/InputFunctions/domainGeometry.h | 7 ++++++- 5 files changed, 11 insertions(+), 6 deletions(-) diff --git a/include/GMGPolar/build_rhs_f.h b/include/GMGPolar/build_rhs_f.h index ff2d5f0a..031d4bb9 100644 --- a/include/GMGPolar/build_rhs_f.h +++ b/include/GMGPolar/build_rhs_f.h @@ -1,6 +1,6 @@ #include "../../include/GMGPolar/gmgpolar.h" -template +template void GMGPolar::discretize_rhs_f(const Level& level, Vector rhs_f) { const PolarGrid& grid = level.grid(); diff --git a/include/GMGPolar/gmgpolar.h b/include/GMGPolar/gmgpolar.h index 81ec2e38..d8dbdf0d 100644 --- a/include/GMGPolar/gmgpolar.h +++ b/include/GMGPolar/gmgpolar.h @@ -25,7 +25,7 @@ class LevelCache; #include "test_cases.h" #include "igmgpolar.h" -template +template class GMGPolar : public IGMGPolar { public: diff --git a/include/GMGPolar/setup.h b/include/GMGPolar/setup.h index 4ac84235..5372aa27 100644 --- a/include/GMGPolar/setup.h +++ b/include/GMGPolar/setup.h @@ -1,5 +1,5 @@ -template +template void GMGPolar::setup() { LIKWID_START("Setup"); diff --git a/include/GMGPolar/writeToVTK.h b/include/GMGPolar/writeToVTK.h index a3cc2b59..2f79d1b7 100644 --- a/include/GMGPolar/writeToVTK.h +++ b/include/GMGPolar/writeToVTK.h @@ -1,6 +1,6 @@ #include "../../include/GMGPolar/gmgpolar.h" -template +template void GMGPolar::writeToVTK(const std::filesystem::path& file_path, const PolarGrid& grid) { const auto filename = file_path.stem().string() + ".vtu"; @@ -59,7 +59,7 @@ void GMGPolar::writeToVTK(const std::filesystem::path& file_path << "\n"; } -template +template void GMGPolar::writeToVTK(const std::filesystem::path& file_path, const Level& level, ConstVector grid_function) { diff --git a/include/InputFunctions/domainGeometry.h b/include/InputFunctions/domainGeometry.h index 7a7f83fb..86408994 100644 --- a/include/InputFunctions/domainGeometry.h +++ b/include/InputFunctions/domainGeometry.h @@ -30,8 +30,11 @@ class DomainGeometry virtual double dFy_dt(double r, double theta) const = 0; }; +namespace concepts +{ + template -concept DomainGeometryConcept = !std::same_as && requires(const T geom, double r, double theta) { +concept DomainGeometry = !std::same_as && requires(const T geom, double r, double theta) { { geom.Fx(r, theta) } -> std::convertible_to; { geom.Fy(r, theta) } -> std::convertible_to; { geom.dFx_dr(r, theta) } -> std::convertible_to; @@ -39,3 +42,5 @@ concept DomainGeometryConcept = !std::same_as && requires(con { geom.dFx_dt(r, theta) } -> std::convertible_to; { geom.dFy_dt(r, theta) } -> std::convertible_to; }; + +} // namespace concepts From 60a822248611dac9fca45b47debe8729a10c88a8 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 17 Dec 2025 15:20:22 +0100 Subject: [PATCH 10/13] Commit missing file --- include/ConfigParser/test_selection.h | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 include/ConfigParser/test_selection.h diff --git a/include/ConfigParser/test_selection.h b/include/ConfigParser/test_selection.h new file mode 100644 index 00000000..891e9dca --- /dev/null +++ b/include/ConfigParser/test_selection.h @@ -0,0 +1,5 @@ +#pragma once +#include +#include "../../include/GMGPolar/test_cases.h" + +using DomainGeometryVariant = std::variant; From a50ee04f47f39ef4c1a8b50244f00a2d36b273e3 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Sat, 20 Dec 2025 16:40:31 +0100 Subject: [PATCH 11/13] Add explanatory comments --- src/ConfigParser/select_test_case.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/ConfigParser/select_test_case.cpp b/src/ConfigParser/select_test_case.cpp index 08e62292..040c964c 100644 --- a/src/ConfigParser/select_test_case.cpp +++ b/src/ConfigParser/select_test_case.cpp @@ -3,11 +3,17 @@ std::unique_ptr ConfigParser::solver() const { + // Create local aliases so the class doesn't need to be captured by the lamda + // These are references, not copies. const PolarGrid& grid = grid_; const DensityProfileCoefficients& density_profile_coefficients = *density_profile_coefficients_; + + // Create a solver specialized to the active domain geometry. return std::visit( [&grid, &density_profile_coefficients](auto const& domain_geometry) { + // Deduce the concrete geometry type using DomainGeomType = std::decay_t; + // Construct the solver specialized for this geometry type. return static_cast>( std::make_unique>(grid, domain_geometry, density_profile_coefficients)); }, From 1ca2ab563317ee33ca70c77c840d81f6168580cd Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Sat, 20 Dec 2025 16:43:03 +0100 Subject: [PATCH 12/13] Create variable and add comment --- src/ConfigParser/select_test_case.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/ConfigParser/select_test_case.cpp b/src/ConfigParser/select_test_case.cpp index 040c964c..a0ca463f 100644 --- a/src/ConfigParser/select_test_case.cpp +++ b/src/ConfigParser/select_test_case.cpp @@ -14,8 +14,10 @@ std::unique_ptr ConfigParser::solver() const // Deduce the concrete geometry type using DomainGeomType = std::decay_t; // Construct the solver specialized for this geometry type. - return static_cast>( - std::make_unique>(grid, domain_geometry, density_profile_coefficients)); + std::unique_ptr solver = + std::make_unique>(grid, domain_geometry, density_profile_coefficients); + // The lambdas must return objects of identical type + return solver; }, *domain_geometry_); } From e00f0380e85ccce96aa3680c94f71b163e6b004a Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Sat, 20 Dec 2025 16:45:51 +0100 Subject: [PATCH 13/13] Clang format --- src/GMGPolar/setup.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/GMGPolar/setup.cpp b/src/GMGPolar/setup.cpp index 38bf262f..dc14afe9 100644 --- a/src/GMGPolar/setup.cpp +++ b/src/GMGPolar/setup.cpp @@ -1,6 +1,5 @@ #include "../../include/GMGPolar/gmgpolar.h" - int IGMGPolar::chooseNumberOfLevels(const PolarGrid& finestGrid) { const int minRadialNodes = 5;