diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8b81b769..6b8cf593 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -75,7 +75,7 @@ to turn off this check")
 endif()
 
 set(GINKGO_CHECKOUT_VERSION
-   "ogl_0600_gko190"
+    "ogl_0600_gko190"
     CACHE STRING "Use specific version of ginkgo")
 
 if(OGL_DEVEL_TOOLS)
@@ -216,13 +216,14 @@ target_sources(
           src/MatrixWrapper/SparsityPattern.cpp
           src/MatrixWrapper/Distributed.cpp
           src/MatrixWrapper/Combination.cpp
+          src/Preconditioner/CoarseSolver.cpp
           src/Repartitioner.cpp
           src/Solver/CG.cpp
-	  src/Solver/PipeCG.cpp
-	  src/Solver/BiCGStab.cpp
-	  src/Solver/GMRES.cpp
-	  src/Solver/Multigrid.cpp
-	  # src/Solver/IR.cpp
+          src/Solver/PipeCG.cpp
+          src/Solver/BiCGStab.cpp
+          src/Solver/GMRES.cpp
+          src/Solver/Multigrid.cpp
+          # src/Solver/IR.cpp
 )
 
 enable_sanitizers(
diff --git a/include/OGL/DevicePersistent/Vector.hpp b/include/OGL/DevicePersistent/Vector.hpp
index c8610f59..7b0c6334 100644
--- a/include/OGL/DevicePersistent/Vector.hpp
+++ b/include/OGL/DevicePersistent/Vector.hpp
@@ -212,6 +212,7 @@ class PersistentVector
         tmp.set_executor(ref_exec);
 
         MPI_Request copy_back_req;
+
         MPI_Iscatterv(tmp.get_data(), repartAllToAll.send_counts.data(),
                       repartAllToAll.send_offsets.data(), MPI_DOUBLE,
                       const_cast<T *>(memory_), repartAllToAll.recv_counts[0],
diff --git a/include/OGL/Preconditioner.hpp b/include/OGL/Preconditioner.hpp
index 2bd3d5f1..c4194f51 100644
--- a/include/OGL/Preconditioner.hpp
+++ b/include/OGL/Preconditioner.hpp
@@ -8,11 +8,19 @@
 
 #include "OGL/DevicePersistent/Base.hpp"
 #include "OGL/MatrixWrapper/Distributed.hpp"
+#include "OGL/Preconditioner/Cholesky.hpp"
+#include "OGL/Preconditioner/ISAI.hpp"
+#include "OGL/Preconditioner/Jacobi.hpp"
+#include "OGL/Preconditioner/LU.hpp"
+#include "OGL/Preconditioner/Multigrid.hpp"
+#include "OGL/Preconditioner/Schwarz.hpp"
 
 #include "fvCFD.H"
 #include "regIOobject.H"
 
 namespace Foam {
+
+
 class Preconditioner {
     using mtx = gko::matrix::Csr<scalar>;
     using bj = gko::preconditioner::Jacobi<scalar, label>;
@@ -52,110 +60,6 @@ class Preconditioner {
           verbose_(verbose)
     {}
 
-    template <typename PrecondFactory>
-    std::shared_ptr<gko::LinOp> wrap_multi_level_schwarz(
-        std::shared_ptr<const gko::LinOp> gkomatrix,
-        std::shared_ptr<gko::Executor> device_exec,
-        std::shared_ptr<PrecondFactory> precond, const dictionary &d,
-        label local_rows) const
-    {
-        using pgm = gko::multigrid::Pgm<scalar, label>;
-        using fc = gko::multigrid::FixedCoarsening<scalar, label>;
-        using solver = gko::solver::Cg<scalar>;
-
-        auto selCoarseRows = d.lookupOrDefault("selCoarseRows", label(5));
-        auto fixedCoarsening =
-            d.lookupOrDefault<Switch>("fixedCoarsening", false);
-        auto coarseWeight = d.lookupOrDefault("coarseWeight", scalar(0.01));
-        auto solveNormC =
-            d.lookupOrDefault("reductionCoarseSolver", label(1e-6));
-        auto maxIterCoarse = d.lookupOrDefault("maxIterCoarse", label(50));
-
-        word msg = "Generate multi level schwarz:\n\tfixedCoarsening " +
-                   std::to_string(fixedCoarsening) + "\n\tselCoarseRows " +
-                   std::to_string(selCoarseRows) + "\n\trelTolCoarse " +
-                   std::to_string(solveNormC) + "\n\tmaxIterCoarse " +
-                   std::to_string(maxIterCoarse) + "\n\tcoarseWeigth" +
-                   std::to_string(coarseWeight);
-        MLOG_0(verbose_, msg)
-
-        auto pre_factory = ras::build().with_local_solver(
-            bj::build().with_skip_sorting(true).with_max_block_size(1u).on(
-                device_exec));
-
-        auto coarse_solver = gko::share(
-            solver::build()
-                .with_preconditioner(pre_factory)
-                .with_criteria(
-                    gko::stop::Iteration::build().with_max_iters(maxIterCoarse),
-                    gko::stop::ResidualNorm<scalar>::build()
-                        .with_reduction_factor(solveNormC))
-                .on(device_exec));
-
-        if (fixedCoarsening) {
-            auto n_rows = local_rows / selCoarseRows;
-            auto sel_rows =
-                gko::array<label>(gko::ReferenceExecutor::create(), n_rows);
-            for (auto i = 0; i < sel_rows.get_size(); i++) {
-                sel_rows.get_data()[i] = selCoarseRows * i;
-            }
-
-            sel_rows.set_executor(device_exec);
-            auto pgm_fac = gko::share(fc::build()
-                                          .with_skip_sorting(true)
-                                          .with_coarse_rows(sel_rows)
-                                          .on(device_exec));
-
-            return gko::share(ras::build()
-                                  .with_local_solver(precond)
-                                  .with_coarse_level(pgm_fac)
-                                  .with_l1_smoother(false)
-                                  .with_coarse_solver(coarse_solver)
-                                  .with_coarse_weight(coarseWeight)
-                                  .on(device_exec)
-                                  ->generate(gkomatrix));
-        } else {
-            auto pgm_fac = gko::share(
-                pgm::build().with_skip_sorting(true).on(device_exec));
-
-            return gko::share(ras::build()
-                                  .with_local_solver(precond)
-                                  .with_coarse_level(pgm_fac)
-                                  .with_l1_smoother(false)
-                                  .with_coarse_solver(coarse_solver)
-                                  .on(device_exec)
-                                  ->generate(gkomatrix));
-        }
-    }
-
-    template <typename PrecondFactory>
-    std::shared_ptr<gko::LinOp> wrap_schwarz(
-        std::shared_ptr<const gko::LinOp> gkomatrix,
-        std::shared_ptr<gko::Executor> device_exec,
-        std::unique_ptr<PrecondFactory> precond) const
-    {
-        auto local = gko::as<RepartDistMatrix>(gkomatrix)->get_local();
-        return gko::share(
-            ras::build()
-                .with_generated_local_solver(precond->generate(local))
-                .on(device_exec)
-                ->generate(
-                    gko::as<RepartDistMatrix>(gkomatrix)->get_dist_mtx()));
-    }
-
-    template <typename PrecondFactory, typename Factorization>
-    std::shared_ptr<gko::LinOp> wrap_schwarz(
-        std::shared_ptr<const gko::LinOp> gkomatrix,
-        std::shared_ptr<gko::Executor> device_exec,
-        std::unique_ptr<PrecondFactory> precond,
-        std::shared_ptr<Factorization> factorization) const
-    {
-        return gko::share(
-            ras::build()
-                .with_generated_local_solver(precond->generate(factorization))
-                .on(device_exec)
-                ->generate(gkomatrix));
-    }
 
     std::shared_ptr<gko::LinOp> init_preconditioner_impl(
         const word name, const dictionary &d,
@@ -167,373 +71,20 @@ class Preconditioner {
             d.lookupOrDefault<Switch>("multiLevelSchwarz", false);
 
         if (name == "BJ") {
-            // TODO for non constant system matrix reuse block pointers
-            label max_block_size(d.lookupOrDefault("maxBlockSize", label(1)));
-            word precision(d.lookupOrDefault("precision", word("double")));
-
-            word msg = "Generate preconditioner " + name + "<" + precision +
-                       "> MaxBlockSize " + std::to_string(max_block_size);
-            MLOG_0(verbose_, msg)
-
-            if (precision == "double") {
-                if (multi_level_schwarz) {
-                    auto pre_factory = gko::share(
-                        dbj::build()
-                            .with_skip_sorting(skip_sorting)
-                            .with_max_block_size(
-                                static_cast<gko::uint32>(max_block_size))
-                            .on(device_exec));
-                    auto gkodistmatrix =
-                        gko::as<RepartDistMatrix>(gkomatrix)->get_dist_matrix();
-                    auto local =
-                        gko::as<RepartDistMatrix>(gkomatrix)->get_local();
-                    auto local_rows = local->get_size()[0];
-                    return wrap_multi_level_schwarz(gkodistmatrix, device_exec,
-                                                    pre_factory, d, local_rows);
-                } else {
-                    auto pre_factory =
-                        dbj::build()
-                            .with_skip_sorting(skip_sorting)
-                            .with_max_block_size(
-                                static_cast<gko::uint32>(max_block_size))
-                            .on(device_exec);
-                    auto gkodistmatrix =
-                        gko::as<RepartDistMatrix>(gkomatrix)->get_dist_matrix();
-                    return wrap_schwarz(gkomatrix, device_exec,
-                                        std::move(pre_factory));
-                }
-            } else {
-                auto pre_factory =
-                    fbj::build()
-                        .with_skip_sorting(skip_sorting)
-                        .with_max_block_size(
-                            static_cast<gko::uint32>(max_block_size))
-                        .on(device_exec);
-                return wrap_schwarz(gkomatrix, device_exec,
-                                    std::move(pre_factory));
-            }
-        }
-        if (name == "ILU") {
-            label iterations(d.lookupOrDefault("iterations", label(0)));
-            word msg = "Generate preconditioner " + name;
-            MLOG_0(verbose_, msg)
-
-            auto factorization_factory =
-                gko::factorization::ParIlu<scalar, label>::build()
-                    .with_skip_sorting(skip_sorting)
-                    .with_iterations(iterations)
-                    .on(device_exec);
-            auto gkodistmatrix =
-                gko::as<RepartDistMatrix>(gkomatrix)->get_dist_matrix();
-            auto factorization = gko::share(factorization_factory->generate(
-                gko::as<gko::experimental::distributed::Matrix<
-                    scalar, label, label>>(gkodistmatrix)
-                    ->get_local_matrix()));
-            auto precond_factory =
-                gko::preconditioner::Ilu<>::build().on(device_exec);
-            return wrap_schwarz(gkodistmatrix, device_exec,
-                                std::move(precond_factory), factorization);
+            return BlockJacobi(device_exec, gkomatrix, d, verbose_).create();
         }
-        if (name == "ILUT") {
-            word msg = "Generate preconditioner " + name;
-            MLOG_0(verbose_, msg)
-            auto gkodistmatrix =
-                gko::as<RepartDistMatrix>(gkomatrix)->get_dist_matrix();
-
-            auto factorization_factory =
-                gko::factorization::ParIlut<scalar, label>::build()
-                    .with_skip_sorting(skip_sorting)
-                    .on(device_exec);
-
-            auto factorization = gko::share(factorization_factory->generate(
-                gko::as<gko::experimental::distributed::Matrix<>>(gkodistmatrix)
-                    ->get_local_matrix()));
-            auto precond_factory =
-                gko::preconditioner::Ilu<>::build().on(device_exec);
-            return wrap_schwarz(gkodistmatrix, device_exec,
-                                std::move(precond_factory), factorization);
+        if (name == "Multigrid") {
+            return Multigrid(device_exec, gkomatrix, d, verbose_)
+                .create(/*db_, exec_handler*/);
         }
-        if (name == "IRILU") {
-            auto trisolve_factory =
-                ir::build()
-                    .with_solver(
-                        bj::build().with_max_block_size(1u).on(device_exec))
-                    .with_criteria(
-                        gko::stop::Iteration::build().with_max_iters(5u).on(
-                            device_exec))
-                    .on(device_exec);
-
-            // Generate an ILU preconditioner factory by setting lower and
-            // upper triangular solver - in this case the previously defined
-            // iterative refinement method.
-            auto precond_factory =
-                gko::preconditioner::Ilu<ir, ir>::build()
-                    .with_l_solver(gko::clone(trisolve_factory))
-                    .with_u_solver(gko::clone(trisolve_factory))
-                    .on(device_exec);
-
-            auto factorization_factory =
-                gko::factorization::Ilu<scalar, label>::build()
-                    .with_skip_sorting(skip_sorting)
-                    .on(device_exec);
-
-            auto factorization = gko::share(factorization_factory->generate(
-                gko::as<gko::experimental::distributed::Matrix<>>(gkomatrix)
-                    ->get_local_matrix()));
-
-            // Use incomplete factors to generate ILU preconditioner
-            return wrap_schwarz(gkomatrix, device_exec,
-                                std::move(precond_factory), factorization);
+        if (name == "ILU") {
+            return LU(device_exec, gkomatrix, d, verbose_).create();
         }
         if (name == "IC") {
-            word msg = "Generate preconditioner " + name;
-            MLOG_0(verbose_, msg)
-
-            auto factorization_factory =
-                gko::factorization::Ic<scalar, label>::build()
-                    .with_skip_sorting(skip_sorting)
-                    .on(device_exec);
-            auto gkodistmatrix =
-                gko::as<RepartDistMatrix>(gkomatrix)->get_dist_matrix();
-            auto factorization = gko::share(factorization_factory->generate(
-                gko::as<gko::experimental::distributed::Matrix<
-                    scalar, label, label>>(gkodistmatrix)
-                    ->get_local_matrix()));
-            auto precond_factory =
-                gko::preconditioner::Ic<>::build().on(device_exec);
-            return wrap_schwarz(gkodistmatrix, device_exec,
-                                std::move(precond_factory), factorization);
-        }
-
-        if (name == "ICT") {
-            bool approx_select(d.lookupOrDefault("approximateSelect", true));
-            word msg = "Generate preconditioner " + name +
-                       " with approximate select " +
-                       std::to_string(approx_select);
-            MLOG_0(verbose_, msg)
-
-            auto factorization_factory =
-                gko::factorization::ParIct<scalar, label>::build()
-                    .with_skip_sorting(skip_sorting)
-                    .on(device_exec);
-
-            auto gkodistmatrix =
-                gko::as<RepartDistMatrix>(gkomatrix)->get_dist_matrix();
-            auto ic_factorization = gko::share(factorization_factory->generate(
-                gko::as<gko::experimental::distributed::Matrix<
-                    scalar, label, label>>(gkodistmatrix)
-                    ->get_local_matrix()));
-
-            auto precond_factory =
-                gko::preconditioner::Ic<>::build().on(device_exec);
-
-            return wrap_schwarz(gkomatrix, device_exec,
-                                std::move(precond_factory), ic_factorization);
+            return Cholesky(device_exec, gkomatrix, d, verbose_).create();
         }
         if (name == "ISAI") {
-            label sparsity_power(d.lookupOrDefault("sparsityPower", label(1)));
-
-            word msg = "Generate preconditioner " + name + " SparsityPower " +
-                       std::to_string(sparsity_power);
-            MLOG_0(verbose_, msg)
-
-            auto pre_factory =
-                gko::preconditioner::Isai<gko::preconditioner::isai_type::spd,
-                                          scalar, label>::build()
-                    .with_skip_sorting(skip_sorting)
-                    .with_sparsity_power(sparsity_power)
-                    .on(device_exec);
-
-            return wrap_schwarz(gkomatrix, device_exec, std::move(pre_factory));
-        }
-        if (name == "GISAI") {
-            label sparsity_power(d.lookupOrDefault("sparsityPower", label(1)));
-
-            word msg = "Generate preconditioner " + name + " SparsityPower " +
-                       std::to_string(sparsity_power);
-            MLOG_0(verbose_, msg)
-
-            auto pre_factory = gko::preconditioner::Isai<
-                                   gko::preconditioner::isai_type::general,
-                                   scalar, label>::build()
-                                   .with_skip_sorting(skip_sorting)
-                                   .with_sparsity_power(sparsity_power)
-                                   .on(device_exec);
-
-            return wrap_schwarz(gkomatrix, device_exec, std::move(pre_factory));
-        }
-        if (name == "Multigrid") {
-            word type = d.lookupOrDefault("type", word("Schwarz"));
-
-            auto maxIterCoarseS(d.lookupOrDefault("maxIterCoarse", label(1)));
-            auto solveNorm = d.lookupOrDefault("relTolCoarse", scalar(1e-6));
-            auto relaxFac = d.lookupOrDefault("relaxationFactor", scalar(0.9));
-            auto cycleName = d.lookupOrDefault("cycle", word("v"));
-            auto maxLevels = d.lookupOrDefault("maxLevels", label(20));
-            auto minRowsC = d.lookupOrDefault("minCoarseRows", label(64000));
-            auto smoother = d.lookupOrDefault("smoother", word("Jacobi"));
-            auto coarseSolver =
-                d.lookupOrDefault("coarseSolver", word("Jacobi"));
-            auto maxIterS = d.lookupOrDefault("maxIterSmoother", label(1));
-
-            gko::solver::multigrid::cycle cycle;
-            if (cycleName == "v") cycle = gko::solver::multigrid::cycle::v;
-            if (cycleName == "w") cycle = gko::solver::multigrid::cycle::w;
-            if (cycleName == "f") cycle = gko::solver::multigrid::cycle::f;
-
-            word msg = "Generate preconditioner: " + name +
-                       "\n\tmaxLevels: " + std::to_string(maxLevels) +
-                       "\n\tminCoarseRows: " + std::to_string(minRowsC) +
-                       "\n\tSmoother: " + smoother +
-                       "\n\trelaxationFactor: " + std::to_string(relaxFac) +
-                       "\n\tmaxIterSmoother: " + std::to_string(maxIterS) +
-                       "\n\tcoarseSolver: " + coarseSolver +
-                       "\n\tmaxIterCoarse: " + std::to_string(maxIterCoarseS) +
-                       "\n\tinnerSolverNorm: " + std::to_string(solveNorm) +
-                       "\n\tcycle: " + cycleName + " type: " + type;
-            MLOG_0(verbose_, msg)
-
-
-            std::shared_ptr<gko::LinOpFactory> bjfac{};
-
-            if (smoother == "Jacobi") {
-                bjfac = bj::build()
-                            .with_max_block_size(1u)
-                            .with_skip_sorting(true)
-                            .on(device_exec);
-            }
-            if (smoother == "SOR") {
-                bjfac = sor::build()
-                            .with_skip_sorting(true)
-                            .with_symmetric(false)
-                            .on(device_exec);
-            }
-            if (smoother == "SSOR") {
-                bjfac = sor::build()
-                            .with_skip_sorting(true)
-                            .with_symmetric(true)
-                            .on(device_exec);
-            }
-            if (bjfac == nullptr) {
-                FatalErrorInFunction << "Unknown smoother: " << smoother
-                                     << "\nValid Choices: Jacobi, SOR, SSOR"
-                                     << abort(FatalError);
-            }
-
-            auto single_it = it::build().with_max_iters(1u);
-            auto coarse_solve_it = gko::stop::Iteration::build().with_max_iters(
-                static_cast<gko::uint32>(maxIterCoarseS));
-            auto coarse_solve_norm =
-                gko::stop::ResidualNorm<scalar>::build().with_reduction_factor(
-                    solveNorm);
-            auto smoother_it = gko::stop::Iteration::build().with_max_iters(
-                static_cast<gko::uint32>(maxIterS));
-
-            auto smoother_gen =
-                type == "Distributed"
-                    ? gko::share(ir::build()
-                                     .with_solver(
-                                         ras::build().with_local_solver(bjfac))
-                                     .with_relaxation_factor(relaxFac)
-                                     .with_criteria(smoother_it)
-                                     .on(device_exec))
-                    : gko::share(ir::build()
-                                     .with_solver(bjfac)
-                                     .with_relaxation_factor(relaxFac)
-                                     .with_criteria(smoother_it)
-                                     .on(device_exec));
-
-            if (type == "Schwarz") {
-                std::shared_ptr<const gko::LinOpFactory> coarsest_solver{};
-                if (coarseSolver == "CG") {
-                    coarsest_solver = gko::share(
-                        cg::build()
-                            .with_preconditioner(bjfac)
-                            .with_criteria(coarse_solve_it, coarse_solve_norm)
-                            .on(device_exec));
-                }
-                if (coarseSolver == "Jacobi") {
-                    coarsest_solver =
-                        gko::share(ir::build()
-                                       .with_solver(bjfac)
-                                       .with_relaxation_factor(relaxFac)
-                                       .with_criteria(coarse_solve_it)
-                                       .on(device_exec));
-                }
-                if (coarsest_solver == nullptr) {
-                    FatalErrorInFunction << "Unknown smoother: " << coarseSolver
-                                         << "\nValid Choices: CG, Jacobi"
-                                         << abort(FatalError);
-                }
-
-                auto pre_factory =
-                    mg::build()
-                        .with_max_levels(static_cast<gko::uint32>(maxLevels))
-                        .with_cycle(cycle)
-                        .with_min_coarse_rows(
-                            static_cast<gko::uint32>(minRowsC))
-                        .with_pre_smoother(smoother_gen)
-                        .with_post_uses_pre(true)
-                        .with_mg_level(
-                            pgm::build().with_deterministic(false).on(
-                                device_exec))
-                        .with_coarsest_solver(coarsest_solver)
-                        .with_criteria(single_it)
-                        .on(device_exec);
-                return wrap_schwarz(gkomatrix, device_exec,
-                                    std::move(pre_factory));
-            }
-
-            if (type == "Distributed") {
-                std::shared_ptr<const gko::LinOpFactory> coarsest_solver{};
-                if (coarseSolver == "CG") {
-                    coarsest_solver = gko::share(
-                        ir::build()
-                            .with_solver(ras::build().with_local_solver(bjfac))
-                            .with_relaxation_factor(relaxFac)
-                            .with_criteria(coarse_solve_it)
-                            .on(device_exec));
-                }
-                if (coarseSolver == "Jacobi") {
-                    coarsest_solver = gko::share(
-                        cg::build()
-                            .with_preconditioner(
-                                ras::build().with_local_solver(bjfac))
-                            .with_criteria(coarse_solve_it, coarse_solve_norm)
-                            .on(device_exec));
-                }
-                if (coarsest_solver == nullptr) {
-                    FatalErrorInFunction << "Unknown smoother: " << coarseSolver
-                                         << "\nValid Choices: CG, Jacobi"
-                                         << abort(FatalError);
-                }
-                auto gkodistmatrix =
-                    gko::as<RepartDistMatrix>(gkomatrix)->get_dist_matrix();
-                auto smoother_gen = gko::share(
-                    ir::build()
-                        .with_solver(ras::build().with_local_solver(bjfac))
-                        .with_relaxation_factor(relaxFac)
-                        .with_criteria(smoother_it)
-                        .on(device_exec));
-                auto ret = gko::share(
-                    gko::solver::Multigrid::build()
-                        .with_max_levels(maxLevels)
-                        .with_mg_level(gko::multigrid::Pgm<scalar>::build()
-                                           .with_deterministic(true))
-                        .with_min_coarse_rows(minRowsC)
-                        .with_coarsest_solver(coarsest_solver)
-                        .with_criteria(it::build().with_max_iters(2u))
-                        .with_smoother_iters(maxIterS)
-                        .with_pre_smoother(smoother_gen)
-                        .with_post_uses_pre(true)
-                        .with_cycle(cycle)
-                        .on(device_exec)
-                        ->generate(gkodistmatrix));
-                return ret;
-            }
-            FatalErrorInFunction << "Unknown Multigrid type: " << type
-                                 << "\nValid Choices: Schwarz, Distributed"
-                                 << abort(FatalError);
+            return ISAI(device_exec, gkomatrix, d, verbose_).create();
         }
         if (name == "none") {
             return {};
@@ -653,4 +204,6 @@ class Preconditioner {
         return generated_precond;
     }
 };
+
+
 }  // namespace Foam
diff --git a/include/OGL/Preconditioner/Cholesky.hpp b/include/OGL/Preconditioner/Cholesky.hpp
new file mode 100644
index 00000000..143a1279
--- /dev/null
+++ b/include/OGL/Preconditioner/Cholesky.hpp
@@ -0,0 +1,103 @@
+// SPDX-FileCopyrightText: 2025 OGL authors
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+#pragma once
+
+#include "OGL/Preconditioner/Schwarz.hpp"
+
+class Cholesky  // : public PreconditionerWrapper
+{
+    using dic = gko::preconditioner::Jacobi<double, label>;
+    using fic = gko::preconditioner::Jacobi<float, label>;
+
+    std::shared_ptr<gko::Executor> exec_;
+    std::shared_ptr<const gko::LinOp> mtx_;
+    const dictionary &d_;
+    const label verbose_;
+    bool skip_sorting_;
+    bool multi_level_schwarz_;
+    word precision_;
+    word factorization_;
+
+
+public:
+    Cholesky(std::shared_ptr<gko::Executor> exec,
+             std::shared_ptr<const gko::LinOp> mtx, const dictionary &d,
+             label verbose)
+        : exec_(exec),
+          mtx_(mtx),
+          d_(d),
+          verbose_(verbose),
+          skip_sorting_(d.lookupOrDefault<Switch>("skipSorting", true)),
+          multi_level_schwarz_(
+              d.lookupOrDefault<Switch>("multiLevelSchwarz", false)),
+          precision_(d.lookupOrDefault("precision", word("double"))),
+          factorization_(d.lookupOrDefault("factorization", word("IC")))
+    {
+        word msg = "Generate " + factorization_ +
+                   "  preconditioner:\n\tprecision: " + precision_;
+        MLOG_0(verbose_, msg)
+    }
+
+
+    std::shared_ptr<gko::LinOp> generate_factorization() const
+    {
+        if (factorization_ == "IC") {
+            auto factorization_factory =
+                gko::factorization::Ic<scalar, label>::build()
+                    .with_skip_sorting(skip_sorting_)
+                    .on(exec_);
+
+            auto gkodistmatrix =
+                gko::as<RepartDistMatrix>(mtx_)->get_dist_matrix();
+            return factorization_factory->generate(
+                gko::as<gko::experimental::distributed::Matrix<scalar, label,
+                                                               label>>(
+                    gkodistmatrix)
+                    ->get_local_matrix());
+        }
+        if (factorization_ == "ParIC") {
+            label iterations = d_.lookupOrDefault("iterations", label(5));
+            auto factorization_factory =
+                gko::factorization::ParIct<scalar, label>::build()
+                    .with_skip_sorting(skip_sorting_)
+                    .with_iterations(iterations)
+                    .on(exec_);
+
+            auto gkodistmatrix =
+                gko::as<RepartDistMatrix>(mtx_)->get_dist_matrix();
+            return factorization_factory->generate(
+                gko::as<gko::experimental::distributed::Matrix<scalar, label,
+                                                               label>>(
+                    gkodistmatrix)
+                    ->get_local_matrix());
+        }
+        if (factorization_ == "ParICT") {
+            label iterations = d_.lookupOrDefault("iterations", label(5));
+            label fillInLimit = d_.lookupOrDefault("fillInLimit", label(2));
+            auto factorization_factory =
+                gko::factorization::ParIct<scalar, label>::build()
+                    .with_skip_sorting(skip_sorting_)
+                    .with_fill_in_limit(fillInLimit)
+                    .with_iterations(iterations)
+                    .on(exec_);
+
+            auto gkodistmatrix =
+                gko::as<RepartDistMatrix>(mtx_)->get_dist_matrix();
+            return factorization_factory->generate(
+                gko::as<gko::experimental::distributed::Matrix<scalar, label,
+                                                               label>>(
+                    gkodistmatrix)
+                    ->get_local_matrix());
+        }
+    }
+
+    virtual std::shared_ptr<gko::LinOp> create()
+    {
+        auto precond_factory = gko::preconditioner::Ic<>::build().on(exec_);
+        return dispatch_schwarz(
+            mtx_, exec_,
+            gko::share(precond_factory->generate(generate_factorization())), d_,
+            verbose_);
+    }
+};
diff --git a/include/OGL/Preconditioner/CoarseSolver.hpp b/include/OGL/Preconditioner/CoarseSolver.hpp
new file mode 100644
index 00000000..0f803a3c
--- /dev/null
+++ b/include/OGL/Preconditioner/CoarseSolver.hpp
@@ -0,0 +1,11 @@
+// SPDX-FileCopyrightText: 2024 OGL authors
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+#pragma once
+
+#include <ginkgo/ginkgo.hpp>
+
+#include "fvCFD.H"
+
+std::shared_ptr<const gko::LinOpFactory> generate_coarse_solver(
+    std::shared_ptr<gko::Executor> exec, const dictionary &d, label verbose);
diff --git a/include/OGL/Preconditioner/ISAI.hpp b/include/OGL/Preconditioner/ISAI.hpp
new file mode 100644
index 00000000..3538836f
--- /dev/null
+++ b/include/OGL/Preconditioner/ISAI.hpp
@@ -0,0 +1,81 @@
+// SPDX-FileCopyrightText: 2025 OGL authors
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+#pragma once
+
+#include "OGL/Preconditioner/Schwarz.hpp"
+
+class ISAI  // : public PreconditionerWrapper
+{
+    using dic = gko::preconditioner::Jacobi<double, label>;
+    using fic = gko::preconditioner::Jacobi<float, label>;
+
+    std::shared_ptr<gko::Executor> exec_;
+    std::shared_ptr<const gko::LinOp> mtx_;
+    const dictionary &d_;
+    const label verbose_;
+    bool skip_sorting_;
+    bool multi_level_schwarz_;
+    word type_;
+    label sparsityPower_;
+
+
+public:
+    ISAI(std::shared_ptr<gko::Executor> exec,
+         std::shared_ptr<const gko::LinOp> mtx, const dictionary &d,
+         label verbose)
+        : exec_(exec),
+          mtx_(mtx),
+          d_(d),
+          verbose_(verbose),
+          skip_sorting_(d.lookupOrDefault<Switch>("skipSorting", true)),
+          multi_level_schwarz_(
+              d.lookupOrDefault<Switch>("multiLevelSchwarz", false)),
+          type_(d.lookupOrDefault("type", word("general"))),
+          sparsityPower_(d.lookupOrDefault("sparsityPower", label(1)))
+    {
+        word msg = "Generate " + type_ + "ISAI" +
+                   "  preconditioner:\n\tsparsityPower: " +
+                   std::to_string(sparsityPower_);
+        MLOG_0(verbose_, msg)
+    }
+
+    auto generate_precond_factory_spd()
+    {
+        return gko::preconditioner::Isai<gko::preconditioner::isai_type::spd,
+                                         scalar, label>::build()
+            .with_skip_sorting(skip_sorting_)
+            .with_sparsity_power(sparsityPower_)
+            .on(exec_);
+    }
+
+    auto generate_precond_factory_general()
+    {
+        return gko::preconditioner::Isai<
+                   gko::preconditioner::isai_type::general, scalar,
+                   label>::build()
+            .with_skip_sorting(skip_sorting_)
+            .with_sparsity_power(sparsityPower_)
+            .on(exec_);
+    }
+
+    virtual std::shared_ptr<gko::LinOp> create()
+    {
+        if (type_ == "SPD") {
+            auto b = generate_precond_factory_spd();
+            return dispatch_schwarz(
+                mtx_, exec_,
+                gko::share(
+                    b->generate(gko::as<RepartDistMatrix>(mtx_)->get_local())),
+                d_, verbose_);
+        }
+        if (type_ == "General") {
+            auto b = generate_precond_factory_general();
+            return dispatch_schwarz(
+                mtx_, exec_,
+                gko::share(
+                    b->generate(gko::as<RepartDistMatrix>(mtx_)->get_local())),
+                d_, verbose_);
+        }
+    }
+};
diff --git a/include/OGL/Preconditioner/Jacobi.hpp b/include/OGL/Preconditioner/Jacobi.hpp
new file mode 100644
index 00000000..ca71027f
--- /dev/null
+++ b/include/OGL/Preconditioner/Jacobi.hpp
@@ -0,0 +1,70 @@
+// SPDX-FileCopyrightText: 2025 OGL authors
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+#pragma once
+
+#include "OGL/Preconditioner/Schwarz.hpp"
+
+class BlockJacobi  // : public PreconditionerWrapper
+{
+    using dbj = gko::preconditioner::Jacobi<double, label>;
+    using fbj = gko::preconditioner::Jacobi<float, label>;
+
+    std::shared_ptr<gko::Executor> exec_;
+    std::shared_ptr<const gko::LinOp> mtx_;
+    const dictionary &d_;
+    const label verbose_;
+    bool skip_sorting_;
+    bool multi_level_schwarz_;
+    label max_block_size_;
+    word precision_;
+
+
+public:
+    BlockJacobi(std::shared_ptr<gko::Executor> exec,
+                std::shared_ptr<const gko::LinOp> mtx, const dictionary &d,
+                label verbose)
+        : exec_(exec),
+          mtx_(mtx),
+          d_(d),
+          verbose_(verbose),
+          skip_sorting_(d.lookupOrDefault<Switch>("skipSorting", true)),
+          multi_level_schwarz_(
+              d.lookupOrDefault<Switch>("multiLevelSchwarz", false)),
+          max_block_size_(d.lookupOrDefault("maxBlockSize", label(1))),
+          precision_(d.lookupOrDefault("precision", word("double")))
+    {
+        word msg = "Generate (Block)-Jacobi preconditioner:\n\tprecision: " +
+                   precision_ + "\n\tmaxBlockSize " +
+                   std::to_string(max_block_size_);
+        MLOG_0(verbose_, msg)
+    }
+
+    virtual std::shared_ptr<gko::LinOp> create()
+    {
+        auto builder = [this](auto b) {
+            return b.with_skip_sorting(skip_sorting_)
+                .with_max_block_size(static_cast<gko::uint32>(max_block_size_))
+                .on(exec_);
+        };
+
+        if (precision_ == "double") {
+            auto b = builder(dbj::build());
+            return dispatch_schwarz(
+                mtx_, exec_,
+                gko::share(
+                    b->generate(gko::as<RepartDistMatrix>(mtx_)->get_local())),
+                d_, verbose_);
+        }
+        if (precision_ == "float") {
+            auto b = builder(fbj::build());
+            return dispatch_schwarz(
+                mtx_, exec_,
+                gko::share(
+                    b->generate(gko::as<RepartDistMatrix>(mtx_)->get_local())),
+                d_, verbose_);
+        }
+
+        return {};
+    }
+};
diff --git a/include/OGL/Preconditioner/LU.hpp b/include/OGL/Preconditioner/LU.hpp
new file mode 100644
index 00000000..6733b828
--- /dev/null
+++ b/include/OGL/Preconditioner/LU.hpp
@@ -0,0 +1,102 @@
+// SPDX-FileCopyrightText: 2025 OGL authors
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+#pragma once
+
+#include "OGL/Preconditioner/Schwarz.hpp"
+
+class LU  // : public PreconditionerWrapper
+{
+    using dic = gko::preconditioner::Jacobi<double, label>;
+    using fic = gko::preconditioner::Jacobi<float, label>;
+
+    std::shared_ptr<gko::Executor> exec_;
+    std::shared_ptr<const gko::LinOp> mtx_;
+    const dictionary &d_;
+    const label verbose_;
+    bool skip_sorting_;
+    bool multi_level_schwarz_;
+    word precision_;
+    word factorization_;
+
+
+public:
+    LU(std::shared_ptr<gko::Executor> exec,
+       std::shared_ptr<const gko::LinOp> mtx, const dictionary &d,
+       label verbose)
+        : exec_(exec),
+          mtx_(mtx),
+          d_(d),
+          verbose_(verbose),
+          skip_sorting_(d.lookupOrDefault<Switch>("skipSorting", true)),
+          multi_level_schwarz_(
+              d.lookupOrDefault<Switch>("multiLevelSchwarz", false)),
+          precision_(d.lookupOrDefault("precision", word("double"))),
+          factorization_(d.lookupOrDefault("factorization", word("ParILU")))
+    {
+        word msg = "Generate " + factorization_ +
+                   "  preconditioner:\n\tprecision: " + precision_;
+        MLOG_0(verbose_, msg)
+    }
+
+
+    std::shared_ptr<gko::LinOp> generate_factorization() const
+    {
+        if (factorization_ == "ILU") {
+            auto factorization_factory =
+                gko::factorization::Ilu<scalar, label>::build()
+                    .with_skip_sorting(skip_sorting_)
+                    .on(exec_);
+
+            auto gkodistmatrix =
+                gko::as<RepartDistMatrix>(mtx_)->get_dist_matrix();
+            return factorization_factory->generate(
+                gko::as<gko::experimental::distributed::Matrix<scalar, label,
+                                                               label>>(
+                    gkodistmatrix)
+                    ->get_local_matrix());
+        }
+        if (factorization_ == "ParILU") {
+            label iterations = d_.lookupOrDefault("iterations", label(5));
+            auto factorization_factory =
+                gko::factorization::ParIlut<scalar, label>::build()
+                    .with_skip_sorting(skip_sorting_)
+                    .with_iterations(iterations)
+                    .on(exec_);
+
+            auto gkodistmatrix =
+                gko::as<RepartDistMatrix>(mtx_)->get_dist_matrix();
+            return factorization_factory->generate(
+                gko::as<gko::experimental::distributed::Matrix<scalar, label,
+                                                               label>>(
+                    gkodistmatrix)
+                    ->get_local_matrix());
+        }
+        if (factorization_ == "ParILUT") {
+            label iterations = d_.lookupOrDefault("iterations", label(5));
+            label fillInLimit = d_.lookupOrDefault("fillInLimit", label(2));
+            auto factorization_factory =
+                gko::factorization::ParIlut<scalar, label>::build()
+                    .with_skip_sorting(skip_sorting_)
+                    .with_fill_in_limit(fillInLimit)
+                    .with_iterations(iterations)
+                    .on(exec_);
+            auto gkodistmatrix =
+                gko::as<RepartDistMatrix>(mtx_)->get_dist_matrix();
+            return factorization_factory->generate(
+                gko::as<gko::experimental::distributed::Matrix<scalar, label,
+                                                               label>>(
+                    gkodistmatrix)
+                    ->get_local_matrix());
+        }
+    }
+
+    virtual std::shared_ptr<gko::LinOp> create()
+    {
+        auto precond_factory = gko::preconditioner::Ilu<>::build().on(exec_);
+        return dispatch_schwarz(
+            mtx_, exec_,
+            gko::share(precond_factory->generate(generate_factorization())), d_,
+            verbose_);
+    }
+};
diff --git a/include/OGL/Preconditioner/Multigrid.hpp b/include/OGL/Preconditioner/Multigrid.hpp
new file mode 100644
index 00000000..35c3979e
--- /dev/null
+++ b/include/OGL/Preconditioner/Multigrid.hpp
@@ -0,0 +1,205 @@
+// SPDX-FileCopyrightText: 2025 OGL authors
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+#pragma once
+
+#include "OGL/Preconditioner/CoarseSolver.hpp"
+#include "OGL/Preconditioner/Schwarz.hpp"
+
+class Multigrid {
+    using ir = gko::solver::Ir<scalar>;
+    using it = gko::stop::Iteration;
+    using ras =
+        gko::experimental::distributed::preconditioner::Schwarz<scalar, label,
+                                                                label>;
+    using sor = gko::preconditioner::Sor<scalar, label>;
+    using dbj = gko::preconditioner::Jacobi<double, label>;
+    using fbj = gko::preconditioner::Jacobi<float, label>;
+    using mg = gko::solver::Multigrid;
+    using pgm = gko::multigrid::Pgm<scalar, label>;
+
+    std::shared_ptr<gko::Executor> exec_;
+
+    std::shared_ptr<const gko::LinOp> mtx_;
+
+    const dictionary &d_;
+
+    const label verbose_;
+
+    bool skip_sorting_;
+
+    bool multi_level_schwarz_;
+
+    word type_;
+    scalar relaxFac_;
+    word cycleName_;
+    label maxLevels_;
+    label minRowsC_;
+    word smoother_;
+    word coarsening_;
+    label maxIterS_;
+
+public:
+    Multigrid(std::shared_ptr<gko::Executor> exec,
+              std::shared_ptr<const gko::LinOp> mtx, const dictionary &d,
+              label verbose)
+        : exec_(exec),
+          mtx_(mtx),
+          d_(d),
+          verbose_(verbose),
+          skip_sorting_(d.lookupOrDefault<Switch>("skipSorting", true)),
+          multi_level_schwarz_(
+              d.lookupOrDefault<Switch>("multiLevelSchwarz", false)),
+          type_(d.lookupOrDefault("type", word("Schwarz"))),
+          relaxFac_(d.lookupOrDefault("relaxationFactor", scalar(0.9))),
+          cycleName_(d.lookupOrDefault("cycle", word("v"))),
+          maxLevels_(d.lookupOrDefault("maxLevels", label(20))),
+          minRowsC_(d.lookupOrDefault("minCoarseRows", label(64000))),
+          smoother_(d.lookupOrDefault("smoother", word("Jacobi"))),
+          coarsening_(d.lookupOrDefault("coarsening", word("GAMG"))),
+          maxIterS_(d.lookupOrDefault("maxIterSmoother", label(1)))
+    {
+        word msg = "\nGenerate Multigrid preconditioner:\n\tmaxLevels:" +
+                   std::to_string(maxLevels_) +
+                   "\n\tminCoarseRows: " + std::to_string(minRowsC_) +
+                   "\n\tSmoother: " + smoother_ +
+                   "\n\trelaxationFactor: " + std::to_string(relaxFac_) +
+                   "\n\tmaxIterSmoother: " + std::to_string(maxIterS_) +
+                   "\n\tcoarsening: " + coarsening_ +
+                   "\n\tcycle: " + cycleName_ + "\n\ttype: " + type_;
+        MLOG_0(verbose_, msg)
+        // FatalErrorInFunction << "Unknown Multigrid type: " << type
+        //                      << "\nValid Choices: Schwarz, Distributed"
+        //                      << abort(FatalError);
+    }
+
+    std::shared_ptr<gko::LinOp> create(
+        /* const objectRegistry &db,
+          const ExecutorHandler &exec_handler*/
+    )
+    {
+        gko::solver::multigrid::cycle cycle;
+        if (cycleName_ == "v") cycle = gko::solver::multigrid::cycle::v;
+        if (cycleName_ == "w") cycle = gko::solver::multigrid::cycle::w;
+        if (cycleName_ == "f") cycle = gko::solver::multigrid::cycle::f;
+
+        std::shared_ptr<gko::LinOpFactory> smootherFac{};
+        if (smoother_ == "Jacobi") {
+            smootherFac =
+                dbj::build().with_max_block_size(1u).with_skip_sorting(true).on(
+                    exec_);
+        }
+        if (smoother_ == "SOR") {
+            smootherFac =
+                sor::build().with_skip_sorting(true).with_symmetric(false).on(
+                    exec_);
+        }
+        if (smoother_ == "SSOR") {
+            smootherFac =
+                sor::build().with_skip_sorting(true).with_symmetric(true).on(
+                    exec_);
+        }
+        if (smootherFac == nullptr) {
+            FatalErrorInFunction << "Unknown smoother: " << smoother_
+                                 << "\nValid Choices: Jacobi, SOR, SSOR"
+                                 << abort(FatalError);
+        }
+
+        // std::shared_ptr<gko::matrix::Csr<scalar, label>> coarseningWeight;
+        // if (coarsening_ == "GAMG") {
+        // }
+
+        // if (coarsening_ == "PGM") {
+        //     // auto repartDistMtx = gko::as<RepartDistMatrix>(mtx_);
+        //     // auto cw = gko::as<gko::matrix::Csr<scalar, label>>(
+        //     //     repartDistMtx->get_local_matrix());
+        //     coarseningWeight = nullptr;
+        //     // std::const_pointer_cast<gko::matrix::Csr<scalar, label>>(cw);
+        // }
+        // // if (coarseningWeight == nullptr) {
+        // //     FatalErrorInFunction << "Unknown coarsening: " << coarsening_
+        // //                          << "\nValid Choices: GAMG, PGM"
+        // //                          << abort(FatalError);
+        // // }
+
+        auto single_it = it::build().with_max_iters(1u);
+        auto smoother_it = gko::stop::Iteration::build().with_max_iters(
+            static_cast<gko::uint32>(maxIterS_));
+
+        auto smoother_gen =
+            type_ == "Distributed"
+                ? gko::share(ir::build()
+                                 .with_solver(ras::build().with_local_solver(
+                                     smootherFac))
+                                 .with_relaxation_factor(relaxFac_)
+                                 .with_criteria(smoother_it)
+                                 .on(exec_))
+                : gko::share(ir::build()
+                                 .with_solver(smootherFac)
+                                 .with_relaxation_factor(relaxFac_)
+                                 .with_criteria(smoother_it)
+                                 .on(exec_));
+
+        if (type_ == "Schwarz") {
+            auto pre_factory = gko::share(
+                mg::build()
+                    .with_max_levels(static_cast<gko::uint32>(maxLevels_))
+                    .with_cycle(cycle)
+                    .with_min_coarse_rows(static_cast<gko::uint32>(minRowsC_))
+                    .with_pre_smoother(smoother_gen)
+                    .with_smoother_iters(maxIterS_)
+                    .with_post_uses_pre(true)
+                    .with_mg_level(
+                        pgm::build()
+                            .with_deterministic(true)
+                            // .with_local_weight_mtx(coarseningWeight)
+                            .on(exec_))
+                    .with_coarsest_solver(
+                        generate_coarse_solver(exec_, d_, verbose_))
+                    .with_criteria(single_it)
+                    .on(exec_));
+            return wrap_schwarz(
+                mtx_, exec_,
+                gko::share(pre_factory->generate(
+                    gko::as<RepartDistMatrix>(mtx_)->get_local())));
+        }
+
+        if (type_ == "Distributed") {
+            // std::shared_ptr<const gko::LinOpFactory> coarsest_solver{};
+            // NOTE does not support distributed currently
+            // if (coarseSolver_ == "Direct") {
+            //     coarsest_solver = gko::share(
+            //         gko::experimental::solver::Direct<scalar, label>::build()
+            //             .on(exec_));
+            // }
+            auto gkodistmatrix =
+                gko::as<RepartDistMatrix>(mtx_)->get_dist_matrix();
+            auto smoother_gen = gko::share(
+                ir::build()
+                    .with_solver(ras::build().with_local_solver(smootherFac))
+                    .with_relaxation_factor(relaxFac_)
+                    .with_criteria(smoother_it)
+                    .on(exec_));
+            auto ret = gko::share(
+                mg::build()
+                    .with_max_levels(maxLevels_)
+                    .with_cycle(cycle)
+                    .with_min_coarse_rows(minRowsC_)
+                    .with_pre_smoother(smoother_gen)
+                    .with_smoother_iters(maxIterS_)
+                    .with_post_uses_pre(true)
+                    .with_mg_level(
+                        pgm::build()
+                            // .with_local_weight_mtx(coarseningWeight)
+                            .with_deterministic(true))
+                    .with_coarsest_solver(
+                        generate_coarse_solver(exec_, d_, verbose_))
+                    .with_criteria(single_it)
+                    .on(exec_)
+                    ->generate(gkodistmatrix));
+            return ret;
+        }
+
+        return {};
+    }
+};
diff --git a/include/OGL/Preconditioner/PreconditionerWrapper.hpp b/include/OGL/Preconditioner/PreconditionerWrapper.hpp
new file mode 100644
index 00000000..939c0e08
--- /dev/null
+++ b/include/OGL/Preconditioner/PreconditionerWrapper.hpp
@@ -0,0 +1,32 @@
+// SPDX-FileCopyrightText: 2024 OGL authors
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+#pragma once
+
+#include <ginkgo/ginkgo.hpp>
+
+namespace Foam {
+
+class PreconditionerWrapper {
+protected:
+    std::shared_ptr<gko::Executor> exec_;
+
+    std::shared_ptr<const gko::LinOp> mtx_;
+
+public:
+    PreconditionerWrapper(std::shared_ptr<gko::Executor> exec,
+                          std::shared_ptr<const gko::LinOp> mtx)
+        : exec_(exec), mtx_(mtx)
+    {}
+
+    /* @brief instantiate a local preconditioner and return LinOp Ptr */
+    virtual std::shared_ptr<gko::LinOp> create();
+
+    /* @brief instantiate a local preconditioner and return LinOp Ptr */
+    // virtual std::shared_ptr<gko::LinOp> create_local();
+
+    // /* @brief instantiate a local preconditioner and return LinOp Ptr */
+    // virtual std::shared_ptr<gko::LinOp> create_multilevel();
+};
+
+}  // namespace Foam
diff --git a/include/OGL/Preconditioner/Schwarz.hpp b/include/OGL/Preconditioner/Schwarz.hpp
new file mode 100644
index 00000000..ac079460
--- /dev/null
+++ b/include/OGL/Preconditioner/Schwarz.hpp
@@ -0,0 +1,133 @@
+// SPDX-FileCopyrightText: 2024 OGL authors
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+#pragma once
+
+#include <ginkgo/ginkgo.hpp>
+#include "OGL/Preconditioner/CoarseSolver.hpp"
+
+namespace Foam {
+template <typename PrecondFactory>
+std::shared_ptr<gko::LinOp> wrap_schwarz(
+    std::shared_ptr<const gko::LinOp> gkomatrix,
+    std::shared_ptr<gko::Executor> device_exec,
+    std::shared_ptr<PrecondFactory> precond)
+{
+    using ras =
+        gko::experimental::distributed::preconditioner::Schwarz<scalar, label,
+                                                                label>;
+    // auto local = gko::as<RepartDistMatrix>(gkomatrix)->get_local();
+    return gko::share(
+        ras::build()
+            .with_generated_local_solver(precond)
+            .on(device_exec)
+            ->generate(gko::as<RepartDistMatrix>(gkomatrix)->get_dist_mtx()));
+}
+
+// template <typename PrecondFactory, typename Factorization>
+// std::shared_ptr<gko::LinOp> wrap_schwarz(
+//     std::shared_ptr<const gko::LinOp> gkomatrix,
+//     std::shared_ptr<gko::Executor> device_exec,
+//     std::shared_ptr<PrecondFactory> precond,
+//     std::shared_ptr<Factorization> factorization)
+// {
+//     using ras =
+//         gko::experimental::distributed::preconditioner::Schwarz<scalar,
+//         label,
+//                                                                 label>;
+//     return gko::share(
+//         ras::build()
+//             .with_generated_local_solver(precond->generate(factorization))
+//             .on(device_exec)
+//             ->generate(gkomatrix));
+// }
+//
+
+template <typename Precond>
+std::shared_ptr<gko::LinOp> wrap_multi_level_schwarz(
+    std::shared_ptr<const gko::LinOp> mtx,
+    std::shared_ptr<gko::Executor> device_exec,
+    std::shared_ptr<Precond> precond, const dictionary &d, label verbose)
+{
+    using pgm = gko::multigrid::Pgm<scalar, label>;
+    using fc = gko::multigrid::FixedCoarsening<scalar, label>;
+    using bj = gko::preconditioner::Jacobi<scalar, label>;
+    using solver = gko::solver::Cg<scalar>;
+    using ras =
+        gko::experimental::distributed::preconditioner::Schwarz<scalar, label,
+                                                                label>;
+
+    auto distmtx = gko::as<RepartDistMatrix>(mtx)->get_dist_matrix();
+    auto local_rows =
+        gko::as<RepartDistMatrix>(mtx)->get_local()->get_size()[0];
+
+    auto coarse_solver = generate_coarse_solver(device_exec, d, verbose);
+    auto coarsening = d.lookupOrDefault<word>("coarsening", word("PGM"));
+    auto coarseWeight = d.lookupOrDefault("coarseWeight", scalar(0.1));
+
+    if (coarsening == "fixed") {
+        auto selCoarseRows = d.lookupOrDefault("selCoarseRows", label(5));
+        auto solveNormC =
+            d.lookupOrDefault("reductionCoarseSolver", label(1e-6));
+        auto maxIterCoarse = d.lookupOrDefault("maxIterCoarse", label(50));
+        word msg =
+            "\nGenerate multi level schwarz:\n\tcoarsening: "
+            "fixed\n\tselCoarseRows: " +
+            std::to_string(selCoarseRows) +
+            "\n\tcoarseWeight: " + std::to_string(coarseWeight);
+        MLOG_0(verbose, msg)
+
+        auto n_rows = local_rows / selCoarseRows;
+        auto sel_rows =
+            gko::array<label>(gko::ReferenceExecutor::create(), n_rows);
+        for (auto i = 0; i < sel_rows.get_size(); i++) {
+            sel_rows.get_data()[i] = selCoarseRows * i;
+        }
+        sel_rows.set_executor(device_exec);
+        auto coarsening_fac = gko::share(
+            fc::build().with_skip_sorting(true).with_coarse_rows(sel_rows).on(
+                device_exec));
+
+        return gko::share(ras::build()
+                              .with_generated_local_solver(precond)
+                              .with_coarse_level(coarsening_fac)
+                              .with_l1_smoother(false)
+                              .with_coarse_solver(coarse_solver)
+                              .with_coarse_weight(coarseWeight)
+                              .on(device_exec)
+                              ->generate(distmtx));
+    }
+    if (coarsening == "PGM") {
+        word msg =
+            "\nGenerate multi level schwarz:\n\tfixedCoarsening: "
+            "PGM\n\tcoarseWeight: " +
+            std::to_string(coarseWeight);
+        MLOG_0(verbose, msg)
+        auto pgm_fac =
+            gko::share(pgm::build().with_skip_sorting(true).on(device_exec));
+        return gko::share(ras::build()
+                              .with_generated_local_solver(precond)
+                              .with_coarse_level(pgm_fac)
+                              .with_l1_smoother(false)
+                              .with_coarse_weight(coarseWeight)
+                              .with_coarse_solver(coarse_solver)
+                              .on(device_exec)
+                              ->generate(distmtx));
+    }
+}
+
+
+template <typename PrecondFactory>
+std::shared_ptr<gko::LinOp> dispatch_schwarz(
+    std::shared_ptr<const gko::LinOp> mtx, std::shared_ptr<gko::Executor> exec,
+    std::shared_ptr<PrecondFactory> precond, const dictionary &d, label verbose)
+{
+    if (d.lookupOrDefault<Switch>("multiLevelSchwarz", false)) {
+        return wrap_multi_level_schwarz(mtx, exec, precond,
+                                        d.subDict("multiLevelConfig"), verbose);
+    } else {
+        return wrap_schwarz(mtx, exec, std::move(precond));
+    }
+}
+
+}  // namespace Foam
diff --git a/include/OGL/StoppingCriterion.hpp b/include/OGL/StoppingCriterion.hpp
index 944acd60..8ca0f18e 100644
--- a/include/OGL/StoppingCriterion.hpp
+++ b/include/OGL/StoppingCriterion.hpp
@@ -220,11 +220,11 @@ class StoppingCriterion {
             }
         }
 
-        word msg = "Creating stopping criterion with minIter " +
-                   std::to_string(minIter) + " frequency " +
-                   std::to_string(frequency) + " prev_solve_iters " +
-                   std::to_string(prev_solve_iters) + " adapt_minIter_  " +
-                   std::to_string(adapt_minIter_) + " prev_rel_cost  ";
+        word msg = "\nCreating stopping criterion\n\tminIter: " +
+                   std::to_string(minIter) +
+                   "\n\tfrequency: " + std::to_string(frequency) +
+                   "\n\tprev_solve_iters: " + std::to_string(prev_solve_iters) +
+                   "\n\tadapt_minIter:  " + std::to_string(adapt_minIter_);
 
         MLOG_0(verbose, msg)
 
diff --git a/include/OGL/lduLduBase.hpp b/include/OGL/lduLduBase.hpp
index 3e997f77..f49cec82 100644
--- a/include/OGL/lduLduBase.hpp
+++ b/include/OGL/lduLduBase.hpp
@@ -308,7 +308,6 @@ class lduLduBase : public OGL_Info,
 
         TIME_WITH_FIELDNAME(verbose_, copy_x_back, this->fieldName(),
                             dist_x.copy_back();)
-
         auto bandwidth_copy_back =
             sizeof(scalar) * psi.size() / delta_t_copy_x_back / 1000.0;
 
@@ -335,7 +334,8 @@ class lduLduBase : public OGL_Info,
             std::to_string(time_per_dof) + std::string(" [ns]") +
             std::string("\n\tTime per iteration and DOF: ") +
             std::to_string(time_per_iter_and_dof) + std::string(" [ns]") +
-            std::string("\n\tRetrieve results bandwidth ");  // +
+            std::string("\n\tRetrieve results bandwidth ");
+
         std::to_string(bandwidth_copy_back) + std::string(" [GByte/s]");
         MLOG_0(verbose_, msg)
 
diff --git a/src/MatrixWrapper/Distributed.cpp b/src/MatrixWrapper/Distributed.cpp
index 6897cf7a..7ce40fc1 100644
--- a/src/MatrixWrapper/Distributed.cpp
+++ b/src/MatrixWrapper/Distributed.cpp
@@ -399,32 +399,14 @@ void update_impl(
                               exec_handler, rank]() {
         // NOTE if symmetric (get it from host_A) we can skip id=0 and wait till
         // id=1 has been copied to use device copy
-        //
         for (auto [id, comm_pattern, data_ptr] : all_to_all_update_data) {
-            // auto start = std::chrono::steady_clock::now();
             auto repartAllToAll =
                 compute_repart_allToall(exec_handler, comm_pattern, rank);
-            // auto end = std::chrono::steady_clock::now();
-            // auto delta_t =
-            // std::chrono::duration_cast<std::chrono::microseconds>(end -
-            // start).count()/1000.0; std::cout << __FILE__ << ":" << "delta t "
-            // << delta_t << " [ms]\n";
 
             auto [length, send_data_ptr] = host_A->get_interface_data(id);
-            // communicate_values(ref_exec, device_exec, repart_comm,
-            // repartAllToAll,
-            //                    send_data_ptr, data_ptr, force_host_buffer);
-            // std::cout << __FILE__ <<
-            //     " Pstream::rank " << Pstream::myProcNo() <<
-            //     " repart_rank() "  << repart_comm->rank() <<
-            //     " send_offsets.back() "  <<
-            //     " id " << id <<
-            //     repartAllToAll.send_offsets.back() << " recv_counts: "  <<
-            //     repartAllToAll.recv_counts << " recv_offsets: "  <<
-            //     repartAllToAll.recv_offsets <<
-            // std::endl;
-
             if (id == 0 && host_A->get_symmetric()) {
+                // if symmetric we can skip id ==0 since it is the same as id==1
+
             } else {
                 MPI_Request request;
                 MPI_Igatherv(send_data_ptr, repartAllToAll.send_offsets.back(),
@@ -441,7 +423,6 @@ void update_impl(
                 auto [zid, zcomm_pattern, zdata_ptr] =
                     all_to_all_update_data[0];
                 // copy recv size data from data_ptr to zdata_ptr
-                //
                 label recv_buffer_size = repartAllToAll.recv_offsets.back();
                 auto l_view = gko::array<scalar>::view(
                     device_exec, recv_buffer_size, data_ptr);
@@ -466,7 +447,6 @@ void update_impl(
         }
     };
 
-
     // perform pairwise communications
     // this update interface data which needs communication
     auto pairwise_communicate = [comm, ref_exec, device_exec,
diff --git a/src/Preconditioner/CoarseSolver.cpp b/src/Preconditioner/CoarseSolver.cpp
new file mode 100644
index 00000000..14dac825
--- /dev/null
+++ b/src/Preconditioner/CoarseSolver.cpp
@@ -0,0 +1,76 @@
+// SPDX-FileCopyrightText: 2024 OGL authors
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#include "OGL/Preconditioner/CoarseSolver.hpp"
+#include "OGL/common.hpp"
+
+std::shared_ptr<const gko::LinOpFactory> generate_coarse_solver(
+    std::shared_ptr<gko::Executor> exec, const dictionary &d, label verbose)
+{
+    using dbj = gko::preconditioner::Jacobi<double, label>;
+    using cg = gko::solver::Cg<scalar>;
+    using bicgstab = gko::solver::Bicgstab<scalar>;
+    using gmres = gko::solver::Gmres<scalar>;
+    using ir = gko::solver::Ir<scalar>;
+    using ras =
+        gko::experimental::distributed::preconditioner::Schwarz<scalar, label,
+                                                                label>;
+
+    auto hasSubDict = d.found("coarseSolverConfig");
+    auto solverDict =
+        (hasSubDict) ? d.subDict("coarseSolverConfig") : dictionary{};
+
+    word solver = solverDict.lookupOrDefault("solver", word("CG"));
+    scalar solveNorm = solverDict.lookupOrDefault("relTol", scalar(1e-6));
+    label maxIter = solverDict.lookupOrDefault("maxIter", label(1));
+
+    word msg = "\nGenerate coarse solver\n\tsolver: " + solver +
+               "\n\trelTolCoarse: " + std::to_string(solveNorm) +
+               "\n\tmaxIter: " + std::to_string(maxIter);
+    MLOG_0(verbose, msg)
+
+    auto solve_it = gko::stop::Iteration::build().with_max_iters(
+        static_cast<gko::uint32>(maxIter));
+    auto solve_red =
+        gko::stop::ResidualNorm<scalar>::build().with_reduction_factor(
+            solveNorm);
+
+    std::shared_ptr<const gko::LinOpFactory> coarsest_solver = {};
+    auto bjfac = ras::build().with_local_solver(
+        dbj::build().with_skip_sorting(true).with_max_block_size(1u).on(exec));
+
+    if (solver == "CG") {
+        coarsest_solver = gko::share(cg::build()
+                                         .with_preconditioner(bjfac)
+                                         .with_criteria(solve_it, solve_red)
+                                         .on(exec));
+    }
+    if (solver == "BiCGStab") {
+        coarsest_solver = gko::share(bicgstab::build()
+                                         .with_preconditioner(bjfac)
+                                         .with_criteria(solve_it, solve_red)
+                                         .on(exec));
+    }
+    if (solver == "GMRES") {
+        coarsest_solver = gko::share(gmres::build()
+                                         .with_preconditioner(bjfac)
+                                         .with_criteria(solve_it, solve_red)
+                                         .on(exec));
+    }
+    if (solver == "Jacobi") {
+        scalar relaxFac{
+            solverDict.lookupOrDefault("relaxationFactor", scalar(0.9))};
+        coarsest_solver = gko::share(ir::build()
+                                         .with_solver(bjfac)
+                                         .with_relaxation_factor(relaxFac)
+                                         .with_criteria(solve_it, solve_red)
+                                         .on(exec));
+    }
+    if (coarsest_solver == nullptr) {
+        FatalErrorInFunction << "Unknown coarse solver: " << solver
+                             << "\nValid choices: CG, BiCGStab, GMRES, Jacobi"
+                             << abort(FatalError);
+    }
+    return coarsest_solver;
+}
diff --git a/src/StoppingCriterion.cpp b/src/StoppingCriterion.cpp
index 07c7c71b..60be435e 100644
--- a/src/StoppingCriterion.cpp
+++ b/src/StoppingCriterion.cpp
@@ -174,6 +174,7 @@ bool StoppingCriterion::OpenFOAMDistStoppingCriterion::check_impl(
                               end_eval - start_eval)
                               .count() /
                           1.0;
+    // std::cout << __FILE__ << "time " << *(parameters_.time) << " [mu s]\n";
     return result;
 }