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"
#include "linearAlgebra/utilities/LinearSolverParameters.hpp"
 
namespace geos
{
 
/*
 * 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( LinearSolverParameters::Block params );
 
  void setupBlock( localIndex const blockIndex,
                   stdVector< DofManager::SubComponent > blockDofs,
                   std::unique_ptr< PreconditionerBase< LAI > > solver,
                   real64 const scaling = 1.0 );
 
  void setupBlock( localIndex const blockIndex,
                   std::vector< DofManager::SubComponent > blockDofs,
                   PreconditionerBase< LAI > * const solver,
                   real64 const scaling = 1.0 );
 
  void setProlongation( localIndex const blockIndex,
                        Matrix const & P );
 
 
  virtual void setup( Matrix const & mat ) override;
 
  virtual void apply( Vector const & src, Vector & dst ) const override;
 
  virtual void clear() override;
 
 
  BlockOperator< Vector, Matrix > const & blocks() const
  {
    return m_matBlocks;
  }
 
private:
 
  void reinitialize( Matrix const & mat );
 
  void applyBlockScaling();
 
  void computeSchurComplement();
 
  LinearSolverParameters::Block m_params;
 
  std::array< stdVector< DofManager::SubComponent >, 2 > m_blockDofs{};
 
  std::array< Matrix const *, 2 > m_prolongators{};
 
  std::array< Matrix, 2 > m_prolongatorsOwned{};
 
  BlockOperator< Vector, Matrix > m_matBlocks;
 
  std::array< PreconditionerBase< LAI > *, 2 > m_solvers{};
 
  std::array< std::unique_ptr< PreconditionerBase< LAI > >, 2 > m_solversOwned{};
 
  std::array< real64, 2 > m_scaling{};
 
  mutable BlockVector< Vector > m_rhs;
 
  mutable BlockVector< Vector > m_sol;
};
 
} //namespace geos
 
#endif //GEOS_LINEARALGEBRA_SOLVERS_BLOCKPRECONDITIONER_HPP_