SetUpSchurImpl.cpp

Go to the documentation of this file.

/** @file SetUpSchurImpl
 * @brief
 * @date 2023-05-13
 *
 * @license{This project is released under the MIT License.}
 *
 */
 
struct SetUpSchurImpl : public EshelbianCore::SetUpSchur {
 
  SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj<Mat> m,
                 SmartPetscObj<Mat> p, EshelbianCore *ep_core_ptr)
      : SetUpSchur(), mField(m_field), M(m), P(p), epCorePtr(ep_core_ptr) {
    if (S) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Is expected that schur matrix is not allocated. This is "
          "possible only is PC is set up twice");
    }
  }
  virtual ~SetUpSchurImpl() { S.reset(); }
 
  MoFEMErrorCode setUp(KSP solver);
  MoFEMErrorCode preProc();
  MoFEMErrorCode postProc();
 
private:
  SmartPetscObj<Mat> M;
  SmartPetscObj<Mat> P;
  SmartPetscObj<Mat> S;
  SmartPetscObj<AO> aoUp;
  SmartPetscObj<IS> isA00;
  SmartPetscObj<IS> isShift;
 
  MoFEM::Interface &mField;
  EshelbianCore *epCorePtr;
};
 
MoFEMErrorCode SetUpSchurImpl::setUp(KSP solver) {
  MoFEMFunctionBegin;
 
  auto get_ents_by_dim = [&](int dim) {
    Range ents;
    CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents);
    return ents;
  };
 
  std::vector<std::string> schur_field_list{
 
      epCorePtr->stretchTensor, 
      
      epCorePtr->bubbleField,
 
      epCorePtr->rotAxis,       
 
      epCorePtr->spatialL2Disp,
 
      epCorePtr->piolaStress
 
  };
  std::vector<std::string> field_list{epCorePtr->piolaStress,
                                      epCorePtr->contactDisp};
 
  auto create_schur_dm = [&](SmartPetscObj<DM> &dm_sub) {
    MoFEMFunctionBegin;
 
    dm_sub = createDM(mField.get_comm(), "DMMOFEM");
    CHKERR DMMoFEMCreateSubDM(dm_sub, epCorePtr->dmElastic, "SUB_SCHUR");
    CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
    CHKERR DMMoFEMSetIsPartitioned(dm_sub, PETSC_TRUE);
    CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->elementVolumeName);
    CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->naturalBcElement);
    CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->contactElement);
 
    auto faces_ptr = boost::make_shared<Range>(get_ents_by_dim(SPACE_DIM - 1));
    std::vector<boost::shared_ptr<Range>> dm_range_list{faces_ptr, nullptr};
 
    if (field_list.size() != dm_range_list.size()) {
      SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
              "Both ranges should have the same size");
    }
 
    int r_idx = 0;
    for (auto f : field_list) {
      MOFEM_LOG("EP", Sev::inform) << "Add schur field: " << f;
      CHKERR DMMoFEMAddSubFieldRow(dm_sub, f, dm_range_list[r_idx]);
      CHKERR DMMoFEMAddSubFieldCol(dm_sub, f, dm_range_list[r_idx]);
      ++r_idx;
    }
    CHKERR DMSetUp(dm_sub);
    MoFEMFunctionReturn(0);
  };
 
  auto create_a00_is = [&](SmartPetscObj<IS> &is_a00) {
    MoFEMFunctionBegin;
    auto vols = get_ents_by_dim(SPACE_DIM);
    std::vector<SmartPetscObj<IS>> is_vec;
    std::vector<Range *> range_list_ptr(schur_field_list.size(), nullptr);
    range_list_ptr.back() = &vols;
 
    int r_idx = 0;
    for (auto f : schur_field_list) {
      SmartPetscObj<IS> is;
      CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
          "ELASTIC_PROBLEM", ROW, f, 0, MAX_DOFS_ON_ENTITY, is,
          range_list_ptr[r_idx]);
      is_vec.push_back(is);
      ++r_idx;
    }
    if (is_vec.size()) {
      is_a00 = is_vec[0];
      for (auto i = 1; i < is_vec.size(); ++i) {
        IS is_union_raw;
        CHKERR ISExpand(is_a00, is_vec[i], &is_union_raw);
        is_a00 = SmartPetscObj<IS>(is_union_raw);
      }
      CHKERR ISSort(is_a00);
    }
    MoFEMFunctionReturn(0);
  };
 
  PC pc;
  CHKERR KSPSetFromOptions(solver);
  CHKERR KSPGetPC(solver, &pc);
  PetscBool is_pcfs = PETSC_FALSE;
  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
  if (is_pcfs) {
 
    SmartPetscObj<DM> schur_dm, a00_dm;
    CHKERR create_schur_dm(schur_dm);
    CHKERR create_a00_is(isA00);
 
    if (S) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Is expected that schur matrix is not allocated. This is "
          "possible only is PC is set up twice");
    }
 
    S = createDMMatrix(schur_dm);
 
    auto set_ops = [&]() {
      MoFEMFunctionBegin;
      auto range_ptr = boost::make_shared<Range>(get_ents_by_dim(SPACE_DIM));
      std::vector<boost::shared_ptr<Range>> ranges_list(schur_field_list.size(),
                                                        nullptr);
      ranges_list.back() = range_ptr;
      std::vector<SmartPetscObj<AO>> ao_list(schur_field_list.size(),
                                             SmartPetscObj<AO>());
      std::vector<SmartPetscObj<Mat>> mat_list(schur_field_list.size(),
                                               SmartPetscObj<Mat>());
      std::vector<bool> symm_list(schur_field_list.size(), false);
 
      auto dm_is = getDMSubData(schur_dm)->getSmartRowIs();
      aoUp = createAOMappingIS(dm_is, PETSC_NULL);
      ao_list.back() = aoUp;
      mat_list.back() = S;
 
      const bool symmetric_system =
          EshelbianCore::gradApperoximator <= MODERATE_ROT &&
          EshelbianCore::rotSelector == SMALL_ROT;
      epCorePtr->elasticFeLhs->getOpPtrVector().push_front(
          new OpSchurAssembleBegin());
      epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
          new OpSchurAssembleEnd<SCHUR_DGESV>(schur_field_list, ranges_list,
                                              ao_list, mat_list, symm_list,
                                              symmetric_system));
      epCorePtr->elasticBcLhs->getOpPtrVector().push_front(
          new OpSchurAssembleBegin());
      epCorePtr->elasticBcLhs->getOpPtrVector().push_back(
          new OpSchurAssembleEnd<SCHUR_DGESV>(schur_field_list, ranges_list,
                                              ao_list, mat_list, symm_list,
                                              symmetric_system));
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_pc = [&]() {
      MoFEMFunctionBegin;
      CHKERR PCFieldSplitSetIS(pc, NULL, isA00);
      CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_ops();
    CHKERR set_pc();
 
  } else {
    epCorePtr->elasticFeLhs->getOpPtrVector().push_front(
        new OpSchurAssembleBegin());
    epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
        new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
    epCorePtr->elasticBcLhs->getOpPtrVector().push_front(
        new OpSchurAssembleBegin());
    epCorePtr->elasticBcLhs->getOpPtrVector().push_back(
        new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode SetUpSchurImpl::preProc() {
  MoFEMFunctionBegin;
  MOFEM_LOG("EP", Sev::verbose) << "Zero Schur";
  if (SetUpSchurImpl::S) {
    CHKERR MatZeroEntries(S);
  }
  CHKERR MatZeroEntries(M);
  CHKERR MatZeroEntries(P);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode SetUpSchurImpl::postProc() {
  MoFEMFunctionBegin;
  if (S) {
    MOFEM_LOG("EP", Sev::verbose) << "Assemble Schur";
    CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
  }
  CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
  CHKERR MatAXPY(P, 1, M, SAME_NONZERO_PATTERN);
 
  MoFEMFunctionReturn(0);
}
 
boost::shared_ptr<EshelbianCore::SetUpSchur>
EshelbianCore::SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field,
                                            SmartPetscObj<Mat> m,
                                            SmartPetscObj<Mat> p,
                                            EshelbianCore *ep_core_ptr) {
  return boost::shared_ptr<SetUpSchur>(
      new SetUpSchurImpl(m_field, m, p, ep_core_ptr));
}