Remove the superclass BoundaryConditions

include/ConfigParser/config_parser.h

@@ -26,7 +26,7 @@ class ConfigParser
     // Test Case
     const DomainGeometryVariant& domainGeometry() const;
     const DensityProfileCoefficientsVariant& densityProfileCoefficients() const;
-    const BoundaryConditions& boundaryConditions() const;
+    const BoundaryConditionsVariant& boundaryConditions() const;
     const SourceTerm& sourceTerm() const;
     const ExactSolution& exactSolution() const;
     std::unique_ptr<IGMGPolar> solver() const;
@@ -62,7 +62,7 @@ class ConfigParser
     // Input Functions
     std::unique_ptr<const DomainGeometryVariant> domain_geometry_;
     std::unique_ptr<const DensityProfileCoefficientsVariant> density_profile_coefficients_;
-    std::unique_ptr<const BoundaryConditions> boundary_conditions_;
+    std::unique_ptr<const BoundaryConditionsVariant> boundary_conditions_;
     std::unique_ptr<const SourceTerm> source_term_;
     std::unique_ptr<const ExactSolution> exact_solution_;
     // General solver output and visualization settings

include/ConfigParser/test_selection.h

@@ -1,9 +1,30 @@
 #pragma once
 #include <variant>
 #include "../../include/GMGPolar/test_cases.h"
+#include "../../include/InputFunctions/BoundaryConditions/cartesianR2_Boundary_CircularGeometry.h"
+#include "../../include/InputFunctions/BoundaryConditions/cartesianR6_Boundary_CzarnyGeometry.h"
+#include "../../include/InputFunctions/BoundaryConditions/polarR6_Boundary_CzarnyGeometry.h"
+#include "../../include/InputFunctions/BoundaryConditions/refined_Boundary_CzarnyGeometry.h"
+#include "../../include/InputFunctions/BoundaryConditions/cartesianR2_Boundary_CzarnyGeometry.h"
+#include "../../include/InputFunctions/BoundaryConditions/cartesianR6_Boundary_ShafranovGeometry.h"
+#include "../../include/InputFunctions/BoundaryConditions/polarR6_Boundary_ShafranovGeometry.h"
+#include "../../include/InputFunctions/BoundaryConditions/refined_Boundary_ShafranovGeometry.h"
+#include "../../include/InputFunctions/BoundaryConditions/cartesianR2_Boundary_ShafranovGeometry.h"
+#include "../../include/InputFunctions/BoundaryConditions/polarR6_Boundary_CircularGeometry.h"
+#include "../../include/InputFunctions/BoundaryConditions/refined_Boundary_CircularGeometry.h"
+#include "../../include/InputFunctions/BoundaryConditions/cartesianR6_Boundary_CircularGeometry.h"
+#include "../../include/InputFunctions/BoundaryConditions/polarR6_Boundary_CulhamGeometry.h"
+#include "../../include/InputFunctions/BoundaryConditions/refined_Boundary_CulhamGeometry.h"
 
 using DomainGeometryVariant = std::variant<CircularGeometry, ShafranovGeometry, CzarnyGeometry, CulhamGeometry>;
 
 using DensityProfileCoefficientsVariant =
     std::variant<PoissonCoefficients, SonnendruckerCoefficients, SonnendruckerGyroCoefficients, ZoniCoefficients,
                  ZoniGyroCoefficients, ZoniShiftedCoefficients, ZoniShiftedGyroCoefficients>;
+
+using BoundaryConditionsVariant = std::variant<
+    CartesianR2_Boundary_CircularGeometry, CartesianR6_Boundary_CzarnyGeometry, PolarR6_Boundary_CzarnyGeometry,
+    Refined_Boundary_CzarnyGeometry, CartesianR2_Boundary_CzarnyGeometry, CartesianR6_Boundary_ShafranovGeometry,
+    PolarR6_Boundary_ShafranovGeometry, Refined_Boundary_ShafranovGeometry, CartesianR2_Boundary_ShafranovGeometry,
+    PolarR6_Boundary_CircularGeometry, Refined_Boundary_CircularGeometry, CartesianR6_Boundary_CircularGeometry,
+    PolarR6_Boundary_CulhamGeometry, Refined_Boundary_CulhamGeometry>;

include/GMGPolar/build_rhs_f.h

@@ -139,3 +139,56 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::discretize_rhs_f(cons
         }
     }
 }
