ConstrainMatrixCtx.hpp

Go to the documentation of this file.

/** \file ConstrainMatrixCtx.hpp
 *
 * Matrix-free projection and lifting operations for a supplied constraint
 * matrix and constraint-space inverse action.
 */
 
#ifndef MOFEM_CONSTRAIN_MATRIX_CTX_HPP
  #define MOFEM_CONSTRAIN_MATRIX_CTX_HPP
 
namespace MoFEM {
 
/**
 * \brief API for matrices that project onto, or lift from, a constraint space.
 * \ingroup projection_matrix
 */
struct ConstrainMatrixCtx {
 
  virtual ~ConstrainMatrixCtx() = default;
 
  /**
   * Set the constraint matrix C and its row-space inverse action atomically.
   *
   * @p inverse_cct represents (C C^T)^-1 on the selected constraint space, or
   * the corresponding embedded generalized inverse for a masked full space.
   * It must have the row layout of C in both dimensions and support MatMult().
   * MatMultTranspose() is additionally required by applyRT() and
   * applyCTC_QTKQ(). Both PETSc objects are retained.
   *
   * All data derived from the previous pair are invalidated. For in-place
   * value changes which preserve layouts and sparsity, update @p inverse_cct
   * in the same transaction and call recalculateCT() and recalculateCTC().
   */
  virtual MoFEMErrorCode setConstrainMatrix(Mat constrain_matrix,
                                            Mat inverse_cct) = 0;
 
  /** Get the constraint matrix C. */
  virtual Mat getConstrainMatrix() const = 0;
 
  /** Get the supplied inverse action for C C^T. */
  virtual Mat getInverseCCT() const = 0;
 
  /**
   * Set the matrix K used by C^T C + Q^T K Q and retain a PETSc reference.
   * Work vectors derived from the previous K are invalidated.
   */
  virtual MoFEMErrorCode setSystemMatrix(Mat system_matrix) = 0;
 
  /** Get the matrix K used by C^T C + Q^T K Q. */
  virtual Mat getSystemMatrix() const = 0;
 
  /**
   * Initialize matrices, vectors, and the scatter used by Q, P, R, and R^T.
   * The operation is idempotent until the C/inverse-action pair is replaced.
   *
   * @param x Vector defining the full problem layout.
   */
  virtual MoFEMErrorCode initializeQorP(Vec x) = 0;
 
  /**
   * Initialize matrices and work vectors used by C^T C + Q^T K Q.
   * The operation is idempotent until C or K is replaced.
   */
  virtual MoFEMErrorCode initializeQTKQ() = 0;
 
  /**
   * Recalculate C^T after an in-place value change in C.
   * The layout and nonzero structure of C must be unchanged. The caller must
   * update the supplied inverse action in the same transaction. Projection
   * MATSHELL objects already created from this context do not receive a PETSc
   * state increment; recreate them before passing the changed operator to a
   * consumer which caches by matrix state.
   */
  virtual MoFEMErrorCode recalculateCT() = 0;
 
  /**
   * Recalculate C^T C after an in-place value change in C.
   * The layout and nonzero structure of C must be unchanged.
   */
  virtual MoFEMErrorCode recalculateCTC() = 0;
 
  /** Apply Q = I - C^T S C for the supplied inverse action S. */
  virtual MoFEMErrorCode applyQ(Vec x, Vec f) = 0;
 
  /** Apply P = C^T S C for the supplied inverse action S. */
  virtual MoFEMErrorCode applyP(Vec x, Vec f) = 0;
 
  /** Apply R = C^T S for the supplied inverse action S. */
  virtual MoFEMErrorCode applyR(Vec x, Vec f) = 0;
 
  /** Apply R^T = S^T C for the supplied inverse action S. */
  virtual MoFEMErrorCode applyRT(Vec x, Vec f) = 0;
 
