diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8b81b769..86f8d3b2 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_gko110"
     CACHE STRING "Use specific version of ginkgo")
 
 if(OGL_DEVEL_TOOLS)
diff --git a/include/OGL/MatrixWrapper/Distributed.hpp b/include/OGL/MatrixWrapper/Distributed.hpp
index bda53b87..7060bca5 100644
--- a/include/OGL/MatrixWrapper/Distributed.hpp
+++ b/include/OGL/MatrixWrapper/Distributed.hpp
@@ -13,7 +13,6 @@
 #include "OGL/MatrixWrapper/HostMatrix.hpp"
 #include "OGL/Repartitioner.hpp"
 
-
 /* The RepartDistMatrix class is a wrapper around Ginkgos distributed Matrix
  * class
  *
@@ -38,20 +37,60 @@ class RepartDistMatrix
     using reorder_map_type =
         std::tuple<std::shared_ptr<gko::array<label>>, scalar *,
                    std::shared_ptr<gko::array<label>>>;
-    using all_to_all_data = std::tuple<label, AllToAllPattern, scalar *>;
+    using all_to_all_data = std::tuple<label, AllToAllPattern, scalar *, label>;
 
     struct pairwise_data {
-        label id;          // original interface id on orig rank
-        label send;        // 0 - send, 1, receive, 2 same_rank
-        label comm_rank;   // other side of communication
-        label length;      // length of the interface to communicate
-        label send_id;     //
-        scalar *recv_ptr;  //
+        label id;           // original interface id on orig rank
+        label send;         // 0 - send, 1, receive, 2 same_rank
+        label comm_rank;    // other side of communication
+        label length;       // length of the interface to communicate
+        label send_id;      //
+        scalar *recv_ptr;   // where to put the received data (linop data)
+        label recv_offset;  // offset to begin of linop
     };
 
     using gko::EnableLinOp<RepartDistMatrix>::convert_to;
     using gko::EnableLinOp<RepartDistMatrix>::move_to;
 
+    /* @brief replaces all data ptrs in all_to_all_data with new_ptr
+     */
+    std::vector<reorder_map_type> update_reorder_map_ptr(
+        std::vector<reorder_map_type> in, scalar *new_ptr) const
+    {
+        std::vector<reorder_map_type> out;
+        auto [map, _, pad] = in[0];  //) {
+        out.emplace_back(map, new_ptr, pad);
+        // }
+        return out;
+    }
+
+    /* @brief replaces all data ptrs in all_to_all_data with new_ptr
+     */
+    std::vector<all_to_all_data> update_all_to_all_recv_ptr(
+        std::vector<all_to_all_data> in, scalar *new_ptr) const
+    {
+        std::vector<all_to_all_data> out;
+        for (auto [id, comm_pattern, data_ptr, offset] : in) {
+            out.emplace_back(id, comm_pattern, new_ptr, offset);
+        }
+        return out;
+    }
+
+    /* @brief replaces all data ptrs in pairwise_data with new_ptr
+     */
+    std::vector<pairwise_data> update_pairwise_recv_ptr(
+        std::vector<pairwise_data> in, scalar *new_ptr) const
+    {
+        std::vector<pairwise_data> out;
+        for (auto [linop_id, mode, comm_rank, length, send_id, recv_ptr,
+                   offset] : in) {
+            out.push_back(pairwise_data{linop_id, mode, comm_rank, length,
+                                        send_id, new_ptr, offset});
+        }
+        return out;
+    }
+
+
     std::shared_ptr<const gko::LinOp> get_dist_matrix() const
     {
         return this->dist_mtx_;
@@ -80,6 +119,33 @@ class RepartDistMatrix
         }
     }
 
