diff --git a/include/ConfigParser/config_parser.h b/include/ConfigParser/config_parser.h
index c8d7f01f..cb4661b7 100644
--- a/include/ConfigParser/config_parser.h
+++ b/include/ConfigParser/config_parser.h
@@ -25,7 +25,7 @@ class ConfigParser
 
     // Test Case
     const DomainGeometryVariant& domainGeometry() const;
-    const DensityProfileCoefficients& densityProfileCoefficients() const;
+    const DensityProfileCoefficientsVariant& densityProfileCoefficients() const;
     const BoundaryConditions& boundaryConditions() const;
     const SourceTerm& sourceTerm() const;
     const ExactSolution& exactSolution() const;
@@ -62,7 +62,7 @@ class ConfigParser
     cmdline::parser parser_;
     // Input Functions
     std::unique_ptr<const DomainGeometryVariant> domain_geometry_;
-    std::unique_ptr<const DensityProfileCoefficients> density_profile_coefficients_;
+    std::unique_ptr<const DensityProfileCoefficientsVariant> density_profile_coefficients_;
     std::unique_ptr<const BoundaryConditions> boundary_conditions_;
     std::unique_ptr<const SourceTerm> source_term_;
     std::unique_ptr<const ExactSolution> exact_solution_;
diff --git a/include/ConfigParser/test_selection.h b/include/ConfigParser/test_selection.h
index 891e9dca..922a67d7 100644
--- a/include/ConfigParser/test_selection.h
+++ b/include/ConfigParser/test_selection.h
@@ -3,3 +3,7 @@
 #include "../../include/GMGPolar/test_cases.h"
 
 using DomainGeometryVariant = std::variant<CircularGeometry, ShafranovGeometry, CzarnyGeometry, CulhamGeometry>;
+
+using DensityProfileCoefficientsVariant =
+    std::variant<PoissonCoefficients, SonnendruckerCoefficients, SonnendruckerGyroCoefficients, ZoniCoefficients,
+                 ZoniGyroCoefficients, ZoniShiftedCoefficients, ZoniShiftedGyroCoefficients>;
diff --git a/include/GMGPolar/build_rhs_f.h b/include/GMGPolar/build_rhs_f.h
index 031d4bb9..6915c0cd 100644
--- a/include/GMGPolar/build_rhs_f.h
+++ b/include/GMGPolar/build_rhs_f.h
@@ -1,7 +1,7 @@
 #include "../../include/GMGPolar/gmgpolar.h"
 
