KrylovSolver.hpp

/*
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 * Copyright (c) 2019, Lawrence Livermore National Security, LLC.
 *
 * Produced at the Lawrence Livermore National Laboratory
 *
 * LLNL-CODE-746361
 *
 * All rights reserved. See COPYRIGHT for details.
 *
 * This file is part of the GEOSX Simulation Framework.
 *
 * GEOSX is a free software; you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License (as published by the
 * Free Software Foundation) version 2.1 dated February 1999.
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 */

#ifndef GEOSX_LINEARALGEBRA_SOLVERS_KRYLOVSOLVER_HPP_
#define GEOSX_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 geosx
{

template< typename Vector > class LinearOperator;
template< typename Vector > class BlockVectorView;

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( LinearOperator< Vector > const & matrix,
                LinearOperator< Vector > const & precond,
                real64 const tolerance,
                localIndex const maxIterations,
                integer const verbosity );

  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();
  }

  LinearSolverResult const & result() const
  {
    return m_result;
  }

  arrayView1d< real64 const > history() const
  {
    return m_residualNorms;
  }

  integer getLogLevel() const
  {
    return m_logLevel;
  }

  virtual string methodName() const = 0;

private:


  template< typename VEC >
  struct VectorStorageHelper
  {
    using type = VEC;

    static VEC createFrom( VEC const & src )
    {
      VEC v;
      v.createWithLocalSize( src.localSize(), src.getComm() );
      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 ).createWithLocalSize( src.block( i ).localSize(), src.block( i ).getComm() );
      }
      return v;
    }
  };


protected:

  using VectorTemp = typename VectorStorageHelper< VECTOR >::type;

  static VectorTemp createTempVector( Vector const & src )
  {
    return VectorStorageHelper< VECTOR >::createFrom( src );
  }

  void logProgress( localIndex const iter, real64 const rnorm ) const
  {
    m_residualNorms[iter] = rnorm;
    GEOSX_LOG_LEVEL_RANK_0( 2, methodName() << " iteration " << iter << ": residual = " << rnorm );
  }

  void logResult() const
  {
    GEOSX_LOG_LEVEL_RANK_0( 1, methodName() << ' ' <<
                            ( m_result.success() ? "converged" : "failed to converge" ) <<
                            " in " << m_result.numIterations << " iterations " <<
                            "(" << m_result.solveTime << " s)" );
  }

  LinearOperator< Vector > const & m_operator;

  LinearOperator< Vector > const & m_precond;

  real64 m_tolerance;

  localIndex m_maxIterations;

  integer m_logLevel;

  mutable LinearSolverResult m_result;

  mutable array1d< real64 > m_residualNorms;
};

} //namespace geosx

#endif //GEOSX_LINEARALGEBRA_SOLVERS_KRYLOVSOLVER_HPP_