diff --git a/.gitmodules b/.gitmodules
index e69de29bb..ac7a68ff2 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "vendor/GMGPolar"]
+	path = vendor/GMGPolar
+	url = https://github.com/SciCompMod/GMGPolar.git
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9128a0207..604939ab6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -21,6 +21,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Add more labels to memory allocations.
 - Add a `NDIdentityInterpolationBuilder` class.
 - Add a new abstract class `IPolarPoissonLikeSolver`.
+- Add class `GMGPolarPoissonLikeSolver` to allow the use of [GMGPolar](https://github.com/SciCompMod/GMGPolar) as a polar Poisson solver.
 
 ### Fixed
 
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7af55032f..c7fd99b96 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -62,6 +62,12 @@ find_package(LAPACKE REQUIRED)
 
 find_package(Ginkgo 1.8 REQUIRED)
 
+set(GMGPOLAR_BUILD_TESTS OFF CACHE BOOL "")
+add_subdirectory(vendor/GMGPolar EXCLUDE_FROM_ALL SYSTEM)
+# Suppress -Werror=unused-variable inherited via CMAKE_CXX_FLAGS_INIT so that GMGPolar warnings don't break the build.
+target_compile_options(GMGPolarLib PRIVATE -Wno-error=unused-variable)
+
+
 ###############################################################################################
 #                              Build libraries and executables
 ###############################################################################################
diff --git a/src/pde_solvers/CMakeLists.txt b/src/pde_solvers/CMakeLists.txt
index 35e9185a0..eef6382dd 100644
--- a/src/pde_solvers/CMakeLists.txt
+++ b/src/pde_solvers/CMakeLists.txt
@@ -10,6 +10,7 @@ target_link_libraries("pde_solvers"
     INTERFACE
         DDC::core
         DDC::fft
+        GMGPolar::GMGPolarLib
         gslx::coord_transformations
         gslx::data_types
         gslx::matrix_tools
diff --git a/src/pde_solvers/gmg_polar_poisson_like_solver.hpp b/src/pde_solvers/gmg_polar_poisson_like_solver.hpp
new file mode 100644
index 000000000..69aab4e52
--- /dev/null
+++ b/src/pde_solvers/gmg_polar_poisson_like_solver.hpp
@@ -0,0 +1,299 @@
+// SPDX-License-Identifier: MIT
+#pragma once
+#include <cmath>
+#include <vector>
+
+#include <GMGPolar/gmgpolar.h>
+
+#include "ddc_alias_inline_functions.hpp"
+#include "ipolar_poisson_like_solver.hpp"
+
+namespace GMGPolarTools {
+
+/**
+ * @brief Wraps a gyselalibxx coordinate mapping to satisfy the GMGPolar DomainGeometry concept.
+ * @tparam ToPhysicalMapping A mapping from (r, theta) curvilinear coordinates to (x, y) Cartesian.
+ */
+template <class ToPhysicalMapping>
+class MappingToDomainGeometry
+{
+    using R = typename ToPhysicalMapping::curvilinear_tag_r;
+    using Theta = typename ToPhysicalMapping::curvilinear_tag_theta;
+
+    using X = typename ToPhysicalMapping::cartesian_tag_x;
+    using Y = typename ToPhysicalMapping::cartesian_tag_y;
+
+    using R_cov = typename R::Dual;
+    using Theta_cov = typename Theta::Dual;
+
+private:
+    ToPhysicalMapping m_to_physical;
+
+public:
+    /// Construct the wrapper class
+    explicit MappingToDomainGeometry(ToPhysicalMapping to_physical) : m_to_physical(to_physical) {}
+
+    /// X(r, theta)
+    double Fx(const double& r, const double& theta) const
+    {
+        return Coord<X>(m_to_physical(Coord<R, Theta>(r, theta)));
+    }
+    /// Y(r, theta)
+    double Fy(const double& r, const double& theta) const
+    {
+        return Coord<Y>(m_to_physical(Coord<R, Theta>(r, theta)));
+    }
+    /// d/dr X(r, theta)
+    double dFx_dr(const double& r, const double& theta) const
+    {
+        return m_to_physical.template jacobian_component<X, R_cov>(Coord<R, Theta>(r, theta));
+    }
+    /// d/dr Y(r, theta)
+    double dFy_dr(const double& r, const double& theta) const
+    {
+        return m_to_physical.template jacobian_component<Y, R_cov>(Coord<R, Theta>(r, theta));
+    }
+    /// d/(d theta) X(r, theta)
+    double dFx_dt(const double& r, const double& theta) const
+    {
+        return m_to_physical.template jacobian_component<X, Theta_cov>(Coord<R, Theta>(r, theta));
+    }
+    /// d/(d theta) Y(r, theta)
+    double dFy_dt(const double& r, const double& theta) const
+    {
+        return m_to_physical.template jacobian_component<Y, Theta_cov>(Coord<R, Theta>(r, theta));
+    }
+};
+
+/**
+ * @brief Homogeneous Dirichlet boundary conditions satisfying the GMGPolar BoundaryConditions concept.
+ */
+class HomogeneousDirichletBoundaryConditions
+{
+public:
+    /// The value of the solution on the boundary
+    static double u_D(const double& r, const double& theta)
+    {
+        return 0.0;
+    }
+    /// The value of the solution on the inner boundary (at r=rmin). Required for the concept, not needed here.
+    static double u_D_Interior(const double& r, const double& theta)
+    {
+        // Only used if DirBC_Interior = true
+        assert(false);
+        return 0.0;
+    }
+};
+
+/**
+ * @brief Wraps gyselalibxx spline-represented coefficients to satisfy the GMGPolar
+ *        DensityProfileCoefficients concept.
+ * @tparam SplineEvaluator_host A 2D host-space spline evaluator for (BSplinesR, BSplinesTheta).
+ * @tparam BSplinesR The radial B-spline type.
+ * @tparam BSplinesTheta The poloidal B-spline type.
+ */
+template <class SplineEvaluator_host, class BSplinesR, class BSplinesTheta>
+class PolarPoissonLikeCoefficients
+{
+    static_assert(
+            std::is_same_v<typename SplineEvaluator_host::memory_space, Kokkos::HostSpace>,
+            "SplineEvaluator_host must operate on Kokkos::HostSpace");
+
+    using R = typename BSplinesR::continuous_dimension_type;
+    using Theta = typename BSplinesTheta::continuous_dimension_type;
+
+    using DConstSplineRTheta_host
+            = DConstField<IdxRange<BSplinesR, BSplinesTheta>, Kokkos::HostSpace>;
+
+private:
+    SplineEvaluator_host m_evaluator;
+    DConstSplineRTheta_host m_coeff_alpha;
+    DConstSplineRTheta_host m_coeff_beta;
+
+public:
+    /// Build th class instance
+    PolarPoissonLikeCoefficients(
+            SplineEvaluator_host evaluator,
+            DConstSplineRTheta_host coeff_alpha,
+            DConstSplineRTheta_host coeff_beta)
+        : m_evaluator(evaluator)
+        , m_coeff_alpha(coeff_alpha)
+        , m_coeff_beta(coeff_beta)
+    {
+    }
+
+    /// The coefficient alpha in the Poisson-like equation
+    double alpha(const double& r, const double& theta) const
+    {
+        return m_evaluator(Coord<R, Theta>(r, theta), m_coeff_alpha);
+    }
+    /// The coefficient beta in the Poisson-like equation
+    double beta(const double& r, const double& theta) const
+    {
+        return m_evaluator(Coord<R, Theta>(r, theta), m_coeff_beta);
+    }
+
+    /// Required for the concept, only used in custom mesh generation (refinement_radius); not needed here.
+    static double getAlphaJump()
+    {
+        assert(false);
+        return 0.0;
+    }
+};
+
+} // namespace GMGPolarTools
+
+/**
+ * @brief A Poisson-like solver using the GMGPolar multigrid library.
+ *
+ * Solves -∇·(α∇φ) + βφ = ρ on a polar domain with homogeneous Dirichlet BCs
+ * at the outer boundary, using an across-the-origin discretisation at r = 0.
+ *
+ * @tparam ToPhysicalMapping    Mapping from (r,θ) to (x,y).
+ * @tparam GridR                Discrete radial grid.
+ * @tparam GridTheta            Discrete poloidal grid.
+ * @tparam BSplinesR            Radial B-spline space.
+ * @tparam BSplinesTheta        Poloidal B-spline space.
+ * @tparam SplineBuilder_host   2D host-space spline builder for (GridR × GridTheta).
+ * @tparam SplineEvaluator_host 2D host-space spline evaluator for (BSplinesR × BSplinesTheta).
+ */
+template <
+        class ToPhysicalMapping,
+        class GridR,
+        class GridTheta,
+        class BSplinesR,
+        class BSplinesTheta,
+        class SplineBuilder_host,
+        class SplineEvaluator_host>
+class GMGPolarPoissonLikeSolver
+    : public IPolarPoissonLikeSolver<IdxRange<GridR, GridTheta>, IdxRange<GridR, GridTheta>>
+{
+    using IdxRangeR = IdxRange<GridR>;
+    using IdxRangeTheta = IdxRange<GridTheta>;
+    using IdxRangeRTheta = IdxRange<GridR, GridTheta>;
+    using IdxRTheta = Idx<GridR, GridTheta>;
+    using IdxR = Idx<GridR>;
+    using IdxTheta = Idx<GridTheta>;
+    using IdxStepRTheta = IdxStep<GridR, GridTheta>;
+
+    using SplineRThetaMem_host = DFieldMem<IdxRange<BSplinesR, BSplinesTheta>, Kokkos::HostSpace>;
+
+    using DomainGeometry = GMGPolarTools::MappingToDomainGeometry<ToPhysicalMapping>;
+    using DensityCoeffs = GMGPolarTools::
+            PolarPoissonLikeCoefficients<SplineEvaluator_host, BSplinesR, BSplinesTheta>;
+
+private:
+    DomainGeometry const m_domain_geom;
+    SplineBuilder_host const& m_builder;
+    SplineEvaluator_host const& m_evaluator;
+    SplineRThetaMem_host m_coeff_alpha;
+    SplineRThetaMem_host m_coeff_beta;
+    DensityCoeffs const m_density_coeffs;
+
+
+public:
+    /**
+     * @brief Construct a GMGPolarPoissonLikeSolver.
+     *
+     * @param[in] to_physical The mapping from the logical to the physical domain.
+     * @param[in] builder A builder to construct the coefficients of the interpolation.
+     * @param[in] evaluator The evaluator for the interpolation.
+     */
+    GMGPolarPoissonLikeSolver(
+            ToPhysicalMapping to_physical,
+            SplineBuilder_host const& builder,
+            SplineEvaluator_host const& evaluator)
+        : m_domain_geom(to_physical)
+        , m_builder(builder)
+        , m_evaluator(evaluator)
+        , m_coeff_alpha(get_spline_idx_range(m_builder))
+        , m_coeff_beta(get_spline_idx_range(m_builder))
+        , m_density_coeffs(
+                  m_evaluator,
+                  get_const_field(m_coeff_alpha),
+                  get_const_field(m_coeff_beta))
+    {
+    }
+
+    /**
+     * @brief Rebuild the internal spline representations of α and β from grid values.
+     * @param[in] alpha Values of α at the grid interpolation points.
+     * @param[in] beta  Values of β at the grid interpolation points.
+     */
+    void update_coefficients(DConstField<IdxRangeRTheta> alpha, DConstField<IdxRangeRTheta> beta)
+            override
+    {
+        auto alpha_host = ddc::create_mirror_view_and_copy(alpha);
+        auto beta_host = ddc::create_mirror_view_and_copy(beta);
+        m_builder(get_field(m_coeff_alpha), get_const_field(alpha_host));
+        m_builder(get_field(m_coeff_beta), get_const_field(beta_host));
+    }
+
+    /**
+     * @brief Solve the Poisson-like equation.
+     *
+     * @param[out] phi The solution @f$\phi@f$ on the grid.
+     * @param[in]  rho The right-hand side @f$\rho@f$ on the grid.
+     */
+    void operator()(DField<IdxRangeRTheta> phi, DConstField<IdxRangeRTheta> rho) const override
+    {
+        // Copy rho to host
+        auto rho_host = ddc::create_mirror_view_and_copy(rho);
+
+        IdxRangeRTheta idx_range = get_idx_range(phi);
+        IdxRangeR idx_range_r(idx_range);
+        IdxRangeTheta idx_range_theta(idx_range);
+        IdxRangeTheta idx_range_theta_with_poloidal_point(
+                idx_range_theta.front(),
+                idx_range_theta.extents() + 1);
+
+        host_t<DFieldMem<IdxRangeR>> r_coords(idx_range_r);
+        host_t<DFieldMem<IdxRangeTheta>> theta_coords(idx_range_theta_with_poloidal_point);
+        ddcHelper::dump_coordinates(Kokkos::DefaultHostExecutionSpace(), get_field(r_coords));
+        ddcHelper::dump_coordinates(Kokkos::DefaultHostExecutionSpace(), get_field(theta_coords));
+
+        gmgpolar::PolarGrid const polar_grid(
+                r_coords.allocation_kokkos_view(),
+                theta_coords.allocation_kokkos_view());
+
+        gmgpolar::GMGPolar<DomainGeometry, DensityCoeffs>
+                solver(polar_grid, m_domain_geom, m_density_coeffs);
+
+        // ----------------------------------------------------------------
+        // Solver parameters
+        solver.DirBC_Interior(false); // Use across-origin discretisation
+        solver.FMG(true);
+        solver.FMG_iterations(3);
+        solver.FMG_cycle(MultigridCycleType::F_CYCLE);
+        solver.extrapolation(ExtrapolationType::IMPLICIT_EXTRAPOLATION);
+        solver.maxLevels(7);
+        solver.preSmoothingSteps(1);
+        solver.postSmoothingSteps(1);
+        solver.multigridCycle(MultigridCycleType::F_CYCLE);
+        solver.maxIterations(150);
+        solver.residualNormType(ResidualNormType::EUCLIDEAN);
+        solver.absoluteTolerance(1e-50);
+        solver.relativeTolerance(1e-6);
+        // ----------------------------------------------------------------
+
+        solver.setup();
+
+        // Source term: maps GMGPolar (i_r, i_theta) indices to rho grid values
+        GMGPolarTools::HomogeneousDirichletBoundaryConditions const bcs;
+        solver.solve(bcs, rho_host.allocation_kokkos_view());
+
+        // Copy solution back to phi
+        auto phi_host = ddc::create_mirror_view(Kokkos::DefaultHostExecutionSpace(), phi);
+        //Kokkos::View<double*, Kokkos::LayoutRight, Kokkos::HostSpace> solution = solver.solution();
+        Kokkos::View<double*> solution = solver.solution();
+
+        ddc::host_for_each(idx_range, [&](IdxRTheta idx) {
+            IdxStepRTheta offset(idx - idx_range.front());
+            int i_r = ddc::select<GridR>(offset);
+            int i_theta = ddc::select<GridTheta>(offset);
+            phi_host(idx) = solution[polar_grid.index(i_r, i_theta)];
+        });
+
+        ddc::parallel_deepcopy(phi, get_const_field(phi_host));
+    }
+};
diff --git a/tests/geometryRTheta/polar_poisson/CMakeLists.txt b/tests/geometryRTheta/polar_poisson/CMakeLists.txt
index 97f42d8bf..0533b49e3 100644
--- a/tests/geometryRTheta/polar_poisson/CMakeLists.txt
+++ b/tests/geometryRTheta/polar_poisson/CMakeLists.txt
@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MIT
 
 foreach(MAPPING_TYPE "CIRCULAR_MAPPING" "CZARNY_MAPPING")
   foreach(SOLUTION "CURVILINEAR_SOLUTION" "CARTESIAN_SOLUTION")
@@ -32,3 +33,19 @@ foreach(MAPPING_TYPE "CIRCULAR_MAPPING" "CZARNY_MAPPING")
     endif()
   endforeach()
 endforeach()
+
+add_executable(gmg_polar_poisson_like_solver_test
+    gmgpolarpoissonsolver.cpp
+    ../../main.cpp
+)
+target_link_libraries(gmg_polar_poisson_like_solver_test
+    PUBLIC
+        GTest::gtest
+	    GTest::gmock
+        DDC::core
+        DDC::splines
+        gslx::geometry_RTheta
+        gslx::pde_solvers
+)
+
+gtest_discover_tests(gmg_polar_poisson_like_solver_test DISCOVERY_MODE PRE_TEST)
diff --git a/tests/geometryRTheta/polar_poisson/gmgpolarpoissonsolver.cpp b/tests/geometryRTheta/polar_poisson/gmgpolarpoissonsolver.cpp
new file mode 100644
index 000000000..a1bad637e
--- /dev/null
+++ b/tests/geometryRTheta/polar_poisson/gmgpolarpoissonsolver.cpp
@@ -0,0 +1,144 @@
+// SPDX-License-Identifier: MIT
+//
+// Tests the GMGPolarPoissonLikeSolver on a circular domain using a manufactured solution.
+//
+// Exact solution:  phi(r, theta) = C * r^6 * (r-1)^6 * cos(3*theta), C = 1e-4 * 2^12
+//
+// With alpha=1, beta=0, the circular mapping gives the standard polar Laplacian:
+//   -Delta phi = -(phi_rr + phi_r/r + phi_tt/r^2) = rho
+//
+// The test passes if the L-inf error is below 1e-2.
+
+#include <cmath>
+#include <iostream>
+
+#include <ddc/ddc.hpp>
+
+#include <gtest/gtest.h>
+
+#include "circular_to_cartesian.hpp"
+#include "ddc_alias_inline_functions.hpp"
+#include "geometry_r_theta.hpp"
+#include "gmg_polar_poisson_like_solver.hpp"
+#include "mesh_builder.hpp"
+#include "spline_definitions_r_theta.hpp"
+
+using Mapping = CircularToCartesian<R, Theta, X, Y>;
+
+using GMGSolver = GMGPolarPoissonLikeSolver<
+        Mapping,
+        GridR,
+        GridTheta,
+        BSplinesR,
+        BSplinesTheta,
+        SplineRThetaBuilder_host,
+        SplineRThetaEvaluatorNullBound_host>;
+
+namespace {
+
+// C = 1e-4 * 2^12
+constexpr double C_MANUF = 1e-4 * 4096.0;
+// Angular mode number
+constexpr int M_MODE = 3;
+
+double phi_exact(double r, double theta)
+{
+    return C_MANUF * std::pow(r, 6) * std::pow(r - 1.0, 6) * std::cos(M_MODE * theta);
+}
+
+// Analytical RHS: rho = -(f'' + f'/r - m^2 * f / r^2) * cos(m*theta)
+// where f(r) = C * r^6 * (r-1)^6
+double rho_exact(double r, double theta)
+{
+    if (r == 0.0) {
+        return 0.0;
+    }
+    const double f = C_MANUF * std::pow(r, 6) * std::pow(r - 1.0, 6);
+    const double fp = 6.0 * C_MANUF * std::pow(r, 5) * std::pow(r - 1.0, 5) * (2.0 * r - 1.0);
+    const double fpp = 6.0 * C_MANUF * std::pow(r, 4) * std::pow(r - 1.0, 4)
+                       * (5.0 * std::pow(2.0 * r - 1.0, 2) + 2.0 * r * (r - 1.0));
+    return -(fpp + fp / r - static_cast<double>(M_MODE * M_MODE) * f / (r * r))
+           * std::cos(M_MODE * theta);
+}
+
+} // namespace
+
+TEST(GMGPolarIntegration, PoissonEquation)
+{
+    CoordR const r_min(1e-5);
+    CoordR const r_max(1.0);
+    IdxStepR const r_ncells(16 - BSDegreeR + 1); // To have a number of grid cells that is 2**N
+
+    CoordTheta const theta_min(0.0);
+    CoordTheta const theta_max(2.0 * M_PI);
+    IdxStepTheta const theta_ncells(64);
+
+    std::vector<CoordR> r_break_points = build_uniform_break_points(r_min, r_max, r_ncells);
+    std::vector<CoordTheta> theta_break_points
+            = build_uniform_break_points(theta_min, theta_max, theta_ncells);
+
+    ddc::init_discrete_space<BSplinesR>(r_break_points);
+    ddc::init_discrete_space<BSplinesTheta>(theta_break_points);
+    ddc::init_discrete_space<GridR>(SplineInterpPointsR::get_sampling<GridR>());
+    ddc::init_discrete_space<GridTheta>(SplineInterpPointsTheta::get_sampling<GridTheta>());
+
+    IdxRangeR const idx_range_r(SplineInterpPointsR::get_domain<GridR>());
+    IdxRangeTheta const idx_range_theta(SplineInterpPointsTheta::get_domain<GridTheta>());
+    IdxRangeRTheta const grid(idx_range_r, idx_range_theta);
+
+    ddc::NullExtrapolationRule bv_r_min;
+    ddc::NullExtrapolationRule bv_r_max;
+    ddc::PeriodicExtrapolationRule<Theta> bv_theta_min;
+    ddc::PeriodicExtrapolationRule<Theta> bv_theta_max;
+
+    SplineRThetaBuilder_host const builder_host(grid);
+    SplineRThetaEvaluatorNullBound_host const
+            evaluator_host(bv_r_min, bv_r_max, bv_theta_min, bv_theta_max);
+
+    const Mapping mapping;
+
+    GMGSolver solver(mapping, builder_host, evaluator_host);
+
+    // Set alpha = 1, beta = 0 everywhere
+    host_t<DFieldMemRTheta> alpha_host(grid);
+    host_t<DFieldMemRTheta> beta_host(grid);
+    ddc::parallel_fill(get_field(alpha_host), 1.0);
+    ddc::parallel_fill(get_field(beta_host), 0.0);
+
+    DFieldMemRTheta alpha_alloc(grid);
+    DFieldMemRTheta beta_alloc(grid);
+    ddc::parallel_deepcopy(get_field(alpha_alloc), get_const_field(alpha_host));
+    ddc::parallel_deepcopy(get_field(beta_alloc), get_const_field(beta_host));
+
+    solver.update_coefficients(get_const_field(alpha_alloc), get_const_field(beta_alloc));
+
+    // Compute the RHS on the grid
+    host_t<DFieldMemRTheta> rho_host(grid);
+    ddc::host_for_each(grid, [&](IdxRTheta irtheta) {
+        double const r = ddc::coordinate(ddc::select<GridR>(irtheta));
+        double const theta = ddc::coordinate(ddc::select<GridTheta>(irtheta));
+        rho_host(irtheta) = rho_exact(r, theta);
+    });
+
+    DFieldMemRTheta rho_alloc(grid);
+    ddc::parallel_deepcopy(get_field(rho_alloc), get_const_field(rho_host));
+
+    // Solve
+    DFieldMemRTheta phi_alloc(grid);
+    ddc::parallel_fill(get_field(phi_alloc), 0.0);
+    solver(get_field(phi_alloc), get_const_field(rho_alloc));
+
+    // Check L-inf error
+    auto phi_result_host = ddc::create_mirror_view_and_copy(get_const_field(phi_alloc));
+    double max_err = 0.0;
+    ddc::host_for_each(grid, [&](IdxRTheta irtheta) {
+        double const r = ddc::coordinate(ddc::select<GridR>(irtheta));
+        double const theta = ddc::coordinate(ddc::select<GridTheta>(irtheta));
+        double const err = std::abs(phi_result_host(irtheta) - phi_exact(r, theta));
+        max_err = std::max(max_err, err);
+    });
+
+    std::cout << "Max L-inf error: " << max_err << std::endl;
+
+    EXPECT_LT(max_err, 1e-6);
+}
diff --git a/vendor/GMGPolar b/vendor/GMGPolar
new file mode 160000
index 000000000..c2da216ee
--- /dev/null
+++ b/vendor/GMGPolar
@@ -0,0 +1 @@
+Subproject commit c2da216eedebdb011ab74acfe4ab2a0039117e0d