218                                                      {
  219  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
  222  dm_field->aO = SmartPetscObj<AO>(ao, true);
  224}
  225 
  227                                                       bool create_sub_matrix,
  228                                                       bool shell_sub_a) {
  229  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
  232  dm_field->coarseningIS.emplace_back(SmartPetscObj<IS>(is, true));
  233  dm_field->shellMatrixCtxPtr.push_back(
new PCMGSubMatrixCtx(
A, is));
 
  234  if (is) {
  235    auto is2 = SmartPetscObj<IS>(is, true);
  236    if (dm_field->aO) {
  239      CHKERR AOApplicationToPetscIS(dm_field->aO, is2);
 
  240    }
  241    if (create_sub_matrix) {
  242      if (shell_sub_a) {
  246        MPI_Comm comm;
  247        CHKERR PetscObjectGetComm((PetscObject)
A, &comm);
 
  248        Mat sub_a_raw;
  250                              &(dm_field->shellMatrixCtxPtr.back()),
  251                              &sub_a_raw);
  252        CHKERR MatShellSetOperation(sub_a_raw, MATOP_MULT,
 
  253                                    (void (*)(void))sub_mat_mult);
  254        CHKERR MatShellSetOperation(sub_a_raw, MATOP_MULT_ADD,
 
  255                                    (void (*)(void))sub_mat_mult_add);
  256        CHKERR MatShellSetOperation(sub_a_raw, MATOP_SOR,
 
  257                                    (void (*)(void))sub_mat_sor);
  258        dm_field->kspOperators.emplace_back(
  259            SmartPetscObj<Mat>(sub_a_raw, false));
  260      } else {
  261        Mat sub_a_raw;
  262#if PETSC_VERSION_GE(3, 8, 0)
  263        CHKERR MatCreateSubMatrix(
A, is2, is2, MAT_INITIAL_MATRIX, &sub_a_raw);
 
  264#else
  265        CHKERR MatGetSubMatrix(
A, is2, is2, MAT_INITIAL_MATRIX, &sub_a_raw);
 
  266#endif
  267        dm_field->kspOperators.emplace_back(
  268            SmartPetscObj<Mat>(sub_a_raw, false));
  269      }
  270    } else {
  271      dm_field->kspOperators.emplace_back(SmartPetscObj<Mat>(
A, 
true));
 
  272    }
  273  } else {
  275  }
  276  PetscInfo(dm, "Push back IS to DMMGViaApproxOrders\n");
  278}
  279 
  281  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
  284  if (dm_field->coarseningIS.back()) {
  285    dm_field->coarseningIS.pop_back();
  286  }
  287  dm_field->kspOperators.pop_back();
  288  PetscInfo(dm, "Pop back IS to DMMGViaApproxOrders\n");
  290}
  291 
  293  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
  296  CHKERR dm_field->destroyCoarseningIS();
 
  297  PetscInfo(dm, "Clear DMs data structures\n");
  299}
  300 
  303  CHKERR DMRegister(sname, DMCreate_MGViaApproxOrders);
 
  305}
  306 
  310}
  311 
  313  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
  315  if (!dm->data) {
  316    dm->data = new DMMGViaApproxOrdersCtx();
  317  } else {
  318    ((DMCtxImpl *)(dm->data))->incrementReference();
  319  }
  321 
  327 
  328  CHKERR DMKSPSetComputeOperators(dm, ksp_set_operators, NULL);
 
  329#if PETSC_VERSION_GE(3, 17, 0)
  330  PetscInfo(dm, "Create DMMGViaApproxOrders reference = %d\n",
  331            ((DMCtxImpl *)(dm->data))->useCount());
  332#else
  333  PetscInfo1(dm, "Create DMMGViaApproxOrders reference = %d\n",
  334             ((DMCtxImpl *)(dm->data))->useCount());
  335#endif
  337}
  338 
  340  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
  343  CHKERR dm_field->destroyCoarseningIS();
 
  346}
  347 
  349 
  350  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
  353 
  354  int leveldown = dm->leveldown;
  355 
  356  if (dm_field->kspOperators.empty()) {
  358  } else {
  359    MPI_Comm comm;
  360    CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
 
  361    if (dm_field->kspOperators.empty()) {
  363              "data inconsistency, operator can not be set");
  364    }
  365    if (static_cast<int>(dm_field->kspOperators.size()) < leveldown) {
  367              "data inconsistency, no IS for that level");
  368    }
  369    *
M = dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown];
 
  370    CHKERR PetscObjectReference((PetscObject)*M);
 
  372  }
  373 
  374#if PETSC_VERSION_GE(3, 17, 0)
  375  PetscInfo(dm, "Create Matrix DMMGViaApproxOrders leveldown = %d\n",
  376             leveldown);
  377#else
  378  PetscInfo1(dm, "Create Matrix DMMGViaApproxOrders leveldown = %d\n",
  379             leveldown);
  380#endif
  381 
  383}
  384 
  386  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
  389  CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
 
  390  CHKERR DMCreate(comm, dmc);
 
  391  (*dmc)->data = dm->data;
  393  CHKERR DMGetType(dm, &type);
 
  394  CHKERR DMSetType(*dmc, type);
 
  395  dm_field->coarsenDM = SmartPetscObj<DM>(*dmc, false);
  396#if PETSC_VERSION_GE(3, 17, 0)
  397  PetscInfo(dm, "Coarsen DMMGViaApproxOrders leveldown = %d\n", dm->leveldown);
  398#else
  399  PetscInfo1(dm, "Coarsen DMMGViaApproxOrders leveldown = %d\n", dm->leveldown);
  400#endif
  402}
  403 
  405                                                       Vec *vec) {
  406  PetscValidHeaderSpecific(dm1, DM_CLASSID, 1);
  407  PetscValidHeaderSpecific(dm2, DM_CLASSID, 1);
  409 
  410  MPI_Comm comm;
  411  CHKERR PetscObjectGetComm((PetscObject)dm1, &comm);
 
  412 
  414 
  415  DM dm_down = dm1;
  416  DM dm_up = dm2;
  417 
  418  int dm_down_leveldown = dm_down->leveldown;
  419  int dm_up_leveldown = dm_up->leveldown;
  420#if PETSC_VERSION_GE(3, 17, 0)
  421  PetscInfo(dm1,
  422            "Create interpolation DMMGViaApproxOrders dm1_leveldown = %d "
  423            "dm2_leveldown = %d\n",
  424            dm_down_leveldown, dm_up_leveldown);
  425#else
  426  PetscInfo2(dm1,
  427             "Create interpolation DMMGViaApproxOrders dm1_leveldown = %d "
  428             "dm2_leveldown = %d\n",
  429             dm_down_leveldown, dm_up_leveldown);
  430#endif
  431 
  432  IS is_down, is_up;
  433  {
  434    
  436    if (static_cast<int>(dm_field->coarseningIS.size()) < dm_down_leveldown) {
  438    }
  439    is_down = dm_field->coarseningIS[dm_field->coarseningIS.size() - 1 -
  440                                     dm_down_leveldown];
  441    CHKERR ISGetSize(is_down, &M);
 
  442    CHKERR ISGetLocalSize(is_down, &
m);
 
  443  }
  444  {
  445    
  447    if (static_cast<int>(dm_field->coarseningIS.size()) < dm_up_leveldown) {
  449    }
  450    is_up =
  451        dm_field
  452            ->coarseningIS[dm_field->coarseningIS.size() - 1 - dm_up_leveldown];
  454    CHKERR ISGetLocalSize(is_up, &
n);
 
  455  }
  456 
  457  
  458  
  459 
  460  CHKERR MatCreate(comm, mat);
 
  462  CHKERR MatSetType(*mat, MATMPIAIJ);
 
  463  CHKERR MatMPIAIJSetPreallocation(*mat, 1, PETSC_NULLPTR, 0, PETSC_NULLPTR);
 
  464 
  465  
  466  PetscLayout rmap, cmap;
  467  CHKERR MatGetLayouts(*mat, &rmap, &cmap);
 
  468  int rstart, rend, cstart, cend;
  469  CHKERR PetscLayoutGetRange(rmap, &rstart, &rend);
 
  470  CHKERR PetscLayoutGetRange(cmap, &cstart, &cend);
 
  471 
  472  const int *row_indices_ptr, *col_indices_ptr;
  473  CHKERR ISGetIndices(is_down, &row_indices_ptr);
 
  474  CHKERR ISGetIndices(is_up, &col_indices_ptr);
 
  475 
  476  map<int, int> idx_map;
  477  for (
int ii = 0; ii < 
m; ii++) {
 
  478    idx_map[row_indices_ptr[ii]] = rstart + ii;
  479  }
  480 
  481  CHKERR MatZeroEntries(*mat);
 
  482  
  483  for (
int jj = 0; jj < 
n; jj++) {
 
  484    map<int, int>::iterator mit = idx_map.find(col_indices_ptr[jj]);
  485    if (mit != idx_map.end()) {
  486      CHKERR MatSetValue(*mat, mit->second, cstart + jj, 1, INSERT_VALUES);
 
  487    }
  488  }
  489 
  490  CHKERR ISRestoreIndices(is_down, &row_indices_ptr);
 
  491  CHKERR ISRestoreIndices(is_up, &col_indices_ptr);
 
  492 
  493  CHKERR MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
 
  494  CHKERR MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);
 
  495 
  496  if (vec != NULL) {
  497    *vec = PETSC_NULLPTR;
  498  }
  499 
  501}
  502 
  504  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
  506  int leveldown = dm->leveldown;
  508  if (dm_field->kspOperators.empty()) {
  510  } else {
  511#if PETSC_VERSION_GE(3, 5, 3)
  513        dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
  515#else
  517        dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
  519#endif
  521  }
  522#if PETSC_VERSION_GE(3, 17, 0)
  523  PetscInfo(dm, "Create global vector DMMGViaApproxOrders leveldown = %d\n",
  524            dm->leveldown);
  525#else
  526  PetscInfo1(dm, "Create global vector DMMGViaApproxOrders leveldown = %d\n",
  527             dm->leveldown);
  528#endif
  530}
  531 
  534  PetscOptionsBegin(PETSC_COMM_WORLD, "",
  535                           "MOFEM Multi-Grid (Orders) pre-conditioner", "none");
  536 
  537  CHKERR PetscOptionsInt(
"-mofem_mg_levels", 
"nb levels of multi-grid solver",
 
  538                         "", 2, &nbLevels, PETSC_NULLPTR);
  539  CHKERR PetscOptionsInt(
"-mofem_mg_coarse_order",
 
  540                         "approximation order of coarse level", "", 2,
  541                         &coarseOrder, PETSC_NULLPTR);
  542  CHKERR PetscOptionsInt(
"-mofem_mg_order_at_last_level", 
"order at last level",
 
  543                         "", 100, &orderAtLastLevel, PETSC_NULLPTR);
  544  CHKERR PetscOptionsInt(
"-mofem_mg_verbose", 
"nb levels of multi-grid solver",
 
  545                         "", 0, &verboseLevel, PETSC_NULLPTR);
  546  PetscBool shell_sub_a = shellSubA ? PETSC_TRUE : PETSC_FALSE;
  547  CHKERR PetscOptionsBool(
"-mofem_mg_shell_a", 
"use shell matrix as sub matrix",
 
  548                          "", shell_sub_a, &shell_sub_a, NULL);
  549  shellSubA = (shellSubA == PETSC_TRUE);
  550 
  551  PetscOptionsEnd();
  553}
  554 
  555MoFEMErrorCode PCMGSetUpViaApproxOrdersCtx::createIsAtLevel(
int kk, IS *is) {
 
  559  
  562  const Problem *problem_ptr;
  564  int order_at_next_level = kk + coarseOrder;
  565  if (kk == nbLevels - 1) {
  566    int first = problem_ptr->getNumeredRowDofsPtr()
  567                    ->get<PetscGlobalIdx_mi_tag>()
  568                    .lower_bound(0)
  569                    ->get()
  570                    ->getPetscGlobalDofIdx();
  571    CHKERR ISCreateStride(PETSC_COMM_WORLD, problem_ptr->getNbLocalDofsRow(),
 
  572                          first, 1, is);
  574    
  575  }
  576  string problem_name = problem_ptr->getName();
  578                                              order_at_next_level, is);
  580}
  581 
  582MoFEMErrorCode PCMGSetUpViaApproxOrdersCtx::destroyIsAtLevel(
int kk, IS *is) {
 
  586}
  587 
  589PCMGSetUpViaApproxOrdersCtx::buildProlongationOperator(bool use_mat_a,
  590                                                       int verb) {
  592  verb = verb > verboseLevel ? verb : verboseLevel;
  593 
  594  MPI_Comm comm;
  595  CHKERR PetscObjectGetComm((PetscObject)dM, &comm);
 
  596 
  597  MOFEM_LOG_C(
"PCMGViaApproximationOrders", Sev::inform, 
"set MG levels %u\n",
 
  598              nbLevels);
  599 
  602 
  603  std::vector<IS> is_vec(nbLevels + 1);
 
  604  std::vector<
int> is_glob_size(nbLevels + 1), is_loc_size(nbLevels + 1);
 
  605 
  606  for (int kk = 0; kk < nbLevels; kk++) {
  607 
  608    
  609    CHKERR createIsAtLevel(kk, &is_vec[kk]);
 
  610    CHKERR ISGetSize(is_vec[kk], &is_glob_size[kk]);
 
  611    CHKERR ISGetLocalSize(is_vec[kk], &is_loc_size[kk]);
 
  612 
  614                "Nb. dofs at level [ %d ] global %u local %d\n", kk,
  615                is_glob_size[kk], is_loc_size[kk]);
  616 
  617    
  618    if (is_glob_size[kk] == 0) {
  620    }
  621  }
  622 
  623  for (int kk = 0; kk != nbLevels; kk++) {
  624    if (kk == nbLevels - 1 && use_mat_a) {
  626                                                     false);
  627    } else {
  628      if (kk > 0) {
  629        
  631                                                       shellSubA);
  632      } else {
  633        
  634        
  636                                                       false);
  637      }
  638    }
  639  }
  640 
  641  for (unsigned int kk = 0; kk < is_vec.size(); kk++) {
  642    CHKERR destroyIsAtLevel(kk, &is_vec[kk]);
 
  643  }
  644 
  647 
  649}
  650 
  651boost::shared_ptr<PCMGSetUpViaApproxOrdersCtx>
  653  return boost::make_shared<PCMGSetUpViaApproxOrdersCtx>(dm, 
A, use_shell_mat);
 
  654}
  655 
  657    PC pc, boost::shared_ptr<PCMGSetUpViaApproxOrdersCtx> ctx, int verb) {
  658 
  659  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  661 
  662  auto core_log = logging::core::get();
  663  core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(),
  664                                            "PCMGViaApproximationOrders"));
  665  LogManager::setLog("PCMGViaApproximationOrders");
  666  MOFEM_LOG_TAG(
"PCMGViaApproximationOrders", 
"PCMGViaApproximationOrders");
 
  667 
  668  MOFEM_LOG(
"PCMGViaApproximationOrders", Sev::verbose)
 
  669      << "Setup PCMGSetUpViaApproxOrders";
  670 
  672  CHKERR ctx->buildProlongationOperator(
true, verb);
 
  673 
  674#if PETSC_VERSION_GE(3, 8, 0)
  675  CHKERR PCMGSetGalerkin(pc, PC_MG_GALERKIN_NONE);
 
  676#else
  677  CHKERR PCMGSetGalerkin(pc, PETSC_FALSE);
 
  678#endif
  679 
  680  CHKERR PCMGSetLevels(pc, ctx->nbLevels, NULL);
 
  681 
  682  MOFEM_LOG(
"PCMGViaApproximationOrders", Sev::noisy)
 
  683      << "Setup PCMGSetUpViaApproxOrders <- end";
  684 
  686}
  687 
  688} 
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_LOG_C(channel, severity, format,...)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
MoFEMErrorCode DMMGViaApproxOrdersPopBackCoarseningIS(DM dm)
Pop IS form MG via approximation orders.
PetscErrorCode DMDestroy_MoFEM(DM dm)
Destroys dm with MoFEM data structure.
MoFEMErrorCode DMMGViaApproxOrdersPushBackCoarseningIS(DM dm, IS is, Mat A, bool create_sub_matrix, bool shell_sub_a)
Push back coarsening level to MG via approximation orders.
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
MoFEMErrorCode DMCoarsen_MGViaApproxOrders(DM dm, MPI_Comm comm, DM *dmc)
Coarsen DM.
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
MoFEMErrorCode DMMGViaApproxOrdersClearCoarseningIS(DM dm)
Clear approximation orders.
MoFEMErrorCode DMCreateMatrix_MGViaApproxOrders(DM dm, Mat *M)
Create matrix for Multi-Grid via approximation orders.
PetscErrorCode DMSetOperators_MoFEM(DM dm)
Set operators for MoFEM dm.
MoFEMErrorCode isCreateProblemOrder(const std::string problem_name, RowColData rc, int min_order, int max_order, IS *is) const
create IS for given order range (collective)
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
const double n
refractive index of diffusive medium
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
MoFEMErrorCode DMCreateInterpolation_MGViaApproxOrders(DM dm1, DM dm2, Mat *mat, Vec *vec)
Create interpolation matrix between data managers dm1 and dm2.
static MoFEMErrorCode ksp_set_operators(KSP ksp, Mat A, Mat B, void *ctx)
MoFEMErrorCode DMCreateGlobalVector_MGViaApproxOrders(DM dm, Vec *g)
Create global vector for DMGViaApproxOrders.
boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > createPCMGSetUpViaApproxOrdersCtx(DM dm, Mat A, bool use_shell_mat)
createPCMGSetUpViaApproxOrdersCtx
MoFEMErrorCode DMCreate_MGViaApproxOrders(DM dm)
Create DM data structure for Multi-Grid via approximation orders.
PetscErrorCode DMDestroy_MGViaApproxOrders(DM dm)
Destroy DM.
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > ctx, int verb)
Function build MG structure.
FTensor::Index< 'M', 3 > M
FTensor::Index< 'm', 3 > m
Deprecated interface functions.
Section manager is used to create indexes and sections.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.