schur_elastic.cpp

Go to the documentation of this file.

/**
 * @file schur_elastic.cpp
 * @example mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp
 *
 * @copyright Anonymous authors (c) 2025 under the MIT license (see LICENSE for
 * details)
 */
 
#include <MoFEM.hpp>
 
using namespace MoFEM;
 
[[maybe_unused]] constexpr int BASE_DIM = 1; //< Dimension of the base functions
 
//! [Define dimension]
constexpr int SPACE_DIM =
    EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
//! [Define dimension]
constexpr AssemblyType A = (SCHUR_ASSEMBLE)
                               ? AssemblyType::BLOCK_SCHUR
                               : AssemblyType::PETSC; //< selected assembly type
 
constexpr IntegrationType I =
    IntegrationType::GAUSS; //< selected integration type
 
//! [Define entities]
using EntData = EntitiesFieldData::EntData;
using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
using BoundaryEle =
    PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
using DomainEleOp = DomainEle::UserDataOperator;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
//! [Define entities]
 
struct DomainBCs {};
struct BoundaryBCs {};
 
using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<A>::LinearForm<I>;
using OpDomainRhsBCs = DomainRhsBCs::OpFlux<DomainBCs, 1, SPACE_DIM>;
using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<A>::LinearForm<I>;
using OpBoundaryRhsBCs = BoundaryRhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<A>::BiLinearForm<I>;
using OpBoundaryLhsBCs = BoundaryLhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
 
template <int DIM> struct PostProcEleByDim;
 
template <> struct PostProcEleByDim<2> {
  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
  using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
};
 
template <> struct PostProcEleByDim<3> {
  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
  using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
};
 
using PostProcEleDomain = PostProcEleByDim<SPACE_DIM>::PostProcEleDomain;
using SideEle = PostProcEleByDim<SPACE_DIM>::SideEle;
using PostProcEleBdy = PostProcEleByDim<SPACE_DIM>::PostProcEleBdy;
 
#include <ElasticExample.hpp>
 
struct ElasticSchurExample : public ElasticExample {
  ElasticSchurExample(MoFEM::Interface &field) : ElasticExample(field) {}
 
  MoFEMErrorCode runProblem();
 
protected:
  MoFEMErrorCode kspSetUpAndSolve(SmartPetscObj<KSP> solver);
};
 
//! [Run problem]
MoFEMErrorCode ElasticSchurExample::runProblem() {
  MoFEMFunctionBegin;
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  CHKERR outputResults();
  CHKERR checkResults();
  MoFEMFunctionReturn(0);
}
//! [Run problem]
 
struct SetUpSchur {
  static boost::shared_ptr<SetUpSchur>
  createSetUpSchur(MoFEM::Interface &m_field);
  virtual MoFEMErrorCode setUp(SmartPetscObj<KSP> solver) = 0;
 
protected:
  SetUpSchur() = default;
};
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
 
  // Initialisation of MoFEM/PETSc and MOAB data structures
  const char param_file[] = "param_file.petsc";
  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
 
  auto core_log = logging::core::get();
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
 
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmSync(), "FieldEvaluator"));
  LogManager::setLog("FieldEvaluator");
  MOFEM_LOG_TAG("FieldEvaluator", "field_eval");
 
  try {
 
    //! [Register MoFEM discrete manager in PETSc]
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
    DMType dm_name_mg = "DMMOFEM_MG";
    CHKERR DMRegister_MGViaApproxOrders(dm_name_mg);
    //! [Register MoFEM discrete manager in PETSc
 
    //! [Create MoAB]
    moab::Core mb_instance;              ///< mesh database
    moab::Interface &moab = mb_instance; ///< mesh database interface
    //! [Create MoAB]
 
    //! [Create MoFEM]
    MoFEM::Core core(moab);           ///< finite element database
    MoFEM::Interface &m_field = core; ///< finite element database interface
    //! [Create MoFEM]
 
    //! [Example]
    ElasticSchurExample ex(m_field);
    CHKERR ex.runProblem();
    //! [Example]
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
}
 
MoFEMErrorCode ElasticSchurExample::kspSetUpAndSolve(SmartPetscObj<KSP> solver) {
  MoFEMFunctionBegin;
 
  MOFEM_LOG_CHANNEL("TIMER");
  MOFEM_LOG_TAG("TIMER", "timer");
 
  if (A == AssemblyType::BLOCK_SCHUR || A == AssemblyType::SCHUR) {
    auto schur_ptr = SetUpSchur::createSetUpSchur(mField);
    CHKERR schur_ptr->setUp(solver);
  }
 
 
  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp";
  CHKERR KSPSetUp(solver);
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp <= Done";
 
  DM dm;
  CHKERR KSPGetDM(solver, &dm);
  auto D = createDMVector(dm);
  auto F = vectorDuplicate(D);
 
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve";
  CHKERR KSPSolve(solver, F, D);
  CHKERR mField.getInterface<TestingInterface>()->checkSolver(
      solver, "elasticity schur");
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve <= Done";
 
  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
  MoFEMFunctionReturn(0);
};
 
