Schur.hpp

Go to the documentation of this file.

/**
 * @file Schur.hpp
 * @brief Assemble Schur complement
 * @date 2023-02-02
 *
 * @copyright Copyright (c) 2023
 *
 * To create nested system of Schur complements, you push sequence of operator,
 * to set up data on entities, and then assemble complements.
 *
 * \note Try separate floating points operations from book keeping. Also, align
 * memory that blocks follow floating point operations.
 *
 */
 
#ifndef __SCHUR_HPP__
#define __SCHUR_HPP__
 
namespace MoFEM {
 
constexpr const char MoFEM_BLOCK_MAT[] = "mofem_block_mat";
 
/**
 * @brief Structure to register events for Schur block assembly and solver
 */
struct SchurEvents {
  static PetscLogEvent MOFEM_EVENT_schurMatSetValues;
  static PetscLogEvent MOFEM_EVENT_opSchurAssembleEnd;
  static PetscLogEvent MOFEM_EVENT_BlockStructureSetValues;
  static PetscLogEvent MOFEM_EVENT_BlockStructureMult;
  static PetscLogEvent MOFEM_EVENT_BlockStructureSolve;
  static PetscLogEvent MOFEM_EVENT_AssembleSchurMat;
  static PetscLogEvent MOFEM_EVENT_zeroRowsAndCols;
  SchurEvents();
};
 
struct SchurElemMats;
struct SchurElemMatsBlock;
struct SchurElemMatsPreconditionedBlock;
 
struct OpSchurAssembleBase : public ForcesAndSourcesCore::UserDataOperator {
 
  OpSchurAssembleBase() = delete;
 
private:
  using UserDataOperator::UserDataOperator;
};
 
OpSchurAssembleBase *createOpSchurAssembleBegin();
 
/**
 * @brief Construct a new Op Schur Assemble End object
 *
 * @param fields_name list of fields (can be empty)
 * @param field_ents list of entities on which schur complement is applied
 * @param schur_aos maps dofs indices from main problem to schur complement
 * @param schur_mat schur matrix
 * @param sym_schur true if schur (matrix) complement is symmetric
 * @param symm_op true if block diagonal is symmetric
 */
OpSchurAssembleBase *
createOpSchurAssembleEnd(std::vector<std::string> fields_name,
                         std::vector<boost::shared_ptr<Range>> field_ents,
                         SmartPetscObj<AO> ao = SmartPetscObj<AO>(),
                         SmartPetscObj<Mat> schur = SmartPetscObj<Mat>(),
                         bool sym_schur = false, bool symm_op = false);
 
/**
 * @brief Create a Op Schur Zero Rows And Cols object
 * 
 * @param marker_ptr 
 * @return OpSchurAssembleBase* 
 */
OpSchurAssembleBase *createOpSchurZeroRowsAndCols(
    boost::shared_ptr<std::vector<unsigned char>> marker_ptr, double diag_val);
 
using SchurFieldPair = std::pair<std::string, std::string>;
 
using SchurFEOpsFEandFields = std::vector<
 
    std::pair<std::string, std::vector<SchurFieldPair>>
 
    >;
 
/**
 * @brief Create block matrix structure
 *
 * For each listed element and each listed pair of fields, one dense block is
 * allocated per pair of (row entity, column entity). Blocks are stored
 * row-major, values are shared by block matrix and by nested Schur matrices
 * created with createSchurNestedMatrixStruture. Block values are allocated on
 * first use.
 *
 * @param dm problem DM
 * @param schur_fe_op_vec list of elements with pairs of (row, column) fields
 * @return boost::shared_ptr<BlockStructure>
 */
boost::shared_ptr<BlockStructure> createBlockMatStructure(
 
    DM dm,                                //< dm
    SchurFEOpsFEandFields schur_fe_op_vec //< block elements
 
);
 
using SchurShellMatData =
    std::pair<SmartPetscObj<Mat>, boost::shared_ptr<BlockStructure>>;
 
/**
 * @brief Create a Schur Mat object
 *
 * @param dm
 * @param data
 * @return std::pair<SmartPetscObj<Mat>, boost::shared_ptr<BlockStructure>>
 */
SchurShellMatData createBlockMat(DM dm, boost::shared_ptr<BlockStructure> data);
 
using NestSchurData = std::tuple<
 
    std::array<SmartPetscObj<Mat>, 4>,
    std::array<boost::shared_ptr<BlockStructure>, 4>,
    boost::shared_ptr<BlockStructure>,
    std::pair<SmartPetscObj<IS>, SmartPetscObj<IS>>
 
