BlockPreconditioner.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2018-2020 Total, S.A
 * Copyright (c) 2019-     GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */

#ifndef GEOSX_LINEARALGEBRA_SOLVERS_BLOCKPRECONDITIONER_HPP_
#define GEOSX_LINEARALGEBRA_SOLVERS_BLOCKPRECONDITIONER_HPP_

#include "linearAlgebra/DofManager.hpp"
#include "linearAlgebra/solvers/PreconditionerBase.hpp"
#include "linearAlgebra/utilities/BlockOperator.hpp"
#include "linearAlgebra/utilities/BlockVector.hpp"

namespace geosx
{

enum class SchurComplementOption
{
  None,                  
  FirstBlockDiagonal,    
  RowsumDiagonalProbing, 
  FirstBlockUserDefined  
};

enum class BlockScalingOption
{
  None,          
  FrobeniusNorm, 
  UserProvided   
};

enum class BlockShapeOption
{
  Diagonal,            
  UpperTriangular,     
  LowerUpperTriangular 
};

/*
 * Since formulas in Doxygen are broken with 1.8.13 and certain versions of ghostscript,
 * keeping this documentation in a separate comment block for now. Should be moved into
 * documentation of BlockPreconditioner class.
 *
 * This class implements a 2x2 block preconditioner of the form:
 * @f$
 * M^{-1} =
 * \begin{bmatrix}
 * I_{1} & -M_{1}^{-1}A_{12} \\
 *       &  I_{2}
 * \end{bmatrix}
 * \begin{bmatrix}
 * M_{1}^{-1} &            \\
 *            & M_{2}^{-1}
 * \end{bmatrix}
 * \begin{bmatrix}
 *  I_{1}            &       \\
 * -A_{21}M_{1}^{-1} & I_{2}
 * \end{bmatrix}
 * @f$
 * with
 * @f$ M_1^{-1} ~= A_{11}^{-1} @f$, @f$ M_2^{-1} ~= S_{22}^{-1} @f$, and
 * @f$ S_{22} ~= A_{22} - A_{21} A_{11}^{-1} A_{12} @f$ being the (optional) approximate Schur complement.
 *
 * The first step (predictor solve with @f$ M_{1} @f$) is optional. Enabling it results in a more powerful
 * preconditioner, but involves two applications of @f$ M_{1}^{-1} @f$ instead of just one on each solve.
 *
 * The user provides individual block solvers @f$ M_{1} @f$ and @f$ M_{2} @f$ as well as
 * a description of the block split of monolithic matrix in terms of DOF components.
 */

template< typename LAI >
class BlockPreconditioner : public PreconditionerBase< LAI >
{
public:

  using Base = PreconditionerBase< LAI >;

  using Vector = typename Base::Vector;

  using Matrix = typename Base::Matrix;

  explicit BlockPreconditioner( BlockShapeOption const shapeOption,
                                SchurComplementOption const schurOption,
                                BlockScalingOption const scalingOption );

  virtual ~BlockPreconditioner() override;

  void setupBlock( localIndex const blockIndex,
                   std::vector< DofManager::SubComponent > blockDofs,
                   std::unique_ptr< PreconditionerBase< LAI > > solver,
                   real64 const scaling = 1.0 );


  using PreconditionerBase< LAI >::compute;

  virtual void compute( Matrix const & mat,
                        DofManager const & dofManager ) override;

  virtual void apply( Vector const & src, Vector & dst ) const override;

  virtual void clear() override;


private:

  void reinitialize( Matrix const & mat, DofManager const & dofManager );

  void applyBlockScaling();

  void computeSchurComplement();

  BlockShapeOption m_shapeOption;

  SchurComplementOption m_schurOption;

  BlockScalingOption m_scalingOption;

  std::array< std::vector< DofManager::SubComponent >, 2 > m_blockDofs;

  std::array< Matrix, 2 > m_restrictors;

  std::array< Matrix, 2 > m_prolongators;

  BlockOperator< Vector, Matrix > m_matBlocks;

  std::array< std::unique_ptr< PreconditionerBase< LAI > >, 2 > m_solvers;

  std::array< real64, 2 > m_scaling;

  mutable BlockVector< Vector > m_rhs;

  mutable BlockVector< Vector > m_sol;
};

} //namespace geosx

#endif //GEOSX_LINEARALGEBRA_SOLVERS_BLOCKPRECONDITIONER_HPP_