KspCtx.cpp

Go to the documentation of this file.

/**
 * @brief Context for PETSc KSP, i.e. linear solver
 * @author Anonymous contributor under MIT license
 *
 */
 
namespace MoFEM {
 
struct KspCtx::KspCtx::KspCtxImpl {
 
  KspCtxImpl(MoFEM::Interface &m_field, std::string problem_name);
 
  virtual ~KspCtxImpl() = default;
 
  /**
   * @return return reference to vector with FEMethod to calculate matrix
   */
  FEMethodsSequence &getSetOperators() { return loopsOperator; }
 
  /**
   * @return return vector to vector with FEMethod to vector
   */
  FEMethodsSequence &getComputeRhs() { return loopsRhs; }
 
  /**
   * The sequence of BasicMethod is executed before residual is calculated. It
   * can be used to setup data structures, e.g. zero global variable which is
   * integrated in domain, e.g. for calculation of strain energy.
   *
   * @return reference to BasicMethod for preprocessing
   */
  BasicMethodsSequence &getPreProcComputeRhs() { return preProcessRhs; }
 
  /**
   * The sequence of BasicMethod is executed after residual is calculated. It
   * can be used to setup data structures, e.g. aggregate data from processors
   * or to apply essential boundary conditions.
   *
   * @return reference to BasicMethod for postprocessing
   */
  BasicMethodsSequence &getPostProcComputeRhs() { return postProcessRhs; }
 
  /**
   * @return reference to BasicMethod for preprocessing
   */
  BasicMethodsSequence &getPreProcSetOperators() { return preProcessOperator; }
 
  /**
   * The sequence of BasicMethod is executed after tangent matrix is calculated.
   * It can be used to setup data structures, e.g. aggregate data from
   * processors or to apply essential boundary conditions.
   *
   * @return reference to BasicMethod for postprocessing
   */
  BasicMethodsSequence &getPostProcSetOperators() {
    return postProcessOperator;
  }
 
  friend PetscErrorCode KspRhs(KSP ksp, Vec f, void *ctx);
  friend PetscErrorCode KspMat(KSP ksp, Mat A, Mat B, void *ctx);
 
  /**
   * @brief Clear loops
   *
   * @return MoFEMErrorCode
   */
  MoFEMErrorCode clearLoops();
 
protected:
  PetscLogEvent MOFEM_EVENT_KspRhs;
  PetscLogEvent MOFEM_EVENT_KspMat;
 
  boost::movelib::unique_ptr<bool> vecAssembleSwitch;
  boost::movelib::unique_ptr<bool> matAssembleSwitch;
 
  MoFEM::Interface &mField;
  moab::Interface &moab;
 
  std::string problemName; ///< Problem name
  MoFEMTypes bH;           ///< If set to MF_EXIST check if element exist
 
