From 677dfacdcddacb36bee86441eaee4031a62efff9 Mon Sep 17 00:00:00 2001
From: Gregor Olenik <go@hpsim.de>
Date: Thu, 7 Aug 2025 15:31:25 +0200
Subject: [PATCH 01/12] WIP gamg coarsening

---
 include/OGL/Preconditioner.hpp | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/include/OGL/Preconditioner.hpp b/include/OGL/Preconditioner.hpp
index 2bd3d5f1..a5d11c2b 100644
--- a/include/OGL/Preconditioner.hpp
+++ b/include/OGL/Preconditioner.hpp
@@ -372,6 +372,7 @@ class Preconditioner {
             auto maxLevels = d.lookupOrDefault("maxLevels", label(20));
             auto minRowsC = d.lookupOrDefault("minCoarseRows", label(64000));
             auto smoother = d.lookupOrDefault("smoother", word("Jacobi"));
+            auto coarsening = d.lookupOrDefault("coarsening", word("GAMG"));
             auto coarseSolver =
                 d.lookupOrDefault("coarseSolver", word("Jacobi"));
             auto maxIterS = d.lookupOrDefault("maxIterSmoother", label(1));
@@ -387,13 +388,13 @@ class Preconditioner {
                        "\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 + " type: " + type;
             MLOG_0(verbose_, msg)
 
-
             std::shared_ptr<gko::LinOpFactory> bjfac{};
 
             if (smoother == "Jacobi") {
@@ -420,6 +421,11 @@ class Preconditioner {
                                      << abort(FatalError);
             }
 
+	    std::shared_ptr<gko::LinOp> coarseningWeight;
+	    if (coarsening=="GAMG"){
+		    coarseningWeight=gkomatrix;
+	    }
+
             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));

From ba7f67a746d769dc8cf7ef6e614c298b02ad5c2d Mon Sep 17 00:00:00 2001
From: Gregor Olenik <go@hpsim.de>
Date: Fri, 15 Aug 2025 16:57:16 +0200
Subject: [PATCH 02/12] wip gamg weights

---
 CMakeLists.txt                 |  2 +-
 include/OGL/Preconditioner.hpp | 29 ++++++++++++++++++++++++++---
 2 files changed, 27 insertions(+), 4 deletions(-)

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/Preconditioner.hpp b/include/OGL/Preconditioner.hpp
index a5d11c2b..6640b167 100644
--- a/include/OGL/Preconditioner.hpp
+++ b/include/OGL/Preconditioner.hpp
@@ -421,9 +421,28 @@ class Preconditioner {
                                      << abort(FatalError);
             }
 
-	    std::shared_ptr<gko::LinOp> coarseningWeight;
+	    std::shared_ptr<gko::matrix::Csr<double, int>> coarseningWeight;
 	    if (coarsening=="GAMG"){
-		    coarseningWeight=gkomatrix;
+
+		auto& fvmesh = db_.template lookupObjectRef<fvMesh>("fvSchemes");
+		// weights[facei] -> column
+		auto weights = 	mag
+		(
+		cmptMultiply
+		(
+		fvmesh.Sf().primitiveField()
+		/sqrt(fvmesh.magSf().primitiveField()),
+		vector(1, 1.01, 1.02)
+		)
+		);
+
+		// here ldu adressing is needed
+		// 1. create row, cols, vals vector
+		// 2. fill with diagonals i=j=row=col
+	        // 3. fill with off-diagonals values from weights and row, cols from ldu
+
+
+		coarseningWeight=nullptr;
 	    }
 
             auto single_it = it::build().with_max_iters(1u);
@@ -481,7 +500,10 @@ class Preconditioner {
                         .with_pre_smoother(smoother_gen)
                         .with_post_uses_pre(true)
                         .with_mg_level(
-                            pgm::build().with_deterministic(false).on(
+                            pgm::build()
+			    .with_deterministic(false)
+			    .with_local_weight_mtx(coarseningWeight)
+			    .on(
                                 device_exec))
                         .with_coarsest_solver(coarsest_solver)
                         .with_criteria(single_it)
@@ -525,6 +547,7 @@ class Preconditioner {
                     gko::solver::Multigrid::build()
                         .with_max_levels(maxLevels)
                         .with_mg_level(gko::multigrid::Pgm<scalar>::build()
+			    .with_local_weight_mtx(coarseningWeight)
                                            .with_deterministic(true))
                         .with_min_coarse_rows(minRowsC)
                         .with_coarsest_solver(coarsest_solver)

From 42b16faca64683f751c7f48225bd9123169a6a3b Mon Sep 17 00:00:00 2001
From: Gregor Olenik <gregor.olenik@web.de>
Date: Fri, 15 Aug 2025 16:58:26 +0200
Subject: [PATCH 03/12] format

---
 include/OGL/Preconditioner.hpp | 56 ++++++++++++++++------------------
 1 file changed, 26 insertions(+), 30 deletions(-)

diff --git a/include/OGL/Preconditioner.hpp b/include/OGL/Preconditioner.hpp
index 6640b167..91022068 100644
--- a/include/OGL/Preconditioner.hpp
+++ b/include/OGL/Preconditioner.hpp
@@ -421,29 +421,25 @@ class Preconditioner {
                                      << abort(FatalError);
             }
 
-	    std::shared_ptr<gko::matrix::Csr<double, int>> coarseningWeight;
-	    if (coarsening=="GAMG"){
-
-		auto& fvmesh = db_.template lookupObjectRef<fvMesh>("fvSchemes");
-		// weights[facei] -> column
-		auto weights = 	mag
-		(
-		cmptMultiply
-		(
-		fvmesh.Sf().primitiveField()
-		/sqrt(fvmesh.magSf().primitiveField()),
-		vector(1, 1.01, 1.02)
-		)
-		);
-
-		// here ldu adressing is needed
-		// 1. create row, cols, vals vector
-		// 2. fill with diagonals i=j=row=col
-	        // 3. fill with off-diagonals values from weights and row, cols from ldu
-
-
-		coarseningWeight=nullptr;
-	    }
+            std::shared_ptr<gko::matrix::Csr<double, int>> coarseningWeight;
+            if (coarsening == "GAMG") {
+                auto &fvmesh =
+                    db_.template lookupObjectRef<fvMesh>("fvSchemes");
+                // weights[facei] -> column
+                auto weights =
+                    mag(cmptMultiply(fvmesh.Sf().primitiveField() /
+                                         sqrt(fvmesh.magSf().primitiveField()),
+                                     vector(1, 1.01, 1.02)));
+
+                // here ldu adressing is needed
+                // 1. create row, cols, vals vector
+                // 2. fill with diagonals i=j=row=col
+                // 3. fill with off-diagonals values from weights and row, cols
+                // from ldu
+
+
+                coarseningWeight = nullptr;
+            }
 
             auto single_it = it::build().with_max_iters(1u);
             auto coarse_solve_it = gko::stop::Iteration::build().with_max_iters(
@@ -501,10 +497,9 @@ class Preconditioner {
                         .with_post_uses_pre(true)
                         .with_mg_level(
                             pgm::build()
-			    .with_deterministic(false)
-			    .with_local_weight_mtx(coarseningWeight)
-			    .on(
-                                device_exec))
+                                .with_deterministic(false)
+                                .with_local_weight_mtx(coarseningWeight)
+                                .on(device_exec))
                         .with_coarsest_solver(coarsest_solver)
                         .with_criteria(single_it)
                         .on(device_exec);
@@ -546,9 +541,10 @@ class Preconditioner {
                 auto ret = gko::share(
                     gko::solver::Multigrid::build()
                         .with_max_levels(maxLevels)
-                        .with_mg_level(gko::multigrid::Pgm<scalar>::build()
-			    .with_local_weight_mtx(coarseningWeight)
-                                           .with_deterministic(true))
+                        .with_mg_level(
+                            gko::multigrid::Pgm<scalar>::build()
+                                .with_local_weight_mtx(coarseningWeight)
+                                .with_deterministic(true))
                         .with_min_coarse_rows(minRowsC)
                         .with_coarsest_solver(coarsest_solver)
                         .with_criteria(it::build().with_max_iters(2u))

From 868b5f6258c545ff5fe717e64b9f9755e7cfefac Mon Sep 17 00:00:00 2001
From: Gregor Olenik <gregor.olenik@web.de>
Date: Mon, 25 Aug 2025 11:34:35 +0200
Subject: [PATCH 04/12] wip pass exec handler

---
 include/OGL/Preconditioner.hpp | 35 ++++++++++++++++++++++++++++++----
 include/OGL/lduLduBase.hpp     |  6 +++---
 2 files changed, 34 insertions(+), 7 deletions(-)

diff --git a/include/OGL/Preconditioner.hpp b/include/OGL/Preconditioner.hpp
index 91022068..ce86f354 100644
--- a/include/OGL/Preconditioner.hpp
+++ b/include/OGL/Preconditioner.hpp
@@ -160,8 +160,9 @@ class Preconditioner {
     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);
@@ -431,6 +432,23 @@ class Preconditioner {
                                          sqrt(fvmesh.magSf().primitiveField()),
                                      vector(1, 1.01, 1.02)));
 
+                //   const std::shared_ptr<HostMatrixWrapper>
+                //   weight_matrix_wrapper
+                //       {std::make_shared<HostMatrixWrapper>(
+                // exec_handler, db_, matrix.diag().size(),
+                // matrix.upper().size(), matrix.symmetric(),
+                // matrix.diag().begin(), matrix.upper().begin(),
+                // matrix.lower().begin(), matrix.lduAddr(), interfaceBouCoeffs,
+                // interfaceIntCoeffs, interfaces, solverControls, fieldName,
+                // verbose_)
+                // };
+
+
+                // std::shared_ptr<RepartDistMatrix> dist_weight_mtx =
+                //     create_distributed(&exec_handler, repartitioner,
+                //                        weigh_matrix_wrapper, "Csr" true,
+                //                        verbose_);
+
                 // here ldu adressing is needed
                 // 1. create row, cols, vals vector
                 // 2. fill with diagonals i=j=row=col
@@ -571,10 +589,19 @@ class Preconditioner {
         return {};
     }
 
+    class PreconditionerWrapper {
+
+        virtual std::shared_ptr<gko::LinOp> create() {
+
+        }
+
+    };
+
     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;
@@ -655,7 +682,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;
@@ -667,7 +694,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);
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(

From f5ee3a14e83a2c03ba8435398d47c93283bc59f8 Mon Sep 17 00:00:00 2001
From: Gregor Olenik <gregor.olenik@web.de>
Date: Tue, 26 Aug 2025 08:46:00 +0200
Subject: [PATCH 05/12] wip move preconditioner to separate implementation
 files

---
 include/OGL/Preconditioner.hpp                | 41 +++-------------
 include/OGL/Preconditioner/Jacobi.hpp         | 47 +++++++++++++++++++
 .../Preconditioner/PreconditionerWrapper.hpp  | 15 ++++++
 include/OGL/Preconditioner/Schwarz.hpp        | 38 +++++++++++++++
 4 files changed, 106 insertions(+), 35 deletions(-)
 create mode 100644 include/OGL/Preconditioner/Jacobi.hpp
 create mode 100644 include/OGL/Preconditioner/PreconditionerWrapper.hpp
 create mode 100644 include/OGL/Preconditioner/Schwarz.hpp

diff --git a/include/OGL/Preconditioner.hpp b/include/OGL/Preconditioner.hpp
index ce86f354..cb4d60b9 100644
--- a/include/OGL/Preconditioner.hpp
+++ b/include/OGL/Preconditioner.hpp
@@ -8,11 +8,15 @@
 
 #include "OGL/DevicePersistent/Base.hpp"
 #include "OGL/MatrixWrapper/Distributed.hpp"
+#include "OGL/Preconditioner/Jacobi.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>;
@@ -128,34 +132,6 @@ class Preconditioner {
         }
     }
 
-    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,
@@ -589,13 +565,6 @@ class Preconditioner {
         return {};
     }
 
-    class PreconditionerWrapper {
-
-        virtual std::shared_ptr<gko::LinOp> create() {
-
-        }
-
-    };
 
     std::shared_ptr<gko::LinOp> init_preconditioner(
         std::shared_ptr<const gko::LinOp> gkomatrix,
@@ -705,4 +674,6 @@ class Preconditioner {
         return generated_precond;
     }
 };
+
+
 }  // namespace Foam
diff --git a/include/OGL/Preconditioner/Jacobi.hpp b/include/OGL/Preconditioner/Jacobi.hpp
new file mode 100644
index 00000000..4e1dbd97
--- /dev/null
+++ b/include/OGL/Preconditioner/Jacobi.hpp
@@ -0,0 +1,47 @@
+// SPDX-FileCopyrightText: 2025 OGL authors
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+#pragma once
+
+#include "PreconditionerWrapper.hpp"
+#include "OGL/Preconditioner/Schwarz.hpp"
+
+class BJ : public PreconditionerWrapper {
+    using dbj = gko::preconditioner::Jacobi<double, label>;
+    using fbj = gko::preconditioner::Jacobi<float, label>;
+
+    const label verbose_;
+    bool skip_sorting_;
+    label max_block_size_;
+    word precision_;
+
+
+public:
+    BJ(const dictionary &d, label verbose)
+        : verbose_(verbose),
+          skip_sorting_(d.lookupOrDefault<Switch>("skipSorting", true)),
+          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_local()
+    {
+        // 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));
+    }
+
+
+};
diff --git a/include/OGL/Preconditioner/PreconditionerWrapper.hpp b/include/OGL/Preconditioner/PreconditionerWrapper.hpp
new file mode 100644
index 00000000..0e92f5ef
--- /dev/null
+++ b/include/OGL/Preconditioner/PreconditionerWrapper.hpp
@@ -0,0 +1,15 @@
+// SPDX-FileCopyrightText: 2024 OGL authors
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+#pragma once
+
+#include <ginkgo/ginkgo.hpp>
+
+class PreconditionerWrapper {
+
+    /* @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();
+};
diff --git a/include/OGL/Preconditioner/Schwarz.hpp b/include/OGL/Preconditioner/Schwarz.hpp
new file mode 100644
index 00000000..fc552c91
--- /dev/null
+++ b/include/OGL/Preconditioner/Schwarz.hpp
@@ -0,0 +1,38 @@
+// 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::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));
+    }
+
+    }

From b8c01acf02d7604f1b8c05aa2a6e92c551a0a0e2 Mon Sep 17 00:00:00 2001
From: Gregor Olenik <gregor.olenik@web.de>
Date: Tue, 26 Aug 2025 10:07:28 +0200
Subject: [PATCH 06/12] wip move preconditioner to separate implementation
 files

---
 include/OGL/Preconditioner.hpp                | 122 +---------------
 include/OGL/Preconditioner/Jacobi.hpp         |  66 ++++++---
 .../Preconditioner/PreconditionerWrapper.hpp  |  21 ++-
 include/OGL/Preconditioner/Schwarz.hpp        | 138 ++++++++++++++----
 4 files changed, 177 insertions(+), 170 deletions(-)

diff --git a/include/OGL/Preconditioner.hpp b/include/OGL/Preconditioner.hpp
index cb4d60b9..a5ef6b58 100644
--- a/include/OGL/Preconditioner.hpp
+++ b/include/OGL/Preconditioner.hpp
@@ -56,82 +56,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));
-        }
-    }
-
 
     std::shared_ptr<gko::LinOp> init_preconditioner_impl(
         const word name, const dictionary &d,
@@ -144,51 +68,7 @@ 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));
-            }
+            return BlockJacobi(device_exec, gkomatrix, d, verbose_).create();
         }
         if (name == "ILU") {
             label iterations(d.lookupOrDefault("iterations", label(0)));
diff --git a/include/OGL/Preconditioner/Jacobi.hpp b/include/OGL/Preconditioner/Jacobi.hpp
index 4e1dbd97..40ea00c7 100644
--- a/include/OGL/Preconditioner/Jacobi.hpp
+++ b/include/OGL/Preconditioner/Jacobi.hpp
@@ -3,45 +3,71 @@
 // SPDX-License-Identifier: GPL-3.0-or-later
 #pragma once
 
-#include "PreconditionerWrapper.hpp"
 #include "OGL/Preconditioner/Schwarz.hpp"
 
-class BJ : public PreconditionerWrapper {
+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:
-    BJ(const dictionary &d, label verbose)
-        : verbose_(verbose),
+    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)
+        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_local()
+    virtual std::shared_ptr<gko::LinOp> create()
     {
-        // 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));
-    }
+        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_);
+        };
+
+        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/PreconditionerWrapper.hpp b/include/OGL/Preconditioner/PreconditionerWrapper.hpp
index 0e92f5ef..939c0e08 100644
--- a/include/OGL/Preconditioner/PreconditionerWrapper.hpp
+++ b/include/OGL/Preconditioner/PreconditionerWrapper.hpp
@@ -5,11 +5,28 @@
 
 #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_local();
+    virtual std::shared_ptr<gko::LinOp> create();
 
     /* @brief instantiate a local preconditioner and return LinOp Ptr */
-    virtual std::shared_ptr<gko::LinOp> create_multilevel();
+    // 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
index fc552c91..63976aa4 100644
--- a/include/OGL/Preconditioner/Schwarz.hpp
+++ b/include/OGL/Preconditioner/Schwarz.hpp
@@ -6,33 +6,117 @@
 #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::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>
+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)
+{
+    using ras =
+        gko::experimental::distributed::preconditioner::Schwarz<scalar, label,
+                                                                label>;
 
-    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));
-    }
+    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)
+{
+    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

From 25653f7fd7e6b72debafd7d45dbb614a3bfb41ea Mon Sep 17 00:00:00 2001
From: Gregor Olenik <gregor.olenik@web.de>
Date: Wed, 27 Aug 2025 10:08:29 +0200
Subject: [PATCH 07/12] wip clean up distribute matrix

---
 include/OGL/MatrixWrapper/Distributed.hpp |  39 +-
 include/OGL/MatrixWrapper/lduMatrix.hpp   |   0
 include/OGL/Preconditioner.hpp            | 523 +++++++---------------
 include/OGL/Preconditioner/Jacobi.hpp     |   9 +-
 include/OGL/Preconditioner/Multigrid.hpp  | 262 +++++++++++
 include/OGL/Preconditioner/Schwarz.hpp    |   4 +-
 include/OGL/Solver/Multigrid.hpp          |   1 -
 src/MatrixWrapper/Distributed.cpp         |  28 +-
 8 files changed, 474 insertions(+), 392 deletions(-)
 create mode 100644 include/OGL/MatrixWrapper/lduMatrix.hpp
 create mode 100644 include/OGL/Preconditioner/Multigrid.hpp

diff --git a/include/OGL/MatrixWrapper/Distributed.hpp b/include/OGL/MatrixWrapper/Distributed.hpp
index bda53b87..7220f16a 100644
--- a/include/OGL/MatrixWrapper/Distributed.hpp
+++ b/include/OGL/MatrixWrapper/Distributed.hpp
@@ -80,6 +80,25 @@ class RepartDistMatrix
         }
     }
 