  /** Apply C^T C + Q^T K Q. */
  virtual MoFEMErrorCode applyCTC_QTKQ(Vec x, Vec f) = 0;
};
 
/**
 * Create a context whose scatter is built from two DMs.
 *
 * @param x_dm DM owning the full vector layout.
 * @param y_dm DM whose column numbering matches the constraint matrix.
 * @param constrain_matrix Constraint matrix C.
 * @param inverse_cct Supplied inverse action for C C^T.
 */
boost::shared_ptr<ConstrainMatrixCtx>
createConstrainMatrixCtx(DM x_dm, DM y_dm, Mat constrain_matrix,
                         Mat inverse_cct);
 
/**
 * Create a context using a supplied scatter.
 *
 * The communicator is obtained from the scatter.
 *
 * @param scatter Scatter between the full and constrained vector layouts.
 * @param constrain_matrix Constraint matrix C.
 * @param inverse_cct Supplied inverse action for C C^T.
 */
boost::shared_ptr<ConstrainMatrixCtx>
createConstrainMatrixCtx(VecScatter scatter, Mat constrain_matrix,
                         Mat inverse_cct);
 
/**
 * Create a KSP-backed shell action for (C C^T)^-1.
 *
 * The returned square MATSHELL retains @p constrain_matrix. It refreshes
 * C C^T and the KSP when values of C change, applies KSPSolve() in MatMult(),
 * and applies KSPSolveTranspose() in MatMultTranspose(). A change of C layout
 * requires a new shell matrix. Updating C through this backend does not
 * refresh C^T cached by a ConstrainMatrixCtx; call recalculateCT() and, when
 * needed, recalculateCTC() on that context in the same update transaction.
 * A zero-dimensional constraint space is applied as the zero action without
 * constructing a KSP.
 *
 * If @p ksp is nullptr, the shell creates and configures a KSP using
 * @p options_prefix. The default gives options such as
 * `-constraint_matrix_ksp_type` and `-constraint_matrix_pc_type`. A supplied
 * KSP remains caller-configured except that every multiply enforces a zero
 * initial guess, as required for PETSc MatMult overwrite semantics.
 *
 * @param constrain_matrix Constraint matrix C.
 * @param ksp Optional caller-configured KSP.
 * @param options_prefix Prefix for an internally created KSP.
 */
SmartPetscObj<Mat>
createCCTInverseKSP(Mat constrain_matrix, KSP ksp = nullptr,
                    const std::string &options_prefix = "constraint_matrix_");
 
/**
 * Create a lifetime-safe shell matrix for Q.
 *
 * The returned PETSc matrix retains @p ctx until the matrix is destroyed.
 * @param ctx Constraint-matrix context.
 * @param x Vector defining the full problem layout.
 */
SmartPetscObj<Mat>
createProjectionMatrixQ(const boost::shared_ptr<ConstrainMatrixCtx> &ctx,
                        Vec x);
 
/** Create a lifetime-safe shell matrix for P. */
SmartPetscObj<Mat>
createProjectionMatrixP(const boost::shared_ptr<ConstrainMatrixCtx> &ctx,
                        Vec x);
 
/** Create a lifetime-safe shell matrix for R. */
SmartPetscObj<Mat>
createProjectionMatrixR(const boost::shared_ptr<ConstrainMatrixCtx> &ctx,
                        Vec x);
 
/** Create a lifetime-safe shell matrix for R^T. */
SmartPetscObj<Mat>
createProjectionMatrixRT(const boost::shared_ptr<ConstrainMatrixCtx> &ctx,
                         Vec x);
 
/** Create a lifetime-safe shell matrix for C^T C + Q^T K Q. */
SmartPetscObj<Mat> createProjectionMatrixCTCPlusQTKQ(
    const boost::shared_ptr<ConstrainMatrixCtx> &ctx, Vec x);
 
} // namespace MoFEM
 
#endif // MOFEM_CONSTRAIN_MATRIX_CTX_HPP
 
/**
 \defgroup projection_matrix Constraint Projection Matrix
 \ingroup mofem_petsc_solvers
**/