  FEMethodsSequence loopsOperator; ///< Sequence of finite elements instances
                                   ///< assembling tangent matrix
  FEMethodsSequence loopsRhs;      ///< Sequence of finite elements
                                   ///< instances assembling residual vector
  BasicMethodsSequence preProcessOperator;  ///< Sequence of methods run before
                                            ///< tangent matrix is assembled
  BasicMethodsSequence postProcessOperator; ///< Sequence of methods run after
                                            ///< tangent matrix is assembled
  BasicMethodsSequence preProcessRhs;       ///< Sequence of methods run before
                                            ///< residual is assembled
  BasicMethodsSequence postProcessRhs;      ///< Sequence of methods run after
                                            ///< residual is assembled
};
 
KspCtx::KspCtx(MoFEM::Interface &m_field, std::string problem_name)
    : kspCtxImpl(boost::movelib::make_unique<KspCtx::KspCtxImpl>(
          m_field, problem_name)) {}
 
KspCtx::~KspCtx() = default;
 
FEMethodsSequence &KspCtx::getSetOperators() {
  return kspCtxImpl->getSetOperators();
}
 
FEMethodsSequence &KspCtx::getComputeRhs() {
  return kspCtxImpl->getComputeRhs();
}
 
BasicMethodsSequence &KspCtx::getPreProcComputeRhs() {
  return kspCtxImpl->getPreProcComputeRhs();
}
 
BasicMethodsSequence &KspCtx::getPostProcComputeRhs() {
  return kspCtxImpl->getPostProcComputeRhs();
}
 
BasicMethodsSequence &KspCtx::getPreProcSetOperators() {
  return kspCtxImpl->getPreProcSetOperators();
}
 
BasicMethodsSequence &KspCtx::getPostProcSetOperators() {
  return kspCtxImpl->getPostProcSetOperators();
}
 
MoFEMErrorCode KspCtx::clearLoops() { return kspCtxImpl->clearLoops(); }
 
KspCtx::KspCtxImpl::KspCtxImpl(MoFEM::Interface &m_field,
                               std::string problem_name)
    : mField(m_field), moab(m_field.get_moab()), problemName(problem_name),
      bH(MF_EXIST) {
  PetscLogEventRegister("LoopKSPRhs", 0, &MOFEM_EVENT_KspRhs);
  PetscLogEventRegister("LoopKSPMat", 0, &MOFEM_EVENT_KspMat);
}
 
MoFEMErrorCode KspCtx::KspCtxImpl::clearLoops() {
  MoFEMFunctionBegin;
  loopsOperator.clear();
  loopsRhs.clear();
  preProcessOperator.clear();
  postProcessOperator.clear();
  preProcessRhs.clear();
  postProcessRhs.clear();
  MoFEMFunctionReturn(0);
}
 
PetscErrorCode KspRhs(KSP ksp, Vec f, void *ctx) {
  MoFEMFunctionBegin;
  KspCtx::KspCtxImpl *ksp_ctx = static_cast<KspCtx *>(ctx)->kspCtxImpl.get();
  PetscLogEventBegin(ksp_ctx->MOFEM_EVENT_KspRhs, 0, 0, 0, 0);
 
  ksp_ctx->vecAssembleSwitch = boost::movelib::make_unique<bool>(true);
 
  auto cache_ptr = boost::make_shared<CacheTuple>();
  CHKERR ksp_ctx->mField.cache_problem_entities(ksp_ctx->problemName,
                                                cache_ptr);
 
  auto set = [&](auto &fe) {
    fe.ksp = ksp;
    fe.ksp_ctx = KspMethod::CTX_SETFUNCTION;
    fe.data_ctx = PetscData::CtxSetF;
    fe.ksp_f = f;
    fe.cacheWeakPtr = cache_ptr;
  };
 
  auto unset = [&](auto &fe) {
    fe.ksp_ctx = KspMethod::CTX_KSPNONE;
    fe.data_ctx = PetscData::CtxSetNone;
  };
 
  // pre-process
  for (auto &bit : ksp_ctx->preProcessRhs) {
    bit->vecAssembleSwitch = boost::move(ksp_ctx->vecAssembleSwitch);
    set(*bit);
    CHKERR ksp_ctx->mField.problem_basic_method_preProcess(ksp_ctx->problemName,
                                                           *bit);
    unset(*bit);
    ksp_ctx->vecAssembleSwitch = boost::move(bit->vecAssembleSwitch);
  }
 
  // operators
  for (auto &lit : ksp_ctx->loopsRhs) {
    lit.second->vecAssembleSwitch = boost::move(ksp_ctx->vecAssembleSwitch);
    set(*lit.second);
    CHKERR ksp_ctx->mField.loop_finite_elements(ksp_ctx->problemName, lit.first,
                                                *(lit.second), nullptr,
                                                ksp_ctx->bH, cache_ptr);
    unset(*lit.second);
    ksp_ctx->vecAssembleSwitch = boost::move(lit.second->vecAssembleSwitch);
  }
 
  // post-process
  for (auto &bit : ksp_ctx->postProcessRhs) {
    bit->vecAssembleSwitch = boost::move(ksp_ctx->vecAssembleSwitch);
    set(*bit);
    CHKERR ksp_ctx->mField.problem_basic_method_postProcess(
        ksp_ctx->problemName, *bit);
    unset(*bit);
    ksp_ctx->vecAssembleSwitch = boost::move(bit->vecAssembleSwitch);
  }
 
  if (*ksp_ctx->vecAssembleSwitch) {
    CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(f);
    CHKERR VecAssemblyEnd(f);
  }
  PetscLogEventEnd(ksp_ctx->MOFEM_EVENT_KspRhs, 0, 0, 0, 0);
  MoFEMFunctionReturn(0);
}
PetscErrorCode KspMat(KSP ksp, Mat A, Mat B, void *ctx) {
  // PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
  MoFEMFunctionBegin;
  KspCtx::KspCtxImpl *ksp_ctx = static_cast<KspCtx *>(ctx)->kspCtxImpl.get();
  PetscLogEventBegin(ksp_ctx->MOFEM_EVENT_KspMat, 0, 0, 0, 0);
 
  ksp_ctx->matAssembleSwitch = boost::movelib::make_unique<bool>(true);
  auto cache_ptr = boost::make_shared<CacheTuple>();
  CHKERR ksp_ctx->mField.cache_problem_entities(ksp_ctx->problemName,
                                                cache_ptr);
 
  auto set = [&](auto &fe) {
    fe.ksp = ksp;
    fe.ksp_A = A;
    fe.ksp_B = B;
    fe.ksp_ctx = KspMethod::CTX_OPERATORS;
    fe.data_ctx = PetscData::CtxSetA | PetscData::CtxSetB;
    fe.cacheWeakPtr = cache_ptr;
  };
 
  auto unset = [&](auto &fe) {
    fe.ksp_ctx = KspMethod::CTX_KSPNONE;
    fe.data_ctx = PetscData::CtxSetNone;
  };
 
  auto ent_data_cache = boost::make_shared<std::vector<EntityCacheDofs>>();
  auto ent_row_cache =
      boost::make_shared<std::vector<EntityCacheNumeredDofs>>();
  auto ent_col_cache =
      boost::make_shared<std::vector<EntityCacheNumeredDofs>>();
 
  // pre-procsess
  for (auto &bit : ksp_ctx->preProcessOperator) {
    bit->matAssembleSwitch = boost::move(ksp_ctx->matAssembleSwitch);
    set(*bit);
    CHKERR ksp_ctx->mField.problem_basic_method_preProcess(ksp_ctx->problemName,
                                                           *bit);
    unset(*bit);
    ksp_ctx->matAssembleSwitch = boost::move(bit->matAssembleSwitch);
  }
 
  // operators
  for (auto &lit : ksp_ctx->loopsOperator) {
    lit.second->matAssembleSwitch = boost::move(ksp_ctx->matAssembleSwitch);
    set(*lit.second);
    CHKERR ksp_ctx->mField.loop_finite_elements(ksp_ctx->problemName, lit.first,
                                                *(lit.second), nullptr,
                                                ksp_ctx->bH, cache_ptr);
    unset(*lit.second);
    ksp_ctx->matAssembleSwitch = boost::move(lit.second->matAssembleSwitch);
  }
 
  // post-process
  for (auto &bit : ksp_ctx->postProcessOperator) {
    bit->matAssembleSwitch = boost::move(ksp_ctx->matAssembleSwitch);
    set(*bit);
    CHKERR ksp_ctx->mField.problem_basic_method_postProcess(
        ksp_ctx->problemName, *bit);
    unset(*bit);
    ksp_ctx->matAssembleSwitch = boost::move(bit->matAssembleSwitch);
  }
 
  if (ksp_ctx->matAssembleSwitch) {
    CHKERR MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
  }
  PetscLogEventEnd(ksp_ctx->MOFEM_EVENT_KspMat, 0, 0, 0, 0);
  MoFEMFunctionReturn(0);
}
 
} // namespace MoFEM