KrylovSolver.hpp

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 TotalEnergies
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_LINEARALGEBRA_SOLVERS_KRYLOVSOLVER_HPP_
#define GEOS_LINEARALGEBRA_SOLVERS_KRYLOVSOLVER_HPP_
 
#include "linearAlgebra/utilities/BlockVectorView.hpp"
#include "linearAlgebra/utilities/BlockVector.hpp"
#include "linearAlgebra/utilities/BlockOperatorView.hpp"
#include "linearAlgebra/utilities/LinearSolverParameters.hpp"
#include "linearAlgebra/utilities/LinearSolverResult.hpp"
 
namespace geos
{
 
template< typename VECTOR >
class KrylovSolver : public LinearOperator< VECTOR >
{
public:
 
  using Base = LinearOperator< VECTOR >;
 
  using Vector = typename Base::Vector;
 
  static std::unique_ptr< KrylovSolver< VECTOR > > create( LinearSolverParameters const & parameters,
                                                           LinearOperator< VECTOR > const & matrix,
                                                           LinearOperator< VECTOR > const & precond );
 
  KrylovSolver( LinearSolverParameters params,
                LinearOperator< Vector > const & matrix,
                LinearOperator< Vector > const & precond );
 
  virtual ~KrylovSolver() override = default;
 
  virtual void solve( Vector const & b, Vector & x ) const = 0;
 
 
  virtual void apply( Vector const & src, Vector & dst ) const override final
  {
    solve( src, dst );
  }
 
  virtual globalIndex numGlobalRows() const override final
  {
    return m_operator.numGlobalRows();
  }
 
  virtual globalIndex numGlobalCols() const override final
  {
    return m_operator.numGlobalCols();
  }
 
  virtual localIndex numLocalRows() const override final
  {
    return m_operator.numLocalRows();
  }
 
  virtual localIndex numLocalCols() const override final
  {
    return m_operator.numLocalCols();
  }
 
  virtual MPI_Comm comm() const override final
  {
    return m_operator.comm();
  }
 
  LinearSolverParameters const & parameters() const
  {
    return m_params;
  }
 
  LinearSolverResult const & result() const
  {
    return m_result;
  }
 
  arrayView1d< real64 const > history() const
  {
    return m_residualNorms;
  }
 
  virtual string methodName() const = 0;
 
private:
 
 
  template< typename VEC >
  struct VectorStorageHelper
  {
    using type = VEC;
 
    static VEC createFrom( VEC const & src )
    {
      VEC v;
      v.create( src.localSize(), src.comm() );
      return v;
    }
  };
 
  template< typename VEC >
  struct VectorStorageHelper< BlockVectorView< VEC > >
  {
    using type = BlockVector< VEC >;
 
    static BlockVector< VEC > createFrom( BlockVectorView< VEC > const & src )
    {
      BlockVector< VEC > v( src.blockSize() );
      for( localIndex i = 0; i < src.blockSize(); ++i )
      {
        v.block( i ).create( src.block( i ).localSize(), src.block( i ).comm() );
      }
      return v;
    }
  };
 
 
protected:
 
  using VectorTemp = typename VectorStorageHelper< VECTOR >::type;
 
  static VectorTemp createTempVector( Vector const & src )
  {
    return VectorStorageHelper< VECTOR >::createFrom( src );
  }
 
  void logProgress() const;
 
  void logResult() const;
 
  LinearSolverParameters m_params;
 
  LinearOperator< Vector > const & m_operator;
 
  LinearOperator< Vector > const & m_precond;
 
  mutable LinearSolverResult m_result;
 
  mutable array1d< real64 > m_residualNorms;
};
 
} //namespace geos
 
#endif //GEOS_LINEARALGEBRA_SOLVERS_KRYLOVSOLVER_HPP_