+    std::shared_ptr<RepartDistMatrix> clone() const
+    {
+        auto new_dist_mtx = gko::share(this->dist_mtx_->clone());
+        auto new_local_ptr = gko::as<gko::matrix::Csr<scalar, label>>(
+            new_dist_mtx->get_local_matrix());
+        const scalar *new_data_ptr = new_local_ptr->get_const_values();
+
+        auto all_to_all_data = update_all_to_all_recv_ptr(
+            this->all_to_all_update_data_, const_cast<scalar *>(new_data_ptr));
+
+        auto new_non_local_ptr = gko::as<gko::matrix::Csr<scalar, label>>(
+            new_dist_mtx->get_non_local_matrix());
+        const scalar *new_non_local_data_ptr =
+            new_non_local_ptr->get_const_values();
+        auto pairwise_update_data = update_pairwise_recv_ptr(
+            this->pairwise_update_data_,
+            const_cast<scalar *>(new_non_local_data_ptr));
+
+        auto new_reorder_map = update_reorder_map_ptr(
+            this->reorder_maps_, const_cast<scalar *>(new_data_ptr));
+
+        return std::make_shared<RepartDistMatrix>(
+            this->get_executor(), this->get_communicator(),
+            this->matrix_format_, new_dist_mtx, this->repartitioner_,
+            this->fuse_, all_to_all_data, pairwise_update_data, new_reorder_map,
+            this->compress_to_global_);
+    }
 
     /**
      * Copy-assigns a CombinationMatrix matrix. Preserves executor, copies
@@ -101,7 +167,6 @@ class RepartDistMatrix
             this->pairwise_update_data_ = other.pairwise_update_data_;
             this->reorder_maps_ = other.reorder_maps_;
             this->compress_to_global_ = other.compress_to_global_;
-            this->linops_ = other.linops_;
         }
         return *this;
     }
@@ -126,7 +191,6 @@ class RepartDistMatrix
                 std::move(other.pairwise_update_data_);
             this->reorder_maps_ = std::move(other.reorder_maps_);
             this->compress_to_global_ = std::move(other.compress_to_global_);
-            this->linops_ = std::move(other.linops_);
         }
         return *this;
     }
@@ -135,6 +199,9 @@ class RepartDistMatrix
     void update(const ExecutorHandler &exec_handler,
                 std::shared_ptr<const HostMatrixWrapper> host_A, label verbose);
 
+    void update(const ExecutorHandler &exec_handler, const scalar *diag_ptr,
+                const scalar *face_ptr, label verbose);
+
     word get_matrix_format() const { return matrix_format_; }
 
     std::shared_ptr<const gko::LinOp> get_dist_mtx() const { return dist_mtx_; }
@@ -147,8 +214,10 @@ class RepartDistMatrix
                      std::vector<all_to_all_data> all_to_all_update_data,
                      std::vector<pairwise_data> pairwise_update_data,
                      std::vector<reorder_map_type> reorder_maps,
-                     std::vector<label> compress_to_global,
-                     std::map<label, scalar *> linops)
+                     std::vector<label> compress_to_global
+                     // ,
+                     // std::map<label, scalar *> linops
+                     )
         : gko::EnableLinOp<RepartDistMatrix>(exec),
           gko::experimental::distributed::DistributedBase(comm),
           fuse_(fuse),
@@ -158,8 +227,7 @@ class RepartDistMatrix
           pairwise_update_data_(pairwise_update_data),
           repartitioner_(repartitioner),
           reorder_maps_(reorder_maps),
-          compress_to_global_(compress_to_global),
-          linops_(linops)
+          compress_to_global_(compress_to_global)
     {
         this->set_size(dist_mtx_->get_size());
     }
@@ -219,18 +287,17 @@ class RepartDistMatrix
 
     std::shared_ptr<dist_mtx> dist_mtx_;
 
+    // id, comm_pattern, data_ptr, offset
     std::vector<all_to_all_data> all_to_all_update_data_;
 
     std::vector<pairwise_data> pairwise_update_data_;
 
     std::shared_ptr<const Repartitioner> repartitioner_;
 
+    // map, data_ptr, pad
     std::vector<reorder_map_type> reorder_maps_;
 
     std::vector<label> compress_to_global_;
-
-    std::map<label, scalar *>
-        linops_;  // map between linop id and shared ptr to linop
 };
 
 std::shared_ptr<const gko::LinOp> get_local(
diff --git a/include/OGL/MatrixWrapper/lduMatrix.hpp b/include/OGL/MatrixWrapper/lduMatrix.hpp
new file mode 100644
index 00000000..e69de29b
diff --git a/include/OGL/Preconditioner.hpp b/include/OGL/Preconditioner.hpp
index 2bd3d5f1..dc61fc71 100644
--- a/include/OGL/Preconditioner.hpp
+++ b/include/OGL/Preconditioner.hpp
@@ -8,17 +8,23 @@
 
 #include "OGL/DevicePersistent/Base.hpp"
 #include "OGL/MatrixWrapper/Distributed.hpp"
+#include "OGL/Preconditioner/Cholesky.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>;
     using fbj = gko::preconditioner::Jacobi<float, label>;
     using dbj = gko::preconditioner::Jacobi<double, label>;
-    using sor = gko::preconditioner::Sor<scalar, label>;
     using ic = gko::preconditioner::Ic<>;
     using ir = gko::solver::Ir<scalar>;
     using it = gko::stop::Iteration;
@@ -52,489 +58,168 @@ 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,
         std::shared_ptr<const gko::LinOp> gkomatrix,
-        std::shared_ptr<gko::Executor> device_exec) const
+        const ExecutorHandler &exec_handler) const
     {
+        auto device_exec = exec_handler.get_device_exec();
         bool skip_sorting = d.lookupOrDefault<Switch>("skipSorting", true);
         bool multi_level_schwarz =
             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);
-        }
-        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);
+            return BlockJacobi(device_exec, gkomatrix, d, verbose_).create();
         }
-        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 == "Multigrid") {
+            return Multigrid(device_exec, gkomatrix, d, verbose_)
+                .create(db_, exec_handler);
         }
         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);
+        if (name == "ILU") {
+            return LU(device_exec, gkomatrix, d, verbose_).create();
         }
+        // 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);
+        // }
+        // 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 == "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 == "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);
+        // }
+        // 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 == "none") {
             return {};
         }
@@ -546,10 +231,12 @@ class Preconditioner {
         return {};
     }
 
+
     std::shared_ptr<gko::LinOp> init_preconditioner(
         std::shared_ptr<const gko::LinOp> gkomatrix,
-        std::shared_ptr<gko::Executor> device_exec) const
+        const ExecutorHandler &exec_handler) const
     {
+        auto device_exec = exec_handler.get_device_exec();
         const word precond_store_name =
             sys_matrix_name_ + "Cached_preconditinoner";
         const fileName path = precond_store_name;
@@ -630,7 +317,7 @@ class Preconditioner {
                 set_next_caching(sys_matrix_name_, db_, caching_period);
 
                 auto generated_precond =
-                    init_preconditioner_impl(name, d, gkomatrix, device_exec);
+                    init_preconditioner_impl(name, d, gkomatrix, exec_handler);
 
                 auto precond_ptr = prev_precond.get_ptr();
                 precond_ptr = generated_precond;
@@ -642,7 +329,7 @@ class Preconditioner {
         cache = get_next_caching(sys_matrix_name_, db_);
 
         auto generated_precond =
-            init_preconditioner_impl(name, d, gkomatrix, device_exec);
+            init_preconditioner_impl(name, d, gkomatrix, exec_handler);
 
         auto po = new DevicePersistentBase<gko::LinOp>(IOobject(path, db_),
                                                        generated_precond);
@@ -653,4 +340,6 @@ class Preconditioner {
         return generated_precond;
     }
 };
+
+
 }  // namespace Foam
diff --git a/include/OGL/Preconditioner/.Multigrid.hpp.swo b/include/OGL/Preconditioner/.Multigrid.hpp.swo
new file mode 100644
index 00000000..04ce7461
Binary files /dev/null and b/include/OGL/Preconditioner/.Multigrid.hpp.swo differ
diff --git a/include/OGL/Preconditioner/Cholesky.hpp b/include/OGL/Preconditioner/Cholesky.hpp
new file mode 100644
index 00000000..8c3571e2
--- /dev/null
+++ b/include/OGL/Preconditioner/Cholesky.hpp
@@ -0,0 +1,72 @@
+// 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_;
+
+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")))
+    {
+        word msg = "Generate incomplete cholesky (IC) preconditioner:\n\tprecision: " +
+                   precision_; 
+        MLOG_0(verbose_, msg)
+    }
+
+    virtual std::shared_ptr<gko::LinOp> create()
+    {
+     auto factorization_factory = gko::share(
+	gko::factorization::Ic<scalar, label>::build()
+	    .with_skip_sorting(skip_sorting_)
+	    .on(exec_));
+     auto gkodistmatrix =
+	 gko::as<RepartDistMatrix>(mtx_)->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::share(
+	 gko::preconditioner::Ic<>::build().on(exec_));
+	
+    return wrap_schwarz(gkodistmatrix, exec_,
+			std::move(precond_factory), factorization);
+
+        // auto wrapper = [this](auto f) {
+        //     if (multi_level_schwarz_) {
+        //         auto distmtx =
+        //             gko::as<RepartDistMatrix>(mtx_)->get_dist_matrix();
+        //         auto local = gko::as<RepartDistMatrix>(mtx_)->get_local();
+        //         auto local_rows = local->get_size()[0];
+        //         return wrap_multi_level_schwarz(distmtx, exec_, f, d_,
+        //                                         local_rows, verbose_);
+        //     } else {
+        //         return wrap_schwarz(mtx_, exec_, std::move(f));
+        //     }
+        // };
+        //
+        // return wrapper(factorization_factory);
+    }
+};
diff --git a/include/OGL/Preconditioner/Jacobi.hpp b/include/OGL/Preconditioner/Jacobi.hpp
new file mode 100644
index 00000000..219c70d6
--- /dev/null
+++ b/include/OGL/Preconditioner/Jacobi.hpp
@@ -0,0 +1,74 @@
+// 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 gko::share(b.with_skip_sorting(skip_sorting_)
+                                  .with_max_block_size(
+                                      static_cast<gko::uint32>(max_block_size_))
+                                  .on(exec_));
+        };
+
+        auto wrapper = [this](auto f) {
+            if (multi_level_schwarz_) {
+                auto distmtx =
+                    gko::as<RepartDistMatrix>(mtx_)->get_dist_matrix();
+                auto local = gko::as<RepartDistMatrix>(mtx_)->get_local();
+                auto local_rows = local->get_size()[0];
+                return wrap_multi_level_schwarz(distmtx, exec_, f, d_,
+                                                local_rows, verbose_);
+            } else {
+                return wrap_schwarz(mtx_, exec_, std::move(f));
+            }
+        };
+
+        if (precision_ == "double") {
+            return wrapper(builder(dbj::build()));
+        }
+        if (precision_ == "float") {
+            return wrapper(builder(fbj::build()));
+        }
+
+        return {};
+    }
+};
diff --git a/include/OGL/Preconditioner/LU.hpp b/include/OGL/Preconditioner/LU.hpp
new file mode 100644
index 00000000..ac36cdac
--- /dev/null
+++ b/include/OGL/Preconditioner/LU.hpp
@@ -0,0 +1,75 @@
+// 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_;
+
+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")))
+    {
+        word msg = "Generate incomplete LU (ILU) preconditioner:\n\tprecision: " +
+                   precision_; 
+        MLOG_0(verbose_, msg)
+    }
+
+    virtual std::shared_ptr<gko::LinOp> create()
+    {
+            label iterations {10};//(d.lookupOrDefault("iterations", label(1)));
+
+	auto precond_factory = gko::share(
+		gko::preconditioner::Ilu<>::build().on(exec_));
+
+        auto wrapper = [this](auto f) {
+            auto factorization_factory =
+                gko::factorization::Ilu<scalar, label>::build()
+                    .with_skip_sorting(skip_sorting_)
+                    //.with_iterations(iterations)
+                    .on(exec_);
+            auto gkodistmatrix =
+                gko::as<RepartDistMatrix>(mtx_)->get_dist_matrix();
+            auto factorization = gko::share(factorization_factory->generate(
+                gko::as<gko::experimental::distributed::Matrix<
+                    scalar, label, label>>(gkodistmatrix)
+                    ->get_local_matrix()));
+            if (multi_level_schwarz_) {
+                // auto distmtx =
+                //     gko::as<RepartDistMatrix>(mtx_)->get_dist_matrix();
+                // auto local = gko::as<RepartDistMatrix>(mtx_)->get_local();
+                // auto local_rows = local->get_size()[0];
+                // return wrap_multi_level_schwarz(distmtx, exec_, f, d_,
+                //                                 local_rows, verbose_);
+            } else {
+                return wrap_schwarz(mtx_, exec_, std::move(f));
+            return wrap_schwarz(gkodistmatrix, exec_,
+                                std::move(f), factorization);
+            }
+        };
+
+	return wrapper(precond_factory);
+
+    }
+};
diff --git a/include/OGL/Preconditioner/Multigrid.hpp b/include/OGL/Preconditioner/Multigrid.hpp
new file mode 100644
index 00000000..b26fabd0
--- /dev/null
+++ b/include/OGL/Preconditioner/Multigrid.hpp
@@ -0,0 +1,309 @@
+// SPDX-FileCopyrightText: 2025 OGL authors
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+#pragma once
+
+#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 cg = gko::solver::Cg<scalar>;
+    using bicgstab = gko::solver::Bicgstab<scalar>;
+    using gmres = gko::solver::Gmres<scalar>;
+    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_;
+    label maxIterCoarseS_;
+    scalar solveNorm_;
+    scalar relaxFac_;
+    word cycleName_;
+    label maxLevels_;
+    label minRowsC_;
+    word smoother_;
+    word coarsening_;
+    word coarseSolver_;
+    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"))),
+          maxIterCoarseS_(d.lookupOrDefault("maxIterCoarse", label(1))),
+          solveNorm_(d.lookupOrDefault("relTolCoarse", scalar(1e-6))),
+          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"))),
+          coarseSolver_(d.lookupOrDefault("coarseSolver", word("Jacobi"))),
+          maxIterS_(d.lookupOrDefault("maxIterSmoother", label(1)))
+    {
+        word msg = "Generate 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\tcoarseSolver: " + coarseSolver_ +
+                   "\n\tmaxIterCoarse: " + std::to_string(maxIterCoarseS_) +
+                   "\n\tinnerSolverNorm: " + std::to_string(solveNorm_) +
+                   "\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::LinOpFactory> bjfac =
+            dbj::build().with_max_block_size(1u).with_skip_sorting(true).on(
+                exec_);
+
+        std::shared_ptr<gko::matrix::Csr<scalar, label>> coarseningWeight;
+        if (coarsening_ == "GAMG") {
+            auto repartDistMtx = gko::as<RepartDistMatrix>(mtx_)->clone();
+            auto &fvmesh = db.template lookupObjectRef<fvMesh>("fvSchemes");
+            auto weights =
+                mag(cmptMultiply(fvmesh.Sf().primitiveField() /
+                                     sqrt(fvmesh.magSf().primitiveField()),
+                                 vector(-1, -1.01, -1.02)))();
+
+            std::vector<scalar> diag(mtx_->get_size()[0], 0.00001);
+
+            repartDistMtx->update(exec_handler, diag.data(), weights.begin(),
+                                  verbose_);
+
+            auto cw = gko::as<gko::matrix::Csr<scalar, label>>(
+                repartDistMtx->get_local_matrix());
+            coarseningWeight =
+                std::const_pointer_cast<gko::matrix::Csr<scalar, label>>(cw);
+
+            // std::cout << " coarsening weight:\n";
+            // for (auto i=0; i<coarseningWeight->get_num_stored_elements();
+            // i++){
+            //   std::cout << i << ":" << coarseningWeight->get_values()[i] <<
+            //   ",";
+            // }
+            // std::cout << "\n";
+        }
+
+        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 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(
+                                     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") {
+            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(exec_));
+            }
+            if (coarseSolver_ == "BiCGStab") {
+                coarsest_solver = gko::share(
+                    bicgstab::build()
+                        .with_preconditioner(bjfac)
+                        .with_criteria(coarse_solve_it, coarse_solve_norm)
+                        .on(exec_));
+            }
+            if (coarseSolver_ == "GMRES") {
+	    std::cout << __FILE__ << __LINE__ << " GMRES 1\n";
+                coarsest_solver = gko::share(
+                    gmres::build()
+                        .with_preconditioner(bjfac)
+                        .with_criteria(coarse_solve_it, coarse_solve_norm)
+                        .on(exec_));
+	    std::cout << __FILE__ << __LINE__ << " GMRES 2 done\n";
+            }
+            if (coarseSolver_ == "Jacobi") {
+                coarsest_solver =
+                    gko::share(ir::build()
+                                   .with_solver(bjfac)
+                                   .with_relaxation_factor(relaxFac_)
+                                   .with_criteria(coarse_solve_it)
+                                   .on(exec_));
+            }
+            if (coarsest_solver == nullptr) {
+                FatalErrorInFunction << "Unknown smoother: " << coarseSolver_
+                                     << "\nValid Choices: CG, Jacobi"
+                                     << abort(FatalError);
+            }
+
+	    std::cout << __FILE__ << __LINE__ << " MG 1\n";
+            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(coarsest_solver)
+                    .with_criteria(single_it)
+                    .on(exec_));
+	    std::cout << __FILE__ << __LINE__ << " MG 2\n";
+            return wrap_schwarz(mtx_, exec_, std::move(pre_factory));
+        }
+
+        if (type_ == "Distributed") {
+            std::shared_ptr<const gko::LinOpFactory> coarsest_solver{};
+            if (coarseSolver_ == "CG") {
+                coarsest_solver = gko::share(
+                    cg::build()
+                        .with_preconditioner(
+                            ras::build().with_local_solver(bjfac))
+                        .with_criteria(coarse_solve_it, coarse_solve_norm)
+                        .on(exec_));
+            }
+            if (coarseSolver_ == "BiCGStab") {
+                coarsest_solver = gko::share(
+                    bicgstab::build()
+                        .with_preconditioner(
+                            ras::build().with_local_solver(bjfac))
+                        .with_criteria(coarse_solve_it, coarse_solve_norm)
+                        .on(exec_));
+            }
+            if (coarseSolver_ == "Jacobi") {
+                coarsest_solver = gko::share(
+                    ir::build()
+                        .with_solver(ras::build().with_local_solver(bjfac))
+                        .with_relaxation_factor(relaxFac_)
+                        .with_criteria(coarse_solve_it)
+                        .on(exec_));
+            }
+            // NOTE does not support distributed currently
+            // if (coarseSolver_ == "Direct") {
+            //     coarsest_solver = gko::share(
+            //         gko::experimental::solver::Direct<scalar, label>::build()
+            //             .on(exec_));
+            // }
+            if (coarsest_solver == nullptr) {
+                FatalErrorInFunction << "Unknown smoother: " << coarseSolver_
+                                     << "\nValid Choices: CG, Jacobi"
+                                     << abort(FatalError);
+            }
+            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(coarsest_solver)
+                    .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..220d5471
--- /dev/null
+++ b/include/OGL/Preconditioner/Schwarz.hpp
@@ -0,0 +1,124 @@
+// SPDX-FileCopyrightText: 2024 OGL authors
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+#pragma once
+
+#include <ginkgo/ginkgo.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>;
+
+    std::cout << __FILE__ << __LINE__ << " MG 3\n";
+    auto local = gko::as<RepartDistMatrix>(gkomatrix)->get_local();
+    std::cout << __FILE__ << __LINE__ << " MG 4\n";
+    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::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 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, 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 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));
+    }
+}
+
+}  // namespace Foam
diff --git a/include/OGL/Solver/Multigrid.hpp b/include/OGL/Solver/Multigrid.hpp
index b76674ad..210867e3 100644
--- a/include/OGL/Solver/Multigrid.hpp
+++ b/include/OGL/Solver/Multigrid.hpp
@@ -24,7 +24,6 @@ class GKOMultigridFactory {
     using ras =
         gko::experimental::distributed::preconditioner::Schwarz<scalar, label,
                                                                 label>;
-
     using dist_vec = gko::experimental::distributed::Vector<scalar>;
     using dist_mtx =
         gko::experimental::distributed::Matrix<scalar, label, label>;
diff --git a/include/OGL/lduLduBase.hpp b/include/OGL/lduLduBase.hpp
index e584adbe..a5c750be 100644
--- a/include/OGL/lduLduBase.hpp
+++ b/include/OGL/lduLduBase.hpp
@@ -268,9 +268,9 @@ class lduLduBase : public OGL_Info,
                 dist_b.get_vector()->scale(dense_scaling.get());)
         }
 
-        TIME_WITH_FIELDNAME(verbose_, init_precond, this->fieldName(),
-                            auto precond = this->init_preconditioner(
-                                dist_A_v, exec_handler_.get_device_exec());)
+        TIME_WITH_FIELDNAME(
+            verbose_, init_precond, this->fieldName(),
+            auto precond = this->init_preconditioner(dist_A_v, exec_handler_);)
 
         bool active = repartitioner->get_repart_size() != 0;
         bool export_system(
diff --git a/src/MatrixWrapper/Distributed.cpp b/src/MatrixWrapper/Distributed.cpp
index 72eaae43..62b57cc6 100644
--- a/src/MatrixWrapper/Distributed.cpp
+++ b/src/MatrixWrapper/Distributed.cpp
@@ -100,7 +100,7 @@ void generate_pairwise_update_data(
             if (fuse) {
                 get_id = (start_offset == 0) ? -1 : 0;
             }
-            recv_ptr = linops[get_id] + offset;
+            recv_ptr = linops[get_id];  // + offset;
             offset += (fuse) ? length : 0;
         }
         // std::cout << __FILE__ << __LINE__ << " set update_data rank " << rank
@@ -109,7 +109,7 @@ void generate_pairwise_update_data(
         //           << send_data_ptr << " recv_data_ptr " << recv_data_ptr
         //           << "\n";
         update_data.push_back(RepartDistMatrix::pairwise_data{
-            linop_id, mode, comm_rank, length, send_id, recv_ptr});
+            linop_id, mode, comm_rank, length, send_id, recv_ptr, offset});
     }
 }
 
@@ -132,9 +132,9 @@ void generate_alltoall_update_data(
 
         // NOTE Probably dont need to store linops[linop-idx] because we can
         // just reuse linops
-        update_data.emplace_back(
-            in->get_id()[i], comm_pattern,
-            (owner) ? linops[linop_idx] + linop_offset_store : nullptr);
+        update_data.emplace_back(in->get_id()[i], comm_pattern,
+                                 (owner) ? linops[linop_idx] : nullptr,
+                                 linop_offset_store);
 
         linop_offset_store += (fuse) ? recv_size : 0;
     }
@@ -260,7 +260,6 @@ void RepartDistMatrix::write(const ExecutorHandler &exec_handler,
     auto non_local = gko::share(
         gko::matrix::Coo<scalar, label>::create(exec_handler.get_ref_exec()));
 
-
     if (fuse_) {
         gko::as<LocalMatrixType>(dist_mtx_->get_local_matrix())
             ->convert_to(local.get());
@@ -375,13 +374,61 @@ void reorder_interface_impl(const ExecutorHandler &exec_handler,
     }
 }
 
+void advanced_update_impl(
+    const ExecutorHandler &exec_handler, const scalar *diag_weights,
+    const scalar *face_weights,
+    std::vector<RepartDistMatrix::all_to_all_data> &all_to_all_update_data,
+    std::vector<RepartDistMatrix::pairwise_data> &pairwise_update_data,
+    std::vector<RepartDistMatrix::reorder_map_type> &reorder_maps, bool fuse,
+    label verbose)
+{
+    auto comm = exec_handler.get_host_comm();
+    // auto repart_comm = exec_handler.get_repart_comm();
+    auto ref_exec = exec_handler.get_ref_exec();
+    auto rank = exec_handler.get_host_rank();
+    auto device_exec = exec_handler.get_device_exec();
+    bool force_host_buffer = exec_handler.get_gko_force_host_buffer();
+    word fieldname = "foo";
+
+    // perform all-to-all updates first
+    auto all_to_all_update = [comm, ref_exec, device_exec,
+                              all_to_all_update_data, force_host_buffer,
+                              face_weights, diag_weights]() {
+        for (auto [id, comm_pattern, data_ptr, offset] :
+             all_to_all_update_data) {
+            // auto [length, send_data_ptr] = host_A->get_interface_data(id);
+            if (id == 2) {
+                communicate_values(ref_exec, device_exec, comm, comm_pattern,
+                                   diag_weights, data_ptr + offset,
+                                   force_host_buffer);
+            } else {
+                communicate_values(ref_exec, device_exec, comm, comm_pattern,
+                                   face_weights, data_ptr + offset,
+                                   force_host_buffer);
+            }
+        }
+    };
+
+    TIME_WITH_FIELDNAME(verbose, perform_all_to_all_update, fieldname,
+                        all_to_all_update(););
+
+    auto reorder_data = [reorder_maps, exec_handler]() {
+        for (auto [reorder_map, data_ptr, pad] : reorder_maps) {
+            reorder_interface_impl(exec_handler, reorder_map, pad, data_ptr);
+        }
+    };
+
+    TIME_WITH_FIELDNAME(verbose, reorder_matrix_data, fieldname,
+                        reorder_data();)
+}
+
 void update_impl(
     const ExecutorHandler &exec_handler,
     std::shared_ptr<const HostMatrixWrapper> host_A,
     std::vector<RepartDistMatrix::all_to_all_data> &all_to_all_update_data,
     std::vector<RepartDistMatrix::pairwise_data> &pairwise_update_data,
     std::vector<RepartDistMatrix::reorder_map_type> &reorder_maps, bool fuse,
-    std::map<label, scalar *> linops, label verbose)
+    label verbose)
 {
     auto comm = exec_handler.get_host_comm();
     // auto repart_comm = exec_handler.get_repart_comm();
@@ -395,10 +442,12 @@ void update_impl(
     auto all_to_all_update = [comm, ref_exec, device_exec,
                               all_to_all_update_data, host_A,
                               force_host_buffer]() {
-        for (auto [id, comm_pattern, data_ptr] : all_to_all_update_data) {
+        for (auto [id, comm_pattern, data_ptr, offset] :
+             all_to_all_update_data) {
             auto [length, send_data_ptr] = host_A->get_interface_data(id);
             communicate_values(ref_exec, device_exec, comm, comm_pattern,
-                               send_data_ptr, data_ptr, force_host_buffer);
+                               send_data_ptr, data_ptr + offset,
+                               force_host_buffer);
         }
     };
 
@@ -419,9 +468,8 @@ void update_impl(
     // this update interface data which needs communication
     auto pairwise_communicate = [comm, ref_exec, device_exec,
                                  pairwise_update_data, host_A, rank,
-                                 &get_send_ptr, &linops, fuse,
-                                 force_host_buffer]() {
-        for (auto [id, mode, comm_rank, length, send_id, recv_ptr] :
+                                 &get_send_ptr, fuse, force_host_buffer]() {
+        for (auto [id, mode, comm_rank, length, send_id, recv_ptr, offset] :
              pairwise_update_data) {
             std::vector<scalar> send_buffer;
             send_buffer.reserve(length);
@@ -437,11 +485,12 @@ void update_impl(
                 if (force_host_buffer) {
                     auto tmp = gko::array<scalar>(ref_exec, length);
                     comm->recv(ref_exec, tmp.get_data(), length, comm_rank, 0);
-                    auto recv_view =
-                        gko::array<scalar>::view(device_exec, length, recv_ptr);
+                    auto recv_view = gko::array<scalar>::view(
+                        device_exec, length, recv_ptr + offset);
                     recv_view = tmp;
                 } else {
-                    comm->recv(device_exec, recv_ptr, length, comm_rank, 0);
+                    comm->recv(device_exec, recv_ptr + offset, length,
+                               comm_rank, 0);
                 }
             }
             if (mode == 2) {
@@ -478,12 +527,21 @@ void RepartDistMatrix::update(const ExecutorHandler &exec_handler,
                               std::shared_ptr<const HostMatrixWrapper> host_A,
                               label verbose)
 {
-    SIMPLE_TIME(verbose, perform_matrix_update,
-                update_impl(exec_handler, host_A, all_to_all_update_data_,
-                            pairwise_update_data_, reorder_maps_, fuse_,
-                            linops_, verbose););
+    SIMPLE_TIME(
+        verbose, perform_matrix_update,
+        update_impl(exec_handler, host_A, all_to_all_update_data_,
+                    pairwise_update_data_, reorder_maps_, fuse_, verbose););
 }
 
+void RepartDistMatrix::update(const ExecutorHandler &exec_handler,
+                              const scalar *diag_ptr, const scalar *face_ptr,
+                              label verbose)
+{
+    SIMPLE_TIME(verbose, perform_matrix_update,
+                advanced_update_impl(
+                    exec_handler, diag_ptr, face_ptr, all_to_all_update_data_,
+                    pairwise_update_data_, reorder_maps_, fuse_, verbose););
+}
 
 template <typename LocalMatrixType,
           typename NonLocalMatrixType = gko::matrix::Coo<scalar, label>>
@@ -496,7 +554,6 @@ std::shared_ptr<RepartDistMatrix> create_impl(
     using dist_mtx =
         gko::experimental::distributed::Matrix<scalar, label, label>;
     using dim_type = gko::dim<2>::dimension_type;
-    label rank = exec_handler.get_host_rank();
     auto exec = exec_handler.get_ref_exec();
     auto host_comm = *exec_handler.get_host_comm().get();
     exec_handler.init_device_comm();
@@ -526,14 +583,13 @@ std::shared_ptr<RepartDistMatrix> create_impl(
     auto device_exec = exec_handler.get_device_exec();
     auto ranks_per_owner = repartitioner->get_ranks_per_gpu();
     bool reparts = ranks_per_owner > 1;
-    // TODO rename this to something useful
-    // this is the part of the update data structure
-    std::map<label, scalar *> linops;
+    std::map<label, scalar *> host_to_linop_map;
     auto [loc_rows, loc_cols, loc_map, loc_ids] =
         (fuse) ? repart_loc_sparsity->get_fused_vecs(false)
                : repart_loc_sparsity->get_vecs(false, reparts);
     auto local_linops = generate_inner_linops<LocalMatrixType>(
-        exec_handler, repart_dim, loc_rows, loc_cols, loc_ids, linops, fuse);
+        exec_handler, repart_dim, loc_rows, loc_cols, loc_ids,
+        host_to_linop_map, fuse);
 
     auto partition = gko::share(
         gko::experimental::distributed::build_partition_from_local_size<label,
@@ -545,6 +601,7 @@ std::shared_ptr<RepartDistMatrix> create_impl(
     auto recv_connections =
         repart_comm_pattern->compute_recv_connections(exec_handler, partition);
 
+    // label rank = exec_handler.get_host_rank();
     // std::cout << __FILE__ << " rank " << rank
     //     << " repart size "  << repartitioner->get_repart_size()
     //     << " recv_connections size "  << recv_connections.get_size()
@@ -576,7 +633,7 @@ std::shared_ptr<RepartDistMatrix> create_impl(
 
     auto non_local_linops = generate_inner_linops<NonLocalMatrixType>(
         exec_handler, repart_non_local_dim, non_loc_rows, non_loc_cols,
-        non_loc_ids, linops, fuse);
+        non_loc_ids, host_to_linop_map, fuse);
 
     auto compress_to_global =
         repart_non_loc_sparsity->compute_to_global_map(fuse);
@@ -591,8 +648,8 @@ std::shared_ptr<RepartDistMatrix> create_impl(
     std::vector<RepartDistMatrix::all_to_all_data> all_to_all_update_data;
     SIMPLE_TIME(verbose, generate_all_to_all_update_data,
                 generate_alltoall_update_data<LocalMatrixType>(
-                    exec_handler, local_sparsity, linops, fuse, owner,
-                    ranks_per_owner, all_to_all_update_data););
+                    exec_handler, local_sparsity, host_to_linop_map, fuse,
+                    owner, ranks_per_owner, all_to_all_update_data););
 
     size_t start_local_offset = 0;
     if (fuse && owner) {
@@ -605,17 +662,18 @@ std::shared_ptr<RepartDistMatrix> create_impl(
     }
 
     std::vector<RepartDistMatrix::pairwise_data> pairwise_update_data;
+    SIMPLE_TIME(verbose, generate_local_pairwise_data,
+                generate_pairwise_update_data<LocalMatrixType>(
+                    exec_handler, host_A,
+                    (!owner) ? local_sparsity : repart_loc_sparsity, fuse,
+                    repartitioner, host_to_linop_map, start_local_offset,
+                    pairwise_update_data););
     SIMPLE_TIME(
-        verbose, generate_local_pairwise_data,
-        generate_pairwise_update_data<LocalMatrixType>(
+        verbose, generate_non_local_pairwise_data,
+        generate_pairwise_update_data<NonLocalMatrixType>(
             exec_handler, host_A,
-            (!owner) ? local_sparsity : repart_loc_sparsity, fuse,
-            repartitioner, linops, start_local_offset, pairwise_update_data););
-    SIMPLE_TIME(verbose, generate_non_local_pairwise_data,
-                generate_pairwise_update_data<NonLocalMatrixType>(
-                    exec_handler, host_A,
-                    (!owner) ? non_local_sparsity : repart_non_loc_sparsity,
-                    fuse, repartitioner, linops, 0, pairwise_update_data););
+            (!owner) ? non_local_sparsity : repart_non_loc_sparsity, fuse,
+            repartitioner, host_to_linop_map, 0, pairwise_update_data););
 
     // sort the pairwise update_data by interface id so that it is
     // consistent across ranks
@@ -671,15 +729,15 @@ std::shared_ptr<RepartDistMatrix> create_impl(
                     exec_handler, non_local_linops, non_loc_map, non_local_pad,
                     reorder_maps););
 
-    SIMPLE_TIME(verbose, perform_matrix_update,
-                update_impl(exec_handler, host_A, all_to_all_update_data,
-                            pairwise_update_data, reorder_maps, fuse, linops,
-                            verbose););
+    SIMPLE_TIME(
+        verbose, perform_matrix_update,
+        update_impl(exec_handler, host_A, all_to_all_update_data,
+                    pairwise_update_data, reorder_maps, fuse, verbose););
 
     return std::make_shared<RepartDistMatrix>(
         device_exec, host_comm, matrix_format, dist_A, repartitioner, fuse,
         all_to_all_update_data, pairwise_update_data, reorder_maps,
-        compress_to_global, linops);
+        compress_to_global);
 }
 
 void write_distributed(const ExecutorHandler &exec_handler, word field_name,
diff --git a/src/StoppingCriterion.cpp b/src/StoppingCriterion.cpp
index 38f9c13d..9dec159c 100644
--- a/src/StoppingCriterion.cpp
+++ b/src/StoppingCriterion.cpp
@@ -150,12 +150,17 @@ bool StoppingCriterion::OpenFOAMDistStoppingCriterion::check_impl(
     if (residual_norm < parameters_.openfoam_absolute_tolerance) {
         result = true;
     }
-    // check if relative tolerance is hit
+    // check if relative convergence tolerance is hit
     if (parameters_.openfoam_relative_tolerance > 0 &&
         residual_norm <
             parameters_.openfoam_relative_tolerance * init_residual) {
         result = true;
     }
+    // check if relative divergence criteria is hit
+    if (parameters_.openfoam_relative_tolerance > 0 &&
+        residual_norm > 20) {
+        result = true;
+    }
 
     if (result) {
         this->set_all_statuses(stoppingId, setFinalized, stop_status);