    >;
 
/**
 * @brief Create structure of nested matrix [A11 A10; A01 A00]
 *
 * Storage contract:
 *
 * - Blocks of A00 (both entities in block DM) of one finite element are
 *   stored in a buffer of that element. Eliminated fields have to be element
 *   local (L2, broken spaces, interior DOFs); an entity shared by two elements
 *   is an error. Elements with the same block pattern (entities, numbers of
 *   DOFs, listed pairs) share the layout of their buffers.
 * - Buffers start with listed blocks only. At the first factorisation each
 *   pattern is analysed on a probe element: entities whose diagonal block is
 *   an admissible pivot are eliminated block by block (block LU in place, fill
 *   only where the pattern creates it), the remaining entities form a dense
 *   core factorised with partial pivoting. Buffers are then relaid in place.
 *   The analysis uses values only, it does not depend on the application.
 * - Block factors of every element are accepted only if the backward error of
 *   a solve with them is small; otherwise the element is factorised dense in
 *   a separate buffer.
 * - Factors are computed by assembleBlockMatSchur or by first MatSolve on the
 *   A00 matrix, and are reused by MatMult, MatSolve and assembleBlockMatSchur
 *   until MatZeroEntries is called. Assembly or MatZeroRowsColumns on
 *   factorised matrix is an error, i.e. each assembly has to start with
 *   MatZeroEntries.
 * - If add_preconditioner_block is true, preconditioner buffers are allocated
 *   only for A00 elements, and factors of (A00 + preconditioner) are stored in
 *   place of preconditioner blocks, A00 blocks are kept intact.
 *   Preconditioner can be assembled only to A00 blocks.
 * - MatSolve is implemented only for the A00 matrix (last matrix in the
 *   array). Processes without A00 elements perform empty solve.
 * - Element matrices not assembled on this process (copies of elements owned
 *   by other processes) are zero and are skipped.
 * - Layout of shared storage is fixed by this function, it can be called only
 *   once for given block structure.
 *
 * Options:
 * - -schur_a00_dense: dense element matrices, no block elimination
 * - -schur_a00_rcond_min (1e-12): smallest reciprocal condition number of
 *   diagonal block admissible as pivot
 * - -schur_a00_pivot_threshold (0, off): block threshold pivoting, smallest
 *   singular value of pivot block relative to its couplings
 * - -schur_a00_backward_error (1e-12): accepted backward error of solve with
 *   block factors (or 100 times the dense backward error on probe)
 *
 * See doc/development/schur.md.
 *
 * @param dms schur dm, and block dm
 * @param block mat A data
 * @param fields_name list of fields
 * @param field_ents list of entities on which schur complement is applied
 * @param add_preconditioner_block add block for preconditioner
 * @return boost::shared_ptr<NestSchurData>
 */
boost::shared_ptr<NestSchurData> createSchurNestedMatrixStruture(
    std::pair<SmartPetscObj<DM>, SmartPetscObj<DM>> dms,
    boost::shared_ptr<BlockStructure> block_mat_data,
 
    std::vector<std::string> fields_name,             //< a00 fields
    std::vector<boost::shared_ptr<Range>> field_ents, //< a00 ents
    bool add_preconditioner_block = false
 
);
 
/**
 * @brief Assemble Schur complement S = A11 - A10 A00^{-1} A01 from block
 * matrix
 *
 * A00 is eliminated element by element. If fields_name and field_ents
 * coincide with A00 of nested matrix (see createSchurNestedMatrixStruture),
 * A00 element matrices are factorised in place and factors are reused by
 * MatSolve. Otherwise element matrices are factorised in scratch memory, and
 * it is an error if A00 is already factorised in place.
 *
 * @param m_field
 * @param B block matrix
 * @param S Schur complement matrix (PETSc matrix)
 * @param fields_name eliminated fields
 * @param field_ents entities of eliminated fields (can be null)
 * @param ao mapping from block matrix indices to S indices
 * @return MoFEMErrorCode
 */
MoFEMErrorCode
assembleBlockMatSchur(MoFEM::Interface &m_field, Mat B, Mat S,
                      std::vector<std::string> fields_name,
                      std::vector<boost::shared_ptr<Range>> field_ents,
                      SmartPetscObj<AO> ao);
 
/**
 * @brief Get block values storage
 *
 * @note After createSchurNestedMatrixStruture it holds blocks outside A00
 * only, A00 blocks are stored in element buffers. Direct access bypasses
 * factorisation guard.
 *
 * @return boost::shared_ptr<std::vector<double>>
 */
boost::shared_ptr<std::vector<double>> getBlockMatStorageMat(Mat B);
 
/**
 * @brief Get preconditioner blocks storage
 *
 * @deprecated Preconditioner blocks are stored in buffers of A00 elements,
 * this function returns an empty vector.
 *
 * @return boost::shared_ptr<std::vector<double>>
 */
boost::shared_ptr<std::vector<double>>
getBlockMatPrconditionerStorageMat(Mat B);
 
/**
 * @brief Switch preconditioner
 * 
 * @param block_mat_data 
 * @return MoFEMErrorCode 
 */
MoFEMErrorCode
schurSwitchPreconditioner(boost::shared_ptr<BlockStructure> block_mat_data);
 
/**
 * @brief Save block matrix as a mesh
 * 
 * @param block_mat_data 
 * @param filename 
 * @return MoFEMErrorCode 
 */
MoFEMErrorCode
schurSaveBlockMesh(boost::shared_ptr<BlockStructure> block_mat_data,
                   std::string filename);
 
/**
 * @brief Create a Mat Diag Blocks object
 *
 * \code {.cpp}
 *
 * auto [nested_mat, nested_data_ptr] = createSchurNestedMatrix(
 *
 *       createSchurNestedMatrixStruture(
 *
 *           {schur_dm, block_dm}, shell_data,
 *
 *           {"TENSOR"}, {nullptr}
 *
 *       )
 *
 *  );
 *
 * \endcode
 *
 * @return Mat
 */
std::pair<SmartPetscObj<Mat>, boost::shared_ptr<NestSchurData>>
createSchurNestedMatrix(boost::shared_ptr<NestSchurData> schur_net_data_ptr);
 
template <>
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr<NestSchurData>);
 