+
+template <concepts::BoundaryConditions BoundaryConditions>
+void IGMGPolar::build_rhs_f(const Level& level, Vector<double> rhs_f, const BoundaryConditions& boundary_conditions,
+                            const SourceTerm& source_term)
+{
+    const PolarGrid& grid = level.grid();
+    assert(std::ssize(rhs_f) == grid.numberOfNodes());
+
+#pragma omp parallel
+    {
+// ----------------------------------------- //
+// Store 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_)) {
+                    rhs_f[grid.index(i_r, i_theta)] = source_term(i_r, i_theta);
+                }
+                else if (i_r == 0 && DirBC_Interior_) {
+                    rhs_f[grid.index(i_r, i_theta)] = boundary_conditions.u_D_Interior(r, theta);
+                }
+                else if (i_r == grid.nr() - 1) {
+                    rhs_f[grid.index(i_r, i_theta)] = boundary_conditions.u_D(r, theta);
+                }
+            }
+        }
+
+// --------------------------------------- //
+// Store rhs values (radial index section) //
+// --------------------------------------- //
+#pragma omp for
+        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_)) {
+                    rhs_f[grid.index(i_r, i_theta)] = source_term(i_r, i_theta);
+                }
+                else if (i_r == 0 && DirBC_Interior_) {
+                    rhs_f[grid.index(i_r, i_theta)] = boundary_conditions.u_D_Interior(r, theta);
+                }
+                else if (i_r == grid.nr() - 1) {
+                    rhs_f[grid.index(i_r, i_theta)] = boundary_conditions.u_D(r, theta);
+                }
+            }
+        }
+    }
+}

include/GMGPolar/igmgpolar.h

@@ -135,6 +135,7 @@ class IGMGPolar
 
     // Solve system with given boundary conditions and source term.
     // Multiple solves with different inputs are supported.
+    template <concepts::BoundaryConditions BoundaryConditions>
     void solve(const BoundaryConditions& boundary_conditions, const SourceTerm& source_term);
 
     /* ---------------------------------------------------------------------- */
@@ -255,6 +256,7 @@ class IGMGPolar
     /* --------------- */
     /* Setup Functions */
     int chooseNumberOfLevels(const PolarGrid& finest_grid);
+    template <concepts::BoundaryConditions BoundaryConditions>
     void build_rhs_f(const Level& level, Vector<double> rhs_f, const BoundaryConditions& boundary_conditions,
                      const SourceTerm& source_term);
     virtual void discretize_rhs_f(const Level& level, Vector<double> rhs_f) = 0;
@@ -326,3 +328,5 @@ class IGMGPolar
     double t_avg_MGC_residual_;
     double t_avg_MGC_directSolver_;
 };
+
+#include "solver.h"

include/GMGPolar/solver.h