struct SetUpSchurImpl : public SetUpSchur {
 
  SetUpSchurImpl(MoFEM::Interface &m_field) : SetUpSchur(), mField(m_field) {
    if (S) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Is expected that schur matrix is not allocated. This is "
          "possible only is if PC is set up twice");
    }
  }
  virtual ~SetUpSchurImpl() = default;
 
  MoFEMErrorCode setUp(SmartPetscObj<KSP> solver);
 
private:
  MoFEMErrorCode setEntities();
  MoFEMErrorCode createSubDM();
  MoFEMErrorCode setOperator();
  MoFEMErrorCode setPC(PC pc);
  MoFEMErrorCode setDiagonalPC(PC pc);
 
  SmartPetscObj<Mat> S;
 
  MoFEM::Interface &mField;
 
  SmartPetscObj<DM> schurDM;
  SmartPetscObj<DM> blockDM;
  Range volEnts;
  Range subEnts;
};
 
MoFEMErrorCode SetUpSchurImpl::setUp(SmartPetscObj<KSP> solver) {
  MoFEMFunctionBegin;
  auto pip = mField.getInterface<PipelineManager>();
  PC pc;
  CHKERR KSPGetPC(solver, &pc);
  PetscBool is_pcfs = PETSC_FALSE;
  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
  if (is_pcfs) {
    if (S) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Is expected that schur matrix is not allocated. This is "
          "possible only is if PC is set up twice");
    }
    CHKERR setEntities();
    CHKERR createSubDM();
 
    // Add data to DM storage
    S = createDMMatrix(schurDM);
    CHKERR MatSetDM(S, PETSC_NULLPTR);
    CHKERR MatSetBlockSize(S, SPACE_DIM);
    CHKERR MatSetOption(S, MAT_SYMMETRIC, PETSC_TRUE);
 
    CHKERR setOperator();
    CHKERR setPC(pc);
 
    if constexpr (A == AssemblyType::BLOCK_SCHUR) {
      // Set DM to use shell block matrix
      DM solver_dm;
      CHKERR KSPGetDM(solver, &solver_dm);
      CHKERR DMSetMatType(solver_dm, MATSHELL);
    }
    CHKERR KSPSetUp(solver);
    CHKERR setDiagonalPC(pc);
 
  } else {
    pip->getOpBoundaryLhsPipeline().push_front(createOpSchurAssembleBegin());
    pip->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd({}, {}));
    pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
    pip->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd({}, {}));
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode SetUpSchurImpl::setEntities() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  CHKERR mField.get_moab().get_entities_by_dimension(simple->getMeshset(),
                                                     SPACE_DIM, volEnts);
  CHKERR mField.get_moab().get_entities_by_handle(simple->getMeshset(),
                                                  subEnts);
  subEnts = subtract(subEnts, volEnts);
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode SetUpSchurImpl::createSubDM() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  auto create_dm = [&](const char *name, auto &ents, auto dm_type) {
    auto dm = createDM(mField.get_comm(), dm_type);
    auto create_dm_imp = [&]() {
      MoFEMFunctionBegin;
      CHKERR DMMoFEMCreateSubDM(dm, simple->getDM(), name);
      CHKERR DMMoFEMSetSquareProblem(dm, PETSC_TRUE);
      CHKERR DMMoFEMAddElement(dm, simple->getDomainFEName());
      auto sub_ents_ptr = boost::make_shared<Range>(ents);
      CHKERR DMMoFEMAddSubFieldRow(dm, "U", sub_ents_ptr);
      CHKERR DMMoFEMAddSubFieldCol(dm, "U", sub_ents_ptr);
      CHKERR DMSetUp(dm);
      MoFEMFunctionReturn(0);
    };
    CHK_THROW_MESSAGE(
        create_dm_imp(),
        "Error in creating schurDM. It is possible that schurDM is "
        "already created");
    return dm;
  };
 
  schurDM = create_dm("SCHUR", subEnts, "DMMOFEM_MG");
  blockDM = create_dm("BLOCK", volEnts, "DMMOFEM");
 
  if constexpr (A == AssemblyType::BLOCK_SCHUR) {
 
    auto get_nested_mat_data = [&]() -> boost::shared_ptr<NestSchurData> {
      auto block_mat_data =
          createBlockMatStructure(simple->getDM(),
 
                                  {{
 
                                      simple->getDomainFEName(),
 
                                      {{"U", "U"}
 
                                      }}}
 
          );
 
      return createSchurNestedMatrixStruture(
 
          {schurDM, blockDM}, block_mat_data,
 
          {"U"}, {boost::make_shared<Range>(volEnts)}
 
      );
    };
 
    auto nested_mat_data = get_nested_mat_data();
 
    CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode SetUpSchurImpl::setOperator() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto pip = mField.getInterface<PipelineManager>();
 
  // Boundary
  auto dm_is = getDMSubData(schurDM)->getSmartRowIs();
  auto ao_up = createAOMappingIS(dm_is, PETSC_NULLPTR);
  pip->getOpBoundaryLhsPipeline().push_front(createOpSchurAssembleBegin());
  pip->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
      {"U"}, {boost::make_shared<Range>(volEnts)}, ao_up, S, true, true));
  // Domain
  pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
  pip->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
      {"U"}, {boost::make_shared<Range>(volEnts)}, ao_up, S, true, true));
 
  auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
  auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
 
  pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
    MoFEMFunctionBegin;
    if (S) {
      CHKERR MatZeroEntries(S);
    }
    MOFEM_LOG("TIMER", Sev::inform) << "Lhs Assemble Begin";
    MoFEMFunctionReturn(0);
  };
 
  post_proc_schur_lhs_ptr->postProcessHook = [this, post_proc_schur_lhs_ptr,
                                              ao_up]() {
    MoFEMFunctionBegin;
    CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
    CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
        mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
    MOFEM_LOG("TIMER", Sev::inform) << "Lhs Assemble End";
    MoFEMFunctionReturn(0);
  };
 
  auto ksp_ctx_ptr = getDMKspCtx(simple->getDM());
  ksp_ctx_ptr->getPreProcSetOperators().push_front(pre_proc_schur_lhs_ptr);
  ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_schur_lhs_ptr);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode SetUpSchurImpl::setPC(PC pc) {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  SmartPetscObj<IS> vol_is;
  mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
      simple->getProblemName(), ROW, "U", 0, SPACE_DIM, vol_is, &volEnts);
  CHKERR PCFieldSplitSetIS(pc, NULL, vol_is);
  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode SetUpSchurImpl::setDiagonalPC(PC pc) {
  MoFEMFunctionBegin;
 
  auto get_pc = [](auto ksp) {
    PC pc_raw;
    CHKERR KSPGetPC(ksp, &pc_raw);
    return SmartPetscObj<PC>(pc_raw, true); // bump reference
  };
 
  if constexpr (A == AssemblyType::BLOCK_SCHUR) {
    auto simple = mField.getInterface<Simple>();
    auto A = createDMBlockMat(simple->getDM());
    auto P = createDMNestSchurMat(simple->getDM());
    CHKERR PCSetOperators(pc, A, P);
 
    KSP *subksp;
    CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULLPTR, &subksp);
    CHKERR setSchurA00MatSolvePC(get_pc(subksp[0]));
 
    auto set_pc_p_mg = [](auto dm, auto pc, auto S) {
      MoFEMFunctionBegin;
      CHKERR PCSetDM(pc, dm);
      PetscBool same = PETSC_FALSE;
      PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
      if (same) {
        MOFEM_LOG("TIMER", Sev::inform) << "Set up MG";
        CHKERR PCMGSetUpViaApproxOrders(
            pc, createPCMGSetUpViaApproxOrdersCtx(dm, S, true));
        CHKERR PCSetFromOptions(pc);
      }
      MoFEMFunctionReturn(0);
    };
 
    auto set_pc_ksp = [&](auto dm, auto pc, auto S) {
      MoFEMFunctionBegin;
      PetscBool same = PETSC_FALSE;
      PetscObjectTypeCompare((PetscObject)pc, PCKSP, &same);
      if (same) {
        MOFEM_LOG("TIMER", Sev::inform) << "Set up inner KSP for PCKSP";
        CHKERR PCSetFromOptions(pc);
        KSP inner_ksp;
        CHKERR PCKSPGetKSP(pc, &inner_ksp);
        CHKERR KSPSetFromOptions(inner_ksp);
        PC inner_pc;
        CHKERR KSPGetPC(inner_ksp, &inner_pc);
        CHKERR PCSetFromOptions(inner_pc);
        CHKERR set_pc_p_mg(dm, inner_pc, S);
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_pc_ksp(schurDM, get_pc(subksp[1]), S);
    CHKERR set_pc_p_mg(schurDM, get_pc(subksp[1]), S);
 
    CHKERR PetscFree(subksp);
  }
  MoFEMFunctionReturn(0);
}
 
boost::shared_ptr<SetUpSchur>
SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field) {
  return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
}