/**
 * @brief Set PC for A00 block
 *
 * @param pc
 * @return MoFEMErrorCode
 */
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj<PC> pc);
 
struct SchurBackendMatSetValuesPtr {
  SchurBackendMatSetValuesPtr() = delete;
  using MatSetValuesPtr = boost::function<MoFEMErrorCode(
      Mat M, const EntitiesFieldData::EntData &row_data,
      const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
      InsertMode iora)>;
  static MatSetValuesPtr matSetValuesPtr; ///< backend assembly function
  static MatSetValuesPtr
      matSetValuesBlockPtr; ///< backend assembly block mat function
  static MatSetValuesPtr
      matSetValuesPreconditionedBlockPtr; ///< backend assembly block
                                          ///< preconditioner mat function
};
 
template <>
MoFEMErrorCode
MatSetValues<SchurElemMats>(Mat M, const EntitiesFieldData::EntData &row_data,
                            const EntitiesFieldData::EntData &col_data,
                            const MatrixDouble &mat, InsertMode iora);
 
template <>
inline MoFEMErrorCode
VecSetValues<SchurElemMats>(Vec V, const EntitiesFieldData::EntData &data,
                            const VectorDouble &nf, InsertMode iora) {
  return VecSetValues<EssentialBcStorage>(V, data, nf, iora);
}
 
template <>
MoFEMErrorCode MatSetValues<AssemblyTypeSelector<SCHUR>>(
    Mat M, const EntitiesFieldData::EntData &row_data,
    const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
    InsertMode iora);
 
template <>
inline MoFEMErrorCode VecSetValues<AssemblyTypeSelector<SCHUR>>(
    Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf,
    InsertMode iora) {
  return VecSetValues<SchurElemMats>(V, data, nf, iora);
}
 
/***
 * @brief Specialization of MatSetValues for BlockStructure
 */
template <>
MoFEMErrorCode
MatSetValues<BlockStructure>(Mat M, const EntitiesFieldData::EntData &row_data,
                            const EntitiesFieldData::EntData &col_data,
                            const MatrixDouble &mat, InsertMode iora);
 
/***
 * @brief Specialisation of MatSetValues for AssemblyTypeSelector<BLOCK_MAT>
 */
template <>
inline MoFEMErrorCode MatSetValues<AssemblyTypeSelector<BLOCK_MAT>>(
    Mat M, const EntitiesFieldData::EntData &row_data,
    const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
    InsertMode iora) {
  return MatSetValues<BlockStructure>(M, row_data, col_data, mat, iora);
}
 
/***
 * @brief Specialisation of VecSetValues for AssemblyTypeSelector<BLOCK_MAT>
 */
template <>
inline MoFEMErrorCode VecSetValues<AssemblyTypeSelector<BLOCK_MAT>>(
    Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf,
    InsertMode iora) {
  return VecSetValues<EssentialBcStorage>(V, data, nf, iora);
}
 
/***
 * @brief Specialization of MatSetValues for SchurElemMatsBlock
 */
