From e172bfc8da1db514985693b18b70c767948b8fd3 Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Wed, 29 Apr 2026 16:18:15 -0500
Subject: [PATCH 01/19] mixed: noop mixed-adapter implementation

---
 opm/simulators/linalg/FlexibleSolver_impl.hpp | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/opm/simulators/linalg/FlexibleSolver_impl.hpp b/opm/simulators/linalg/FlexibleSolver_impl.hpp
index f37fd2a2111..d8506ce7d50 100644
--- a/opm/simulators/linalg/FlexibleSolver_impl.hpp
+++ b/opm/simulators/linalg/FlexibleSolver_impl.hpp
@@ -32,8 +32,9 @@
 #include <opm/simulators/linalg/PreconditionerFactoryGPUIncludeWrapper.hpp>
 #include <opm/simulators/linalg/is_gpu_operator.hpp>
 
-#if HAVE_AVX2_EXTENSION
+#if 1 //HAVE_AVX2_EXTENSION
 #include <opm/simulators/linalg/mixed/wrapper.hpp>
+#include <opm/simulators/linalg/mixed/adapter.hpp>
 #endif
 
 #include <dune/common/fmatrix.hh>
@@ -185,7 +186,7 @@ namespace Dune
                                                                             tol, // desired residual reduction factor
                                                                             maxiter, // maximum number of iterations
                                                                             verbosity);
-#if HAVE_AVX2_EXTENSION
+#if 1 //HAVE_AVX2_EXTENSION
           } else if (solver_type == "mixed-bicgstab") {
               if constexpr (Opm::is_gpu_operator_v<Operator>) {
                 OPM_THROW(std::invalid_argument, "mixed-bicgstab solver not supported for GPU operatorsg");
@@ -202,6 +203,17 @@ namespace Dune
                                                                             use_mixed_dilu
                                                                         );
             }
+            // mixed-solver adapter starts here
+          } else if (solver_type == "mixed-adapter") {
+
+                linsolver_ = std::make_shared<Dune::MixedAdapter<Operator,VectorType>>(*linearoperator_for_solver_,
+                                                                            *scalarproduct_,
+                                                                            *preconditioner_,
+                                                                            tol, // desired residual reduction factor
+                                                                            maxiter, // maximum number of iterations
+                                                                            verbosity);
+
+
 #endif
         } else if (solver_type == "loopsolver") {
             linsolver_ = std::make_shared<Dune::LoopSolver<VectorType>>(*linearoperator_for_solver_,

From 22aff658ad2147101755d4d51b033206a7756ba7 Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Thu, 30 Apr 2026 12:22:37 -0500
Subject: [PATCH 02/19] mixed: dummy mixed-precision adapter

---
 opm/simulators/linalg/FlexibleSolver_impl.hpp | 10 +--
 opm/simulators/linalg/mixed/adapter.hpp       | 68 +++++++++++++++++++
 opm/simulators/linalg/setupPropertyTree.cpp   |  3 +-
 3 files changed, 75 insertions(+), 6 deletions(-)
 create mode 100644 opm/simulators/linalg/mixed/adapter.hpp

diff --git a/opm/simulators/linalg/FlexibleSolver_impl.hpp b/opm/simulators/linalg/FlexibleSolver_impl.hpp
index d8506ce7d50..c1bae3b1377 100644
--- a/opm/simulators/linalg/FlexibleSolver_impl.hpp
+++ b/opm/simulators/linalg/FlexibleSolver_impl.hpp
@@ -204,11 +204,11 @@ namespace Dune
                                                                         );
             }
             // mixed-solver adapter starts here
-          } else if (solver_type == "mixed-adapter") {
-
-                linsolver_ = std::make_shared<Dune::MixedAdapter<Operator,VectorType>>(*linearoperator_for_solver_,
-                                                                            *scalarproduct_,
-                                                                            *preconditioner_,
+          } else if (solver_type == "mixed-precision") {
+                printf("made it!\n"); getchar();
+                linsolver_ = std::make_shared<Dune::MixedAdapter<Operator,VectorType>>(linearoperator_for_solver_,
+                                                                            scalarproduct_,
+                                                                            preconditioner_,
                                                                             tol, // desired residual reduction factor
                                                                             maxiter, // maximum number of iterations
                                                                             verbosity);
diff --git a/opm/simulators/linalg/mixed/adapter.hpp b/opm/simulators/linalg/mixed/adapter.hpp
new file mode 100644
index 00000000000..257e9e0d260
--- /dev/null
+++ b/opm/simulators/linalg/mixed/adapter.hpp
@@ -0,0 +1,68 @@
+#ifndef OPM_MIXED_ADAPTER_HEADER_INCLUDED
+#define OPM_MIXED_ADAPTER_HEADER_INCLUDED
+
+//#include <opm/simulators/flow/BlackoilModelParameters.hpp>
+//#include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
+#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
+//#include <dune/istl/solver.hh>
+
+namespace Dune
+{
+template <class Operator, class X>
+class MixedAdapter : public InverseOperator<X,X>
+{
+
+    public:
+    using AbstractPrecondType = Dune::PreconditionerWithUpdate<X, X>;
+    using AbstractScalarProductType = Dune::ScalarProduct<X>;
+
+
+    MixedAdapter(Operator *op,
+                 //Dune::ScalarProduct<X> &sp,
+                 std::shared_ptr<AbstractScalarProductType> sp,
+                 std::shared_ptr<AbstractPrecondType> prec,
+                 const double& tol,
+                 const int& maxiter,
+                 const int& verbosity)
+    {
+        //operator_ = op;
+        //prec_     = prec;
+        //product_  = sp;
+
+        //initialize bicgstab solver from Dune
+        solver_   = std::make_shared<Dune::BiCGSTABSolver<X>>(*op,
+                                                              *sp,
+                                                              *prec,
+                                                              tol, // desired residual reduction factor
+                                                              maxiter, // maximum number of iterations
+                                                              verbosity);
+    }
+
+    virtual void apply (X& x, X& b, InverseOperatorResult& res) override
+    {
+        //apply bicgstab solver from Dune
+        solver_->apply(x,b,res);
+    }
+
+    virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res) override
+    {
+        x=0;
+        b=0;
+        res.reduction = reduction;
+        OPM_THROW(std::invalid_argument, "MixedAdapter::apply(...) not implemented yet.");
+    }
+
+
+    private:
+    using AbstractSolverType = Dune::InverseOperator<X, X>;
+
+    //Operator* operator_;
+    //std::shared_ptr<AbstractPrecondType> prec_;
+    //std::shared_ptr<AbstractScalarProductType> product_;
+    std::shared_ptr<AbstractSolverType> solver_;
+
+};
+
+}
+
+#endif // OPM_MIXED_ADAPTER_HEADER_INCLUDED
diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp
index 6833ca20aac..236ca0528ae 100644
--- a/opm/simulators/linalg/setupPropertyTree.cpp
+++ b/opm/simulators/linalg/setupPropertyTree.cpp
@@ -285,7 +285,8 @@ std::string getSolverString(const FlowLinearSolverParameters& p)
     }
     else
     {
-        return {"bicgstab"};
+        //return {"bicgstab"};
+        return {"mixed-precision"};
     }
 }
 

From 274a354816548403cdccb102bbbd5bc9cfbccb73 Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Fri, 1 May 2026 17:10:17 -0500
Subject: [PATCH 03/19] mixed: working adapter with basic mixed spmv

---
 opm/simulators/linalg/FlexibleSolver_impl.hpp |   6 +-
 opm/simulators/linalg/mixed/adapter.hpp       | 138 ++++++++++++++++--
 2 files changed, 128 insertions(+), 16 deletions(-)