+    std::shared_ptr<const gko::LinOp> clone() const
+    {
+
+        // TODO needs to update pairwise_update_data and all_to_all_update_data
+        // and linops
+
+        // std::map<label, scalar *> linops;
+        // this->linops_;
+        this->all_to_all_update_data_;
+        this->pairwise_update_data_;
+        return std::make_shared<RepartDistMatrix>(
+            this->get_executor(), this->get_communicator(),
+            this->matrix_format_, this->dist_mtx_->clone(),
+            this->repartitioner_, this->fuse_, this->all_to_all_update_data_,
+            this->pairwise_update_data_, this->reorder_maps_,
+            this->compress_to_global_
+            // , this->linops_
+                                                  );
+    }
 
     /**
      * Copy-assigns a CombinationMatrix matrix. Preserves executor, copies
@@ -101,7 +120,7 @@ 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_;
+            // this->linops_ = other.linops_;
         }
         return *this;
     }
@@ -126,7 +145,7 @@ 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_);
+            // this->linops_ = std::move(other.linops_);
         }
         return *this;
     }
@@ -147,8 +166,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 +179,8 @@ 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) //,
+          // linops_(linops)
     {
         this->set_size(dist_mtx_->get_size());
     }
@@ -229,8 +250,10 @@ class RepartDistMatrix
 
     std::vector<label> compress_to_global_;
 
-    std::map<label, scalar *>
-        linops_;  // map between linop id and shared ptr to linop
+    // map between host interface id and pointer to begin of block of  data
+    // before reordering
+    // std::map<label, scalar *>
+    //     linops_;
 };
 
 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 a5ef6b58..4ee761e0 100644
--- a/include/OGL/Preconditioner.hpp
+++ b/include/OGL/Preconditioner.hpp
@@ -9,6 +9,7 @@
 #include "OGL/DevicePersistent/Base.hpp"
 #include "OGL/MatrixWrapper/Distributed.hpp"
 #include "OGL/Preconditioner/Jacobi.hpp"
+#include "OGL/Preconditioner/Multigrid.hpp"
 #include "OGL/Preconditioner/Schwarz.hpp"
 
 #include "fvCFD.H"
@@ -22,7 +23,6 @@ class Preconditioner {
     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;
@@ -70,370 +70,167 @@ class Preconditioner {
         if (name == "BJ") {
             return BlockJacobi(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 == "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);
-        }
-        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 coarsening = d.lookupOrDefault("coarsening", word("GAMG"));
-            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\tcoarsening: " + coarsening +
-                       "\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);
-            }
-
-            std::shared_ptr<gko::matrix::Csr<double, int>> coarseningWeight;
-            if (coarsening == "GAMG") {
-                auto &fvmesh =
-                    db_.template lookupObjectRef<fvMesh>("fvSchemes");
-                // weights[facei] -> column
-                auto weights =
-                    mag(cmptMultiply(fvmesh.Sf().primitiveField() /
-                                         sqrt(fvmesh.magSf().primitiveField()),
-                                     vector(1, 1.01, 1.02)));
-
-                //   const std::shared_ptr<HostMatrixWrapper>
-                //   weight_matrix_wrapper
-                //       {std::make_shared<HostMatrixWrapper>(
-                // exec_handler, db_, matrix.diag().size(),
-                // matrix.upper().size(), matrix.symmetric(),
-                // matrix.diag().begin(), matrix.upper().begin(),
-                // matrix.lower().begin(), matrix.lduAddr(), interfaceBouCoeffs,
-                // interfaceIntCoeffs, interfaces, solverControls, fieldName,
-                // verbose_)
-                // };
-
-
-                // std::shared_ptr<RepartDistMatrix> dist_weight_mtx =
-                //     create_distributed(&exec_handler, repartitioner,
-                //                        weigh_matrix_wrapper, "Csr" true,
-                //                        verbose_);
-
-                // here ldu adressing is needed
-                // 1. create row, cols, vals vector
-                // 2. fill with diagonals i=j=row=col
-                // 3. fill with off-diagonals values from weights and row, cols
-                // from ldu
-
-
-                coarseningWeight = nullptr;
-            }
-
-            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)
-                                .with_local_weight_mtx(coarseningWeight)
-                                .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_local_weight_mtx(coarseningWeight)
-                                .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 Multigrid(device_exec, gkomatrix, d, verbose_)
+                .create(db_, exec_handler);
         }
+        // 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 == "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);
+        // }
+        // 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 {};
         }
diff --git a/include/OGL/Preconditioner/Jacobi.hpp b/include/OGL/Preconditioner/Jacobi.hpp
index 40ea00c7..219c70d6 100644
--- a/include/OGL/Preconditioner/Jacobi.hpp
+++ b/include/OGL/Preconditioner/Jacobi.hpp
@@ -12,7 +12,7 @@ class BlockJacobi  // : public PreconditionerWrapper
 
     std::shared_ptr<gko::Executor> exec_;
     std::shared_ptr<const gko::LinOp> mtx_;
-    const dictionary &d;
+    const dictionary &d_;
     const label verbose_;
     bool skip_sorting_;
     bool multi_level_schwarz_;
@@ -43,9 +43,10 @@ class BlockJacobi  // : public PreconditionerWrapper
     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_);
+            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) {
diff --git a/include/OGL/Preconditioner/Multigrid.hpp b/include/OGL/Preconditioner/Multigrid.hpp
new file mode 100644
index 00000000..f9c6f227
--- /dev/null
+++ b/include/OGL/Preconditioner/Multigrid.hpp
@@ -0,0 +1,262 @@
+// SPDX-FileCopyrightText: 2025 OGL authors
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+#pragma once
+
+#include "OGL/Preconditioner/Schwarz.hpp"
+
+class Multigrid {
+    using sor = gko::preconditioner::Sor<scalar, label>;
+    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_;
+
+    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);
+    }
+
+    virtual 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> bjfac{};
+
+        if (smoother_ == "Jacobi") {
+            bjfac =
+                dbj::build().with_max_block_size(1u).with_skip_sorting(true).on(
+                    exec_);
+        }
+        if (smoother_ == "SOR") {
+            bjfac =
+                sor::build().with_skip_sorting(true).with_symmetric(false).on(
+                    exec_);
+        }
+        if (smoother_ == "SSOR") {
+            bjfac =
+                sor::build().with_skip_sorting(true).with_symmetric(true).on(
+                    exec_);
+        }
+        if (bjfac == nullptr) {
+            FatalErrorInFunction << "Unknown smoother: " << smoother_
+                                 << "\nValid Choices: Jacobi, SOR, SSOR"
+                                 << abort(FatalError);
+        }
+
+        std::shared_ptr<gko::matrix::Csr<double, int>> coarseningWeight;
+        if (coarsening_ == "GAMG") {
+            auto repartDistMtx = gko::as<RepartDistMatrix>(mtx_)->clone();
+            auto &fvmesh = db.template lookupObjectRef<fvMesh>("fvSchemes");
+
+            // // weights[facei] -> column
+            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], 1);
+
+            // const std::shared_ptr<HostMatrixWrapper>
+            //   weight_matrix_wrapper
+            //       {std::make_shared<HostMatrixWrapper>(
+            // exec_handler, db, diag.size(),
+            // weights.size(), true,
+            // diag.begin(), weights.begin(),
+            // weights.begin(), matrix.lduAddr(), interfaceBouCoeffs,
+            // interfaceIntCoeffs, interfaces, solverControls, fieldName,
+            // verbose_)
+            // };
+
+
+            // std::shared_ptr<RepartDistMatrix> dist_weight_mtx =
+            //     create_distributed(&exec_handler, repartitioner,
+            //                        weigh_matrix_wrapper, "Csr" true,
+            //                        verbose_);
+
+            // here ldu adressing is needed
+            // 1. create row, cols, vals vector
+            // 2. fill with diagonals i=j=row=col
+            // 3. fill with off-diagonals values from weights and row, cols
+            // from ldu
+
+
+            coarseningWeight = nullptr;
+        }
+
+        // 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)
+        //                                .with_local_weight_mtx(coarseningWeight)
+        //                                .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_local_weight_mtx(coarseningWeight)
+        //                                .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;
+
+        return {};
+    }
+};
diff --git a/include/OGL/Preconditioner/Schwarz.hpp b/include/OGL/Preconditioner/Schwarz.hpp
index 63976aa4..ef291283 100644
--- a/include/OGL/Preconditioner/Schwarz.hpp
+++ b/include/OGL/Preconditioner/Schwarz.hpp
@@ -10,7 +10,7 @@ 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)
+    std::shared_ptr<PrecondFactory> precond)
 {
     using ras =
         gko::experimental::distributed::preconditioner::Schwarz<scalar, label,
@@ -28,7 +28,7 @@ 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<PrecondFactory> precond,
     std::shared_ptr<Factorization> factorization)
 {
     using ras =
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/src/MatrixWrapper/Distributed.cpp b/src/MatrixWrapper/Distributed.cpp
index 72eaae43..5c5fffc6 100644
--- a/src/MatrixWrapper/Distributed.cpp
+++ b/src/MatrixWrapper/Distributed.cpp
@@ -381,7 +381,7 @@ void update_impl(
     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();
@@ -419,7 +419,7 @@ 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,
+                                 &get_send_ptr, fuse,
                                  force_host_buffer]() {
         for (auto [id, mode, comm_rank, length, send_id, recv_ptr] :
              pairwise_update_data) {
@@ -481,7 +481,7 @@ void RepartDistMatrix::update(const ExecutorHandler &exec_handler,
     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););
+                            verbose););
 }
 
 
@@ -496,7 +496,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 +525,12 @@ 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 +542,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 +574,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,7 +589,7 @@ 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,
+                    exec_handler, local_sparsity, host_to_linop_map, fuse, owner,
                     ranks_per_owner, all_to_all_update_data););
 
     size_t start_local_offset = 0;
@@ -610,12 +608,12 @@ std::shared_ptr<RepartDistMatrix> create_impl(
         generate_pairwise_update_data<LocalMatrixType>(
             exec_handler, host_A,
             (!owner) ? local_sparsity : repart_loc_sparsity, fuse,
-            repartitioner, linops, start_local_offset, pairwise_update_data););
+            repartitioner, host_to_linop_map, 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););
+                    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
@@ -673,13 +671,15 @@ std::shared_ptr<RepartDistMatrix> create_impl(
 
     SIMPLE_TIME(verbose, perform_matrix_update,
                 update_impl(exec_handler, host_A, all_to_all_update_data,
-                            pairwise_update_data, reorder_maps, fuse, linops,
+                            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
+        // , host_to_linop_map
+                                              );
 }
 
 void write_distributed(const ExecutorHandler &exec_handler, word field_name,

From 17c539c82e6a6667147c6660934c186211547299 Mon Sep 17 00:00:00 2001
From: Gregor Olenik <gregor.olenik@web.de>
Date: Mon, 1 Sep 2025 11:44:25 +0200
Subject: [PATCH 08/12] update

---
 include/OGL/MatrixWrapper/Distributed.hpp | 109 +++++---
 include/OGL/Preconditioner/Multigrid.hpp  | 305 ++++++++++++----------
 src/MatrixWrapper/Distributed.cpp         | 135 +++++++---
 3 files changed, 339 insertions(+), 210 deletions(-)

diff --git a/include/OGL/MatrixWrapper/Distributed.hpp b/include/OGL/MatrixWrapper/Distributed.hpp
index 7220f16a..9098f647 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,24 +119,33 @@ class RepartDistMatrix
         }
     }
 
-    std::shared_ptr<const gko::LinOp> clone() const
+    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)
+        );
 
-        // TODO needs to update pairwise_update_data and all_to_all_update_data
-        // and linops
-
-        // std::map<label, scalar *> linops;
-        // this->linops_;
-        this->all_to_all_update_data_;
-        this->pairwise_update_data_;
         return std::make_shared<RepartDistMatrix>(
             this->get_executor(), this->get_communicator(),
-            this->matrix_format_, this->dist_mtx_->clone(),
-            this->repartitioner_, this->fuse_, this->all_to_all_update_data_,
-            this->pairwise_update_data_, this->reorder_maps_,
-            this->compress_to_global_
-            // , this->linops_
-                                                  );
+            this->matrix_format_, new_dist_mtx, this->repartitioner_,
+            this->fuse_, all_to_all_data, pairwise_update_data,
+            new_reorder_map, this->compress_to_global_);
     }
 
     /**
@@ -120,7 +168,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;
     }
@@ -145,7 +192,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;
     }
@@ -154,6 +200,11 @@ 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_; }
@@ -179,8 +230,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());
     }
@@ -240,20 +290,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_;
-
-    // map between host interface id and pointer to begin of block of  data
-    // before reordering
-    // std::map<label, scalar *>
-    //     linops_;
 };
 
 std::shared_ptr<const gko::LinOp> get_local(
diff --git a/include/OGL/Preconditioner/Multigrid.hpp b/include/OGL/Preconditioner/Multigrid.hpp
index f9c6f227..1f367239 100644
--- a/include/OGL/Preconditioner/Multigrid.hpp
+++ b/include/OGL/Preconditioner/Multigrid.hpp
@@ -6,9 +6,17 @@
 #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 mg = gko::solver::Multigrid;
+    using pgm = gko::multigrid::Pgm<scalar, label>;
 
     std::shared_ptr<gko::Executor> exec_;
 
@@ -74,188 +82,201 @@ class Multigrid {
         //                      << abort(FatalError);
     }
 
-    virtual std::shared_ptr<gko::LinOp> create(
-        const objectRegistry &db, const ExecutorHandler &exec_handler)
+    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> bjfac{};
-
+        std::shared_ptr<gko::LinOpFactory> smootherFac{};
         if (smoother_ == "Jacobi") {
-            bjfac =
+            smootherFac =
                 dbj::build().with_max_block_size(1u).with_skip_sorting(true).on(
                     exec_);
         }
         if (smoother_ == "SOR") {
-            bjfac =
+            smootherFac =
                 sor::build().with_skip_sorting(true).with_symmetric(false).on(
                     exec_);
         }
         if (smoother_ == "SSOR") {
-            bjfac =
+            smootherFac =
                 sor::build().with_skip_sorting(true).with_symmetric(true).on(
                     exec_);
         }
-        if (bjfac == nullptr) {
+        if (smootherFac == nullptr) {
             FatalErrorInFunction << "Unknown smoother: " << smoother_
                                  << "\nValid Choices: Jacobi, SOR, SSOR"
                                  << abort(FatalError);
         }
 
-        std::shared_ptr<gko::matrix::Csr<double, int>> coarseningWeight;
+        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");
-
-            // // weights[facei] -> column
             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], 1);
+                                 vector(-1, -1.01, -1.02)))();
 
-            // const std::shared_ptr<HostMatrixWrapper>
-            //   weight_matrix_wrapper
-            //       {std::make_shared<HostMatrixWrapper>(
-            // exec_handler, db, diag.size(),
-            // weights.size(), true,
-            // diag.begin(), weights.begin(),
-            // weights.begin(), matrix.lduAddr(), interfaceBouCoeffs,
-            // interfaceIntCoeffs, interfaces, solverControls, fieldName,
-            // verbose_)
-            // };
+            std::vector<scalar> diag(mtx_->get_size()[0], 0.00001);
 
+            repartDistMtx->update(
+                                  exec_handler,
+                                  diag.data(),
+                                  weights.begin(),
+                                  verbose_);
 
-            // std::shared_ptr<RepartDistMatrix> dist_weight_mtx =
-            //     create_distributed(&exec_handler, repartitioner,
-            //                        weigh_matrix_wrapper, "Csr" true,
-            //                        verbose_);
-
-            // here ldu adressing is needed
-            // 1. create row, cols, vals vector
-            // 2. fill with diagonals i=j=row=col
-            // 3. fill with off-diagonals values from weights and row, cols
-            // from ldu
+            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";
+        }
 
-            coarseningWeight = nullptr;
+        if (coarsening_ == "PGM")
+        {
+            auto repartDistMtx = gko::as<RepartDistMatrix>(mtx_);
+            auto cw = gko::as<gko::matrix::Csr<scalar, label>>(
+                repartDistMtx->get_local_matrix());
+            coarseningWeight =
+                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 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));
+        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(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);
-        //     }
+        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_ == "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);
+            }
 
-        //     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)
-        //                                .with_local_weight_mtx(coarseningWeight)
-        //                                .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));
-        // }
+            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(false)
+                                       .with_local_weight_mtx(coarseningWeight)
+                                       .on(exec_))
+                    .with_coarsest_solver(coarsest_solver)
+                    .with_criteria(single_it)
+                    .on(exec_));
+            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(
-        //             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_local_weight_mtx(coarseningWeight)
-        //                                .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;
+        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_ == "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(false))
+                    .with_coarsest_solver(coarsest_solver)
+                    .with_criteria(single_it)
+                    .on(exec_)
+                    ->generate(gkodistmatrix));
+            return ret;
+        }
 
         return {};
     }
diff --git a/src/MatrixWrapper/Distributed.cpp b/src/MatrixWrapper/Distributed.cpp
index 5c5fffc6..8dfefe5f 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,6 +374,56 @@ 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,
@@ -395,10 +444,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 +470,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, 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 +487,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 +529,22 @@ 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_,
-                            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>>
@@ -530,7 +591,8 @@ std::shared_ptr<RepartDistMatrix> create_impl(
         (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, host_to_linop_map, 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,
@@ -589,8 +651,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, host_to_linop_map, 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) {
@@ -603,17 +665,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, host_to_linop_map, 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, host_to_linop_map, 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
@@ -669,17 +732,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,
-                            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
-        // , host_to_linop_map
-                                              );
+        compress_to_global);
 }
 
 void write_distributed(const ExecutorHandler &exec_handler, word field_name,

From a541605ec2b1f0da73247a1469acf76fa2177991 Mon Sep 17 00:00:00 2001
From: Gregor Olenik <gregor.olenik@web.de>
Date: Mon, 1 Sep 2025 11:45:56 +0200
Subject: [PATCH 09/12] update

---
 include/OGL/MatrixWrapper/Distributed.hpp | 15 ++++------
 include/OGL/Preconditioner/Multigrid.hpp  | 30 +++++++++----------
 src/MatrixWrapper/Distributed.cpp         | 35 +++++++++++------------
 3 files changed, 36 insertions(+), 44 deletions(-)

diff --git a/include/OGL/MatrixWrapper/Distributed.hpp b/include/OGL/MatrixWrapper/Distributed.hpp
index 9098f647..7060bca5 100644
--- a/include/OGL/MatrixWrapper/Distributed.hpp
+++ b/include/OGL/MatrixWrapper/Distributed.hpp
@@ -58,7 +58,7 @@ class RepartDistMatrix
         std::vector<reorder_map_type> in, scalar *new_ptr) const
     {
         std::vector<reorder_map_type> out;
-        auto [map, _, pad] = in[0]; //) {
+        auto [map, _, pad] = in[0];  //) {
         out.emplace_back(map, new_ptr, pad);
         // }
         return out;
@@ -138,14 +138,13 @@ class RepartDistMatrix
             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)
-        );
+            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_);
+            this->fuse_, all_to_all_data, pairwise_update_data, new_reorder_map,
+            this->compress_to_global_);
     }
 
     /**
@@ -200,10 +199,8 @@ 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);
+    void update(const ExecutorHandler &exec_handler, const scalar *diag_ptr,
+                const scalar *face_ptr, label verbose);
 
     word get_matrix_format() const { return matrix_format_; }
 
diff --git a/include/OGL/Preconditioner/Multigrid.hpp b/include/OGL/Preconditioner/Multigrid.hpp
index 1f367239..6a03fa04 100644
--- a/include/OGL/Preconditioner/Multigrid.hpp
+++ b/include/OGL/Preconditioner/Multigrid.hpp
@@ -113,8 +113,8 @@ class Multigrid {
         }
 
         std::shared_ptr<gko::LinOpFactory> bjfac =
-                dbj::build().with_max_block_size(1u).with_skip_sorting(true).on(
-                    exec_);
+            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") {
@@ -127,10 +127,7 @@ class Multigrid {
 
             std::vector<scalar> diag(mtx_->get_size()[0], 0.00001);
 
-            repartDistMtx->update(
-                                  exec_handler,
-                                  diag.data(),
-                                  weights.begin(),
+            repartDistMtx->update(exec_handler, diag.data(), weights.begin(),
                                   verbose_);
 
             auto cw = gko::as<gko::matrix::Csr<scalar, label>>(
@@ -139,14 +136,15 @@ class Multigrid {
                 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] << ",";
+            // for (auto i=0; i<coarseningWeight->get_num_stored_elements();
+            // i++){
+            //   std::cout << i << ":" << coarseningWeight->get_values()[i] <<
+            //   ",";
             // }
             // std::cout << "\n";
         }
 
-        if (coarsening_ == "PGM")
-        {
+        if (coarsening_ == "PGM") {
             auto repartDistMtx = gko::as<RepartDistMatrix>(mtx_);
             auto cw = gko::as<gko::matrix::Csr<scalar, label>>(
                 repartDistMtx->get_local_matrix());
@@ -170,12 +168,12 @@ class Multigrid {
 
         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(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_)
diff --git a/src/MatrixWrapper/Distributed.cpp b/src/MatrixWrapper/Distributed.cpp
index 8dfefe5f..62b57cc6 100644
--- a/src/MatrixWrapper/Distributed.cpp
+++ b/src/MatrixWrapper/Distributed.cpp
@@ -375,9 +375,8 @@ void reorder_interface_impl(const ExecutorHandler &exec_handler,
 }
 
 void advanced_update_impl(
-    const ExecutorHandler &exec_handler,
-    const scalar* diag_weights,
-    const scalar* face_weights,
+    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,
@@ -393,20 +392,19 @@ void advanced_update_impl(
 
     // 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]() {
+                              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);
+            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);
             }
         }
     };
@@ -536,13 +534,12 @@ void RepartDistMatrix::update(const ExecutorHandler &exec_handler,
 }
 
 void RepartDistMatrix::update(const ExecutorHandler &exec_handler,
-    const scalar* diag_ptr,
-    const scalar* face_ptr,
+                              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_,
+    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););
 }
 

From 595d110a072a7136162a5ff4cf12711232f26e0c Mon Sep 17 00:00:00 2001
From: Gregor Olenik <go@hpsim.de>
Date: Fri, 12 Sep 2025 12:04:55 +0200
Subject: [PATCH 10/12] add divergence check

---
 src/StoppingCriterion.cpp | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

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);

From 5f646312173cc87e7b18330eb1328c9455c066a3 Mon Sep 17 00:00:00 2001
From: Gregor Olenik <go@hpsim.de>
Date: Fri, 12 Sep 2025 12:06:03 +0200
Subject: [PATCH 11/12] add further coarse solvers and move implementation to
 specific files

---
 include/OGL/Preconditioner/.Multigrid.hpp.swo | Bin 0 -> 16384 bytes
 include/OGL/Preconditioner/Cholesky.hpp       |  72 +++++++++++++++++
 include/OGL/Preconditioner/LU.hpp             |  75 ++++++++++++++++++
 include/OGL/Preconditioner/Multigrid.hpp      |  52 +++++++++---
 include/OGL/Preconditioner/Schwarz.hpp        |   2 +
 5 files changed, 189 insertions(+), 12 deletions(-)
 create mode 100644 include/OGL/Preconditioner/.Multigrid.hpp.swo
 create mode 100644 include/OGL/Preconditioner/Cholesky.hpp
 create mode 100644 include/OGL/Preconditioner/LU.hpp

diff --git a/include/OGL/Preconditioner/.Multigrid.hpp.swo b/include/OGL/Preconditioner/.Multigrid.hpp.swo
new file mode 100644
index 0000000000000000000000000000000000000000..04ce74617f1f6b00d19bcde354df6a290130db60
GIT binary patch
literal 16384
zcmeHOO>87b6|Q8L1VR!}qJ@Y^Q87rXwmi1S_GZQQ+KwV`))>68mDwZ^Hfr^BmuK4U
z>6%pc?0B;z0f|!p34z3qIB-Y~91<y)NI{AcMM@-2AW(!LClZ9XaL55}!17g9|4qBy
zSRf8a)zaslu6p(Ay{}&VcDvVkZ1WQR{QMHZ^+`fb|LC>BshzvOK6!GPu&2&+me0zJ
z=4-GYEuK5qnGd<=hCUy<K~!7dS*>uvcN2HP<33y9{cs_Oyl{-)+Z$U8JA!#U@`EJc
z5fclS#$gilMc~g5Mx$DlCebip7&tZqcaY_!#f}vJT<f#+?z<j6c1a^>7%&VN1`Gp+
z0mFb{z%XDKFbw?vGmwl<lIP&Y&8aI-r1zUD_Z#Vz?HImV$^U1Xzn!MPQprD&KHr<(
zZ&dQrM{_d_7zPXjh5^HXVZbn87%&VN1`Gp+0mFb{;D5k?>l5-I^u8lE{CNML?f>8X
z1R=i$UIcy!`~V1nec&s=B5*tK_d5xB8DKyg`0x%wUIo4b>;g-`B5*5k3-I3Ug!~?O
z6L<r75qJUk7Vtc<59|R?01p5ba2N28j}!7|;5Wb<z|Vm10@s130RcP!{OvYEeh0h-
zyafCNcow)HI0gLbV}yJU*aU6{e*RHHz6P8J{(37RuK;_%eZYIS5b`eYZQu%U3i#tk
z2>CIv1AGSf!%0HE4{QP_fIs1&<EOx7;6C6U;M2fwalrCC@J--x;7&kt{&sfI^53*A
zO|M{=THZA;WDmPT=Frtu+M2Xxs7j)%w7<`7JLch4CTx3{eOqDMD&Mvz=ShR2tH)P#
zwNTVpr8Z0pU(f5#Y&hbM??d8*K6Jpmal*wH+{h1^pr;uUoPmCtQM%&9hqnV)Fy9#^
zVnvy<6-3*kHR=iGCd|^hcyEt+$s?>E#EG~@Py4+#Q5v36f2AD@3-a&bYh{JHJuVXK
zB0gye9}<F#nZo&oQmre|V{RDubbY`BkHt3adV`58o^I0%(u|-tPMF^!S(Q*#RIwRd
z97mq~Mux(G)}zsW#3vEWzm}oNwGp#vOXkX-xt2>?+Z#^DLi_i4%$znI#(hmzR!@3|
zc~=r07Em**8%hV&ELE9V6|<@ar*OKYNeKoKzU1kupnI5$p)_ljs<uhEBR!4!mWoN4
z>&(#E+}f(9p?k2IFtN^E5i9o&>6wz~HY)R!7HyZdshr|8n4#q&*j+9spPNdhrgev$
zCj<1qIvO&WM#BXQ-2+!P{zcbA52)3ogpl8;sbWiNS)&s{q}hLjPY|*-EjbJ#nYVsa
z#g&GK$6M?w3*#yI2^YQ$O$=GB4TkVRUoA!VXzC*|TZU?`$A(QRfXd@Z&a36)8#TT^
z9&L*YY}dsIXtlD;)csbQhHj6A*5b@exob>CN1dLs0F+1{tkxG@kN1KWRJBWCT||&`
zv?du)xt@&&FKtw{#94#Y>Bus34fU4f%~huZnMG^oyqrDL>2xX%<@DG3ljr6-mZ}Ph
zNvqBfm79BliD;V6JlN~XR=W$-#P~|l6;@~dTt!ha3>~5AG7mK)g_cD&cebJ?jVLYI
zJ&KZID9dS)Ay<wQImXr&ml-)9>BLcGMUI|_!hNMX2@-FxhOn2@zv_Eg`8M5k!<bE*
zpv0U|_E9I+Gn8YMLOI?KMox?wP!4~GSu5(3f?X{K%M&rKRFut0+2k}ugb&qI3B%;T
zvGAvqXj8QWu;kxTvZR6PE=|)q%H3LKRu>9-RZw0opr$s}$$p?FURNy?bYtc7u2t2f
zbYx4FHDzF%jK;k%@NB)rP-6|v#-MZ>RpXx?B=iJzI9*LQdeT+WS$|I@Q?=P(tOfJB
z(jm80yy$;AQFjrV_pGVZqMKCo79rg&TdczPy!N-pc}VMms+>O+RB480U%QaS%2d`=
z!=QR{?4~_RxTf_o(Q%A;r=$KbTU9HU<d@2}z0)7AX#UcvWKGkg%M9hLrZg)Wu@KM0
zr!3|j5?k-5TBp?Q_VzM;m1)t|BmT;6$lavQ693;8^*KV7)s;n$KdM@E@>=(0bNr{5
zYWCD`sYCI{4n|A_Ll$A>ZQFjaYPRhv73mN+ZI|x6x>bZ5(^Vo>>*pZJ+=~-FvTYgJ
zYRy${qaf!XYs$^cb=c~KZXB0ZE!XkQZThJ@ZPk41p4uqCbZ<d)B;5Ycr9AQ&!Q#KW
zbK%Q#Ti7W?F`L`;;dHPY$Ytcl&eq)Wd}oe}xe&`F0pVQ!Vi2-*KDs7?{vf%`4icL#
zb(Zd@sGYjwWWYsC<o<sZ`{s85x&Oa6-G|Ek|7*Y=une34KEVF|P2dGU0AB(=1-y@a
z{mZ~}z&3CL`}Nm>Zvb833~&<o3+ng)cpvyB@C!iJ`!--M!+>GHFkl!k3>XFs1BL;^
zfMMX63_PQ^CB^#+!IH6v=yhbMU5(t?6uI&#mB3Axm1*T`Zmpl5yjs|8$=%h^#bN!y
zOzHhbS$>{J>gdI(o~~%2d7L{998Uz;q6Drbx6dJSf}~Me9nr@L&X?pXE4;0d2mkqV
zebKT_>2F=P&rmgvG~r3$S?N(1sS&SB(kj&dH#-w1w^I(zgi?#-S<K@qvEC$Ic-Ys+
zXckTlY`ZrOLf>kVA8~j(hQAJWvz9tAtF=~Z&(D7(!LgMeC^gm~n3I*WvL6nby{Qlw
zbaX>F0tfw8euCUoNF6WfH<L|buv&AZ9LQb_9P)bDujdb2S}x-{-@40pf91o<qGgL{
zJt~|l@tU2rH%Dx7*zD}8morv=g4X1B8s9jC$Q>&OT$-J=5MJv{y_wKn1vp8n3o1oQ
LR#Y}_g`Ix^z7*hE

literal 0
HcmV?d00001

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/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
index 6a03fa04..b26fabd0 100644
--- a/include/OGL/Preconditioner/Multigrid.hpp
+++ b/include/OGL/Preconditioner/Multigrid.hpp
@@ -15,6 +15,8 @@ class Multigrid {
     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>;
 
@@ -145,17 +147,17 @@ class Multigrid {
         }
 
         if (coarsening_ == "PGM") {
-            auto repartDistMtx = gko::as<RepartDistMatrix>(mtx_);
-            auto cw = gko::as<gko::matrix::Csr<scalar, label>>(
-                repartDistMtx->get_local_matrix());
-            coarseningWeight =
-                std::const_pointer_cast<gko::matrix::Csr<scalar, label>>(cw);
-        }
-        if (coarseningWeight == nullptr) {
-            FatalErrorInFunction << "Unknown coarsening: " << coarsening_
-                                 << "\nValid Choices: GAMG, PGM"
-                                 << abort(FatalError);
+            // 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(
@@ -189,6 +191,22 @@ class Multigrid {
                         .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()
@@ -203,6 +221,7 @@ class Multigrid {
                                      << abort(FatalError);
             }
 
+	    std::cout << __FILE__ << __LINE__ << " MG 1\n";
             auto pre_factory = gko::share(
                 mg::build()
                     .with_max_levels(static_cast<gko::uint32>(maxLevels_))
@@ -212,12 +231,13 @@ class Multigrid {
                     .with_smoother_iters(maxIterS_)
                     .with_post_uses_pre(true)
                     .with_mg_level(pgm::build()
-                                       .with_deterministic(false)
+                                       .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));
         }
 
@@ -231,6 +251,14 @@ class Multigrid {
                         .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()
@@ -268,7 +296,7 @@ class Multigrid {
                     .with_post_uses_pre(true)
                     .with_mg_level(pgm::build()
                                        .with_local_weight_mtx(coarseningWeight)
-                                       .with_deterministic(false))
+                                       .with_deterministic(true))
                     .with_coarsest_solver(coarsest_solver)
                     .with_criteria(single_it)
                     .on(exec_)
diff --git a/include/OGL/Preconditioner/Schwarz.hpp b/include/OGL/Preconditioner/Schwarz.hpp
index ef291283..220d5471 100644
--- a/include/OGL/Preconditioner/Schwarz.hpp
+++ b/include/OGL/Preconditioner/Schwarz.hpp
@@ -16,7 +16,9 @@ std::shared_ptr<gko::LinOp> wrap_schwarz(
         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))

From fd38924d40de1a4b956db6cb12f10627c7df5c52 Mon Sep 17 00:00:00 2001
From: Gregor Olenik <go@hpsim.de>
Date: Fri, 12 Sep 2025 12:06:18 +0200
Subject: [PATCH 12/12] add further coarse solvers and move implementation to
 specific files

---
 include/OGL/Preconditioner.hpp | 29 +++++++++--------------------
 1 file changed, 9 insertions(+), 20 deletions(-)

diff --git a/include/OGL/Preconditioner.hpp b/include/OGL/Preconditioner.hpp
index 4ee761e0..dc61fc71 100644
--- a/include/OGL/Preconditioner.hpp
+++ b/include/OGL/Preconditioner.hpp
@@ -8,7 +8,9 @@
 
 #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"
 
@@ -74,11 +76,17 @@ class Preconditioner {
             return Multigrid(device_exec, gkomatrix, d, verbose_)
                 .create(db_, exec_handler);
         }
+        if (name == "IC") {
+            return Cholesky(device_exec, gkomatrix, d, verbose_).create();
+        }
+        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)
@@ -146,25 +154,6 @@ class Preconditioner {
         //     return wrap_schwarz(gkomatrix, device_exec,
         //                         std::move(precond_factory), factorization);
         // }
-        // 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));