Avoid recreating linear solver for every TEMP solve.

opm/simulators/flow/GenericTemperatureModel.hpp

@@ -36,6 +36,7 @@
 
 #include <opm/models/blackoil/blackoilmodel.hh>
 
+#include <opm/simulators/linalg/FlexibleSolver.hpp>
 #include <opm/simulators/linalg/matrixblock.hh>
 
 #include <cstddef>
@@ -79,6 +80,7 @@ protected:
      */
     void doInit(std::size_t numGridDof);
 
+    void setupLinearSolver(const EnergyMatrix& M);
     bool linearSolve_(const EnergyMatrix& M, EnergyVector& x, EnergyVector& b);
 
     const GridView& gridView_;
@@ -91,6 +93,14 @@ protected:
     std::vector<Scalar> energy_rates_;
     bool doTemp_{false};
     Scalar maxTempChange_{5.0};
+
+    using AbstractSolverType = Dune::InverseOperator<EnergyVector, EnergyVector>;
+    using AbstractOperatorType = Dune::AssembledLinearOperator<EnergyMatrix, EnergyVector, EnergyVector>;
+    using AbstractPreconditionerType = Dune::PreconditionerWithUpdate<EnergyVector, EnergyVector>;
+
+    std::unique_ptr<AbstractSolverType> linear_solver_;
+    std::unique_ptr<AbstractOperatorType> op_;
+    AbstractPreconditionerType* pre_ = nullptr;
 };
 
 } // namespace Opm

opm/simulators/flow/GenericTemperatureModel_impl.hpp

@@ -125,13 +125,11 @@ doInit(std::size_t numGridDof)
 }
 
 template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
-bool GenericTemperatureModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
-linearSolve_(const EnergyMatrix& M, EnergyVector& x, EnergyVector& b)
+void GenericTemperatureModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
+setupLinearSolver(const EnergyMatrix& matrix)
 {
-    x = 0.0;
     Scalar tolerance = 1e-2;
     int maxIter = 100;
-
     int verbosity = 0;
     PropertyTree prm;
     prm.put("maxiter", maxIter);
@@ -141,42 +139,43 @@ linearSolve_(const EnergyMatrix& M, EnergyVector& x, EnergyVector& b)
     prm.put("preconditioner.type", std::string("ParOverILU0"));
 
 #if HAVE_MPI
-    if(gridView_.grid().comm().size() > 1)
+    if (gridView_.grid().comm().size() > 1)
     {
         auto [energyOperator, solver] =
-            createParallelFlexibleSolver<EnergyVector>(gridView_.grid(), M, prm);
-        (void) energyOperator;
-
-        Dune::InverseOperatorResult result;
-        solver->apply(x, b, result);
-
-        // return the result of the solver
-        return result.converged;
+            createParallelFlexibleSolver<EnergyVector>(gridView_.grid(), matrix, prm);
+        op_ = std::move(energyOperator);
+        pre_ = &solver->preconditioner();
+        linear_solver_ = std::move(solver);
     }
     else
 #endif
     {
-        using EnergySolver = Dune::BiCGSTABSolver<EnergyVector>;
-        using EnergyOperator = Dune::MatrixAdapter<EnergyMatrix,EnergyVector,EnergyVector>;
-        using EnergyScalarProduct = Dune::SeqScalarProduct<EnergyVector>;
-        using EnergyPreconditioner = Dune::SeqILU< EnergyMatrix,EnergyVector,EnergyVector>;
-
-        EnergyOperator energyOperator(M);
-        EnergyScalarProduct energyScalarProduct;
-        EnergyPreconditioner energyPreconditioner(M, 0, 1); // results in ILU0
-
-        EnergySolver solver (energyOperator, energyScalarProduct,
-                                  energyPreconditioner, tolerance, maxIter,
-                                  verbosity);
-
-        Dune::InverseOperatorResult result;
-        solver.apply(x, b, result);
-
-        // return the result of the solver
-        return result.converged;
+        using SeqOperator = Dune::MatrixAdapter<EnergyMatrix, EnergyVector, EnergyVector>;
+        using SeqSolver = Dune::FlexibleSolver<SeqOperator>;
+        auto seqOp = std::make_unique<SeqOperator>(matrix);
+        auto dummyWeights = [](){ return EnergyVector(); };
+        auto seqSolver = std::make_unique<SeqSolver>(*seqOp, prm, dummyWeights, 0);
+        op_ = std::move(seqOp);
+        pre_ = &seqSolver->preconditioner();
+        linear_solver_ = std::move(seqSolver);
     }
 }
 
+template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
+bool GenericTemperatureModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
+linearSolve_(const EnergyMatrix& M, EnergyVector& x, EnergyVector& b)
+{
+    if (!linear_solver_) {
+        setupLinearSolver(M);
+    } else {
+        pre_->update();
+    }
+
+    Dune::InverseOperatorResult result;
+    linear_solver_->apply(x, b, result);
+    return result.converged;
+}
+
 } // namespace Opm
 
 #endif // OPM_GENERIC_TEMPERATURE_MODEL_IMPL_HPP