v0.14.0
PCMGSetUpViaApproxOrders.cpp
Go to the documentation of this file.
1 /** \file PCMGSetUpViaApproxOrders.cpp
2  * \brief implementation of multi-grid solver for p- adaptivity
3  */
4 
5 
6 #undef PETSC_VERSION_RELEASE
7 #define PETSC_VERSION_RELEASE 1
8 
9 #if PETSC_VERSION_GE(3, 6, 0)
10 #include <petsc/private/petscimpl.h>
11 #else
12 #include <petsc-private/petscimpl.h>
13 #endif
14 
15 #if PETSC_VERSION_GE(3, 6, 0)
16 #include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/
17 // #include <petsc/private/vecimpl.h> /*I "petscdm.h" I*/
18 #else
19 #include <petsc-private/dmimpl.h> /*I "petscdm.h" I*/
20 #include <petsc-private/vecimpl.h> /*I "petscdm.h" I*/
21 #endif
22 
23 #include <DMMoFEM.hpp>
24 #include <DMCtxImpl.hpp>
25 #include <PCMGSetUpViaApproxOrders.hpp>
26 
27 namespace MoFEM {
28 
30 
31  PCMGSubMatrixCtx(Mat a, IS is);
32  virtual ~PCMGSubMatrixCtx() = default;
33 
34  template <InsertMode MODE>
35  friend MoFEMErrorCode sub_mat_mult_generic(Mat a, Vec x, Vec f);
36  friend MoFEMErrorCode sub_mat_sor(Mat mat, Vec b, PetscReal omega,
37  MatSORType flag, PetscReal shift,
38  PetscInt its, PetscInt lits, Vec x);
39 
41 
47 
48  PetscLogEvent MOFEM_EVENT_mult;
49  PetscLogEvent MOFEM_EVENT_sor;
50 };
51 
52 PCMGSubMatrixCtx::PCMGSubMatrixCtx(Mat a, IS is) : A(a, true), iS(is, true) {
53  PetscLogEventRegister("PCMGSubMatrixCtx_mult", 0, &MOFEM_EVENT_mult);
54  PetscLogEventRegister("PCMGSubMatrixCtx_sor", 0, &MOFEM_EVENT_sor);
55 }
56 
59  if (sCat.use_count() == 0) {
60  auto [x_tmp, f_tmp] = matCreateVecs(A);
61  X = x_tmp;
62  F = f_tmp;
63  sCat = createVecScatter(X, iS, x, PETSC_NULL);
64  }
66 }
67 
68 /**
69  * \brief Structure for DM for multi-grid via approximation orders
70  * \ingroup dm
71  */
73 
74  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
75  MoFEM::UnknownInterface **iface) const;
76 
78  virtual ~DMMGViaApproxOrdersCtx();
79 
81 
84  std::vector<SmartPetscObj<IS>> coarseningIS; ///< Coarsening IS
85  std::vector<SmartPetscObj<Mat>> kspOperators; ///< Get KSP operators
86  boost::ptr_vector<PCMGSubMatrixCtx>
87  shellMatrixCtxPtr; ///< Shell sub-matrix context
88 };
89 
90 /* \brief Set data structures of MG pre-conditioner via approximation orders
91  */
93 
94  PCMGSetUpViaApproxOrdersCtx(DM dm, Mat a, bool shell_sub_a)
95  : dM(dm), A(a), nbLevels(2), coarseOrder(2), orderAtLastLevel(1000),
96  shellSubA(shell_sub_a), verboseLevel(0) {}
97 
98  virtual ~PCMGSetUpViaApproxOrdersCtx() = default;
99 
100  /**
101  * \brief get options from line command
102  * @return error code
103  */
104  virtual MoFEMErrorCode getOptions();
105 
106  /**
107  * \brief Set IS for levels
108  * @param kk level
109  * @param is pointer to IS
110  * @return error code
111  */
112  virtual MoFEMErrorCode createIsAtLevel(int kk, IS *is);
113 
114  /**
115  * \brief Destroy IS if internally created
116  * @param kk level
117  * @param is pointer to is
118  * @return error code
119  */
120  virtual MoFEMErrorCode destroyIsAtLevel(int kk, IS *is);
121 
122  /**
123  * \brief Set up data structures for MG
124  * @param pc MG pre-conditioner
125  * <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCMG.html>
126  * @param verb verbosity level
127  * @return error code
128  */
129  virtual MoFEMErrorCode buildProlongationOperator(bool use_mat_a,
130  int verb = 0);
131 
132  DM dM; ///< Distributed mesh manager
133  Mat A; ///< Matrix at fine level
134 
135  int nbLevels; ///< number of multi-grid levels
136  int coarseOrder; ///< approximation order of coarse level
137  int orderAtLastLevel; ///< set maximal evaluated order
138 
139  bool shellSubA;
141 };
142 
143 template <InsertMode MODE>
145  void *void_ctx;
147  CHKERR MatShellGetContext(a, &void_ctx);
148  PCMGSubMatrixCtx *ctx = (PCMGSubMatrixCtx *)void_ctx;
149  PetscLogEventBegin(ctx->MOFEM_EVENT_mult, 0, 0, 0, 0);
150  CHKERR ctx->initData(x);
151  CHKERR VecScatterBegin(ctx->sCat, x, ctx->X, INSERT_VALUES, SCATTER_REVERSE);
152  CHKERR VecScatterEnd(ctx->sCat, x, ctx->X, INSERT_VALUES, SCATTER_REVERSE);
153  CHKERR MatMult(ctx->A, ctx->X, ctx->F);
154  CHKERR VecScatterBegin(ctx->sCat, ctx->F, f, MODE, SCATTER_FORWARD);
155  CHKERR VecScatterEnd(ctx->sCat, ctx->F, f, MODE, SCATTER_FORWARD);
156  PetscLogEventEnd(ctx->MOFEM_EVENT_mult, 0, 0, 0, 0);
158 }
159 
161  return sub_mat_mult_generic<INSERT_VALUES>(a, x, f);
162 }
163 
165  return sub_mat_mult_generic<ADD_VALUES>(a, x, f);
166 }
167 
168 MoFEMErrorCode sub_mat_sor(Mat mat, Vec b, PetscReal omega, MatSORType flag,
169  PetscReal shift, PetscInt its, PetscInt lits,
170  Vec x) {
171 
172  //FIXME: that is crap implementation of SOR
173 
174  void *void_ctx;
176  CHKERR MatShellGetContext(mat, &void_ctx);
177  PCMGSubMatrixCtx *ctx = (PCMGSubMatrixCtx *)void_ctx;
178  PetscLogEventBegin(ctx->MOFEM_EVENT_sor, 0, 0, 0, 0);
179  CHKERR ctx->initData(x);
180  CHKERR VecScatterBegin(ctx->sCat, b, ctx->X, INSERT_VALUES, SCATTER_REVERSE);
181  CHKERR VecScatterEnd(ctx->sCat, b, ctx->X, INSERT_VALUES, SCATTER_REVERSE);
182  CHKERR MatSOR(ctx->A, ctx->X, omega, flag, shift, its, lits, ctx->F);
183  CHKERR VecScatterBegin(ctx->sCat, ctx->F, x, INSERT_VALUES, SCATTER_FORWARD);
184  CHKERR VecScatterEnd(ctx->sCat, ctx->F, x, INSERT_VALUES, SCATTER_FORWARD);
185  PetscLogEventEnd(ctx->MOFEM_EVENT_sor, 0, 0, 0, 0);
187 }
188 
190 
193  CHKERRABORT(PETSC_COMM_WORLD, ierr);
194 }
195 
198  coarsenDM.reset();
199  aO.reset();
200  coarseningIS.clear();
201  kspOperators.clear();
202  shellMatrixCtxPtr.clear();
204 }
205 
207 DMMGViaApproxOrdersCtx::query_interface(boost::typeindex::type_index type_index,
208  MoFEM::UnknownInterface **iface) const {
209  *iface = static_cast<DMMGViaApproxOrdersCtx *>(
210  const_cast<DMMGViaApproxOrdersCtx *>(this));
211  return 0;
212 }
213 
214 #define GET_DM_FIELD(DM) \
215  auto dm_field = \
216  static_cast<DMCtx *>(DM->data)->getInterface<DMMGViaApproxOrdersCtx>(); \
217  NOT_USED(dm_field)
218 
220  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
222  GET_DM_FIELD(dm);
223  dm_field->aO = SmartPetscObj<AO>(ao, true);
225 }
226 
228  bool create_sub_matrix,
229  bool shell_sub_a) {
230  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
232  GET_DM_FIELD(dm);
233  dm_field->coarseningIS.emplace_back(SmartPetscObj<IS>(is, true));
234  dm_field->shellMatrixCtxPtr.push_back(new PCMGSubMatrixCtx(A, is));
235  if (is) {
236  auto is2 = SmartPetscObj<IS>(is, true);
237  if (dm_field->aO) {
238  is2 = isDuplicate(is);
239  CHKERR ISCopy(is, is2);
240  CHKERR AOApplicationToPetscIS(dm_field->aO, is2);
241  }
242  if (create_sub_matrix) {
243  if (shell_sub_a) {
244  int n, N;
245  CHKERR ISGetSize(is, &N);
246  CHKERR ISGetLocalSize(is, &n);
247  MPI_Comm comm;
248  CHKERR PetscObjectGetComm((PetscObject)A, &comm);
249  Mat sub_a_raw;
250  CHKERR MatCreateShell(comm, n, n, N, N,
251  &(dm_field->shellMatrixCtxPtr.back()),
252  &sub_a_raw);
253  CHKERR MatShellSetOperation(sub_a_raw, MATOP_MULT,
254  (void (*)(void))sub_mat_mult);
255  CHKERR MatShellSetOperation(sub_a_raw, MATOP_MULT_ADD,
256  (void (*)(void))sub_mat_mult_add);
257  CHKERR MatShellSetOperation(sub_a_raw, MATOP_SOR,
258  (void (*)(void))sub_mat_sor);
259  dm_field->kspOperators.emplace_back(
260  SmartPetscObj<Mat>(sub_a_raw, false));
261  } else {
262  Mat sub_a_raw;
263 #if PETSC_VERSION_GE(3, 8, 0)
264  CHKERR MatCreateSubMatrix(A, is2, is2, MAT_INITIAL_MATRIX, &sub_a_raw);
265 #else
266  CHKERR MatGetSubMatrix(A, is2, is2, MAT_INITIAL_MATRIX, &sub_a_raw);
267 #endif
268  dm_field->kspOperators.emplace_back(
269  SmartPetscObj<Mat>(sub_a_raw, false));
270  }
271  } else {
272  dm_field->kspOperators.emplace_back(SmartPetscObj<Mat>(A, true));
273  }
274  } else {
275  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
276  }
277  PetscInfo(dm, "Push back IS to DMMGViaApproxOrders\n");
279 }
280 
282  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
284  GET_DM_FIELD(dm);
285  if (dm_field->coarseningIS.back()) {
286  dm_field->coarseningIS.pop_back();
287  }
288  dm_field->kspOperators.pop_back();
289  PetscInfo(dm, "Pop back IS to DMMGViaApproxOrders\n");
291 }
292 
294  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
296  GET_DM_FIELD(dm);
297  CHKERR dm_field->destroyCoarseningIS();
298  PetscInfo(dm, "Clear DMs data structures\n");
300 }
301 
304  CHKERR DMRegister(sname, DMCreate_MGViaApproxOrders);
306 }
307 
308 static MoFEMErrorCode ksp_set_operators(KSP ksp, Mat A, Mat B, void *ctx) {
311 }
312 
314  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
316  if (!dm->data) {
317  dm->data = new DMMGViaApproxOrdersCtx();
318  } else {
319  ((DMCtxImpl *)(dm->data))->incrementReference();
320  }
322 
323  dm->ops->creatematrix = DMCreateMatrix_MGViaApproxOrders;
324  dm->ops->createglobalvector = DMCreateGlobalVector_MGViaApproxOrders;
325  dm->ops->coarsen = DMCoarsen_MGViaApproxOrders;
326  dm->ops->createinterpolation = DMCreateInterpolation_MGViaApproxOrders;
327  dm->ops->destroy = DMDestroy_MGViaApproxOrders;
328 
329  CHKERR DMKSPSetComputeOperators(dm, ksp_set_operators, NULL);
330  PetscInfo1(dm, "Create DMMGViaApproxOrders reference = %d\n",
331  ((DMCtxImpl *)(dm->data))->useCount());
333 }
334 
335 PetscErrorCode DMDestroy_MGViaApproxOrders(DM dm) {
336  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
337  GET_DM_FIELD(dm);
339  CHKERR dm_field->destroyCoarseningIS();
342 }
343 
345 
346  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
348  GET_DM_FIELD(dm);
349 
350  int leveldown = dm->leveldown;
351 
352  if (dm_field->kspOperators.empty()) {
354  } else {
355  MPI_Comm comm;
356  CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
357  if (dm_field->kspOperators.empty()) {
358  SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
359  "data inconsistency, operator can not be set");
360  }
361  if (static_cast<int>(dm_field->kspOperators.size()) < leveldown) {
362  SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
363  "data inconsistency, no IS for that level");
364  }
365  *M = dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown];
366  CHKERR PetscObjectReference((PetscObject)*M);
367  CHKERR MatSetDM(*M, dm);
368  }
369 
370  PetscInfo1(dm, "Create Matrix DMMGViaApproxOrders leveldown = %d\n",
371  leveldown);
372 
374 }
375 
376 MoFEMErrorCode DMCoarsen_MGViaApproxOrders(DM dm, MPI_Comm comm, DM *dmc) {
377  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
379  GET_DM_FIELD(dm);
380  CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
381  CHKERR DMCreate(comm, dmc);
382  (*dmc)->data = dm->data;
383  DMType type;
384  CHKERR DMGetType(dm, &type);
385  CHKERR DMSetType(*dmc, type);
386  dm_field->coarsenDM = SmartPetscObj<DM>(*dmc, false);
387  PetscInfo1(dm, "Coarsen DMMGViaApproxOrders leveldown = %d\n", dm->leveldown);
389 }
390 
392  Vec *vec) {
393  PetscValidHeaderSpecific(dm1, DM_CLASSID, 1);
394  PetscValidHeaderSpecific(dm2, DM_CLASSID, 1);
396 
397  MPI_Comm comm;
398  CHKERR PetscObjectGetComm((PetscObject)dm1, &comm);
399 
400  int m, n, M, N;
401 
402  DM dm_down = dm1;
403  DM dm_up = dm2;
404 
405  int dm_down_leveldown = dm_down->leveldown;
406  int dm_up_leveldown = dm_up->leveldown;
407 
408  PetscInfo2(dm1,
409  "Create interpolation DMMGViaApproxOrders dm1_leveldown = %d "
410  "dm2_leveldown = %d\n",
411  dm_down_leveldown, dm_up_leveldown);
412 
413  IS is_down, is_up;
414  {
415  // Coarser mesh
416  GET_DM_FIELD(dm_down);
417  if (static_cast<int>(dm_field->coarseningIS.size()) < dm_down_leveldown) {
418  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
419  }
420  is_down = dm_field->coarseningIS[dm_field->coarseningIS.size() - 1 -
421  dm_down_leveldown];
422  CHKERR ISGetSize(is_down, &M);
423  CHKERR ISGetLocalSize(is_down, &m);
424  }
425  {
426  // Finer mesh
427  GET_DM_FIELD(dm_up);
428  if (static_cast<int>(dm_field->coarseningIS.size()) < dm_up_leveldown) {
429  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
430  }
431  is_up =
432  dm_field
433  ->coarseningIS[dm_field->coarseningIS.size() - 1 - dm_up_leveldown];
434  CHKERR ISGetSize(is_up, &N);
435  CHKERR ISGetLocalSize(is_up, &n);
436  }
437 
438  // is_dow rows
439  // is_up columns
440 
441  CHKERR MatCreate(comm, mat);
442  CHKERR MatSetSizes(*mat, m, n, M, N);
443  CHKERR MatSetType(*mat, MATMPIAIJ);
444  CHKERR MatMPIAIJSetPreallocation(*mat, 1, PETSC_NULL, 0, PETSC_NULL);
445 
446  // get matrix layout
447  PetscLayout rmap, cmap;
448  CHKERR MatGetLayouts(*mat, &rmap, &cmap);
449  int rstart, rend, cstart, cend;
450  CHKERR PetscLayoutGetRange(rmap, &rstart, &rend);
451  CHKERR PetscLayoutGetRange(cmap, &cstart, &cend);
452 
453  const int *row_indices_ptr, *col_indices_ptr;
454  CHKERR ISGetIndices(is_down, &row_indices_ptr);
455  CHKERR ISGetIndices(is_up, &col_indices_ptr);
456 
457  map<int, int> idx_map;
458  for (int ii = 0; ii < m; ii++) {
459  idx_map[row_indices_ptr[ii]] = rstart + ii;
460  }
461 
462  CHKERR MatZeroEntries(*mat);
463  // FIXME: Use MatCreateMPIAIJWithArrays and set array directly
464  for (int jj = 0; jj < n; jj++) {
465  map<int, int>::iterator mit = idx_map.find(col_indices_ptr[jj]);
466  if (mit != idx_map.end()) {
467  CHKERR MatSetValue(*mat, mit->second, cstart + jj, 1, INSERT_VALUES);
468  }
469  }
470 
471  CHKERR ISRestoreIndices(is_down, &row_indices_ptr);
472  CHKERR ISRestoreIndices(is_up, &col_indices_ptr);
473 
474  CHKERR MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
475  CHKERR MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);
476 
477  if (vec != NULL) {
478  *vec = PETSC_NULL;
479  }
480 
482 }
483 
485  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
487  int leveldown = dm->leveldown;
488  GET_DM_FIELD(dm);
489  if (dm_field->kspOperators.empty()) {
491  } else {
492 #if PETSC_VERSION_GE(3, 5, 3)
493  CHKERR MatCreateVecs(
494  dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
495  g, NULL);
496 #else
497  CHKERR MatGetVecs(
498  dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
499  g, NULL);
500 #endif
501  CHKERR VecSetDM(*g, dm);
502  }
503  PetscInfo1(dm, "Create global vector DMMGViaApproxOrders leveldown = %d\n",
504  dm->leveldown);
506 }
507 
510  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
511  "MOFEM Multi-Grid (Orders) pre-conditioner", "none");
512  CHKERRG(ierr);
513 
514  CHKERR PetscOptionsInt("-mofem_mg_levels", "nb levels of multi-grid solver",
515  "", 2, &nbLevels, PETSC_NULL);
516  CHKERR PetscOptionsInt("-mofem_mg_coarse_order",
517  "approximation order of coarse level", "", 2,
518  &coarseOrder, PETSC_NULL);
519  CHKERR PetscOptionsInt("-mofem_mg_order_at_last_level", "order at last level",
520  "", 100, &orderAtLastLevel, PETSC_NULL);
521  CHKERR PetscOptionsInt("-mofem_mg_verbose", "nb levels of multi-grid solver",
522  "", 0, &verboseLevel, PETSC_NULL);
523  PetscBool shell_sub_a = shellSubA ? PETSC_TRUE : PETSC_FALSE;
524  CHKERR PetscOptionsBool("-mofem_mg_shell_a", "use shell matrix as sub matrix",
525  "", shell_sub_a, &shell_sub_a, NULL);
526  shellSubA = (shellSubA == PETSC_TRUE);
527 
528  ierr = PetscOptionsEnd();
529  CHKERRG(ierr);
531 }
532 
534  MoFEM::Interface *m_field_ptr;
535  MoFEM::ISManager *is_manager_ptr;
537  // if is last level, take all remaining orders dofs, if any left
538  CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
539  CHKERR m_field_ptr->getInterface(is_manager_ptr);
540  const Problem *problem_ptr;
541  CHKERR DMMoFEMGetProblemPtr(dM, &problem_ptr);
542  int order_at_next_level = kk + coarseOrder;
543  if (kk == nbLevels - 1) {
544  int first = problem_ptr->getNumeredRowDofsPtr()
545  ->get<PetscGlobalIdx_mi_tag>()
546  .lower_bound(0)
547  ->get()
548  ->getPetscGlobalDofIdx();
549  CHKERR ISCreateStride(PETSC_COMM_WORLD, problem_ptr->getNbLocalDofsRow(),
550  first, 1, is);
552  // order_at_next_level = orderAtLastLevel;
553  }
554  string problem_name = problem_ptr->getName();
555  CHKERR is_manager_ptr->isCreateProblemOrder(problem_name, ROW, 0,
556  order_at_next_level, is);
558 }
559 
562  CHKERR ISDestroy(is);
564 }
565 
568  int verb) {
570  verb = verb > verboseLevel ? verb : verboseLevel;
571 
572  MPI_Comm comm;
573  CHKERR PetscObjectGetComm((PetscObject)dM, &comm);
574 
575  MOFEM_LOG_C("PCMGViaApproximationOrders", Sev::inform, "set MG levels %u\n",
576  nbLevels);
577 
578  MOFEM_LOG_CHANNEL("SYNC");
579  MOFEM_LOG_TAG("SYNC", "PCMGViaApproximationOrders")
580 
581  std::vector<IS> is_vec(nbLevels + 1);
582  std::vector<int> is_glob_size(nbLevels + 1), is_loc_size(nbLevels + 1);
583 
584  for (int kk = 0; kk < nbLevels; kk++) {
585 
586  // get indices up to up to give approximation order
587  CHKERR createIsAtLevel(kk, &is_vec[kk]);
588  CHKERR ISGetSize(is_vec[kk], &is_glob_size[kk]);
589  CHKERR ISGetLocalSize(is_vec[kk], &is_loc_size[kk]);
590 
591  MOFEM_LOG_C("SYNC", Sev::inform,
592  "Nb. dofs at level [ %d ] global %u local %d\n", kk,
593  is_glob_size[kk], is_loc_size[kk]);
594 
595  // if no dofs on level kk finish here
596  if (is_glob_size[kk] == 0) {
597  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "no dofs at level");
598  }
599  }
600 
601  for (int kk = 0; kk != nbLevels; kk++) {
602  if (kk == nbLevels - 1 && use_mat_a) {
604  false);
605  } else {
606  if (kk > 0) {
607  // Not coarse level
609  shellSubA);
610  } else {
611  // Coarse lave is compressed matrix allowing for factorization when
612  // needed
614  false);
615  }
616  }
617  }
618 
619  for (unsigned int kk = 0; kk < is_vec.size(); kk++) {
620  CHKERR destroyIsAtLevel(kk, &is_vec[kk]);
621  }
622 
623  MOFEM_LOG_SEVERITY_SYNC(comm, Sev::inform);
624  MOFEM_LOG_CHANNEL("SYNC");
625 
627 }
628 
629 boost::shared_ptr<PCMGSetUpViaApproxOrdersCtx>
630 createPCMGSetUpViaApproxOrdersCtx(DM dm, Mat A, bool use_shell_mat) {
631  return boost::make_shared<PCMGSetUpViaApproxOrdersCtx>(dm, A, use_shell_mat);
632 }
633 
635  PC pc, boost::shared_ptr<PCMGSetUpViaApproxOrdersCtx> ctx, int verb) {
636 
637  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
639 
640  auto core_log = logging::core::get();
642  "PCMGViaApproximationOrders"));
643  LogManager::setLog("PCMGViaApproximationOrders");
644  MOFEM_LOG_TAG("PCMGViaApproximationOrders", "PCMGViaApproximationOrders");
645 
646  MOFEM_LOG("PCMGViaApproximationOrders", Sev::verbose)
647  << "Setup PCMGSetUpViaApproxOrders";
648 
649  CHKERR ctx->getOptions();
650  CHKERR ctx->buildProlongationOperator(true, verb);
651 
652 #if PETSC_VERSION_GE(3, 8, 0)
653  CHKERR PCMGSetGalerkin(pc, PC_MG_GALERKIN_NONE);
654 #else
655  CHKERR PCMGSetGalerkin(pc, PETSC_FALSE);
656 #endif
657 
658  CHKERR PCMGSetLevels(pc, ctx->nbLevels, NULL);
659 
660  MOFEM_LOG("PCMGViaApproximationOrders", Sev::noisy)
661  << "Setup PCMGSetUpViaApproxOrders <- end";
662 
664 }
665 
666 } // namespace MoFEM
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
g
constexpr double g
Definition: shallow_wave.cpp:63
omega
constexpr double omega
Save field DOFS on vertices/tags.
Definition: dynamic_first_order_con_law.cpp:93
MoFEM::PCMGSubMatrixCtx::initData
MoFEMErrorCode initData(Vec x)
Definition: PCMGSetUpViaApproxOrders.cpp:57
MoFEM::PCMGSubMatrixCtx::iS
SmartPetscObj< IS > iS
Definition: PCMGSetUpViaApproxOrders.cpp:43
MoFEM::PCMGSetUpViaApproxOrdersCtx::buildProlongationOperator
virtual MoFEMErrorCode buildProlongationOperator(bool use_mat_a, int verb=0)
Set up data structures for MG.
Definition: PCMGSetUpViaApproxOrders.cpp:567
MoFEM::DMMGViaApproxOrdersCtx::kspOperators
std::vector< SmartPetscObj< Mat > > kspOperators
Get KSP operators.
Definition: PCMGSetUpViaApproxOrders.cpp:85
MOFEM_LOG_SEVERITY_SYNC
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
MoFEM::PCMGSubMatrixCtx::MOFEM_EVENT_mult
PetscLogEvent MOFEM_EVENT_mult
Definition: PCMGSetUpViaApproxOrders.cpp:48
MoFEM::PCMGSetUpViaApproxOrdersCtx::nbLevels
int nbLevels
number of multi-grid levels
Definition: PCMGSetUpViaApproxOrders.cpp:135
MoFEM::createPCMGSetUpViaApproxOrdersCtx
boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > createPCMGSetUpViaApproxOrdersCtx(DM dm, Mat A, bool use_shell_mat)
createPCMGSetUpViaApproxOrdersCtx
Definition: PCMGSetUpViaApproxOrders.cpp:630
MoFEM::PCMGSubMatrixCtx::PCMGSubMatrixCtx
PCMGSubMatrixCtx(Mat a, IS is)
Definition: PCMGSetUpViaApproxOrders.cpp:52
MoFEM::DMMGViaApproxOrdersCtx
Structure for DM for multi-grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:72
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
DMMGViaApproxOrdersCtx
Structure for DM for multi-grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.hpp:22
MoFEM::PCMGSetUpViaApproxOrdersCtx::dM
DM dM
Distributed mesh manager.
Definition: PCMGSetUpViaApproxOrders.cpp:132
MoFEM::DMCreateGlobalVector_MGViaApproxOrders
MoFEMErrorCode DMCreateGlobalVector_MGViaApproxOrders(DM dm, Vec *g)
Create global vector for DMGViaApproxOrders.
Definition: PCMGSetUpViaApproxOrders.cpp:484
MoFEM::sub_mat_sor
MoFEMErrorCode sub_mat_sor(Mat mat, Vec b, PetscReal omega, MatSORType flag, PetscReal shift, PetscInt its, PetscInt lits, Vec x)
Definition: PCMGSetUpViaApproxOrders.cpp:168
DMMoFEM.hpp
Discrete manager interface for MoFEM.
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::DMMGViaApproxOrdersCtx::aO
SmartPetscObj< AO > aO
Definition: PCMGSetUpViaApproxOrders.cpp:83
DMCtxImpl.hpp
Implementation of DM context. You should not use it directly.
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
MoFEM::sub_mat_mult
MoFEMErrorCode sub_mat_mult(Mat a, Vec x, Vec f)
Definition: PCMGSetUpViaApproxOrders.cpp:160
GET_DM_FIELD
#define GET_DM_FIELD(DM)
Definition: PCMGSetUpViaApproxOrders.cpp:214
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::PCMGSetUpViaApproxOrdersCtx
Definition: PCMGSetUpViaApproxOrders.cpp:92
ROW
@ ROW
Definition: definitions.h:136
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1167
MoFEM::PCMGSubMatrixCtx::~PCMGSubMatrixCtx
virtual ~PCMGSubMatrixCtx()=default
Definition: PCMGSetUpViaApproxOrders.cpp:41
MoFEM::DMDestroy_MGViaApproxOrders
PetscErrorCode DMDestroy_MGViaApproxOrders(DM dm)
Destroy DM.
Definition: PCMGSetUpViaApproxOrders.cpp:335
MoFEM::DMMGViaApproxOrdersSetAO
MoFEMErrorCode DMMGViaApproxOrdersSetAO(DM dm, AO ao)
Set DM ordering.
Definition: PCMGSetUpViaApproxOrders.cpp:219
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
MoFEM::DMCtxImpl
Definition: DMCtxImpl.hpp:10
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1197
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::DMMGViaApproxOrdersCtx::~DMMGViaApproxOrdersCtx
virtual ~DMMGViaApproxOrdersCtx()
Definition: PCMGSetUpViaApproxOrders.cpp:191
MoFEM::matCreateVecs
auto matCreateVecs(Mat mat)
Definition: PetscSmartObj.hpp:390
MoFEM::DMMGViaApproxOrdersCtx::DMMGViaApproxOrdersCtx
DMMGViaApproxOrdersCtx()
Definition: PCMGSetUpViaApproxOrders.cpp:189
MoFEM::SmartPetscObj::use_count
int use_count() const
Definition: PetscSmartObj.hpp:109
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::DMDestroy_MoFEM
PetscErrorCode DMDestroy_MoFEM(DM dm)
Destroys dm with MoFEM data structure.
Definition: DMMoFEM.cpp:79
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
convert.type
type
Definition: convert.py:64
MoFEM::ksp_set_operators
static MoFEMErrorCode ksp_set_operators(KSP ksp, Mat A, Mat B, void *ctx)
Definition: PCMGSetUpViaApproxOrders.cpp:308
MoFEM::DMRegister_MGViaApproxOrders
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:302
MoFEM::PCMGSetUpViaApproxOrdersCtx::A
Mat A
Matrix at fine level.
Definition: PCMGSetUpViaApproxOrders.cpp:133
MoFEM::PCMGSetUpViaApproxOrdersCtx::getOptions
virtual MoFEMErrorCode getOptions()
get options from line command
Definition: PCMGSetUpViaApproxOrders.cpp:508
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
MoFEM::PCMGSetUpViaApproxOrdersCtx::PCMGSetUpViaApproxOrdersCtx
PCMGSetUpViaApproxOrdersCtx(DM dm, Mat a, bool shell_sub_a)
Definition: PCMGSetUpViaApproxOrders.cpp:94
MoFEM::PCMGSubMatrixCtx::X
SmartPetscObj< Vec > X
Definition: PCMGSetUpViaApproxOrders.cpp:45
MoFEM::Problem::getNumeredRowDofsPtr
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Definition: ProblemsMultiIndices.hpp:82
MoFEM::sub_mat_mult_generic
MoFEMErrorCode sub_mat_mult_generic(Mat a, Vec x, Vec f)
Definition: PCMGSetUpViaApproxOrders.cpp:144
MoFEM::PCMGSubMatrixCtx::sub_mat_mult_generic
friend MoFEMErrorCode sub_mat_mult_generic(Mat a, Vec x, Vec f)
Definition: PCMGSetUpViaApproxOrders.cpp:144
MoFEM::sub_mat_mult_add
MoFEMErrorCode sub_mat_mult_add(Mat a, Vec x, Vec f)
Definition: PCMGSetUpViaApproxOrders.cpp:164
MoFEM::DMMGViaApproxOrdersClearCoarseningIS
MoFEMErrorCode DMMGViaApproxOrdersClearCoarseningIS(DM dm)
Clear approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:293
MoFEM::PetscGlobalIdx_mi_tag
Definition: TagMultiIndices.hpp:38
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::Problem::getNbLocalDofsRow
DofIdx getNbLocalDofsRow() const
Definition: ProblemsMultiIndices.hpp:378
MoFEM::isDuplicate
auto isDuplicate(IS is)
Definition: PetscSmartObj.hpp:396
MoFEM::PCMGSubMatrixCtx::A
SmartPetscObj< Mat > A
Definition: PCMGSetUpViaApproxOrders.cpp:42
MoFEM::DMMGViaApproxOrdersPopBackCoarseningIS
MoFEMErrorCode DMMGViaApproxOrdersPopBackCoarseningIS(DM dm)
Pop IS form MG via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:281
MoFEM::PCMGSetUpViaApproxOrdersCtx::~PCMGSetUpViaApproxOrdersCtx
virtual ~PCMGSetUpViaApproxOrdersCtx()=default
MoFEM::DMSetOperators_MoFEM
PetscErrorCode DMSetOperators_MoFEM(DM dm)
Set operators for MoFEM dm.
Definition: DMMoFEM.cpp:49
MoFEM::PCMGSubMatrixCtx::sCat
SmartPetscObj< VecScatter > sCat
Definition: PCMGSetUpViaApproxOrders.cpp:44
MoFEM::DMCreateInterpolation_MGViaApproxOrders
MoFEMErrorCode DMCreateInterpolation_MGViaApproxOrders(DM dm1, DM dm2, Mat *mat, Vec *vec)
Create interpolation matrix between data managers dm1 and dm2.
Definition: PCMGSetUpViaApproxOrders.cpp:391
MoFEM::PCMGSetUpViaApproxOrdersCtx::verboseLevel
int verboseLevel
Definition: PCMGSetUpViaApproxOrders.cpp:140
MoFEM::DMMGViaApproxOrdersCtx::shellMatrixCtxPtr
boost::ptr_vector< PCMGSubMatrixCtx > shellMatrixCtxPtr
Shell sub-matrix context.
Definition: PCMGSetUpViaApproxOrders.cpp:87
convert.n
n
Definition: convert.py:82
N
const int N
Definition: speed_test.cpp:3
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
MoFEM::PCMGSetUpViaApproxOrdersCtx::shellSubA
bool shellSubA
Definition: PCMGSetUpViaApproxOrders.cpp:139
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::DMoFEMGetInterfacePtr
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMoFEM.cpp:414
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
MoFEM::PCMGSetUpViaApproxOrdersCtx::orderAtLastLevel
int orderAtLastLevel
set maximal evaluated order
Definition: PCMGSetUpViaApproxOrders.cpp:137
MoFEM::createVecScatter
auto createVecScatter(Vec x, IS ix, Vec y, IS iy)
Create a Vec Scatter object.
Definition: PetscSmartObj.hpp:361
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::DMCreate_MGViaApproxOrders
MoFEMErrorCode DMCreate_MGViaApproxOrders(DM dm)
Create DM data structure for Multi-Grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:313
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:426
MoFEM::DMMGViaApproxOrdersCtx::coarsenDM
SmartPetscObj< DM > coarsenDM
Definition: PCMGSetUpViaApproxOrders.cpp:82
MoFEM::DMMGViaApproxOrdersPushBackCoarseningIS
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.
Definition: PCMGSetUpViaApproxOrders.cpp:227
MoFEM::DMMGViaApproxOrdersCtx::destroyCoarseningIS
MoFEMErrorCode destroyCoarseningIS()
Definition: PCMGSetUpViaApproxOrders.cpp:196
MoFEM::PCMGSetUpViaApproxOrdersCtx::createIsAtLevel
virtual MoFEMErrorCode createIsAtLevel(int kk, IS *is)
Set IS for levels.
Definition: PCMGSetUpViaApproxOrders.cpp:533
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::DMMGViaApproxOrdersCtx::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, MoFEM::UnknownInterface **iface) const
Definition: PCMGSetUpViaApproxOrders.cpp:207
MoFEM::PCMGSetUpViaApproxOrdersCtx::destroyIsAtLevel
virtual MoFEMErrorCode destroyIsAtLevel(int kk, IS *is)
Destroy IS if internally created.
Definition: PCMGSetUpViaApproxOrders.cpp:560
MoFEM::PCMGSubMatrixCtx::F
SmartPetscObj< Vec > F
Definition: PCMGSetUpViaApproxOrders.cpp:46
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEM::PCMGSubMatrixCtx::MOFEM_EVENT_sor
PetscLogEvent MOFEM_EVENT_sor
Definition: PCMGSetUpViaApproxOrders.cpp:49
MoFEM::PCMGSetUpViaApproxOrders
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > ctx, int verb)
Function build MG structure.
Definition: PCMGSetUpViaApproxOrders.cpp:634
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< Mat >
MoFEM::ISManager::isCreateProblemOrder
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)
Definition: ISManager.cpp:204
MoFEM::Problem::getName
auto getName() const
Definition: ProblemsMultiIndices.hpp:372
MoFEM::PCMGSetUpViaApproxOrdersCtx::coarseOrder
int coarseOrder
approximation order of coarse level
Definition: PCMGSetUpViaApproxOrders.cpp:136
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEM::DMCreateMatrix_MGViaApproxOrders
MoFEMErrorCode DMCreateMatrix_MGViaApproxOrders(DM dm, Mat *M)
Create matrix for Multi-Grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:344
MoFEM::DMMGViaApproxOrdersCtx::coarseningIS
std::vector< SmartPetscObj< IS > > coarseningIS
Coarsening IS.
Definition: PCMGSetUpViaApproxOrders.cpp:84
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
MoFEM::PCMGSubMatrixCtx
Definition: PCMGSetUpViaApproxOrders.cpp:29
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::PCMGSubMatrixCtx::sub_mat_sor
friend MoFEMErrorCode sub_mat_sor(Mat mat, Vec b, PetscReal omega, MatSORType flag, PetscReal shift, PetscInt its, PetscInt lits, Vec x)
Definition: PCMGSetUpViaApproxOrders.cpp:168
MoFEM::DMCoarsen_MGViaApproxOrders
MoFEMErrorCode DMCoarsen_MGViaApproxOrders(DM dm, MPI_Comm comm, DM *dmc)
Coarsen DM.
Definition: PCMGSetUpViaApproxOrders.cpp:376