@@ -0,0 +1,179 @@
+
+// =============================================================================
+//   Main Solver Routine
+// =============================================================================
+template <concepts::BoundaryConditions BoundaryConditions>
+void IGMGPolar::solve(const BoundaryConditions& boundary_conditions, const SourceTerm& source_term)
+{
+    auto start_setup_rhs = std::chrono::high_resolution_clock::now();
+
+    /* ------------------------------------- */
+    /* Build rhs_f on Level 0 (finest Level) */
+    /* ------------------------------------- */
+    build_rhs_f(levels_[0], levels_[0].rhs(), boundary_conditions, source_term);
+
+    /* ---------------- */
+    /* Discretize rhs_f */
+    /* ---------------- */
+    int initial_rhs_f_levels = FMG_ ? number_of_levels_ : (extrapolation_ == ExtrapolationType::NONE ? 1 : 2);
+    // Loop through the levels, injecting and discretizing rhs
+    for (int level_depth = 0; level_depth < initial_rhs_f_levels; ++level_depth) {
+        Level& current_level = levels_[level_depth];
+        // Inject rhs if there is a next level
+        if (level_depth + 1 < initial_rhs_f_levels) {
+            Level& next_level = levels_[level_depth + 1];
+            injection(level_depth, next_level.rhs(), current_level.rhs());
+        }
+        // Discretize the rhs for the current level
+        discretize_rhs_f(current_level, current_level.rhs());
+    }
+
+    auto end_setup_rhs = std::chrono::high_resolution_clock::now();
+    t_setup_rhs_       = std::chrono::duration<double>(end_setup_rhs - start_setup_rhs).count();
+
+    LIKWID_START("Solve");
+    auto start_solve = std::chrono::high_resolution_clock::now();
+
+    // Clear solve-phase timings
+    resetSolvePhaseTimings();
+
+    /* ---------------------------- */
+    /* Initialize starting solution */
+    /* ---------------------------- */
+    auto start_initial_approximation = std::chrono::high_resolution_clock::now();
+
+    initializeSolution();
+
+    auto end_initial_approximation = std::chrono::high_resolution_clock::now();
+    t_solve_initial_approximation_ =
+        std::chrono::duration<double>(end_initial_approximation - start_initial_approximation).count();
+
+    // These times are included in the initial approximation and don't count towards the multigrid cyclces.
+    resetAvgMultigridCycleTimings();
+
+    /* --------------------------------------- */
+    /* Start Solver at finest level (depth 0)  */
+    /* --------------------------------------- */
+    Level& level = levels_[0];
+
+    number_of_iterations_                 = 0;
+    double initial_residual_norm          = 1.0;
+    double current_residual_norm          = 1.0;
+    double current_relative_residual_norm = 1.0;
+
+    printIterationHeader(exact_solution_);
+
+    while (number_of_iterations_ < max_iterations_) {
+        /* ---------------------------------------------- */
+        /* Test solution against exact solution if given. */
+        /* ---------------------------------------------- */
+        LIKWID_STOP("Solver");
+        auto start_check_exact_error = std::chrono::high_resolution_clock::now();
+
+        if (exact_solution_ != nullptr)
+            evaluateExactError(level, *exact_solution_);
+
+        auto end_check_exact_error = std::chrono::high_resolution_clock::now();
+        t_check_exact_error_ += std::chrono::duration<double>(end_check_exact_error - start_check_exact_error).count();
+        LIKWID_START("Solver");
+
+        /* ---------------------------- */
+        /* Compute convergence criteria */
+        /* ---------------------------- */
+        auto start_check_convergence = std::chrono::high_resolution_clock::now();
+
+        if (absolute_tolerance_.has_value() || relative_tolerance_.has_value()) {
+            updateResidualNorms(level, number_of_iterations_, initial_residual_norm, current_residual_norm,
+                                current_relative_residual_norm);
+        }
+
+        auto end_check_convergence = std::chrono::high_resolution_clock::now();
+        t_check_convergence_ += std::chrono::duration<double>(end_check_convergence - start_check_convergence).count();
+
+        printIterationInfo(number_of_iterations_, current_residual_norm, current_relative_residual_norm,
+                           exact_solution_);
+
+        if (converged(current_residual_norm, current_relative_residual_norm))
+            break;
+
+        /* ----------------------- */
+        /* Perform Multigrid Cycle */
+        /* ----------------------- */
+        auto start_solve_multigrid_iterations = std::chrono::high_resolution_clock::now();
+
+        switch (multigrid_cycle_) {
+        case MultigridCycleType::V_CYCLE:
+            if (extrapolation_ == ExtrapolationType::NONE) {
+                multigrid_V_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual());
+            }
+            else {
+                implicitlyExtrapolatedMultigrid_V_Cycle(level.level_depth(), level.solution(), level.rhs(),
+                                                        level.residual());
+            }
+            break;
+        case MultigridCycleType::W_CYCLE:
+            if (extrapolation_ == ExtrapolationType::NONE) {
+                multigrid_W_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual());
+            }
+            else {
+                implicitlyExtrapolatedMultigrid_W_Cycle(level.level_depth(), level.solution(), level.rhs(),
+                                                        level.residual());
+            }
+            break;
+        case MultigridCycleType::F_CYCLE:
+            if (extrapolation_ == ExtrapolationType::NONE) {
+                multigrid_F_Cycle(level.level_depth(), level.solution(), level.rhs(), level.residual());
+            }
+            else {
+                implicitlyExtrapolatedMultigrid_F_Cycle(level.level_depth(), level.solution(), level.rhs(),
+                                                        level.residual());
+            }
+            break;
+        default:
+            throw std::invalid_argument("Unknown MultigridCycleType");
+        }
+        number_of_iterations_++;
+
+        auto end_solve_multigrid_iterations = std::chrono::high_resolution_clock::now();
+        t_solve_multigrid_iterations_ +=
+            std::chrono::duration<double>(end_solve_multigrid_iterations - start_solve_multigrid_iterations).count();
+    }
+
+    /* ---------------------- */
+    /* Post-solution analysis */
+    /* ---------------------- */
+    if (number_of_iterations_ > 0) {
+        /* --------------------------------------------- */
+        /* Compute the average Multigrid Iteration times */
+        /* --------------------------------------------- */
+        t_avg_MGC_total_ = t_solve_multigrid_iterations_ / number_of_iterations_;
+        t_avg_MGC_preSmoothing_ /= number_of_iterations_;
+        t_avg_MGC_postSmoothing_ /= number_of_iterations_;
+        t_avg_MGC_residual_ /= number_of_iterations_;
+        t_avg_MGC_directSolver_ /= number_of_iterations_;
+
+        /* -------------------------------- */
+        /* Compute the reduction factor rho */
+        /* -------------------------------- */
+        mean_residual_reduction_factor_ =
+            std::pow(current_residual_norm / initial_residual_norm, 1.0 / number_of_iterations_);
+
+        if (verbose_ > 0) {
+            std::cout << "------------------------------\n";
+            std::cout << "Total Iterations: " << number_of_iterations_ << "\n";
+            std::cout << "Reduction Factor: ρ = " << mean_residual_reduction_factor_ << "\n";
+        }
+    }
+
+    auto end_solve = std::chrono::high_resolution_clock::now();
+    t_solve_total_ = std::chrono::duration<double>(end_solve - start_solve).count() - t_check_exact_error_;
+    LIKWID_STOP("Solve");
+
+    if (paraview_) {
+        writeToVTK("output_solution", level, level.solution());
+        if (exact_solution_ != nullptr) {
+            computeExactError(level, level.solution(), level.residual(), *exact_solution_);
+            writeToVTK("output_error", level, level.residual());
+        }
+    }
+}

