diff --git a/include/GMGPolar/igmgpolar.h b/include/GMGPolar/igmgpolar.h
index 0ba8861f..e44b56e9 100644
--- a/include/GMGPolar/igmgpolar.h
+++ b/include/GMGPolar/igmgpolar.h
@@ -12,8 +12,6 @@ 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"
@@ -22,7 +20,6 @@ class LevelCache;
 #include "../LinearAlgebra/Vector/vector_operations.h"
 #include "../PolarGrid/polargrid.h"
 #include "../Definitions/global_definitions.h"
-#include "test_cases.h"
 
 class IGMGPolar
 {
diff --git a/src/main.cpp b/src/main.cpp
index b1a573ba..452a1aeb 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -6,6 +6,7 @@
 int main(int argc, char* argv[])
 {
     Kokkos::ScopeGuard kokkos_scope(argc, argv);
+
     // Initialize LIKWID markers if enabled
     LIKWID_INIT();
 
@@ -13,69 +14,77 @@ int main(int argc, char* argv[])
     ConfigParser parser;
     parser.parse(argc, argv);
 
-    // Create GMGPolar solver
-    std::unique_ptr<IGMGPolar> 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
-
-    // --- Parallelization and threading settings --- //
-    solver->maxOpenMPThreads(parser.maxOpenMPThreads()); // Maximum OpenMP threads to use
-    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(
-        parser.cacheDensityProfileCoefficients()); // Cache density profile coefficients: alpha, beta
-    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
-
-    // --- Preconditioned Conjugate Gradient settings --- //
-    solver->PCG(parser.PCG()); // Preconditioned Conjugate Gradient mode on/off
-    solver->PCG_FMG(parser.PCG_FMG()); // Use FMG as preconditioner for PCG
-    solver->PCG_FMG_iterations(parser.PCG_FMG_iterations()); // FMG iterations for PCG preconditioner
-    solver->PCG_FMG_cycle(parser.PCG_FMG_cycle()); // FMG cycle type for PCG preconditioner
-    solver->PCG_MG_iterations(parser.PCG_MG_iterations()); // Multigrid iterations for PCG preconditioner
-    solver->PCG_MG_cycle(parser.PCG_MG_cycle()); // Multigrid cycle type for PCG preconditioner
-
-    // --- 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
-
-    // --- Finalize solver setup --- //
-    solver->setup(); // (allocates internal data, prepares operators, etc.)
-
-    // --- Provide optional exact solution --- //
-    solver->setSolution(&parser.exactSolution());
-    // --- Solve Phase --- //
     std::visit(
-        [&](auto const& boundary_condition) {
-            solver->solve(boundary_condition, parser.sourceTerm());
+        [&parser](auto const& domain_geometry, auto const& density_profile_coeffs, auto const& boundary_condition) {
+            // Get the types of the domain geometry and the density profile coefficients
+            using DG = std::decay_t<decltype(domain_geometry)>;
+            using DC = std::decay_t<decltype(density_profile_coeffs)>;
+
+            // Create GMGPolar solver for the selected geometry and coefficient types
+            GMGPolar<DG, DC> solver(parser.grid(), domain_geometry, density_profile_coeffs);
+
+            // --- General solver output and visualization settings --- //
+            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
+            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(
+                parser.cacheDensityProfileCoefficients()); // Cache density profile coefficients: alpha, beta
+            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
+
+            // --- Preconditioned Conjugate Gradient settings --- //
+            solver.PCG(parser.PCG()); // Preconditioned Conjugate Gradient mode on/off
+            solver.PCG_FMG(parser.PCG_FMG()); // Use FMG as preconditioner for PCG
+            solver.PCG_FMG_iterations(parser.PCG_FMG_iterations()); // FMG iterations for PCG preconditioner
+            solver.PCG_FMG_cycle(parser.PCG_FMG_cycle()); // FMG cycle type for PCG preconditioner
+            solver.PCG_MG_iterations(parser.PCG_MG_iterations()); // Multigrid iterations for PCG preconditioner
+            solver.PCG_MG_cycle(parser.PCG_MG_cycle()); // Multigrid cycle type for PCG preconditioner
+
+            // --- 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
+
+            // --- Finalize solver setup --- //
+            solver.setup(); // (allocates internal data, prepares operators, etc.)
+
+            // --- Provide optional exact solution --- //
+            solver.setSolution(&parser.exactSolution());
+
+            // --- Solve Phase --- //
+            solver.solve(boundary_condition, parser.sourceTerm());
+
+            // --- Retrieve solution and associated grid --- //
+            Vector<double> solution = solver.solution();
+            const PolarGrid& grid   = solver.grid();
+
+            // Print timing statistics for each solver phase
+            solver.printTimings();
         },
-        parser.boundaryConditions());
-
-    // --- Retrieve solution and associated grid --- //
-    Vector<double> solution = solver->solution();
-    const PolarGrid& grid   = solver->grid();
+        parser.domainGeometry(), parser.densityProfileCoefficients(), parser.boundaryConditions());
 
     // Finalize LIKWID performance markers
     LIKWID_CLOSE();
 
-    // Print timing statistics for each solver phase
-    solver->printTimings();
-
     return 0;
 }
diff --git a/tests/GMGPolar/convergence_order.cpp b/tests/GMGPolar/convergence_order.cpp
index a75bbbb5..7aa4199b 100644
--- a/tests/GMGPolar/convergence_order.cpp
+++ b/tests/GMGPolar/convergence_order.cpp
@@ -5,6 +5,7 @@
 #include <vector>
 #include <random>
 
+#include "../../include/ConfigParser/test_selection.h"
 #include "../../include/GMGPolar/gmgpolar.h"
 
 template <class T>
diff --git a/tests/GMGPolar/pcg_tests.cpp b/tests/GMGPolar/pcg_tests.cpp
index e1e799fb..a68045bd 100644
--- a/tests/GMGPolar/pcg_tests.cpp
+++ b/tests/GMGPolar/pcg_tests.cpp
@@ -6,6 +6,7 @@
 
 // Including the necessary header from the project
 #include "../../include/GMGPolar/gmgpolar.h"
+#include "../../include/ConfigParser/test_selection.h"
 
 template <class T>
 class PCGTestCase;
diff --git a/tests/GMGPolar/solve_tests.cpp b/tests/GMGPolar/solve_tests.cpp
index dbf55f6c..e6f69c1c 100644
--- a/tests/GMGPolar/solve_tests.cpp
+++ b/tests/GMGPolar/solve_tests.cpp
@@ -5,6 +5,7 @@
 #include <random>
 
 // Including the necessary header from the project
+#include "../../include/ConfigParser/test_selection.h"
 #include "../../include/GMGPolar/gmgpolar.h"
 
 template <class T>