diff --git a/opm/simulators/linalg/FlexibleSolver_impl.hpp b/opm/simulators/linalg/FlexibleSolver_impl.hpp
index c1bae3b1377..c4217497b86 100644
--- a/opm/simulators/linalg/FlexibleSolver_impl.hpp
+++ b/opm/simulators/linalg/FlexibleSolver_impl.hpp
@@ -205,13 +205,13 @@ namespace Dune
             }
             // mixed-solver adapter starts here
           } else if (solver_type == "mixed-precision") {
-                printf("made it!\n"); getchar();
-                linsolver_ = std::make_shared<Dune::MixedAdapter<Operator,VectorType>>(linearoperator_for_solver_,
+                linsolver_ = std::make_shared<Dune::MixedAdapter<Comm,Operator,VectorType>>(linearoperator_for_solver_,
                                                                             scalarproduct_,
                                                                             preconditioner_,
                                                                             tol, // desired residual reduction factor
                                                                             maxiter, // maximum number of iterations
-                                                                            verbosity);
+                                                                            verbosity,
+                                                                            comm);
 
 
 #endif
diff --git a/opm/simulators/linalg/mixed/adapter.hpp b/opm/simulators/linalg/mixed/adapter.hpp
index 257e9e0d260..b955a888409 100644
--- a/opm/simulators/linalg/mixed/adapter.hpp
+++ b/opm/simulators/linalg/mixed/adapter.hpp
@@ -6,31 +6,108 @@
 #include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
 //#include <dune/istl/solver.hh>
 
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/istl/schwarz.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/owneroverlapcopy.hh>
+
+#include <dune/istl/matrixindexset.hh>
+
+/*
+// Type definitions
+typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>> BCRSMat;
+typedef Dune::BlockVector<Dune::FieldVector<double,1>> BVector;
+typedef Dune::OwnerOverlapCopyCommunication<int,int> Communication;
+
+//Assume parallel_A is the assembled BCRSMatrix and comm is the communicator
+//1. Create the parallel operator
+typedef Dune::OverlappingSchwarzOperator<BCRSMat,BVector,BVector,Communication> Operator;
+Operator op(parallel_A,comm);
+
+//2. Create the scalar product
+typedef Dune::OverlappingSchwarzScalarProduct<BVector,Communication> ScalarProduct;
+ScalarProduct sp(comm);
+
+//3. Create the parallel preconditioner (e.g., usingSSOR)
+typedef Dune::SeqSSOR<BCRSMat,BVector,BVector> Smoother;
+typedef Dune::BlockPreconditioner<BVector,BVector,Communication,Smoother> ParPrec;
+Smoother smoother(parallel_A,1,1.0);
+ParPrec pprec(smoother,comm);
+*/
+
 namespace Dune
 {
-template <class Operator, class X>
-class MixedAdapter : public InverseOperator<X,X>
+template <class Comm, class Operator, class Vector>
+class MixedAdapter:public InverseOperator<Vector, Vector>
 {
 
     public:
-    using AbstractPrecondType = Dune::PreconditionerWithUpdate<X, X>;
-    using AbstractScalarProductType = Dune::ScalarProduct<X>;
 
+    using AbstractPrecondType = Dune::PreconditionerWithUpdate<Vector,Vector>;
+    using AbstractScalarProductType = Dune::ScalarProduct<Vector>;
+
+    using typename InverseOperator<Vector, Vector>::domain_type;
+    static constexpr auto block_size = domain_type::block_type::dimension;
+
+    using SingleMatrixType  = Dune::BCRSMatrix<Opm::MatrixBlock<float, block_size, block_size>>;
+#if HAVE_MPI
+    using MixedOperatorType = Dune::OverlappingSchwarzOperator<SingleMatrixType, Vector, Vector, Comm>;
+#else // HAVE_MPI
+    using MixedOperatorType = Dune::MatrixAdapter<SingleMatrixType, Vector, Vector>;
+#endif // HAVE_MPI
 
     MixedAdapter(Operator *op,
-                 //Dune::ScalarProduct<X> &sp,
                  std::shared_ptr<AbstractScalarProductType> sp,
                  std::shared_ptr<AbstractPrecondType> prec,
                  const double& tol,
                  const int& maxiter,
-                 const int& verbosity)
+                 const int& verbosity,
+                 const Comm &comm)
     {
-        //operator_ = op;
-        //prec_     = prec;
-        //product_  = sp;
+
+        auto &A = op->getmat();
+        int nrows = A.N();
+        //int nnz = A.nonzeroes();
+        //int b = A[0][0].N();
+
+        single_matrix_ = std::make_shared<SingleMatrixType>();
+        auto &B = *single_matrix_;
+
+        //copy sparsity pattern
+        int irow=0;
+        int icol=0;
+        MatrixIndexSet sparsity(nrows,nrows);
+        for(auto row=A.begin();row!=A.end();row++)
+        {
+            for(unsigned int i=0;i<row->getsize();i++)
+            {
+                icol = row->getindexptr()[i];
+                sparsity.add(irow,icol);
+            }
+            irow++;
+        }
+        sparsity.exportIdx(B);
+        //B.compress();
+
+        //initializemixedoperator
+        double_operator_ = op;
+
+#if HAVE_MPI
+        mixed_operator_ = std::make_shared<MixedOperatorType>(B,comm);
+#else // HAVE_MPI
+        mixed_operator_ = std::make_shared<MixedOperatorType>(B);
+#endif // HAVE_MPI
+
+        //pointerstodataarrays
+        //double_data_=&A[0][0][0][0];
+        //single_data_=&B[0][0][0][0];
 
         //initialize bicgstab solver from Dune
-        solver_   = std::make_shared<Dune::BiCGSTABSolver<X>>(*op,
+        solver_ = std::make_shared<Dune::BiCGSTABSolver<Vector>>(
+                                                              //*op,
+                                                              *mixed_operator_,
                                                               *sp,
                                                               *prec,
                                                               tol, // desired residual reduction factor
@@ -38,13 +115,40 @@ class MixedAdapter : public InverseOperator<X,X>
                                                               verbosity);
     }
 
-    virtual void apply (X& x, X& b, InverseOperatorResult& res) override
+    virtual void apply(Vector &x, Vector &b, InverseOperatorResult &res) override
     {
+        //demote jacobian to single precision
+        auto &A = double_operator_->getmat();
+        auto &B = *single_matrix_;
+
+        double const *dbl = &A[0][0][0][0];
+        float *flt        = &B[0][0][0][0];
+        int N = A.nonzeroes()*block_size*block_size;
+        for(int k=0;k<N;k++) flt[k] = dbl[k];
+/*
+        //demote jacobian to single precision
+        auto &A = double_operator_->getmat();
+        auto &B = *single_matrix_;
+
+        int irow=0;
+        int icol=0;
+        for(auto row=A.begin();row!=A.end();row++)
+        {
+            for(unsigned int i=0;i<row->getsize();i++)
+            {
+                icol = row->getindexptr()[i];
+                for(int ii=0;ii<block_size;ii++)
+                for(int jj=0;jj<block_size;jj++)
+                B[irow][icol][ii][jj]=A[irow][icol][ii][jj];
+            }
+            irow++;
+        }
+*/
         //apply bicgstab solver from Dune
         solver_->apply(x,b,res);
     }
 
-    virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res) override
+    virtual void apply(Vector &x, Vector &b, double reduction, InverseOperatorResult &res) override
     {
         x=0;
         b=0;
@@ -52,14 +156,22 @@ class MixedAdapter : public InverseOperator<X,X>
         OPM_THROW(std::invalid_argument, "MixedAdapter::apply(...) not implemented yet.");
     }
 
+    virtual Dune::SolverCategory::Category category() const override{return Dune::SolverCategory::overlapping;};
 
     private:
-    using AbstractSolverType = Dune::InverseOperator<X, X>;
+    using AbstractSolverType = Dune::InverseOperator<Vector,Vector>;
+
 
     //Operator* operator_;
     //std::shared_ptr<AbstractPrecondType> prec_;
     //std::shared_ptr<AbstractScalarProductType> product_;
     std::shared_ptr<AbstractSolverType> solver_;
+    std::shared_ptr<MixedOperatorType> mixed_operator_;
+    std::shared_ptr<SingleMatrixType> single_matrix_;
+    Operator *double_operator_;
+    //SingleMatrixType *single_matrix_;
+    double const *double_data_;
+    //float *single_data_;
 
 };
 

From e2898949b2d6e4908566be003facd9c1be7f956a Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Mon, 4 May 2026 09:58:59 -0500
Subject: [PATCH 04/19] mixed: removing stale code

---
 opm/simulators/linalg/mixed/adapter.hpp | 49 +------------------------
 1 file changed, 2 insertions(+), 47 deletions(-)

diff --git a/opm/simulators/linalg/mixed/adapter.hpp b/opm/simulators/linalg/mixed/adapter.hpp
index b955a888409..399e71a629c 100644
--- a/opm/simulators/linalg/mixed/adapter.hpp
+++ b/opm/simulators/linalg/mixed/adapter.hpp
@@ -15,27 +15,6 @@
 
 #include <dune/istl/matrixindexset.hh>
 
-/*
-// Type definitions
-typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>> BCRSMat;
-typedef Dune::BlockVector<Dune::FieldVector<double,1>> BVector;
-typedef Dune::OwnerOverlapCopyCommunication<int,int> Communication;
-
-//Assume parallel_A is the assembled BCRSMatrix and comm is the communicator
-//1. Create the parallel operator
-typedef Dune::OverlappingSchwarzOperator<BCRSMat,BVector,BVector,Communication> Operator;
-Operator op(parallel_A,comm);
-
-//2. Create the scalar product
-typedef Dune::OverlappingSchwarzScalarProduct<BVector,Communication> ScalarProduct;
-ScalarProduct sp(comm);
-
-//3. Create the parallel preconditioner (e.g., usingSSOR)
-typedef Dune::SeqSSOR<BCRSMat,BVector,BVector> Smoother;
-typedef Dune::BlockPreconditioner<BVector,BVector,Communication,Smoother> ParPrec;
-Smoother smoother(parallel_A,1,1.0);
-ParPrec pprec(smoother,comm);
-*/
 
 namespace Dune
 {
@@ -52,11 +31,8 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
     static constexpr auto block_size = domain_type::block_type::dimension;
 
     using SingleMatrixType  = Dune::BCRSMatrix<Opm::MatrixBlock<float, block_size, block_size>>;
-#if HAVE_MPI
     using MixedOperatorType = Dune::OverlappingSchwarzOperator<SingleMatrixType, Vector, Vector, Comm>;
-#else // HAVE_MPI
-    using MixedOperatorType = Dune::MatrixAdapter<SingleMatrixType, Vector, Vector>;
-#endif // HAVE_MPI
+    //using MixedOperatorType = Dune::MatrixAdapter<SingleMatrixType, Vector, Vector>;
 
     MixedAdapter(Operator *op,
                  std::shared_ptr<AbstractScalarProductType> sp,
@@ -94,11 +70,8 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
         //initializemixedoperator
         double_operator_ = op;
 
-#if HAVE_MPI
         mixed_operator_ = std::make_shared<MixedOperatorType>(B,comm);
-#else // HAVE_MPI
-        mixed_operator_ = std::make_shared<MixedOperatorType>(B);
-#endif // HAVE_MPI
+        //mixed_operator_ = std::make_shared<MixedOperatorType>(B);
 
         //pointerstodataarrays
         //double_data_=&A[0][0][0][0];
@@ -125,25 +98,7 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
         float *flt        = &B[0][0][0][0];
         int N = A.nonzeroes()*block_size*block_size;
         for(int k=0;k<N;k++) flt[k] = dbl[k];
-/*
-        //demote jacobian to single precision
-        auto &A = double_operator_->getmat();
-        auto &B = *single_matrix_;
 
-        int irow=0;
-        int icol=0;
-        for(auto row=A.begin();row!=A.end();row++)
-        {
-            for(unsigned int i=0;i<row->getsize();i++)
-            {
-                icol = row->getindexptr()[i];
-                for(int ii=0;ii<block_size;ii++)
-                for(int jj=0;jj<block_size;jj++)
-                B[irow][icol][ii][jj]=A[irow][icol][ii][jj];
-            }
-            irow++;
-        }
-*/
         //apply bicgstab solver from Dune
         solver_->apply(x,b,res);
     }

From 6a5bd2805f27bc0b8ab79771eb972a4ffe150fb4 Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Mon, 4 May 2026 13:26:03 -0500
Subject: [PATCH 05/19] mixed: wrap mixed bsr matrix to interface with operator

---
 opm/simulators/linalg/mixed/MatrixWrapper.hpp | 73 +++++++++++++++++++
 opm/simulators/linalg/mixed/adapter.hpp       | 51 ++++++-------
 opm/simulators/linalg/mixed/bsr.c             | 40 ++++++++++
 opm/simulators/linalg/mixed/bsr.h             | 13 ++++
 4 files changed, 147 insertions(+), 30 deletions(-)
 create mode 100644 opm/simulators/linalg/mixed/MatrixWrapper.hpp

diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp
new file mode 100644
index 00000000000..08befeaa755
--- /dev/null
+++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp
@@ -0,0 +1,73 @@
+#ifndef OPM_MIXED_MATRIX_HEADER_INCLUDED
+#define OPM_MIXED_MATRIX_HEADER_INCLUDED
+
+
+#include <opm/simulators/linalg/mixed/bsr.h>
+
+
+template <class Vector, int b>
+class MatrixWrapper
+{
+    public:
+    virtual void mv(const Vector& x, Vector& y) const;
+    virtual void umv(const Vector& x, Vector& y) const;
+    virtual void usmv(double alpha, const Vector& x, Vector& y) const;
+
+
+    MatrixWrapper(int nrows, int nnz)
+    {
+        nnz_=nnz;
+        M_ = bsr_alloc();
+        bsr_init(M_, nrows, nnz, b);
+
+    }
+
+    void update(double const *data);
+
+    int *rowptr(){return M_->rowptr;}
+    int *colidx(){return M_->colidx;}
+
+    private:
+    int nnz_;
+    bsr_matrix  *M_;
+};
+
+template <class Vector, int b>
+void MatrixWrapper<Vector,b>::
+mv(const Vector& x, Vector& y) const
+{
+    bsr_vmspmv3(M_, &x[0][0], &y[0][0]);
+}
+
+template <class Vector, int b>
+void MatrixWrapper<Vector,b>::
+umv(const Vector& x, Vector& y) const
+{
+    bsr_vmspumv3(M_, &x[0][0], &y[0][0], 1.0);
+}
+
+template <class Vector, int b>
+void MatrixWrapper<Vector,b>::
+usmv(double alpha, const Vector& x, Vector& y) const
+{
+    bsr_vmspumv3(M_, &x[0][0], &y[0][0], alpha);
+}
+
+template <class Vector, int b>
+void MatrixWrapper<Vector,b>::
+update(double const *data)
+{
+    // transpose each dense block to make them column-major
+    int bb=b*b;
+    double B[bb];
+    for(int k=0;k<nnz_;k++)
+    {
+        for(int i=0;i<b;i++) for(int j=0;j<b;j++) B[b*j+i] = data[bb*k + b*i + j];
+        for(int i=0;i<bb;i++) M_->dbl[bb*k + i] = B[i];
+    }
+
+    // downcast to allow mixed precision
+    bsr_downcast(M_);
+}
+
+#endif // OPM_MIXED_MATRIX_HEADER_INCLUDED
diff --git a/opm/simulators/linalg/mixed/adapter.hpp b/opm/simulators/linalg/mixed/adapter.hpp
index 399e71a629c..735747373c7 100644
--- a/opm/simulators/linalg/mixed/adapter.hpp
+++ b/opm/simulators/linalg/mixed/adapter.hpp
@@ -15,6 +15,7 @@
 
 #include <dune/istl/matrixindexset.hh>
 
+#include <opm/simulators/linalg/mixed/MatrixWrapper.hpp>
 
 namespace Dune
 {
@@ -30,7 +31,8 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
     using typename InverseOperator<Vector, Vector>::domain_type;
     static constexpr auto block_size = domain_type::block_type::dimension;
 
-    using SingleMatrixType  = Dune::BCRSMatrix<Opm::MatrixBlock<float, block_size, block_size>>;
+    //using SingleMatrixType  = Dune::BCRSMatrix<Opm::MatrixBlock<float, block_size, block_size>>;
+    using SingleMatrixType  = MatrixWrapper<Vector, block_size>;
     using MixedOperatorType = Dune::OverlappingSchwarzOperator<SingleMatrixType, Vector, Vector, Comm>;
     //using MixedOperatorType = Dune::MatrixAdapter<SingleMatrixType, Vector, Vector>;
 
@@ -45,27 +47,29 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
 
         auto &A = op->getmat();
         int nrows = A.N();
-        //int nnz = A.nonzeroes();
+        int nnz = A.nonzeroes();
         //int b = A[0][0].N();
 
-        single_matrix_ = std::make_shared<SingleMatrixType>();
+        double_data_ = &A[0][0][0][0];
+        single_matrix_ = std::make_shared<SingleMatrixType>(nrows,nnz);
         auto &B = *single_matrix_;
 
-        //copy sparsity pattern
-        int irow=0;
-        int icol=0;
-        MatrixIndexSet sparsity(nrows,nrows);
-        for(auto row=A.begin();row!=A.end();row++)
+        // copy sparsity pattern
+        int *rows = single_matrix_->rowptr();
+        int *cols = single_matrix_->colidx();
+
+        int irow = 0;
+        int icol = 0;
+        rows[0]  = 0;
+        for(auto row=A.begin(); row!=A.end(); row++)
         {
-            for(unsigned int i=0;i<row->getsize();i++)
+            for(unsigned int i=0; i<row->getsize(); i++)
             {
-                icol = row->getindexptr()[i];
-                sparsity.add(irow,icol);
+                cols[icol++] = row->getindexptr()[i];
             }
+            rows[irow+1]     = rows[irow]+row->getsize();
             irow++;
         }
-        sparsity.exportIdx(B);
-        //B.compress();
 
         //initializemixedoperator
         double_operator_ = op;
@@ -73,10 +77,6 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
         mixed_operator_ = std::make_shared<MixedOperatorType>(B,comm);
         //mixed_operator_ = std::make_shared<MixedOperatorType>(B);
 
-        //pointerstodataarrays
-        //double_data_=&A[0][0][0][0];
-        //single_data_=&B[0][0][0][0];
-
         //initialize bicgstab solver from Dune
         solver_ = std::make_shared<Dune::BiCGSTABSolver<Vector>>(
                                                               //*op,
@@ -91,13 +91,9 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
     virtual void apply(Vector &x, Vector &b, InverseOperatorResult &res) override
     {
         //demote jacobian to single precision
-        auto &A = double_operator_->getmat();
-        auto &B = *single_matrix_;
-
-        double const *dbl = &A[0][0][0][0];
-        float *flt        = &B[0][0][0][0];
-        int N = A.nonzeroes()*block_size*block_size;
-        for(int k=0;k<N;k++) flt[k] = dbl[k];
+        //auto &A = double_operator_->getmat();
+        single_matrix_->update(double_data_);
+        //single_matrix_->update(&A[0][0][0][0]);
 
         //apply bicgstab solver from Dune
         solver_->apply(x,b,res);
@@ -117,16 +113,11 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
     using AbstractSolverType = Dune::InverseOperator<Vector,Vector>;
 
 
-    //Operator* operator_;
-    //std::shared_ptr<AbstractPrecondType> prec_;
-    //std::shared_ptr<AbstractScalarProductType> product_;
+    Operator *double_operator_;
     std::shared_ptr<AbstractSolverType> solver_;
     std::shared_ptr<MixedOperatorType> mixed_operator_;
     std::shared_ptr<SingleMatrixType> single_matrix_;
-    Operator *double_operator_;
-    //SingleMatrixType *single_matrix_;
     double const *double_data_;
-    //float *single_data_;
 
 };
 
diff --git a/opm/simulators/linalg/mixed/bsr.c b/opm/simulators/linalg/mixed/bsr.c
index 4900671f7fd..19df0d9163a 100644
--- a/opm/simulators/linalg/mixed/bsr.c
+++ b/opm/simulators/linalg/mixed/bsr.c
@@ -107,6 +107,46 @@ void bsr_vmspmv3(bsr_matrix *A, const double *x, double *y)
     }
 }
 
+void bsr_vmspumv3(bsr_matrix *A, const double *x, double *y, double alpha)
+{
+    int nrows = A->nrows;
+    int *rowptr=A->rowptr;
+    int *colidx=A->colidx;
+    const float *data=A->flt;
+
+    const int b=3;
+
+    __m256d valpha = _mm256_set1_pd(alpha);
+
+    __m256d mm_zeros =_mm256_setzero_pd();
+    for(int i=0;i<nrows;i++)
+    {
+        __m256d vA[3];
+        for(int k=0;k<3;k++) vA[k] = mm_zeros;
+        for(int k=rowptr[i];k<rowptr[i+1];k++)
+        {
+            const float *AA=data+9*k;
+
+            int j = colidx[k];
+            __m256d vx = _mm256_loadu_pd(x+b*j);
+
+            vA[0] += _mm256_cvtps_pd(_mm_loadu_ps(AA+0))*_mm256_permute4x64_pd(vx,0b00000000); //0b01010101
+            vA[1] += _mm256_cvtps_pd(_mm_loadu_ps(AA+3))*_mm256_permute4x64_pd(vx,0b01010101); //0b01010101
+            vA[2] += _mm256_cvtps_pd(_mm_loadu_ps(AA+6))*_mm256_permute4x64_pd(vx,0b10101010); //0b01010101
+        }
+
+        // sum over columns
+        __m256d vy, vz;
+        vz = vA[0] + vA[1] + vA[2];
+
+        double *y_i = y+b*i;
+        vy = _mm256_loadu_pd(y_i);       // optional blend to keep
+        vz = vy + valpha*vz;
+        vz =_mm256_blend_pd(vy,vz,0x7);  // 4th element unchanged
+        _mm256_storeu_pd(y_i,vz);
+    }
+}
+
 void bsr_vdspmv3(bsr_matrix *A, const double *x, double *y)
 {
     int nrows = A->nrows;
diff --git a/opm/simulators/linalg/mixed/bsr.h b/opm/simulators/linalg/mixed/bsr.h
index 11e72787eaf..f093ab5f557 100644
--- a/opm/simulators/linalg/mixed/bsr.h
+++ b/opm/simulators/linalg/mixed/bsr.h
@@ -64,6 +64,19 @@ void bsr_init(bsr_matrix *A, int nrows, int nnz, int b);
  * @param x Pointer to input vector.
  * @param y Pointer to output vector.
  */
+void bsr_vmspumv3(bsr_matrix *A, const double *x, double *y, double alpha);
+
+/**
+ * @brief Sparse matrix-vector multiplication in mixed precision.
+ *
+ * @note Function is specialized for 3x3 block-sparse matrices.
+ * @note Function uses AVX2 intrinsics.
+ *
+ * @param A Pointer to bsr matrix.
+ * @param x Pointer to input vector.
+ * @param y Pointer to output vector.
+ */
+
 void bsr_vmspmv3(bsr_matrix *A, const double *x, double *y);
 
 /**

From 901c794c6027fe0b3886a7e3ba32d5c3b57617c7 Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Mon, 4 May 2026 17:35:16 -0500
Subject: [PATCH 06/19] mixed: unify parallel and sequential

---
 opm/simulators/linalg/mixed/adapter.hpp | 33 ++++++++++++++++++++++---
 1 file changed, 30 insertions(+), 3 deletions(-)

diff --git a/opm/simulators/linalg/mixed/adapter.hpp b/opm/simulators/linalg/mixed/adapter.hpp
index 735747373c7..dad995ef670 100644
--- a/opm/simulators/linalg/mixed/adapter.hpp
+++ b/opm/simulators/linalg/mixed/adapter.hpp
@@ -17,8 +17,26 @@
 
 #include <opm/simulators/linalg/mixed/MatrixWrapper.hpp>
 
+
+
 namespace Dune
 {
+
+
+template <class Matrix, class Vector, class Comm>
+struct MixedOperator
+{
+    using type = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Comm>;
+    //static std::shared_ptr<type> mixed_operator(Matrix &M, Comm comm) {return std:make_shared<type>(M,comm);}
+};
+
+template <class Matrix, class Vector>
+struct MixedOperator<Matrix, Vector, Dune::Amg::SequentialInformation>
+{
+    using type = Dune::MatrixAdapter<Matrix, Vector, Vector>;
+};
+
+
 template <class Comm, class Operator, class Vector>
 class MixedAdapter:public InverseOperator<Vector, Vector>
 {
@@ -33,9 +51,12 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
 
     //using SingleMatrixType  = Dune::BCRSMatrix<Opm::MatrixBlock<float, block_size, block_size>>;
     using SingleMatrixType  = MatrixWrapper<Vector, block_size>;
-    using MixedOperatorType = Dune::OverlappingSchwarzOperator<SingleMatrixType, Vector, Vector, Comm>;
+    //using MixedOperatorType = Dune::OverlappingSchwarzOperator<SingleMatrixType, Vector, Vector, Comm>;
+    using MixedOperatorType = MixedOperator<SingleMatrixType, Vector, Comm>::type;
     //using MixedOperatorType = Dune::MatrixAdapter<SingleMatrixType, Vector, Vector>;
 
+
+
     MixedAdapter(Operator *op,
                  std::shared_ptr<AbstractScalarProductType> sp,
                  std::shared_ptr<AbstractPrecondType> prec,
@@ -73,9 +94,15 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
 
         //initializemixedoperator
         double_operator_ = op;
+        if constexpr (std::is_same_v<Comm, Dune::Amg::SequentialInformation>)
+        {
+            mixed_operator_ = std::make_shared<MixedOperatorType>(B);
+        }
+        else
+        {
+            mixed_operator_ = std::make_shared<MixedOperatorType>(B,comm);
+        }
 
-        mixed_operator_ = std::make_shared<MixedOperatorType>(B,comm);
-        //mixed_operator_ = std::make_shared<MixedOperatorType>(B);
 
         //initialize bicgstab solver from Dune
         solver_ = std::make_shared<Dune::BiCGSTABSolver<Vector>>(

From dfe5a79faef53d0bb269b47a9df5687921939e48 Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Tue, 5 May 2026 12:38:16 -0500
Subject: [PATCH 07/19] mixed: wrap mixed ilu0 preconditioner

---
 .../linalg/StandardPreconditioners_mpi.hpp    |  6 ++
 .../linalg/StandardPreconditioners_serial.hpp |  8 ++
 .../linalg/mixed/PreconditionerWrapper.hpp    | 89 +++++++++++++++++++
 opm/simulators/linalg/setupPropertyTree.cpp   | 19 ++++
 opm/simulators/linalg/setupPropertyTree.hpp   |  1 +
 5 files changed, 123 insertions(+)
 create mode 100644 opm/simulators/linalg/mixed/PreconditionerWrapper.hpp

diff --git a/opm/simulators/linalg/StandardPreconditioners_mpi.hpp b/opm/simulators/linalg/StandardPreconditioners_mpi.hpp
index d2dcf9887e1..266229c3a4b 100644
--- a/opm/simulators/linalg/StandardPreconditioners_mpi.hpp
+++ b/opm/simulators/linalg/StandardPreconditioners_mpi.hpp
@@ -33,6 +33,8 @@
 #include <memory>
 #include <type_traits>
 
+#include <opm/simulators/linalg/mixed/PreconditionerWrapper.hpp>
+
 namespace Opm {
 
 
@@ -163,6 +165,10 @@ struct StandardPreconditioners
             DUNE_UNUSED_PARAMETER(prm);
             return wrapBlockPreconditioner<MultithreadDILU<M, V, V>>(comm, op.getmat());
         });
+        F::addCreator("haugen-ilu0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
+            DUNE_UNUSED_PARAMETER(prm);
+            return wrapBlockPreconditioner<MixedPreconditioner<M,V,V>>(comm, op.getmat());
+        });
         F::addCreator("jac", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
             const int n = prm.get<int>("repeats", 1);
             const double w = prm.get<double>("relaxation", 1.0);
diff --git a/opm/simulators/linalg/StandardPreconditioners_serial.hpp b/opm/simulators/linalg/StandardPreconditioners_serial.hpp
index e18c007ef56..38aacf7e9cf 100644
--- a/opm/simulators/linalg/StandardPreconditioners_serial.hpp
+++ b/opm/simulators/linalg/StandardPreconditioners_serial.hpp
@@ -33,6 +33,10 @@
 #include <memory>
 #include <type_traits>
 
+
+#include <opm/simulators/linalg/mixed/PreconditionerWrapper.hpp>
+
+
 namespace Opm {
 
 template <class X, class Y>
@@ -88,6 +92,10 @@ struct StandardPreconditioners<Operator, Dune::Amg::SequentialInformation, typen
             DUNE_UNUSED_PARAMETER(prm);
             return std::make_shared<MultithreadDILU<M, V, V>>(op.getmat());
         });
+        F::addCreator("haugen-ilu0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
+            DUNE_UNUSED_PARAMETER(prm);
+            return std::make_shared<MixedPreconditioner<M,V,V>>(op.getmat());
+        });
         F::addCreator("mixed-ilu0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
             DUNE_UNUSED_PARAMETER(prm);
             DUNE_UNUSED_PARAMETER(op);
diff --git a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp
new file mode 100644
index 00000000000..a37dff04460
--- /dev/null
+++ b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp
@@ -0,0 +1,89 @@
+#ifndef OPM_MIXED_PREC_HEADER_INCLUDED
+#define OPM_MIXED_PREC_HEADER_INCLUDED
+
+#include <opm/simulators/linalg/mixed/prec.h>
+
+namespace Opm {
+
+template <class M, class X, class Y>
+class MixedPreconditioner : public Dune::PreconditionerWithUpdate<X, Y>
+{
+    public:
+
+    using domain_type = X;
+    static constexpr auto block_size = domain_type::block_type::dimension;
+
+    MixedPreconditioner(const M& A)
+    {
+        int nrows = A.N();
+        int nnz = A.nonzeroes();
+
+        double_data_  = &A[0][0][0][0];
+        prec_         = prec_alloc();
+        mixed_matrix_ = bsr_alloc();
+        bsr_init(mixed_matrix_, nrows, nnz, block_size);
+
+        // copy sparsity pattern
+        int *rows = mixed_matrix_->rowptr;
+        int *cols = mixed_matrix_->colidx;
+
+        int irow = 0;
+        int icol = 0;
+        rows[0]  = 0;
+        for(auto row=A.begin(); row!=A.end(); row++)
+        {
+            for(unsigned int i=0; i<row->getsize(); i++)
+            {
+                cols[icol++] = row->getindexptr()[i];
+            }
+            rows[irow+1]     = rows[irow]+row->getsize();
+            irow++;
+        }
+
+        prec_init(prec_,mixed_matrix_);
+        nnz_ = nnz;
+
+        update();
+    };
+    virtual void update() override;
+    virtual bool hasPerfectUpdate() const override {return true;}
+    virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& y) override {};
+    virtual void post ([[maybe_unused]] X& x) override {};
+    virtual void apply ([[maybe_unused]] X& x, [[maybe_unused]] const Y& y) override;
+    virtual Dune::SolverCategory::Category category() const override { return Dune::SolverCategory::sequential; };
+
+    private:
+    double const *double_data_;
+    bsr_matrix *mixed_matrix_;
+    prec_t *prec_;
+    int nnz_;
+};
+
+template <class M, class X, class Y>
+void MixedPreconditioner<M,X,Y>::
+update ()
+{
+    // transpose each dense block to make them column-major
+    int b = block_size;
+    int bb=b*b;
+    double B[bb];
+    for(int k=0;k<nnz_;k++)
+    {
+        for(int i=0;i<b;i++) for(int j=0;j<b;j++) B[b*j+i] = double_data_[bb*k + b*i + j];
+        for(int i=0;i<bb;i++) mixed_matrix_->dbl[bb*k + i] = B[i];
+    }
+
+    prec_ilu0_factorize(prec_, mixed_matrix_);
+    prec_downcast(prec_);
+}
+
+template <class M, class X, class Y>
+void MixedPreconditioner<M,X,Y>::
+apply ([[maybe_unused]] X& x, [[maybe_unused]] const Y& y)
+{
+    x=y;
+    prec_mapply3c(prec_,&x[0][0]);
+}
+
+}
+#endif // OPM_MIXED_PREC_HEADER_INCLUDED
diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp
index 236ca0528ae..c38b0f83669 100644
--- a/opm/simulators/linalg/setupPropertyTree.cpp
+++ b/opm/simulators/linalg/setupPropertyTree.cpp
@@ -238,6 +238,11 @@ setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters
         return setupILU(conf, p);
     }
 
+    // mixed-precision ILU0
+    if (conf == "haugen-ilu0") {
+        return setupHaugenILU(conf, p);
+    }
+
     // mixed-precision ILU0
     if (conf == "mixed-ilu0") {
         return setupMixedILU(conf, p);
@@ -409,6 +414,20 @@ setupILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParamet
     return prm;
 }
 
+PropertyTree
+setupHaugenILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
+{
+    using namespace std::string_literals;
+    PropertyTree prm;
+    prm.put("tol", p.linear_solver_reduction_);
+    prm.put("maxiter", p.linear_solver_maxiter_);
+    prm.put("verbosity", p.linear_solver_verbosity_);
+    prm.put("solver", "mixed-precision"s);
+    prm.put("preconditioner.type", "haugen-ilu0"s);
+    return prm;
+}
+
+
 PropertyTree
 setupMixedILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
 {
diff --git a/opm/simulators/linalg/setupPropertyTree.hpp b/opm/simulators/linalg/setupPropertyTree.hpp
index a424a573259..4c0edc32bf9 100644
--- a/opm/simulators/linalg/setupPropertyTree.hpp
+++ b/opm/simulators/linalg/setupPropertyTree.hpp
@@ -39,6 +39,7 @@ PropertyTree setupCPR(const std::string& conf, const FlowLinearSolverParameters&
 PropertyTree setupAMG(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupILU(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupMixedILU(const std::string& conf, const FlowLinearSolverParameters& p);
+PropertyTree setupHaugenILU(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupDILU(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupMixedDILU(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupUMFPack(const std::string& conf, const FlowLinearSolverParameters& p);

From b8c1e6838afb955ec03533e0f059dc993b968eda Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Tue, 5 May 2026 17:06:41 -0500
Subject: [PATCH 08/19] mixed: add mixed dilu option

---
 .../linalg/StandardPreconditioners_mpi.hpp     |  4 ++++
 .../linalg/StandardPreconditioners_serial.hpp  |  4 ++++
 .../linalg/mixed/PreconditionerWrapper.hpp     |  6 ++++--
 opm/simulators/linalg/setupPropertyTree.cpp    | 18 ++++++++++++++++++
 opm/simulators/linalg/setupPropertyTree.hpp    |  1 +
 5 files changed, 31 insertions(+), 2 deletions(-)

diff --git a/opm/simulators/linalg/StandardPreconditioners_mpi.hpp b/opm/simulators/linalg/StandardPreconditioners_mpi.hpp
index 266229c3a4b..3b775746710 100644
--- a/opm/simulators/linalg/StandardPreconditioners_mpi.hpp
+++ b/opm/simulators/linalg/StandardPreconditioners_mpi.hpp
@@ -169,6 +169,10 @@ struct StandardPreconditioners
             DUNE_UNUSED_PARAMETER(prm);
             return wrapBlockPreconditioner<MixedPreconditioner<M,V,V>>(comm, op.getmat());
         });
+        F::addCreator("haugen-dilu", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
+            DUNE_UNUSED_PARAMETER(prm);
+            return wrapBlockPreconditioner<MixedPreconditioner<M,V,V>>(comm, op.getmat(), true);
+        });
         F::addCreator("jac", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
             const int n = prm.get<int>("repeats", 1);
             const double w = prm.get<double>("relaxation", 1.0);
diff --git a/opm/simulators/linalg/StandardPreconditioners_serial.hpp b/opm/simulators/linalg/StandardPreconditioners_serial.hpp
index 38aacf7e9cf..c93f46eac65 100644
--- a/opm/simulators/linalg/StandardPreconditioners_serial.hpp
+++ b/opm/simulators/linalg/StandardPreconditioners_serial.hpp
@@ -96,6 +96,10 @@ struct StandardPreconditioners<Operator, Dune::Amg::SequentialInformation, typen
             DUNE_UNUSED_PARAMETER(prm);
             return std::make_shared<MixedPreconditioner<M,V,V>>(op.getmat());
         });
+        F::addCreator("haugen-dilu", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
+            DUNE_UNUSED_PARAMETER(prm);
+            return std::make_shared<MixedPreconditioner<M,V,V>>(op.getmat(),true);
+        });
         F::addCreator("mixed-ilu0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
             DUNE_UNUSED_PARAMETER(prm);
             DUNE_UNUSED_PARAMETER(op);
diff --git a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp
index a37dff04460..5c7744017a7 100644
--- a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp
+++ b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp
@@ -13,7 +13,7 @@ class MixedPreconditioner : public Dune::PreconditionerWithUpdate<X, Y>
     using domain_type = X;
     static constexpr auto block_size = domain_type::block_type::dimension;
 
-    MixedPreconditioner(const M& A)
+    MixedPreconditioner(const M& A, bool use_dilu = false)
     {
         int nrows = A.N();
         int nnz = A.nonzeroes();
@@ -42,6 +42,7 @@ class MixedPreconditioner : public Dune::PreconditionerWithUpdate<X, Y>
 
         prec_init(prec_,mixed_matrix_);
         nnz_ = nnz;
+        use_dilu_ = use_dilu;
 
         update();
     };
@@ -53,6 +54,7 @@ class MixedPreconditioner : public Dune::PreconditionerWithUpdate<X, Y>
     virtual Dune::SolverCategory::Category category() const override { return Dune::SolverCategory::sequential; };
 
     private:
+    bool use_dilu_;
     double const *double_data_;
     bsr_matrix *mixed_matrix_;
     prec_t *prec_;
@@ -73,7 +75,7 @@ update ()
         for(int i=0;i<bb;i++) mixed_matrix_->dbl[bb*k + i] = B[i];
     }
 
-    prec_ilu0_factorize(prec_, mixed_matrix_);
+    use_dilu_ ? prec_dilu_factorize(prec_, mixed_matrix_) : prec_ilu0_factorize(prec_, mixed_matrix_); // choose dilu or ilu0
     prec_downcast(prec_);
 }
 
diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp
index c38b0f83669..bce22e76729 100644
--- a/opm/simulators/linalg/setupPropertyTree.cpp
+++ b/opm/simulators/linalg/setupPropertyTree.cpp
@@ -243,6 +243,11 @@ setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters
         return setupHaugenILU(conf, p);
     }
 
+    // mixed-precision ILU0
+    if (conf == "haugen-dilu") {
+        return setupHaugenDILU(conf, p);
+    }
+
     // mixed-precision ILU0
     if (conf == "mixed-ilu0") {
         return setupMixedILU(conf, p);
@@ -427,6 +432,19 @@ setupHaugenILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverP
     return prm;
 }
 
+PropertyTree
+setupHaugenDILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
+{
+    using namespace std::string_literals;
+    PropertyTree prm;
+    prm.put("tol", p.linear_solver_reduction_);
+    prm.put("maxiter", p.linear_solver_maxiter_);
+    prm.put("verbosity", p.linear_solver_verbosity_);
+    prm.put("solver", "mixed-precision"s);
+    prm.put("preconditioner.type", "haugen-dilu"s);
+    return prm;
+}
+
 
 PropertyTree
 setupMixedILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
diff --git a/opm/simulators/linalg/setupPropertyTree.hpp b/opm/simulators/linalg/setupPropertyTree.hpp
index 4c0edc32bf9..ca43172cc99 100644
--- a/opm/simulators/linalg/setupPropertyTree.hpp
+++ b/opm/simulators/linalg/setupPropertyTree.hpp
@@ -42,6 +42,7 @@ PropertyTree setupMixedILU(const std::string& conf, const FlowLinearSolverParame
 PropertyTree setupHaugenILU(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupDILU(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupMixedDILU(const std::string& conf, const FlowLinearSolverParameters& p);
+PropertyTree setupHaugenDILU(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupUMFPack(const std::string& conf, const FlowLinearSolverParameters& p);
 
 } // namespace Opm

From 3a4eefac48c1cefd1863500b257413ef407c880e Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Thu, 7 May 2026 08:47:15 -0500
Subject: [PATCH 09/19] mixed: optimized sequential scalar product

---
 opm/simulators/linalg/mixed/adapter.hpp | 82 ++++++++++++++++---------
 1 file changed, 54 insertions(+), 28 deletions(-)

diff --git a/opm/simulators/linalg/mixed/adapter.hpp b/opm/simulators/linalg/mixed/adapter.hpp
index dad995ef670..d8ee07191f4 100644
--- a/opm/simulators/linalg/mixed/adapter.hpp
+++ b/opm/simulators/linalg/mixed/adapter.hpp
@@ -1,10 +1,7 @@
 #ifndef OPM_MIXED_ADAPTER_HEADER_INCLUDED
 #define OPM_MIXED_ADAPTER_HEADER_INCLUDED
 
-//#include <opm/simulators/flow/BlackoilModelParameters.hpp>
-//#include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
 #include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
-//#include <dune/istl/solver.hh>
 
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
@@ -15,19 +12,57 @@
 
 #include <dune/istl/matrixindexset.hh>
 
+#include <opm/simulators/linalg/mixed/Vector64b.hpp>
 #include <opm/simulators/linalg/mixed/MatrixWrapper.hpp>
 
 
 
 namespace Dune
 {
+#include <dune/istl/solvers.hh>
+
+
+template<class Vector>
+class SeqOptmizedProduct : public Dune::SeqScalarProduct<Vector>
+{
+public:
+
+    static constexpr auto block_size = Vector::block_type::dimension;
 
+    // Compute the dot product <x, y>
+    virtual double dot(const Vector& vx, const Vector& vy) const override
+    {
+        double const *x = &vx[0][0];
+        double const *y = &vy[0][0];
+
+        int NN = block_size*vx.N();
+
+        // unroll loop in multiples of 8
+        int n=NN/8;
+        int N=8*n;
+        double agg[8];
+        for(int i=0;i<8;i++) agg[i]=0.0;
+        for(int i=0;i<N;i+=8) for(int j=0;j<8;j++) agg[j]+=x[i+j]*y[i+j];
+        for(int j=0;j<4;j++) agg[j]+=agg[j+4];
+        for(int j=0;j<2;j++) agg[j]+=agg[j+2];
+        for(int j=0;j<1;j++) agg[j]+=agg[j+1];
+
+        // trailing end
+        for(int j=N;j<NN;j++) agg[0]+=x[j]*y[j];
+
+        return agg[0];
+    }
+
+    // Compute the norm ||x||
+    virtual double norm(const Vector& x) const override {
+        return std::sqrt(this->dot(x, x));
+    }
+};
 
 template <class Matrix, class Vector, class Comm>
 struct MixedOperator
 {
     using type = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Comm>;
-    //static std::shared_ptr<type> mixed_operator(Matrix &M, Comm comm) {return std:make_shared<type>(M,comm);}
 };
 
 template <class Matrix, class Vector>
@@ -40,22 +75,15 @@ struct MixedOperator<Matrix, Vector, Dune::Amg::SequentialInformation>
 template <class Comm, class Operator, class Vector>
 class MixedAdapter:public InverseOperator<Vector, Vector>
 {
-
     public:
 
     using AbstractPrecondType = Dune::PreconditionerWithUpdate<Vector,Vector>;
     using AbstractScalarProductType = Dune::ScalarProduct<Vector>;
 
-    using typename InverseOperator<Vector, Vector>::domain_type;
-    static constexpr auto block_size = domain_type::block_type::dimension;
-
-    //using SingleMatrixType  = Dune::BCRSMatrix<Opm::MatrixBlock<float, block_size, block_size>>;
-    using SingleMatrixType  = MatrixWrapper<Vector, block_size>;
-    //using MixedOperatorType = Dune::OverlappingSchwarzOperator<SingleMatrixType, Vector, Vector, Comm>;
-    using MixedOperatorType = MixedOperator<SingleMatrixType, Vector, Comm>::type;
-    //using MixedOperatorType = Dune::MatrixAdapter<SingleMatrixType, Vector, Vector>;
-
-
+    static constexpr auto block_size = Vector::block_type::dimension;
+    using MixedMatrixType     = MatrixWrapper<Vector, block_size>;
+    using MixedOperatorType    = MixedOperator<MixedMatrixType, Vector, Comm>::type;
+    using OptimizedProductType = SeqOptmizedProduct<Vector>;
 
     MixedAdapter(Operator *op,
                  std::shared_ptr<AbstractScalarProductType> sp,
@@ -69,15 +97,14 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
         auto &A = op->getmat();
         int nrows = A.N();
         int nnz = A.nonzeroes();
-        //int b = A[0][0].N();
 
         double_data_ = &A[0][0][0][0];
-        single_matrix_ = std::make_shared<SingleMatrixType>(nrows,nnz);
-        auto &B = *single_matrix_;
+        mixed_matrix_ = std::make_shared<MixedMatrixType>(nrows,nnz);
+        auto &B = *mixed_matrix_;
 
         // copy sparsity pattern
-        int *rows = single_matrix_->rowptr();
-        int *cols = single_matrix_->colidx();
+        int *rows = mixed_matrix_->rowptr();
+        int *cols = mixed_matrix_->colidx();
 
         int irow = 0;
         int icol = 0;
@@ -92,23 +119,23 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
             irow++;
         }
 
-        //initializemixedoperator
+        //initialize mixed operator
         double_operator_ = op;
         if constexpr (std::is_same_v<Comm, Dune::Amg::SequentialInformation>)
         {
             mixed_operator_ = std::make_shared<MixedOperatorType>(B);
+            scalar_product_ = std::make_shared<OptimizedProductType>();
         }
         else
         {
             mixed_operator_ = std::make_shared<MixedOperatorType>(B,comm);
+            scalar_product_ = sp;
         }
 
-
         //initialize bicgstab solver from Dune
         solver_ = std::make_shared<Dune::BiCGSTABSolver<Vector>>(
-                                                              //*op,
                                                               *mixed_operator_,
-                                                              *sp,
+                                                              *scalar_product_,
                                                               *prec,
                                                               tol, // desired residual reduction factor
                                                               maxiter, // maximum number of iterations
@@ -118,9 +145,7 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
     virtual void apply(Vector &x, Vector &b, InverseOperatorResult &res) override
     {
         //demote jacobian to single precision
-        //auto &A = double_operator_->getmat();
-        single_matrix_->update(double_data_);
-        //single_matrix_->update(&A[0][0][0][0]);
+        mixed_matrix_->update(double_data_);
 
         //apply bicgstab solver from Dune
         solver_->apply(x,b,res);
@@ -143,7 +168,8 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
     Operator *double_operator_;
     std::shared_ptr<AbstractSolverType> solver_;
     std::shared_ptr<MixedOperatorType> mixed_operator_;
-    std::shared_ptr<SingleMatrixType> single_matrix_;
+    std::shared_ptr<MixedMatrixType> mixed_matrix_;
+    std::shared_ptr<AbstractScalarProductType> scalar_product_;
     double const *double_data_;
 
 };

From 058bf202fec76cb69d21a8e50bff0af22d6178ae Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Fri, 8 May 2026 17:25:25 -0500
Subject: [PATCH 10/19] mixed: add destructors and documentation

---
 opm/simulators/linalg/mixed/MatrixWrapper.hpp | 18 +++++--
 .../linalg/mixed/PreconditionerWrapper.hpp    | 28 +++++++++-
 opm/simulators/linalg/mixed/adapter.hpp       | 52 +++++++++++++++----
 opm/simulators/linalg/setupPropertyTree.cpp   |  3 +-
 .../satfunc/OilPhaseConsistencyChecks.cpp     |  2 +-
 5 files changed, 84 insertions(+), 19 deletions(-)

diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp
index 08befeaa755..10e273a5bd5 100644
--- a/opm/simulators/linalg/mixed/MatrixWrapper.hpp
+++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp
@@ -5,6 +5,10 @@
 #include <opm/simulators/linalg/mixed/bsr.h>
 
 
+//! @brief Wraps c-implementation of mixed-precision matrix
+//!
+//! @tparam Vector the block-vector used by linear operator
+//! @tparam b block size
 template <class Vector, int b>
 class MatrixWrapper
 {
@@ -13,15 +17,20 @@ class MatrixWrapper
     virtual void umv(const Vector& x, Vector& y) const;
     virtual void usmv(double alpha, const Vector& x, Vector& y) const;
 
-
+    //! @brief constructor
+    //!
+    //! @param nrows number of block rows
+    //! @param nnz number of nonzero blocks
     MatrixWrapper(int nrows, int nnz)
     {
         nnz_=nnz;
         M_ = bsr_alloc();
         bsr_init(M_, nrows, nnz, b);
-
     }
 
+    //! @brief destructor
+    ~MatrixWrapper() {bsr_free(M_);}
+
     void update(double const *data);
 
     int *rowptr(){return M_->rowptr;}
@@ -36,6 +45,7 @@ template <class Vector, int b>
 void MatrixWrapper<Vector,b>::
 mv(const Vector& x, Vector& y) const
 {
+    // mixed-precision block spmv (y = M.x)
     bsr_vmspmv3(M_, &x[0][0], &y[0][0]);
 }
 
@@ -43,6 +53,7 @@ template <class Vector, int b>
 void MatrixWrapper<Vector,b>::
 umv(const Vector& x, Vector& y) const
 {
+    // mixed-precision block spmv with update (y += M.x)
     bsr_vmspumv3(M_, &x[0][0], &y[0][0], 1.0);
 }
 
@@ -50,6 +61,7 @@ template <class Vector, int b>
 void MatrixWrapper<Vector,b>::
 usmv(double alpha, const Vector& x, Vector& y) const
 {
+    // scaled mixed-precision block spmv with update (y += alpha *  M.x)
     bsr_vmspumv3(M_, &x[0][0], &y[0][0], alpha);
 }
 
@@ -66,7 +78,7 @@ update(double const *data)
         for(int i=0;i<bb;i++) M_->dbl[bb*k + i] = B[i];
     }
 
-    // downcast to allow mixed precision
+    // downcast to single precision
     bsr_downcast(M_);
 }
 
diff --git a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp
index 5c7744017a7..2f91f2bc1b1 100644
--- a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp
+++ b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp
@@ -5,6 +5,12 @@
 
 namespace Opm {
 
+//! @brief Wraps c-implementation of mixed-precision preconditioner
+//!
+//! @tparam Vector the block-vector used by linear operator
+//! @tparam M block-sparse matrix type
+//! @tparam X block-vector input type
+//! @tparam Y block vector output type
 template <class M, class X, class Y>
 class MixedPreconditioner : public Dune::PreconditionerWithUpdate<X, Y>
 {
@@ -13,17 +19,23 @@ class MixedPreconditioner : public Dune::PreconditionerWithUpdate<X, Y>
     using domain_type = X;
     static constexpr auto block_size = domain_type::block_type::dimension;
 
+    //! @brief constructor
+    //!
+    //! @param A block-sparse matrix
+    //! @param use_dilu toggle between dilu or ilu0 factorization
     MixedPreconditioner(const M& A, bool use_dilu = false)
     {
+        // Access double precision matrix data
         int nrows = A.N();
         int nnz = A.nonzeroes();
-
         double_data_  = &A[0][0][0][0];
         prec_         = prec_alloc();
+
+        // allocate and initialize mixed-precision matrix
         mixed_matrix_ = bsr_alloc();
         bsr_init(mixed_matrix_, nrows, nnz, block_size);
 
-        // copy sparsity pattern
+        // copy sparsity pattern from double preccision matrix
         int *rows = mixed_matrix_->rowptr;
         int *cols = mixed_matrix_->colidx;
 
@@ -40,12 +52,24 @@ class MixedPreconditioner : public Dune::PreconditionerWithUpdate<X, Y>
             irow++;
         }
 
+        // allocate and initialize preconditioner
         prec_init(prec_,mixed_matrix_);
+
+        // attribute initialization
         nnz_ = nnz;
         use_dilu_ = use_dilu;
 
+        // perform matrix factorization
         update();
     };
+
+    //! @brief destructor
+    ~MixedPreconditioner()
+    {
+        bsr_free(mixed_matrix_);
+        prec_free(prec_);
+    }
+
     virtual void update() override;
     virtual bool hasPerfectUpdate() const override {return true;}
     virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& y) override {};
diff --git a/opm/simulators/linalg/mixed/adapter.hpp b/opm/simulators/linalg/mixed/adapter.hpp
index d8ee07191f4..0b4380f3d93 100644
--- a/opm/simulators/linalg/mixed/adapter.hpp
+++ b/opm/simulators/linalg/mixed/adapter.hpp
@@ -12,29 +12,34 @@
 
 #include <dune/istl/matrixindexset.hh>
 
-#include <opm/simulators/linalg/mixed/Vector64b.hpp>
 #include <opm/simulators/linalg/mixed/MatrixWrapper.hpp>
 
 
 
 namespace Dune
 {
-#include <dune/istl/solvers.hh>
+//#include <dune/istl/solvers.hh>
 
 
+//! @brief Optimized sequential scalar product.
+//!
+//! @tparam Vector block-vector class with data stored as contiguous double array
 template<class Vector>
 class SeqOptmizedProduct : public Dune::SeqScalarProduct<Vector>
 {
 public:
 
+    // extract block size
     static constexpr auto block_size = Vector::block_type::dimension;
 
-    // Compute the dot product <x, y>
+    // Compute the dot product <vx, vy>
     virtual double dot(const Vector& vx, const Vector& vy) const override
     {
+        // access underlying data
         double const *x = &vx[0][0];
         double const *y = &vy[0][0];
 
+        // total array length
         int NN = block_size*vx.N();
 
         // unroll loop in multiples of 8
@@ -59,12 +64,21 @@ class SeqOptmizedProduct : public Dune::SeqScalarProduct<Vector>
     }
 };
 
+//! @brief Generalized mixed precision operator interface
+//!
+//! @tparam Matrix the block-matrix used by linear operator
+//! @tparam Vector the block-vector used by linear operator
+//! @tparam Comm the communicator used by linear operator
 template <class Matrix, class Vector, class Comm>
 struct MixedOperator
 {
     using type = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Comm>;
 };
 
+//! @brief Generalized mixed precision operator interface
+//!
+//! @tparam Matrix the block-matrix used by linear operator
+//! @tparam Vector the block-vector used by linear operator
 template <class Matrix, class Vector>
 struct MixedOperator<Matrix, Vector, Dune::Amg::SequentialInformation>
 {
@@ -72,6 +86,11 @@ struct MixedOperator<Matrix, Vector, Dune::Amg::SequentialInformation>
 };
 
 
+//! @brief Wraps mixed precision
+//!
+//! @tparam Comm the communicator passed to FlexibleLinearSolver
+//! @tparam Operator the linear operator passed to FlexibleLinearSolver
+//! @tparam Vector the block vector passed to FlexibleLinearSolver
 template <class Comm, class Operator, class Vector>
 class MixedAdapter:public InverseOperator<Vector, Vector>
 {
@@ -81,10 +100,19 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
     using AbstractScalarProductType = Dune::ScalarProduct<Vector>;
 
     static constexpr auto block_size = Vector::block_type::dimension;
-    using MixedMatrixType     = MatrixWrapper<Vector, block_size>;
+    using MixedMatrixType      = MatrixWrapper<Vector, block_size>;
     using MixedOperatorType    = MixedOperator<MixedMatrixType, Vector, Comm>::type;
     using OptimizedProductType = SeqOptmizedProduct<Vector>;
 
+    //! @brief constructor
+    //!
+    //! @param op the linear operator (assumed double precision)
+    //! @param sp the scalar product
+    //! @param prec the preconditioner to use
+    //! @param reduction the reduction factor passed to the iterative solver
+    //! @param maxit maximum number of iterations for the linear solver
+    //! @param verbose verbosity level
+    //! @param comm the communication object.
     MixedAdapter(Operator *op,
                  std::shared_ptr<AbstractScalarProductType> sp,
                  std::shared_ptr<AbstractPrecondType> prec,
@@ -94,15 +122,17 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
                  const Comm &comm)
     {
 
+        // Access matrix data from double precision operator
         auto &A = op->getmat();
         int nrows = A.N();
         int nnz = A.nonzeroes();
-
         double_data_ = &A[0][0][0][0];
+
+        //allocate mixed matrix
         mixed_matrix_ = std::make_shared<MixedMatrixType>(nrows,nnz);
-        auto &B = *mixed_matrix_;
+        //auto &B = *mixed_matrix_;
 
-        // copy sparsity pattern
+        // copy sparsity pattern from double precision matrix
         int *rows = mixed_matrix_->rowptr();
         int *cols = mixed_matrix_->colidx();
 
@@ -119,16 +149,16 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
             irow++;
         }
 
-        //initialize mixed operator
+        //initialize mixed operator and optimized scalar product
         double_operator_ = op;
         if constexpr (std::is_same_v<Comm, Dune::Amg::SequentialInformation>)
         {
-            mixed_operator_ = std::make_shared<MixedOperatorType>(B);
+            mixed_operator_ = std::make_shared<MixedOperatorType>(*mixed_matrix_);
             scalar_product_ = std::make_shared<OptimizedProductType>();
         }
         else
         {
-            mixed_operator_ = std::make_shared<MixedOperatorType>(B,comm);
+            mixed_operator_ = std::make_shared<MixedOperatorType>(*mixed_matrix_,comm);
             scalar_product_ = sp;
         }
 
@@ -144,7 +174,7 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
 
     virtual void apply(Vector &x, Vector &b, InverseOperatorResult &res) override
     {
-        //demote jacobian to single precision
+        //transpose dense blocks and demote to single precision
         mixed_matrix_->update(double_data_);
 
         //apply bicgstab solver from Dune
diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp
index bce22e76729..f2e0d3cee69 100644
--- a/opm/simulators/linalg/setupPropertyTree.cpp
+++ b/opm/simulators/linalg/setupPropertyTree.cpp
@@ -295,8 +295,7 @@ std::string getSolverString(const FlowLinearSolverParameters& p)
     }
     else
     {
-        //return {"bicgstab"};
-        return {"mixed-precision"};
+        return {"bicgstab"};
     }
 }
 
diff --git a/opm/simulators/utils/satfunc/OilPhaseConsistencyChecks.cpp b/opm/simulators/utils/satfunc/OilPhaseConsistencyChecks.cpp
index df45ad042a8..8a2ee6d65f1 100644
--- a/opm/simulators/utils/satfunc/OilPhaseConsistencyChecks.cpp
+++ b/opm/simulators/utils/satfunc/OilPhaseConsistencyChecks.cpp
@@ -31,7 +31,7 @@ template <typename Scalar>
 void Opm::Satfunc::PhaseChecks::Oil::SOcr_GO<Scalar>::
 testImpl(const EclEpsScalingPointsInfo<Scalar>& endPoints)
 {
-    this->sogcr_ = endPoints.Sogcr;
+    this->sogcr_ = endPoints.Sogcr; //0.0
 
     if (! std::isfinite(this->sogcr_)) {
         this->setViolated();

From 343b5374745b1c21891ff88b9498887b30350037 Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Tue, 12 May 2026 08:52:14 -0500
Subject: [PATCH 11/19] mixed: renaming aesthetics

---
 opm/simulators/linalg/FlexibleSolver_impl.hpp | 10 +++---
 .../linalg/StandardPreconditioners_mpi.hpp    |  4 +--
 .../linalg/StandardPreconditioners_serial.hpp |  8 ++---
 .../mixed/{adapter.hpp => SolverAdapter.hpp}  |  0
 opm/simulators/linalg/setupPropertyTree.cpp   | 32 +++++++++----------
 opm/simulators/linalg/setupPropertyTree.hpp   |  4 +--
 6 files changed, 28 insertions(+), 30 deletions(-)
 rename opm/simulators/linalg/mixed/{adapter.hpp => SolverAdapter.hpp} (100%)

diff --git a/opm/simulators/linalg/FlexibleSolver_impl.hpp b/opm/simulators/linalg/FlexibleSolver_impl.hpp
index c4217497b86..51a6cdc3238 100644
--- a/opm/simulators/linalg/FlexibleSolver_impl.hpp
+++ b/opm/simulators/linalg/FlexibleSolver_impl.hpp
@@ -32,9 +32,9 @@
 #include <opm/simulators/linalg/PreconditionerFactoryGPUIncludeWrapper.hpp>
 #include <opm/simulators/linalg/is_gpu_operator.hpp>
 
-#if 1 //HAVE_AVX2_EXTENSION
+#if HAVE_AVX2_EXTENSION
 #include <opm/simulators/linalg/mixed/wrapper.hpp>
-#include <opm/simulators/linalg/mixed/adapter.hpp>
+#include <opm/simulators/linalg/mixed/SolverAdapter.hpp>
 #endif
 
 #include <dune/common/fmatrix.hh>
@@ -186,7 +186,7 @@ namespace Dune
                                                                             tol, // desired residual reduction factor
                                                                             maxiter, // maximum number of iterations
                                                                             verbosity);
-#if 1 //HAVE_AVX2_EXTENSION
+#if HAVE_AVX2_EXTENSION
           } else if (solver_type == "mixed-bicgstab") {
               if constexpr (Opm::is_gpu_operator_v<Operator>) {
                 OPM_THROW(std::invalid_argument, "mixed-bicgstab solver not supported for GPU operatorsg");
@@ -194,7 +194,7 @@ namespace Dune
                 OPM_THROW(std::invalid_argument, "mixed-bicgstab solver not supported for single precision.");
             } else {
                 const std::string prec_type = prm.get<std::string>("preconditioner.type", "error");
-                bool use_mixed_dilu= (prec_type=="mixed-dilu");
+                bool use_mixed_dilu= (prec_type=="legacy-mixed-dilu");
                 using MatrixType = decltype(linearoperator_for_solver_->getmat());
                 linsolver_ = std::make_shared<Dune::MixedSolver<VectorType,MatrixType>>(
                                                                             linearoperator_for_solver_->getmat(),
@@ -212,8 +212,6 @@ namespace Dune
                                                                             maxiter, // maximum number of iterations
                                                                             verbosity,
                                                                             comm);
-
-
 #endif
         } else if (solver_type == "loopsolver") {
             linsolver_ = std::make_shared<Dune::LoopSolver<VectorType>>(*linearoperator_for_solver_,
diff --git a/opm/simulators/linalg/StandardPreconditioners_mpi.hpp b/opm/simulators/linalg/StandardPreconditioners_mpi.hpp
index 3b775746710..b71c0938a02 100644
--- a/opm/simulators/linalg/StandardPreconditioners_mpi.hpp
+++ b/opm/simulators/linalg/StandardPreconditioners_mpi.hpp
@@ -165,11 +165,11 @@ struct StandardPreconditioners
             DUNE_UNUSED_PARAMETER(prm);
             return wrapBlockPreconditioner<MultithreadDILU<M, V, V>>(comm, op.getmat());
         });
-        F::addCreator("haugen-ilu0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
+        F::addCreator("mixed-ilu0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
             DUNE_UNUSED_PARAMETER(prm);
             return wrapBlockPreconditioner<MixedPreconditioner<M,V,V>>(comm, op.getmat());
         });
-        F::addCreator("haugen-dilu", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
+        F::addCreator("mixed-dilu", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
             DUNE_UNUSED_PARAMETER(prm);
             return wrapBlockPreconditioner<MixedPreconditioner<M,V,V>>(comm, op.getmat(), true);
         });
diff --git a/opm/simulators/linalg/StandardPreconditioners_serial.hpp b/opm/simulators/linalg/StandardPreconditioners_serial.hpp
index c93f46eac65..d5940731888 100644
--- a/opm/simulators/linalg/StandardPreconditioners_serial.hpp
+++ b/opm/simulators/linalg/StandardPreconditioners_serial.hpp
@@ -92,20 +92,20 @@ struct StandardPreconditioners<Operator, Dune::Amg::SequentialInformation, typen
             DUNE_UNUSED_PARAMETER(prm);
             return std::make_shared<MultithreadDILU<M, V, V>>(op.getmat());
         });
-        F::addCreator("haugen-ilu0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
+        F::addCreator("mixed-ilu0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
             DUNE_UNUSED_PARAMETER(prm);
             return std::make_shared<MixedPreconditioner<M,V,V>>(op.getmat());
         });
-        F::addCreator("haugen-dilu", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
+        F::addCreator("mixed-dilu", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
             DUNE_UNUSED_PARAMETER(prm);
             return std::make_shared<MixedPreconditioner<M,V,V>>(op.getmat(),true);
         });
-        F::addCreator("mixed-ilu0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
+        F::addCreator("legacy-mixed-ilu0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
             DUNE_UNUSED_PARAMETER(prm);
             DUNE_UNUSED_PARAMETER(op);
             return std::make_shared<TrivialPreconditioner<V,V>>();
         });
-        F::addCreator("mixed-dilu", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
+        F::addCreator("legacy-mixed-dilu", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
             DUNE_UNUSED_PARAMETER(prm);
             DUNE_UNUSED_PARAMETER(op);
             return std::make_shared<TrivialPreconditioner<V,V>>();
diff --git a/opm/simulators/linalg/mixed/adapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp
similarity index 100%
rename from opm/simulators/linalg/mixed/adapter.hpp
rename to opm/simulators/linalg/mixed/SolverAdapter.hpp
diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp
index f2e0d3cee69..0af2e084501 100644
--- a/opm/simulators/linalg/setupPropertyTree.cpp
+++ b/opm/simulators/linalg/setupPropertyTree.cpp
@@ -239,23 +239,23 @@ setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters
     }
 
     // mixed-precision ILU0
-    if (conf == "haugen-ilu0") {
-        return setupHaugenILU(conf, p);
+    if (conf == "mixed-ilu0") {
+        return setupMixedILU(conf, p);
     }
 
     // mixed-precision ILU0
-    if (conf == "haugen-dilu") {
-        return setupHaugenDILU(conf, p);
+    if (conf == "mixed-dilu") {
+        return setupMixedDILU(conf, p);
     }
 
     // mixed-precision ILU0
-    if (conf == "mixed-ilu0") {
-        return setupMixedILU(conf, p);
+    if (conf == "legacy-mixed-ilu0") {
+        return setupLegacyMixedILU(conf, p);
     }
 
     // mixed-precision ILU0
-    if (conf == "mixed-dilu") {
-        return setupMixedDILU(conf, p);
+    if (conf == "legacy-mixed-dilu") {
+        return setupLegacyMixedDILU(conf, p);
     }
 
     if (conf == "dilu") {
@@ -419,7 +419,7 @@ setupILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParamet
 }
 
 PropertyTree
-setupHaugenILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
+setupMixedILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
 {
     using namespace std::string_literals;
     PropertyTree prm;
@@ -427,12 +427,12 @@ setupHaugenILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverP
     prm.put("maxiter", p.linear_solver_maxiter_);
     prm.put("verbosity", p.linear_solver_verbosity_);
     prm.put("solver", "mixed-precision"s);
-    prm.put("preconditioner.type", "haugen-ilu0"s);
+    prm.put("preconditioner.type", "mixed-ilu0"s);
     return prm;
 }
 
 PropertyTree
-setupHaugenDILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
+setupMixedDILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
 {
     using namespace std::string_literals;
     PropertyTree prm;
@@ -440,13 +440,13 @@ setupHaugenDILU([[maybe_unused]] const std::string& conf, const FlowLinearSolver
     prm.put("maxiter", p.linear_solver_maxiter_);
     prm.put("verbosity", p.linear_solver_verbosity_);
     prm.put("solver", "mixed-precision"s);
-    prm.put("preconditioner.type", "haugen-dilu"s);
+    prm.put("preconditioner.type", "mixed-dilu"s);
     return prm;
 }
 
 
 PropertyTree
-setupMixedILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
+setupLegacyMixedILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
 {
     using namespace std::string_literals;
     PropertyTree prm;
@@ -454,12 +454,12 @@ setupMixedILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverPa
     prm.put("maxiter", p.linear_solver_maxiter_);
     prm.put("verbosity", p.linear_solver_verbosity_);
     prm.put("solver", "mixed-bicgstab"s);
-    prm.put("preconditioner.type", "mixed-ilu0"s);
+    prm.put("preconditioner.type", "legacy-mixed-ilu0"s);
     return prm;
 }
 
 PropertyTree
-setupMixedDILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
+setupLegacyMixedDILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
 {
     using namespace std::string_literals;
     PropertyTree prm;
@@ -467,7 +467,7 @@ setupMixedDILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverP
     prm.put("maxiter", p.linear_solver_maxiter_);
     prm.put("verbosity", p.linear_solver_verbosity_);
     prm.put("solver", "mixed-bicgstab"s);
-    prm.put("preconditioner.type", "mixed-dilu"s);
+    prm.put("preconditioner.type", "legacy-mixed-dilu"s);
     return prm;
 }
 
diff --git a/opm/simulators/linalg/setupPropertyTree.hpp b/opm/simulators/linalg/setupPropertyTree.hpp
index ca43172cc99..0c5b206ce3e 100644
--- a/opm/simulators/linalg/setupPropertyTree.hpp
+++ b/opm/simulators/linalg/setupPropertyTree.hpp
@@ -38,11 +38,11 @@ PropertyTree setupCPRW(const std::string& conf, const FlowLinearSolverParameters
 PropertyTree setupCPR(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupAMG(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupILU(const std::string& conf, const FlowLinearSolverParameters& p);
+PropertyTree setupLegacyMixedILU(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupMixedILU(const std::string& conf, const FlowLinearSolverParameters& p);
-PropertyTree setupHaugenILU(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupDILU(const std::string& conf, const FlowLinearSolverParameters& p);
+PropertyTree setupLegacyMixedDILU(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupMixedDILU(const std::string& conf, const FlowLinearSolverParameters& p);
-PropertyTree setupHaugenDILU(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupUMFPack(const std::string& conf, const FlowLinearSolverParameters& p);
 
 } // namespace Opm

From 82e44bdf1b3acde6e5ca7876b40bd05c3a7fb743 Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Tue, 12 May 2026 09:05:10 -0500
Subject: [PATCH 12/19] mixed: update list of public header files

---
 CMakeLists_files.cmake | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 9a27bb9ee04..c26cec0887d 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -1301,6 +1301,9 @@ if (HAVE_AVX2_EXTENSION)
     opm/simulators/linalg/mixed/prec.h
     opm/simulators/linalg/mixed/vec.h
     opm/simulators/linalg/mixed/wrapper.hpp
+    opm/simulators/linalg/mixed/MatrixWrapper.hpp
+    opm/simulators/linalg/mixed/PreconditionerWrapper.hpp
+    opm/simulators/linalg/mixed/SolverAdapter.hpp
   )
 endif()
 

From 10123856b8643ac437078ea95747ecf37d49beb5 Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Tue, 12 May 2026 09:37:10 -0500
Subject: [PATCH 13/19] mixed: update README file

---
 opm/simulators/linalg/mixed/README.md           | 17 +++++++++++------
 .../utils/satfunc/OilPhaseConsistencyChecks.cpp |  2 +-
 2 files changed, 12 insertions(+), 7 deletions(-)

diff --git a/opm/simulators/linalg/mixed/README.md b/opm/simulators/linalg/mixed/README.md
index 1e27c55ea8d..44ae8572bb1 100644
--- a/opm/simulators/linalg/mixed/README.md
+++ b/opm/simulators/linalg/mixed/README.md
@@ -4,16 +4,18 @@ and a highly optimized mixed-precision implementation of ILU0 and DILU precondit
 bicgstab. Hopefully, this will inspire the exploration of mixed-precision algorithms
 in OPM.
 
-The initial implementation is specialized for 3x3 block-sparse matrices due to their
-importance in reservoir simulation and restricted to serial runs. Extending the work
-to parallel runs is a matter of extending OPM's parallel infrastructure and should be
-relatively straight-forward.
+The initial implementations are specialized for 3x3 block-sparse matrices due to their
+importance in reservoir simulation. The original implementation only works in serial.
+The parallel implementation wraps the optimized mixed-precision matrix-vector multiplication
+and ILU0/DILU preconditioners into building blocks for the ISTL-based implementation of
+bicgstab. This sacrifices some performance for improved modularity. Extending the work to
+block-sparse matrices of arbitrary block size is work in progress.
 
 The mixed-precision solver is selected by the command-line options `--linear-solver=mixed-ilu0`
 or `--linear-solver=mixed-dilu`. The command-line option `--matrix-add-well-contributions=true`
 must also be set as the mixed-precision solver operates directly on block-sparse matrices, not
-on linear operators as other OPM solvers do. For convenience, a wrapper similar to the one below
-can be used.
+on linear operators as other OPM solvers do. For convenience, a wrapper similar to the one
+below can be used.
 
 ``` bash
 OMP_NUM_THREADS=1 mpirun -np 1 --map-by numa --bind-to core build/bin/flow \
@@ -23,3 +25,6 @@ OMP_NUM_THREADS=1 mpirun -np 1 --map-by numa --bind-to core build/bin/flow \
     --linear-solver-max-iter=1024 \
     $@
 ```
+
+To invoke the original serial implementation add the `legacy-` prefix to the mixed-precision
+linear solver options.
diff --git a/opm/simulators/utils/satfunc/OilPhaseConsistencyChecks.cpp b/opm/simulators/utils/satfunc/OilPhaseConsistencyChecks.cpp
index 8a2ee6d65f1..df45ad042a8 100644
--- a/opm/simulators/utils/satfunc/OilPhaseConsistencyChecks.cpp
+++ b/opm/simulators/utils/satfunc/OilPhaseConsistencyChecks.cpp
@@ -31,7 +31,7 @@ template <typename Scalar>
 void Opm::Satfunc::PhaseChecks::Oil::SOcr_GO<Scalar>::
 testImpl(const EclEpsScalingPointsInfo<Scalar>& endPoints)
 {
-    this->sogcr_ = endPoints.Sogcr; //0.0
+    this->sogcr_ = endPoints.Sogcr;
 
     if (! std::isfinite(this->sogcr_)) {
         this->setViolated();

From 960fb43be06e0b03d6cf662549cf49f505a1567b Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Mon, 1 Jun 2026 12:14:40 -0500
Subject: [PATCH 14/19] mixed: minor updates to documentation of bsr matrix

---
 opm/simulators/linalg/mixed/bsr.c | 10 +++++-----
 opm/simulators/linalg/mixed/bsr.h |  2 ++
 2 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/opm/simulators/linalg/mixed/bsr.c b/opm/simulators/linalg/mixed/bsr.c
index 19df0d9163a..ba23e68e01f 100644
--- a/opm/simulators/linalg/mixed/bsr.c
+++ b/opm/simulators/linalg/mixed/bsr.c
@@ -91,7 +91,7 @@ void bsr_vmspmv3(bsr_matrix *A, const double *x, double *y)
             int j = colidx[k];
             __m256d vx = _mm256_loadu_pd(x+b*j);
 
-            vA[0] += _mm256_cvtps_pd(_mm_loadu_ps(AA+0))*_mm256_permute4x64_pd(vx,0x00);
+            vA[0] += _mm256_cvtps_pd(_mm_loadu_ps(AA+0))*_mm256_permute4x64_pd(vx,0x00); // 0b00000000
             vA[1] += _mm256_cvtps_pd(_mm_loadu_ps(AA+3))*_mm256_permute4x64_pd(vx,0x55); // 0b01010101
             vA[2] += _mm256_cvtps_pd(_mm_loadu_ps(AA+6))*_mm256_permute4x64_pd(vx,0xAA); // 0b10101010
         }
@@ -130,9 +130,9 @@ void bsr_vmspumv3(bsr_matrix *A, const double *x, double *y, double alpha)
             int j = colidx[k];
             __m256d vx = _mm256_loadu_pd(x+b*j);
 
-            vA[0] += _mm256_cvtps_pd(_mm_loadu_ps(AA+0))*_mm256_permute4x64_pd(vx,0b00000000); //0b01010101
-            vA[1] += _mm256_cvtps_pd(_mm_loadu_ps(AA+3))*_mm256_permute4x64_pd(vx,0b01010101); //0b01010101
-            vA[2] += _mm256_cvtps_pd(_mm_loadu_ps(AA+6))*_mm256_permute4x64_pd(vx,0b10101010); //0b01010101
+            vA[0] += _mm256_cvtps_pd(_mm_loadu_ps(AA+0))*_mm256_permute4x64_pd(vx,0x00); // 0b00000000
+            vA[1] += _mm256_cvtps_pd(_mm_loadu_ps(AA+3))*_mm256_permute4x64_pd(vx,0x55); // 0b01010101
+            vA[2] += _mm256_cvtps_pd(_mm_loadu_ps(AA+6))*_mm256_permute4x64_pd(vx,0xAA); // 0b10101010
         }
 
         // sum over columns
@@ -168,7 +168,7 @@ void bsr_vdspmv3(bsr_matrix *A, const double *x, double *y)
             int j = colidx[k];
             __m256d vx = _mm256_loadu_pd(x+b*j);
 
-            vA[0] += _mm256_loadu_pd(AA+0)*_mm256_permute4x64_pd(vx,0x00);
+            vA[0] += _mm256_loadu_pd(AA+0)*_mm256_permute4x64_pd(vx,0x00); // 0b00000000
             vA[1] += _mm256_loadu_pd(AA+3)*_mm256_permute4x64_pd(vx,0x55); // 0b01010101
             vA[2] += _mm256_loadu_pd(AA+6)*_mm256_permute4x64_pd(vx,0xAA); // 0b10101010
         }
diff --git a/opm/simulators/linalg/mixed/bsr.h b/opm/simulators/linalg/mixed/bsr.h
index f093ab5f557..a23b7bbc914 100644
--- a/opm/simulators/linalg/mixed/bsr.h
+++ b/opm/simulators/linalg/mixed/bsr.h
@@ -70,6 +70,7 @@ void bsr_vmspumv3(bsr_matrix *A, const double *x, double *y, double alpha);
  * @brief Sparse matrix-vector multiplication in mixed precision.
  *
  * @note Function is specialized for 3x3 block-sparse matrices.
+ * @note The dense 3x3 blocks are assumed to be column-major
  * @note Function uses AVX2 intrinsics.
  *
  * @param A Pointer to bsr matrix.
@@ -83,6 +84,7 @@ void bsr_vmspmv3(bsr_matrix *A, const double *x, double *y);
  * @brief Sparse matrix-vector multiplication in double precision.
  *
  * @note Function is specialized for 3x3 block-sparse matrices.
+ * @note The dense 3x3 blocks are assumed to be column-major
  * @note Function uses AVX2 intrinsics.
  *
  * @param A Pointer to bsr matrix.

From 8cefef12dcd7e6399be24073b13754397ad1a222 Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Mon, 1 Jun 2026 12:54:52 -0500
Subject: [PATCH 15/19] mixed: minor updates to MixedWrapper.hpp

---
 opm/simulators/linalg/mixed/MatrixWrapper.hpp | 21 ++++++++++++-------
 opm/simulators/linalg/mixed/SolverAdapter.hpp |  2 +-
 2 files changed, 14 insertions(+), 9 deletions(-)

diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp
index 10e273a5bd5..c27085ccfba 100644
--- a/opm/simulators/linalg/mixed/MatrixWrapper.hpp
+++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp
@@ -4,13 +4,18 @@
 
 #include <opm/simulators/linalg/mixed/bsr.h>
 
+namespace Opm
+{
 
 //! @brief Wraps c-implementation of mixed-precision matrix
 //!
+//! @note matrix is stored in single precision, but all
+//!       operations are performed in double-precision
+//!
 //! @tparam Vector the block-vector used by linear operator
 //! @tparam b block size
 template <class Vector, int b>
-class MatrixWrapper
+class MixedMatrixWrapper
 {
     public:
     virtual void mv(const Vector& x, Vector& y) const;
@@ -21,7 +26,7 @@ class MatrixWrapper
     //!
     //! @param nrows number of block rows
     //! @param nnz number of nonzero blocks
-    MatrixWrapper(int nrows, int nnz)
+    MixedMatrixWrapper(int nrows, int nnz)
     {
         nnz_=nnz;
         M_ = bsr_alloc();
@@ -29,7 +34,7 @@ class MatrixWrapper
     }
 
     //! @brief destructor
-    ~MatrixWrapper() {bsr_free(M_);}
+    ~MixedMatrixWrapper() {bsr_free(M_);}
 
     void update(double const *data);
 
@@ -42,7 +47,7 @@ class MatrixWrapper
 };
 
 template <class Vector, int b>
-void MatrixWrapper<Vector,b>::
+void MixedMatrixWrapper<Vector,b>::
 mv(const Vector& x, Vector& y) const
 {
     // mixed-precision block spmv (y = M.x)
@@ -50,7 +55,7 @@ mv(const Vector& x, Vector& y) const
 }
 
 template <class Vector, int b>
-void MatrixWrapper<Vector,b>::
+void MixedMatrixWrapper<Vector,b>::
 umv(const Vector& x, Vector& y) const
 {
     // mixed-precision block spmv with update (y += M.x)
@@ -58,7 +63,7 @@ umv(const Vector& x, Vector& y) const
 }
 
 template <class Vector, int b>
-void MatrixWrapper<Vector,b>::
+void MixedMatrixWrapper<Vector,b>::
 usmv(double alpha, const Vector& x, Vector& y) const
 {
     // scaled mixed-precision block spmv with update (y += alpha *  M.x)
@@ -66,7 +71,7 @@ usmv(double alpha, const Vector& x, Vector& y) const
 }
 
 template <class Vector, int b>
-void MatrixWrapper<Vector,b>::
+void MixedMatrixWrapper<Vector,b>::
 update(double const *data)
 {
     // transpose each dense block to make them column-major
@@ -81,5 +86,5 @@ update(double const *data)
     // downcast to single precision
     bsr_downcast(M_);
 }
-
+} // namespace Opm
 #endif // OPM_MIXED_MATRIX_HEADER_INCLUDED
diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp
index 0b4380f3d93..94a87ffbf20 100644
--- a/opm/simulators/linalg/mixed/SolverAdapter.hpp
+++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp
@@ -100,7 +100,7 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
     using AbstractScalarProductType = Dune::ScalarProduct<Vector>;
 
     static constexpr auto block_size = Vector::block_type::dimension;
-    using MixedMatrixType      = MatrixWrapper<Vector, block_size>;
+    using MixedMatrixType      = Opm::MixedMatrixWrapper<Vector, block_size>;
     using MixedOperatorType    = MixedOperator<MixedMatrixType, Vector, Comm>::type;
     using OptimizedProductType = SeqOptmizedProduct<Vector>;
 

From fcb5abb716b9738dddb35eae3d5e0ecc65b7144b Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Mon, 1 Jun 2026 13:18:43 -0500
Subject: [PATCH 16/19] mixed: extract sparsity pattern using iterators

---
 opm/simulators/linalg/mixed/PreconditionerWrapper.hpp | 6 +++---
 opm/simulators/linalg/mixed/SolverAdapter.hpp         | 6 +++---
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp
index 2f91f2bc1b1..7e7818a66b5 100644
--- a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp
+++ b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp
@@ -44,11 +44,11 @@ class MixedPreconditioner : public Dune::PreconditionerWithUpdate<X, Y>
         rows[0]  = 0;
         for(auto row=A.begin(); row!=A.end(); row++)
         {
-            for(unsigned int i=0; i<row->getsize(); i++)
+            for(auto col = row->begin(); col != row->end(); ++col)
             {
-                cols[icol++] = row->getindexptr()[i];
+                cols[icol++] = col.index();
             }
-            rows[irow+1]     = rows[irow]+row->getsize();
+            rows[irow+1] = icol;
             irow++;
         }
 
diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp
index 94a87ffbf20..01b34782769 100644
--- a/opm/simulators/linalg/mixed/SolverAdapter.hpp
+++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp
@@ -141,11 +141,11 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
         rows[0]  = 0;
         for(auto row=A.begin(); row!=A.end(); row++)
         {
-            for(unsigned int i=0; i<row->getsize(); i++)
+            for(auto col = row->begin(); col != row->end(); ++col)
             {
-                cols[icol++] = row->getindexptr()[i];
+                cols[icol++] = col.index();
             }
-            rows[irow+1]     = rows[irow]+row->getsize();
+            rows[irow+1] = icol;
             irow++;
         }
 

From a2d7a809764e296cb678dffb22a0a6f7c2335424 Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Thu, 2 Jul 2026 12:35:29 -0500
Subject: [PATCH 17/19] mixed: remove block size as explicit template parameter

---
 opm/simulators/linalg/mixed/MatrixWrapper.hpp | 27 +++++++++++--------
 opm/simulators/linalg/mixed/SolverAdapter.hpp |  2 +-
 2 files changed, 17 insertions(+), 12 deletions(-)

diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp
index c27085ccfba..f13ff28a740 100644
--- a/opm/simulators/linalg/mixed/MatrixWrapper.hpp
+++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp
@@ -14,10 +14,14 @@ namespace Opm
 //!
 //! @tparam Vector the block-vector used by linear operator
 //! @tparam b block size
-template <class Vector, int b>
+template <class Vector>
 class MixedMatrixWrapper
 {
     public:
+
+    // extract block size
+    static constexpr auto block_size = Vector::block_type::dimension;
+
     virtual void mv(const Vector& x, Vector& y) const;
     virtual void umv(const Vector& x, Vector& y) const;
     virtual void usmv(double alpha, const Vector& x, Vector& y) const;
@@ -30,7 +34,7 @@ class MixedMatrixWrapper
     {
         nnz_=nnz;
         M_ = bsr_alloc();
-        bsr_init(M_, nrows, nnz, b);
+        bsr_init(M_, nrows, nnz, block_size);
     }
 
     //! @brief destructor
@@ -46,36 +50,37 @@ class MixedMatrixWrapper
     bsr_matrix  *M_;
 };
 
-template <class Vector, int b>
-void MixedMatrixWrapper<Vector,b>::
+template <class Vector>
+void MixedMatrixWrapper<Vector>::
 mv(const Vector& x, Vector& y) const
 {
     // mixed-precision block spmv (y = M.x)
     bsr_vmspmv3(M_, &x[0][0], &y[0][0]);
 }
 
-template <class Vector, int b>
-void MixedMatrixWrapper<Vector,b>::
+template <class Vector>
+void MixedMatrixWrapper<Vector>::
 umv(const Vector& x, Vector& y) const
 {
     // mixed-precision block spmv with update (y += M.x)
     bsr_vmspumv3(M_, &x[0][0], &y[0][0], 1.0);
 }
 
-template <class Vector, int b>
-void MixedMatrixWrapper<Vector,b>::
+template <class Vector>
+void MixedMatrixWrapper<Vector>::
 usmv(double alpha, const Vector& x, Vector& y) const
 {
     // scaled mixed-precision block spmv with update (y += alpha *  M.x)
     bsr_vmspumv3(M_, &x[0][0], &y[0][0], alpha);
 }
 
-template <class Vector, int b>
-void MixedMatrixWrapper<Vector,b>::
+template <class Vector>
+void MixedMatrixWrapper<Vector>::
 update(double const *data)
 {
     // transpose each dense block to make them column-major
-    int bb=b*b;
+    int const b = block_size;
+    int const bb=b*b;
     double B[bb];
     for(int k=0;k<nnz_;k++)
     {
diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp
index 01b34782769..4c3f3a74d9e 100644
--- a/opm/simulators/linalg/mixed/SolverAdapter.hpp
+++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp
@@ -100,7 +100,7 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
     using AbstractScalarProductType = Dune::ScalarProduct<Vector>;
 
     static constexpr auto block_size = Vector::block_type::dimension;
-    using MixedMatrixType      = Opm::MixedMatrixWrapper<Vector, block_size>;
+    using MixedMatrixType      = Opm::MixedMatrixWrapper<Vector>;
     using MixedOperatorType    = MixedOperator<MixedMatrixType, Vector, Comm>::type;
     using OptimizedProductType = SeqOptmizedProduct<Vector>;
 

From 5db7457274ab79949c1161b9af8100ec15ae3a11 Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Thu, 2 Jul 2026 12:46:21 -0500
Subject: [PATCH 18/19] mixed: renamed MixedAdapter to MixedBiCGSTABSolver

---
 opm/simulators/linalg/FlexibleSolver_impl.hpp | 2 +-
 opm/simulators/linalg/mixed/SolverAdapter.hpp | 6 +++---
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/opm/simulators/linalg/FlexibleSolver_impl.hpp b/opm/simulators/linalg/FlexibleSolver_impl.hpp
index 51a6cdc3238..b9e85a593b6 100644
--- a/opm/simulators/linalg/FlexibleSolver_impl.hpp
+++ b/opm/simulators/linalg/FlexibleSolver_impl.hpp
@@ -205,7 +205,7 @@ namespace Dune
             }
             // mixed-solver adapter starts here
           } else if (solver_type == "mixed-precision") {
-                linsolver_ = std::make_shared<Dune::MixedAdapter<Comm,Operator,VectorType>>(linearoperator_for_solver_,
+                linsolver_ = std::make_shared<Dune::MixedBiCGSTABSolver<Comm,Operator,VectorType>>(linearoperator_for_solver_,
                                                                             scalarproduct_,
                                                                             preconditioner_,
                                                                             tol, // desired residual reduction factor
diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp
index 4c3f3a74d9e..f4fd2237fb6 100644
--- a/opm/simulators/linalg/mixed/SolverAdapter.hpp
+++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp
@@ -92,7 +92,7 @@ struct MixedOperator<Matrix, Vector, Dune::Amg::SequentialInformation>
 //! @tparam Operator the linear operator passed to FlexibleLinearSolver
 //! @tparam Vector the block vector passed to FlexibleLinearSolver
 template <class Comm, class Operator, class Vector>
-class MixedAdapter:public InverseOperator<Vector, Vector>
+class MixedBiCGSTABSolver:public InverseOperator<Vector, Vector>
 {
     public:
 
@@ -113,7 +113,7 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
     //! @param maxit maximum number of iterations for the linear solver
     //! @param verbose verbosity level
     //! @param comm the communication object.
-    MixedAdapter(Operator *op,
+    MixedBiCGSTABSolver(Operator *op,
                  std::shared_ptr<AbstractScalarProductType> sp,
                  std::shared_ptr<AbstractPrecondType> prec,
                  const double& tol,
@@ -186,7 +186,7 @@ class MixedAdapter:public InverseOperator<Vector, Vector>
         x=0;
         b=0;
         res.reduction = reduction;
-        OPM_THROW(std::invalid_argument, "MixedAdapter::apply(...) not implemented yet.");
+        OPM_THROW(std::invalid_argument, "MixedBiCGSTABSolver::apply(...) not implemented yet.");
     }
 
     virtual Dune::SolverCategory::Category category() const override{return Dune::SolverCategory::overlapping;};

From da0b6c53d8201c3c86a14eabcd3cac3c032b75ce Mon Sep 17 00:00:00 2001
From: Kjetil B Haugen <kjetil@haugenlabs.com>
Date: Thu, 2 Jul 2026 13:19:18 -0500
Subject: [PATCH 19/19] mixed: throw error for unsupported matrix operators

---
 opm/simulators/linalg/mixed/MatrixWrapper.hpp |  2 ++
 opm/simulators/linalg/mixed/SolverAdapter.hpp | 14 ++++++++++++++
 2 files changed, 16 insertions(+)

diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp
index f13ff28a740..1364fb3a6fd 100644
--- a/opm/simulators/linalg/mixed/MatrixWrapper.hpp
+++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp
@@ -32,6 +32,8 @@ class MixedMatrixWrapper
     //! @param nnz number of nonzero blocks
     MixedMatrixWrapper(int nrows, int nnz)
     {
+        if constexpr(block_size!=3) OPM_THROW(std::invalid_argument, "MixedMatrixWrapper only supports block size == 3! \n");
+
         nnz_=nnz;
         M_ = bsr_alloc();
         bsr_init(M_, nrows, nnz, block_size);
diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp
index f4fd2237fb6..ccff6ff94e6 100644
--- a/opm/simulators/linalg/mixed/SolverAdapter.hpp
+++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp
@@ -149,6 +149,20 @@ class MixedBiCGSTABSolver:public InverseOperator<Vector, Vector>
             irow++;
         }
 
+        // The following skeleton is in preparation for better support for various MatrixAdapter. For now, it simply throws an error if
+        // the operator provided is not supported.
+        using MatrixType = std::remove_const_t<std::remove_reference_t<decltype(op->getmat())>>;
+        if constexpr (std::is_same_v<std::remove_pointer_t<Operator>, Dune::MatrixAdapter<MatrixType, Vector, Vector>>)
+        {
+        }
+        else if constexpr (std::is_same_v<std::remove_pointer_t<Operator>, Opm::GhostLastMatrixAdapter<MatrixType, Vector, Vector, Comm>>)
+        {
+        }
+        else
+        {
+            OPM_THROW(std::invalid_argument, "MixedBiCGSTABSolver only supports Dune::MatrixAdapter and Opm::GhostLastMatrixAdapter\n");
+        }
+
         //initialize mixed operator and optimized scalar product
         double_operator_ = op;
         if constexpr (std::is_same_v<Comm, Dune::Amg::SequentialInformation>)