feat: Reuse LS factorization

src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp

@@ -109,6 +109,7 @@ struct LinearSolverParameters
     integer replaceTinyPivot = 1;     ///< Whether to replace tiny pivots by sqrt(epsilon)*norm(A)
     integer iterativeRefine = 1;      ///< Whether to perform iterative refinement
     integer parallel = 1;             ///< Whether to use a parallel solver (instead of a serial one)
+    integer reuseFactorization = 0;   ///< Whether to reuse the LU factorization or not
   }
   direct;                             ///< direct solver parameter struct
 

src/coreComponents/physicsSolvers/LinearSolverParameters.cpp

@@ -86,6 +86,11 @@ LinearSolverParametersInput::LinearSolverParametersInput( string const & name,
     setInputFlag( InputFlags::OPTIONAL ).
     setDescription( "Whether to use a parallel solver (instead of a serial one)" );
 
+  registerWrapper( viewKeyStruct::reuseFactorizationString(), &m_parameters.direct.reuseFactorization ).
+    setApplyDefaultValue( m_parameters.direct.reuseFactorization ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setDescription( "Whether to reuse the LU factorization" );  
+
   registerWrapper( viewKeyStruct::krylovMaxIterString(), &m_parameters.krylov.maxIterations ).
     setApplyDefaultValue( m_parameters.krylov.maxIterations ).
     setInputFlag( InputFlags::OPTIONAL ).

src/coreComponents/physicsSolvers/LinearSolverParameters.hpp

@@ -92,6 +92,8 @@ class LinearSolverParametersInput : public dataRepository::Group
     static constexpr char const * directIterRefString() { return "directIterRef"; }
     /// direct solver parallelism key
     static constexpr char const * directParallelString() { return "directParallel"; }
+    /// reuse factorization key
+    static constexpr char const * reuseFactorizationString() { return "reuseFactorization"; }
 
     /// Krylov max iterations key
     static constexpr char const * krylovMaxIterString() { return "krylovMaxIter"; }

src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp

@@ -1312,16 +1312,21 @@ void PhysicsSolverBase::solveLinearSystem( DofManager const & dofManager,
 
   if( params.solverType == LinearSolverParameters::SolverType::direct || !m_precond )
   {
-    std::unique_ptr< LinearSolverBase< LAInterface > > solver = LAInterface::createSolver( params );
+    if (!m_linearSolver)
+    { 
+      m_linearSolver = LAInterface::createSolver( params );
+      m_linearSolver->setup( matrix );
+    }
+    if ( params.solverType != LinearSolverParameters::SolverType::direct || !params.direct.reuseFactorization )
     {
       Timer timer_setup( m_timers["linear solver setup"] );
-      solver->setup( matrix );
+      m_linearSolver->setup( matrix );
     }
     {
       Timer timer_setup( m_timers["linear solver solve"] );
-      solver->solve( rhs, solution );
+      m_linearSolver->solve( rhs, solution );
     }
-    m_linearSolverResult = solver->result();
+    m_linearSolverResult = m_linearSolver->result();
   }
   else
   {

src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp

@@ -1030,6 +1030,9 @@ class PhysicsSolverBase : public ExecutableGroup
   /// Local system matrix and rhs
   CRSMatrix< real64, globalIndex > m_localMatrix;
 
+  /// Custom linear solver for the "native" solver type
+  std::unique_ptr< LinearSolverBase< LAInterface > > m_linearSolver;
+
   /// Custom preconditioner for the "native" iterative solver
   std::unique_ptr< PreconditionerBase< LAInterface > > m_precond;
 

src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp

@@ -530,6 +530,8 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::setupSystem( DomainPartition &
   solution.create( numLocalRows, MPI_COMM_GEOS );
 
   setUpDflux_dApertureMatrix( domain, dofManager, localMatrix );
+
+  dofManager.printFieldInfo();
 }
 
 template< typename POROMECHANICS_SOLVER >
@@ -690,6 +692,7 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::assembleSystem( real64 const t
                                                                              [&]( localIndex const,
                                                                                   SurfaceElementSubRegion & subRegion )
       {
+        std::cout << "Num of frac elements: " << subRegion.size() << std::endl;
         flowSolver()->accumulationAssemblyLaunch( dofManager,
                                                   subRegion,
                                                   localMatrix,