template <>
MoFEMErrorCode
MatSetValues<SchurElemMatsBlock>(Mat M,
                                 const EntitiesFieldData::EntData &row_data,
                                 const EntitiesFieldData::EntData &col_data,
                                 const MatrixDouble &mat, InsertMode iora);
 
/***
 * @brief Specialization of VecSetValues for SchurElemMatsBlock
 */
template <>
inline MoFEMErrorCode
VecSetValues<SchurElemMatsBlock>(Vec V, const EntitiesFieldData::EntData &data,
                                 const VectorDouble &nf, InsertMode iora) {
  return VecSetValues<EssentialBcStorage>(V, data, nf, iora);
}
 
/***
 * @brief Specialisation of MatSetValues for AssemblyTypeSelector<BLOCK_SCHUR>
 */
template <>
inline MoFEMErrorCode MatSetValues<AssemblyTypeSelector<BLOCK_SCHUR>>(
    Mat M, const EntitiesFieldData::EntData &row_data,
    const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
    InsertMode iora) {
  return MatSetValues<SchurElemMatsBlock>(M, row_data, col_data, mat, iora);
}
 
/***
 * @brief Specialisation of VecSetValues for AssemblyTypeSelector<BLOCK_SCHUR>
 */
template <>
inline MoFEMErrorCode VecSetValues<AssemblyTypeSelector<BLOCK_SCHUR>>(
    Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf,
    InsertMode iora) {
  return VecSetValues<EssentialBcStorage>(V, data, nf, iora);
}
 
/***
 * @brief Specialization of MatSetValues for SchurElemMatsPreconditionedBlock
 */
template <>
MoFEMErrorCode MatSetValues<SchurElemMatsPreconditionedBlock>(
    Mat M, const EntitiesFieldData::EntData &row_data,
    const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
    InsertMode iora);
 
/***
 * @brief Specialisation of MatSetValues for
 * AssemblyTypeSelector<BLOCK_PRECONDITIONER_SCHUR>
 */
template <>
inline MoFEMErrorCode
MatSetValues<AssemblyTypeSelector<BLOCK_PRECONDITIONER_SCHUR>>(
    Mat M, const EntitiesFieldData::EntData &row_data,
    const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
    InsertMode iora) {
  return MatSetValues<SchurElemMatsPreconditionedBlock>(M, row_data, col_data,
                                                        mat, iora);
}
 
/***
 * @brief Specialisation of VecSetValues for AssemblyTypeSelector<BLOCK_SCHUR>
 */
template <>
inline MoFEMErrorCode
VecSetValues<AssemblyTypeSelector<BLOCK_PRECONDITIONER_SCHUR>>(
    Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf,
    InsertMode iora) {
  return ::VecSetValues(V, data.getIndices().size(),
                        &*data.getIndices().begin(), &*nf.begin(), iora);
}
 
/**
 * @deprecated do not use
 */
DEPRECATED OpSchurAssembleBase *createOpSchurAssembleEnd(
    std::vector<std::string> fields_name,
    std::vector<boost::shared_ptr<Range>> field_ents,
    std::vector<SmartPetscObj<AO>> sequence_of_aos,
    std::vector<SmartPetscObj<Mat>> sequence_of_mats,
    std::vector<bool> sym_schur, bool symm_op,
    boost::shared_ptr<BlockStructure> diag_blocks = nullptr);
 
/**
 * @deprecated do not use
 */
DEPRECATED OpSchurAssembleBase *createOpSchurAssembleEnd(
    std::vector<std::string> fields_name,
    std::vector<boost::shared_ptr<Range>> field_ents,
    std::vector<SmartPetscObj<AO>> sequence_of_aos,
    std::vector<SmartPetscObj<Mat>> sequence_of_mats,
    std::vector<bool> sym_schur, std::vector<double> diag_eps, bool symm_op,
    boost::shared_ptr<BlockStructure> diag_blocks = nullptr);
 
/** @deprecated do not use, use createSchurNestedMatrixStruture instead */
DEPRECATED inline boost::shared_ptr<NestSchurData>
getNestSchurData(std::pair<SmartPetscObj<DM>, SmartPetscObj<DM>> dms,
                 boost::shared_ptr<BlockStructure> block_mat_data,
 
                 std::vector<std::string> fields_name, //< a00 fields
                 std::vector<boost::shared_ptr<Range>> field_ents, //< a00 ents
                 bool add_preconditioner_block = false
 
) {
  return createSchurNestedMatrixStruture(dms, block_mat_data, fields_name,
                                         field_ents, add_preconditioner_block);
}
 
} // namespace MoFEM
 
#endif //__SCHUR_HPP__