include/InputFunctions/BoundaryConditions/cartesianR2_Boundary_CircularGeometry.h

@@ -4,16 +4,16 @@
 
 #include "../boundaryConditions.h"
 
-class CartesianR2_Boundary_CircularGeometry : public BoundaryConditions
+class CartesianR2_Boundary_CircularGeometry
 {
 public:
-    CartesianR2_Boundary_CircularGeometry() = default;
     explicit CartesianR2_Boundary_CircularGeometry(double Rmax);
-    virtual ~CartesianR2_Boundary_CircularGeometry() = default;
 
-    double u_D(double r, double theta) const override;
-    double u_D_Interior(double r, double theta) const override;
+    double u_D(double r, double theta) const;
+    double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax = 1.3;
 };
+
+static_assert(concepts::BoundaryConditions<CartesianR2_Boundary_CircularGeometry>);

include/InputFunctions/BoundaryConditions/cartesianR2_Boundary_CzarnyGeometry.h

@@ -4,17 +4,15 @@
 
 #include "../boundaryConditions.h"
 
-class CartesianR2_Boundary_CzarnyGeometry : public BoundaryConditions
+class CartesianR2_Boundary_CzarnyGeometry
 {
 public:
     explicit CartesianR2_Boundary_CzarnyGeometry();
     explicit CartesianR2_Boundary_CzarnyGeometry(double Rmax, double inverse_aspect_ratio_epsilon,
                                                  double ellipticity_e);
 
-    virtual ~CartesianR2_Boundary_CzarnyGeometry() = default;
-
-    double u_D(double r, double theta) const override;
-    double u_D_Interior(double r, double theta) const override;
+    double u_D(double r, double theta) const;
+    double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax                         = 1.3;
@@ -24,3 +22,5 @@ class CartesianR2_Boundary_CzarnyGeometry : public BoundaryConditions
     void initializeGeometry();
     double factor_xi;
 };
+
+static_assert(concepts::BoundaryConditions<CartesianR2_Boundary_CzarnyGeometry>);

include/InputFunctions/BoundaryConditions/cartesianR2_Boundary_ShafranovGeometry.h

@@ -4,18 +4,18 @@
 
 #include "../boundaryConditions.h"
 
-class CartesianR2_Boundary_ShafranovGeometry : public BoundaryConditions
+class CartesianR2_Boundary_ShafranovGeometry
 {
 public:
-    CartesianR2_Boundary_ShafranovGeometry() = default;
     explicit CartesianR2_Boundary_ShafranovGeometry(double Rmax, double elongation_kappa, double shift_delta);
-    virtual ~CartesianR2_Boundary_ShafranovGeometry() = default;
 
-    double u_D(double r, double theta) const override;
-    double u_D_Interior(double r, double theta) const override;
+    double u_D(double r, double theta) const;
+    double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax             = 1.3;
     const double elongation_kappa = 0.3;
     const double shift_delta      = 0.2;
 };
+
+static_assert(concepts::BoundaryConditions<CartesianR2_Boundary_ShafranovGeometry>);