BlockPreconditioner.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * 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_BLOCKPRECONDITIONER_HPP_
#define GEOS_LINEARALGEBRA_SOLVERS_BLOCKPRECONDITIONER_HPP_
 
#include "linearAlgebra/DofManager.hpp"
#include "linearAlgebra/common/PreconditionerBase.hpp"
#include "linearAlgebra/utilities/BlockOperator.hpp"
#include "linearAlgebra/utilities/BlockVector.hpp"
 
namespace geos
{
 
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 >::setup;
 
  virtual void setup( Matrix const & mat ) 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 geos
 
#endif //GEOS_LINEARALGEBRA_SOLVERS_BLOCKPRECONDITIONER_HPP_