solve: mixed precision solver as flexible solver

opm/simulators/flow/BlackoilModel_impl.hpp

@@ -541,16 +541,16 @@ solveJacobianSystem(BVector& x)
             int *rows = bsr_jacobian_->rowptr;
             int *cols = bsr_jacobian_->colidx;
 
-            int irow=0;
-            int icol=0;
-            rows[0]=0;
+            int irow = 0;
+            int icol = 0;
+            rows[0]  = 0;
             for(auto row=jacobian.begin(); row!=jacobian.end(); row++)
             {
-                for(unsigned int i=0;i<row->getsize();i++)
+                for(unsigned int i=0; i<row->getsize(); i++)
                 {
-                    cols[icol++]=row->getindexptr()[i];
+                    cols[icol++] = row->getindexptr()[i];
                 }
-                rows[irow+1]= rows[irow]+row->getsize();
+                rows[irow+1]     = rows[irow]+row->getsize();
                 irow++;
             }
 
@@ -576,8 +576,8 @@ solveJacobianSystem(BVector& x)
 
 
 
-        //linSolver.prepare(jacobian, residual); // what happens here?
-
+        linSolver.prepare(jacobian, residual); // what happens here?
+/*
 
         // transpose each dense block to make them column-major
         double M[9];
@@ -593,13 +593,13 @@ solveJacobianSystem(BVector& x)
         }
 
         //bsr_downcast(bsr_jacobian_);
-
+*/
 
 
         linear_solve_setup_time_ = perfTimer.stop();
 
 
-        //linSolver.setResidual(residual);
+        linSolver.setResidual(residual);
 
 
 
@@ -616,7 +616,7 @@ solveJacobianSystem(BVector& x)
         // copy residual before calling solve
         BVector r(residual);
 */
-
+/*
         //solve Ax=b, where A=jacobian and b=r
         int count = bslv_pbicgstab3(slv_memory_, bsr_jacobian_, &residual[0][0], &x[0][0]);
         linSolver.set_iterations(count);
@@ -624,8 +624,8 @@ solveJacobianSystem(BVector& x)
 //        vec_copy(&residual[0][0], slv_memory_->dtmp[3],n);
         //bslv_info(slv_memory_,count);
         if(count >= slv_memory_->max_iter) printf("------------------------------> %d\n",count);
-
-        //linSolver.solve(x);
+*/
+        linSolver.solve(x);
 /*
         //compute r-Ax, where A=jacobian
         jacobian.mmv(x,r);

opm/simulators/linalg/FlexibleSolver_impl.hpp

@@ -48,6 +48,113 @@
 #endif
 #endif
 
+#include <opm/simulators/linalg/HaugenLabs.h>
+
+namespace Dune
+{
+template <class X, class M>
+class MixedSolver : public InverseOperator<X,X>
+{
+    public:
+
+    MixedSolver()
+    {
+        printf("MixedSolver::MixedSolver()\n");getchar();
+    }
+
+    MixedSolver(double tol, int maxiter)
+    {
+        printf("tol=%.4e\n",tol);
+        printf("maxiter=%d\n",maxiter);
+        printf("MixedSolver::MixedSolver()\n");
+        getchar();
+    }
+
+    MixedSolver(M A, double tol, int maxiter)
+    {
+        int nrows = A.N();
+        int nnz   = A.nonzeroes();
+        int b     = A[0][0].N();
+
+        // create jacobian matrix object and allocate various arrays
+        jacobian_ = bsr_new();
+        bsr_init(jacobian_,nrows,nnz,b);
+
+        // initialize sparsity pattern
+        int *rows = jacobian_->rowptr;
+        int *cols = jacobian_->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++;
+        }
+
+        // print select information
+        bsr_info(jacobian_);
+        bsr_sparsity(jacobian_,"jacobian");
+
+        // allocate and intialize solver memory
+        mem_ =bslv_new();
+        bslv_init(mem_, tol, maxiter, jacobian_);
+
+        //pointer to nonzero blocks
+        data_ = &A[0][0][0][0];
+
+    }
+
+    virtual void apply (X& x, X& b, InverseOperatorResult& res) override
+    {
+        // transpose each dense block to make them column-major
+        double B[9];
+        for(int k=0;k<jacobian_->nnz;k++)
+        {
+            for(int i=0;i<3;i++) for(int j=0;j<3;j++) B[3*j+i] = data_[9*k + 3*i + j];
+            for(int i=0;i<9;i++)
+            {
+                jacobian_->dbl[9*k + i] = B[i];
+                //jacobian_->flt[9*k + i] = B[i];
+            }
+        }
+
+        // downcast to allow mixed precision
+        bsr_downcast(jacobian_);
+
+        // solve linear system
+        int count = bslv_pbicgstab3(mem_, jacobian_, &b[0][0], &x[0][0]);
+        //bslv_info(mem_,count);
+
+        // return convergence information
+        res.converged  = 1;
+        res.reduction  = mem_->e[count];
+        res.iterations = count;
+        //getchar();
+    }
+
+    virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res) override
+    {
+        x=0;
+        b=0;
+        res.reduction = reduction;
+        OPM_THROW(std::invalid_argument, "MixedSolver::apply(4) not implemented yet.");
+
+    }
+
+    private:
+    bsr_matrix  *jacobian_;
+    bslv_memory *mem_;
+    double const *data_;
+};
+
+}
+
 namespace Dune
 {
     /// Create a sequential solver.
@@ -181,6 +288,12 @@ namespace Dune
                                                                             tol, // desired residual reduction factor
                                                                             maxiter, // maximum number of iterations
                                                                             verbosity);
+        } else if (solver_type == "mixed-bicgstab") {
+            using MatrixType = decltype(linearoperator_for_solver_->getmat());
+            linsolver_ = std::make_shared<Dune::MixedSolver<VectorType,MatrixType>>(
+                                                                            linearoperator_for_solver_->getmat(),
+                                                                            tol,
+                                                                            maxiter);
         } else if (solver_type == "loopsolver") {
             linsolver_ = std::make_shared<Dune::LoopSolver<VectorType>>(*linearoperator_for_solver_,
                                                                         *scalarproduct_,
@@ -197,7 +310,6 @@ namespace Dune
                                                                                   restart,
                                                                                   maxiter, // maximum number of iterations
                                                                                   verbosity);
-
         } else {
             if constexpr (!Opm::is_gpu_operator_v<Operator>) {
                 if (solver_type == "flexgmres") {

opm/simulators/linalg/setupPropertyTree.cpp

@@ -241,6 +241,11 @@ setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters
         return setupILU(conf, p);
     }
 
+    // mixed-precision ILU0
+    if (conf == "mixed-ilu0") {
+        return setupMixedILU(conf, p);
+    }
+
     if (conf == "dilu") {
         return setupDILU(conf, p);
     }
@@ -394,6 +399,19 @@ setupILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParamet
     return prm;
 }
 
+PropertyTree
+setupMixedILU([[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-bicgstab"s);
+    return prm;
+}
+
+
 PropertyTree
 setupDILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
 {

opm/simulators/linalg/setupPropertyTree.hpp

@@ -37,6 +37,7 @@ 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 setupMixedILU(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupDILU(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupUMFPack(const std::string& conf, const FlowLinearSolverParameters& p);
 

runme

@@ -1,5 +1,5 @@
 OMP_NUM_THREADS=1 mpirun -np 1 --map-by numa --bind-to core build/bin/flow \
-    --linear-solver=ilu0 \
+    --linear-solver=mixed-ilu0 \
     --linear-solver-reduction=1e-3 \
     --linear-solver-max-iter=200 \
     --matrix-add-well-contributions=true \