-template <concepts::DomainGeometry DomainGeometry>
-void GMGPolar<DomainGeometry>::discretize_rhs_f(const Level& level, Vector<double> rhs_f)
+template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
+void GMGPolar<DomainGeometry, DensityProfileCoefficients>::discretize_rhs_f(const Level& level, Vector<double> rhs_f)
 {
     const PolarGrid& grid = level.grid();
     assert(rhs_f.size() == static_cast<uint>(grid.numberOfNodes()));
diff --git a/include/GMGPolar/gmgpolar.h b/include/GMGPolar/gmgpolar.h
index d8dbdf0d..85bc4d26 100644
--- a/include/GMGPolar/gmgpolar.h
+++ b/include/GMGPolar/gmgpolar.h
@@ -1,31 +1,16 @@
 #pragma once
 
-#include <chrono>
 #include <filesystem>
 #include <iostream>
-#include <memory>
 #include <omp.h>
-#include <optional>
 #include <utility>
 
-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"
+
 #include "igmgpolar.h"
 
-template <concepts::DomainGeometry DomainGeometry>
+template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
 class GMGPolar : public IGMGPolar
 {
 public:
@@ -41,8 +26,9 @@ class GMGPolar : public IGMGPolar
     // - density_profile_coefficients: Coefficients \alpha and \beta defining the PDE.
     GMGPolar(const PolarGrid& grid, const DomainGeometry& domain_geometry,
              const DensityProfileCoefficients& density_profile_coefficients)
-        : IGMGPolar(grid, density_profile_coefficients)
+        : IGMGPolar(grid)
         , domain_geometry_(domain_geometry)
+        , density_profile_coefficients_(density_profile_coefficients)
     {
     }
 
@@ -57,6 +43,7 @@ class GMGPolar : public IGMGPolar
     /* Grid Configuration & Input Functions */
     /* ------------------------------------ */
     const DomainGeometry& domain_geometry_;
+    const DensityProfileCoefficients& density_profile_coefficients_;
 
     /* --------------- */
     /* Setup Functions */
diff --git a/include/GMGPolar/igmgpolar.h b/include/GMGPolar/igmgpolar.h
index 7d72864e..7d4d8692 100644
--- a/include/GMGPolar/igmgpolar.h
+++ b/include/GMGPolar/igmgpolar.h
@@ -35,8 +35,7 @@ class IGMGPolar
     // 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);
+    IGMGPolar(const PolarGrid& grid);
 
     /* ---------------------------------------------------------------------- */
     /* General output & visualization                                         */
@@ -194,7 +193,6 @@ class IGMGPolar
     /* Grid Configuration & Input Functions */
     /* ------------------------------------ */
     const PolarGrid& grid_;
-    const DensityProfileCoefficients& density_profile_coefficients_;
     const ExactSolution* exact_solution_; // Optional exact solution for validation
 
     /* ------------------ */
diff --git a/include/GMGPolar/setup.h b/include/GMGPolar/setup.h
index 5372aa27..9a1042b1 100644
--- a/include/GMGPolar/setup.h
+++ b/include/GMGPolar/setup.h
@@ -1,6 +1,6 @@
 
-template <concepts::DomainGeometry DomainGeometry>
-void GMGPolar<DomainGeometry>::setup()
+template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
+void GMGPolar<DomainGeometry, DensityProfileCoefficients>::setup()
 {
     LIKWID_START("Setup");
     auto start_setup = std::chrono::high_resolution_clock::now();
diff --git a/include/GMGPolar/writeToVTK.h b/include/GMGPolar/writeToVTK.h
index 2f79d1b7..690b9efc 100644
--- a/include/GMGPolar/writeToVTK.h
+++ b/include/GMGPolar/writeToVTK.h
@@ -1,7 +1,8 @@
 #include "../../include/GMGPolar/gmgpolar.h"
 
-template <concepts::DomainGeometry DomainGeometry>
-void GMGPolar<DomainGeometry>::writeToVTK(const std::filesystem::path& file_path, const PolarGrid& grid)
+template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
+void GMGPolar<DomainGeometry, DensityProfileCoefficients>::writeToVTK(const std::filesystem::path& file_path,
+                                                                      const PolarGrid& grid)
 {
     const auto filename = file_path.stem().string() + ".vtu";
 
@@ -59,9 +60,10 @@ void GMGPolar<DomainGeometry>::writeToVTK(const std::filesystem::path& file_path
          << "</VTKFile>\n";
 }
 
-template <concepts::DomainGeometry DomainGeometry>
-void GMGPolar<DomainGeometry>::writeToVTK(const std::filesystem::path& file_path, const Level& level,
-                                          ConstVector<double> grid_function)
+template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
+void GMGPolar<DomainGeometry, DensityProfileCoefficients>::writeToVTK(const std::filesystem::path& file_path,
+                                                                      const Level& level,
+                                                                      ConstVector<double> grid_function)
 {
     const PolarGrid& grid         = level.grid();
     const LevelCache& level_cache = level.levelCache();
diff --git a/include/InputFunctions/densityProfileCoefficients.h b/include/InputFunctions/densityProfileCoefficients.h
index 7cc1339b..c7b08353 100644
--- a/include/InputFunctions/densityProfileCoefficients.h
+++ b/include/InputFunctions/densityProfileCoefficients.h
@@ -11,4 +11,17 @@ class DensityProfileCoefficients
 
     // Only used in custom mesh generation -> refinement_radius
     virtual double getAlphaJump() const = 0;
-};
\ No newline at end of file
+};
+
+namespace concepts
+{
+
+template <typename T>
+concept DensityProfileCoefficients =
+    !std::same_as<T, DensityProfileCoefficients> && requires(const T coeffs, double r, double theta) {
+        { coeffs.alpha(r, theta) } -> std::convertible_to<double>;
+        { coeffs.beta(r, theta) } -> std::convertible_to<double>;
+        { coeffs.getAlphaJump() } -> std::convertible_to<double>;
+    };
+
+} // namespace concepts
diff --git a/src/ConfigParser/config_parser.cpp b/src/ConfigParser/config_parser.cpp
index dd0227f9..e61ddb51 100644
--- a/src/ConfigParser/config_parser.cpp
+++ b/src/ConfigParser/config_parser.cpp
@@ -378,7 +378,7 @@ const DomainGeometryVariant& ConfigParser::domainGeometry() const
     return *domain_geometry_.get();
 }
 
-const DensityProfileCoefficients& ConfigParser::densityProfileCoefficients() const
+const DensityProfileCoefficientsVariant& ConfigParser::densityProfileCoefficients() const
 {
     return *density_profile_coefficients_.get();
 }
diff --git a/src/ConfigParser/select_test_case.cpp b/src/ConfigParser/select_test_case.cpp
index 08e62292..4477f338 100644
--- a/src/ConfigParser/select_test_case.cpp
+++ b/src/ConfigParser/select_test_case.cpp
@@ -3,15 +3,17 @@
 
 std::unique_ptr<IGMGPolar> ConfigParser::solver() const
 {
-    const PolarGrid& grid                                          = grid_;
-    const DensityProfileCoefficients& density_profile_coefficients = *density_profile_coefficients_;
+    const PolarGrid& grid = grid_;
     return std::visit(
-        [&grid, &density_profile_coefficients](auto const& domain_geometry) {
-            using DomainGeomType = std::decay_t<decltype(domain_geometry)>;
+        [&grid](auto const& domain_geometry, auto const& density_profile_coefficients) {
+            using DomainGeomType                 = std::decay_t<decltype(domain_geometry)>;
+            using DensityProfileCoefficientsType = std::decay_t<decltype(density_profile_coefficients)>;
+
             return static_cast<std::unique_ptr<IGMGPolar>>(
-                std::make_unique<GMGPolar<DomainGeomType>>(grid, domain_geometry, density_profile_coefficients));
+                std::make_unique<GMGPolar<DomainGeomType, DensityProfileCoefficientsType>>(
+                    grid, domain_geometry, density_profile_coefficients));
         },
-        *domain_geometry_);
+        *domain_geometry_, *density_profile_coefficients_);
 }
 
 void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType problem_type, AlphaCoeff alpha_type,
@@ -44,16 +46,19 @@ void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType proble
     /* Density Profile Coefficients */
     switch (alpha_type) {
     case AlphaCoeff::POISSON:
-        density_profile_coefficients_ = std::make_unique<PoissonCoefficients>(Rmax, alpha_jump);
+        density_profile_coefficients_ =
+            std::make_unique<DensityProfileCoefficientsVariant>(PoissonCoefficients(Rmax, alpha_jump));
         break;
 
     case AlphaCoeff::SONNENDRUCKER:
         switch (beta_type) {
         case BetaCoeff::ZERO:
-            density_profile_coefficients_ = std::make_unique<SonnendruckerCoefficients>(Rmax, alpha_jump);
+            density_profile_coefficients_ =
+                std::make_unique<DensityProfileCoefficientsVariant>(SonnendruckerCoefficients(Rmax, alpha_jump));
             break;
         case BetaCoeff::ALPHA_INVERSE:
-            density_profile_coefficients_ = std::make_unique<SonnendruckerGyroCoefficients>(Rmax, alpha_jump);
+            density_profile_coefficients_ =
+                std::make_unique<DensityProfileCoefficientsVariant>(SonnendruckerGyroCoefficients(Rmax, alpha_jump));
             break;
         default:
             throw std::runtime_error("Invalid beta.\n");
@@ -63,10 +68,12 @@ void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType proble
     case AlphaCoeff::ZONI:
         switch (beta_type) {
         case BetaCoeff::ZERO:
-            density_profile_coefficients_ = std::make_unique<ZoniCoefficients>(Rmax, alpha_jump);
+            density_profile_coefficients_ =
+                std::make_unique<DensityProfileCoefficientsVariant>(ZoniCoefficients(Rmax, alpha_jump));
             break;
         case BetaCoeff::ALPHA_INVERSE:
-            density_profile_coefficients_ = std::make_unique<ZoniGyroCoefficients>(Rmax, alpha_jump);
+            density_profile_coefficients_ =
+                std::make_unique<DensityProfileCoefficientsVariant>(ZoniGyroCoefficients(Rmax, alpha_jump));
             break;
         default:
             throw std::runtime_error("Invalid beta.\n");
@@ -76,10 +83,12 @@ void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType proble
     case AlphaCoeff::ZONI_SHIFTED:
         switch (beta_type) {
         case BetaCoeff::ZERO:
-            density_profile_coefficients_ = std::make_unique<ZoniShiftedCoefficients>(Rmax, alpha_jump);
+            density_profile_coefficients_ =
+                std::make_unique<DensityProfileCoefficientsVariant>(ZoniShiftedCoefficients(Rmax, alpha_jump));
             break;
         case BetaCoeff::ALPHA_INVERSE:
-            density_profile_coefficients_ = std::make_unique<ZoniShiftedGyroCoefficients>(Rmax, alpha_jump);
+            density_profile_coefficients_ =
+                std::make_unique<DensityProfileCoefficientsVariant>(ZoniShiftedGyroCoefficients(Rmax, alpha_jump));
             break;
         default:
             throw std::runtime_error("Invalid beta.\n");
diff --git a/src/GMGPolar/gmgpolar.cpp b/src/GMGPolar/gmgpolar.cpp
index 4c71c680..9f9ed52d 100644
--- a/src/GMGPolar/gmgpolar.cpp
+++ b/src/GMGPolar/gmgpolar.cpp
@@ -3,9 +3,8 @@
 /* ---------------------------------------------------------------------- */
 /* Constructor & Initialization                                           */
 /* ---------------------------------------------------------------------- */
-IGMGPolar::IGMGPolar(const PolarGrid& grid, const DensityProfileCoefficients& density_profile_coefficients)
+IGMGPolar::IGMGPolar(const PolarGrid& grid)
     : grid_(grid)
-    , density_profile_coefficients_(density_profile_coefficients)
     , exact_solution_(nullptr)
     // General solver output and visualization settings
     , verbose_(0)
diff --git a/tests/ConfigParser/config_parser.cpp b/tests/ConfigParser/config_parser.cpp
index 3e1cf89b..3d501aae 100644
--- a/tests/ConfigParser/config_parser.cpp
+++ b/tests/ConfigParser/config_parser.cpp
@@ -198,6 +198,10 @@ TEST_P(ConfigParserTest, ParseAllGeometryAndProblemCombinations)
     EXPECT_DOUBLE_EQ(parser.absoluteTolerance().value(), absoluteTolerance);
     ASSERT_TRUE(parser.relativeTolerance().has_value());
     EXPECT_DOUBLE_EQ(parser.relativeTolerance().value(), relativeTolerance);
+
+    // Solver
+    std::unique_ptr<IGMGPolar> solver = parser.solver();
+    EXPECT_EQ(&solver->grid(), &parser.grid());
 }
 
 // Define test cases covering all combinations
diff --git a/tests/GMGPolar/convergence_order.cpp b/tests/GMGPolar/convergence_order.cpp
index 145d21b0..fa2ae679 100644
--- a/tests/GMGPolar/convergence_order.cpp
+++ b/tests/GMGPolar/convergence_order.cpp
@@ -152,11 +152,11 @@ std::vector<double> refine(std::vector<double> const& original_points)
     return refined;
 }
 
-std::tuple<double, double> get_gmgpolar_error(PolarGrid const& grid, CzarnyGeometry const& domain_geometry,
-                                              DensityProfileCoefficients const& coefficients,
-                                              BoundaryConditions const& boundary_conditions,
-                                              SourceTerm const& source_term, ExactSolution const& solution,
-                                              ExtrapolationType extrapolation)
+template <class DensityProfileCoefficients>
+std::tuple<double, double>
+get_gmgpolar_error(PolarGrid const& grid, CzarnyGeometry const& domain_geometry,
+                   DensityProfileCoefficients const& coefficients, BoundaryConditions const& boundary_conditions,
+                   SourceTerm const& source_term, ExactSolution const& solution, ExtrapolationType extrapolation)
 {
     GMGPolar gmgpolar(grid, domain_geometry, coefficients);
     gmgpolar.setSolution(&solution);