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 
7 #include <MoFEM.hpp>
8 
9 #include <UnknownInterface.hpp>
10 
11 using namespace MoFEM;
12 
14 
15 #undef PETSC_VERSION_RELEASE
16 #define PETSC_VERSION_RELEASE 1
17 
18 #if PETSC_VERSION_GE(3, 6, 0)
19 #include <petsc/private/petscimpl.h>
20 #else
21 #include <petsc-private/petscimpl.h>
22 #endif
23 
24 #if PETSC_VERSION_GE(3, 6, 0)
25 #include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/
26 // #include <petsc/private/vecimpl.h> /*I "petscdm.h" I*/
27 #else
28 #include <petsc-private/dmimpl.h> /*I "petscdm.h" I*/
29 #include <petsc-private/vecimpl.h> /*I "petscdm.h" I*/
30 #endif
31 
32 PCMGSubMatrixCtx::PCMGSubMatrixCtx(Mat a, IS is) : A(a), iS(is) {
33  // Increase reference of petsc object (works like shared_ptr but unique for
34  // PETSc)
35  ierr = PetscObjectReference((PetscObject)A);
36  CHKERRABORT(PETSC_COMM_WORLD, ierr);
37  ierr = PetscObjectReference((PetscObject)iS);
38  CHKERRABORT(PETSC_COMM_WORLD, ierr);
39 }
40 
42  ierr = MatDestroy(&A);
43  CHKERRABORT(PETSC_COMM_WORLD, ierr);
44  ierr = ISDestroy(&iS);
45  CHKERRABORT(PETSC_COMM_WORLD, ierr);
46 }
47 
50  : PCMGSubMatrixCtx(a, is), isInitisalised(false) {
51  PetscLogEventRegister("PCMGSubMatrixCtx_mult", 0, &MOFEM_EVENT_mult);
52  PetscLogEventRegister("PCMGSubMatrixCtx_sor", 0, &MOFEM_EVENT_sor);
53  }
55  if (isInitisalised) {
56  ierr = VecScatterDestroy(&sCat);
57  CHKERRABORT(PETSC_COMM_WORLD, ierr);
58  ierr = VecDestroy(&X);
59  CHKERRABORT(PETSC_COMM_WORLD, ierr);
60  ierr = VecDestroy(&F);
61  CHKERRABORT(PETSC_COMM_WORLD, ierr);
62  }
63  }
64  template <InsertMode MODE>
65  friend MoFEMErrorCode sub_mat_mult_generic(Mat a, Vec x, Vec f);
66  friend MoFEMErrorCode sub_mat_sor(Mat mat, Vec b, PetscReal omega,
67  MatSORType flag, PetscReal shift,
68  PetscInt its, PetscInt lits, Vec x);
69 
70 public:
73  if (!isInitisalised) {
74  CHKERR MatCreateVecs(A, &X, &F);
75  CHKERR VecScatterCreate(X, iS, x, PETSC_NULL, &sCat);
76  CHKERR VecZeroEntries(X);
77  CHKERR VecZeroEntries(F);
78  isInitisalised = true;
79  }
81  }
82  PetscLogEvent MOFEM_EVENT_mult;
83  PetscLogEvent MOFEM_EVENT_sor;
85 };
86 
87 template <InsertMode MODE>
89  void *void_ctx;
91  CHKERR MatShellGetContext(a, &void_ctx);
93  if (!ctx->isInitisalised) {
94  CHKERR ctx->initData(x);
95  }
96  PetscLogEventBegin(ctx->MOFEM_EVENT_mult, 0, 0, 0, 0);
97  CHKERR VecScatterBegin(ctx->sCat, x, ctx->X, INSERT_VALUES, SCATTER_REVERSE);
98  CHKERR VecScatterEnd(ctx->sCat, x, ctx->X, INSERT_VALUES, SCATTER_REVERSE);
99  CHKERR MatMult(ctx->A, ctx->X, ctx->F);
100  CHKERR VecScatterBegin(ctx->sCat, ctx->F, f, MODE, SCATTER_FORWARD);
101  CHKERR VecScatterEnd(ctx->sCat, ctx->F, f, MODE, SCATTER_FORWARD);
102  PetscLogEventEnd(ctx->MOFEM_EVENT_mult, 0, 0, 0, 0);
104 }
105 
107  return sub_mat_mult_generic<INSERT_VALUES>(a, x, f);
108 }
109 
111  return sub_mat_mult_generic<ADD_VALUES>(a, x, f);
112 }
113 
114 MoFEMErrorCode sub_mat_sor(Mat mat, Vec b, PetscReal omega, MatSORType flag,
115  PetscReal shift, PetscInt its, PetscInt lits,
116  Vec x) {
117  void *void_ctx;
119  CHKERR MatShellGetContext(mat, &void_ctx);
121  if (!ctx->isInitisalised) {
122  CHKERR ctx->initData(x);
123  }
124  PetscLogEventBegin(ctx->MOFEM_EVENT_sor, 0, 0, 0, 0);
125  CHKERR VecScatterBegin(ctx->sCat, b, ctx->X, INSERT_VALUES, SCATTER_REVERSE);
126  CHKERR VecScatterEnd(ctx->sCat, b, ctx->X, INSERT_VALUES, SCATTER_REVERSE);
127  CHKERR MatSOR(ctx->A, ctx->X, omega, flag, shift, its, lits, ctx->F);
128  CHKERR VecScatterBegin(ctx->sCat, ctx->F, x, INSERT_VALUES, SCATTER_FORWARD);
129  CHKERR VecScatterEnd(ctx->sCat, ctx->F, x, INSERT_VALUES, SCATTER_FORWARD);
130  PetscLogEventEnd(ctx->MOFEM_EVENT_sor, 0, 0, 0, 0);
132 }
133 
135  : MoFEM::DMCtx(), aO(PETSC_NULL) {
136  // std::cerr << "create dm\n";
137 }
140  CHKERRABORT(PETSC_COMM_WORLD, ierr);
141 }
142 
145  for (unsigned int ii = 0; ii < coarseningIS.size(); ii++) {
146  if (coarseningIS[ii])
147  CHKERR ISDestroy(&coarseningIS[ii]);
148  }
149  for (unsigned int ii = 0; ii < kspOperators.size(); ii++) {
150  if (kspOperators[ii])
151  CHKERR MatDestroy(&kspOperators[ii]);
152  }
153  if (aO) {
154  CHKERR AODestroy(&aO);
155  }
156  coarseningIS.clear();
157  kspOperators.clear();
158  shellMatrixCtxPtr.clear();
160 }
161 
163 DMMGViaApproxOrdersCtx::query_interface(boost::typeindex::type_index type_index,
164  MoFEM::UnknownInterface **iface) const {
165  *iface = static_cast<DMMGViaApproxOrdersCtx *>(
166  const_cast<DMMGViaApproxOrdersCtx *>(this));
167  return 0;
168 }
169 
170 #define GET_DM_FIELD(DM) \
171  auto dm_field = \
172  static_cast<DMCtx *>(DM->data)->getInterface<DMMGViaApproxOrdersCtx>(); \
173  NOT_USED(dm_field)
174 
176  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
178  GET_DM_FIELD(dm);
179  *ctx = dm_field;
181 }
182 
184  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
186  GET_DM_FIELD(dm);
187  if (dm_field->aO) {
188  // std::cerr << dm_field->aO << std::endl;
189  CHKERR AODestroy(&dm_field->aO);
190  // std::cerr << "destroy ao when adding\n";
191  }
192  dm_field->aO = ao;
193  CHKERR PetscObjectReference((PetscObject)ao);
194  // std::cerr << "add ao\n";
196 }
197 
199  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
201  GET_DM_FIELD(dm);
202  *size = dm_field->coarseningIS.size();
204 }
205 
207  Mat *subA,
208  bool create_sub_matrix,
209  bool shell_sub_a) {
210  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
212  GET_DM_FIELD(dm);
213  dm_field->coarseningIS.push_back(is);
214  dm_field->shellMatrixCtxPtr.push_back(new PCMGSubMatrixCtx_private(A, is));
215  if (is) {
216  CHKERR PetscObjectReference((PetscObject)is);
217  }
218  if (is) {
219  IS is2 = is;
220  if (dm_field->aO) {
221  CHKERR ISDuplicate(is, &is2);
222  CHKERR ISCopy(is, is2);
223  CHKERR AOApplicationToPetscIS(dm_field->aO, is2);
224  }
225  if (create_sub_matrix) {
226  if (shell_sub_a) {
227  int n, N;
228  CHKERR ISGetSize(is, &N);
229  CHKERR ISGetLocalSize(is, &n);
230  MPI_Comm comm;
231  CHKERR PetscObjectGetComm((PetscObject)A, &comm);
232  CHKERR MatCreateShell(comm, n, n, N, N,
233  &(dm_field->shellMatrixCtxPtr.back()), subA);
234  CHKERR MatShellSetOperation(*subA, MATOP_MULT,
235  (void (*)(void))sub_mat_mult);
236  CHKERR MatShellSetOperation(*subA, MATOP_MULT_ADD,
237  (void (*)(void))sub_mat_mult_add);
238  CHKERR MatShellSetOperation(*subA, MATOP_SOR,
239  (void (*)(void))sub_mat_sor);
240  } else {
241 #if PETSC_VERSION_GE(3, 8, 0)
242  CHKERR MatCreateSubMatrix(A, is2, is2, MAT_INITIAL_MATRIX, subA);
243 #else
244  CHKERR MatGetSubMatrix(A, is2, is2, MAT_INITIAL_MATRIX, subA);
245 #endif
246  }
247  }
248  if (dm_field->aO) {
249  CHKERR ISDestroy(&is2);
250  }
251  dm_field->kspOperators.push_back(*subA);
252  CHKERR PetscObjectReference((PetscObject)(*subA));
253  } else {
254  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
255  }
256  PetscInfo(dm, "Push back IS to DMMGViaApproxOrders\n");
258 }
259 
261  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
263  GET_DM_FIELD(dm);
264  if (dm_field->coarseningIS.back()) {
265  CHKERR ISDestroy(&dm_field->coarseningIS.back());
266  dm_field->coarseningIS.pop_back();
267  }
268  if (dm_field->kspOperators.back()) {
269  CHKERR MatDestroy(&dm_field->kspOperators.back());
270  }
271  dm_field->kspOperators.pop_back();
272  PetscInfo(dm, "Pop back IS to DMMGViaApproxOrders\n");
274 }
275 
277  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
279  GET_DM_FIELD(dm);
280  CHKERR dm_field->destroyCoarseningIS();
281  PetscInfo(dm, "Clear DMs data structures\n");
283 }
284 
286  int nb_elems, Mat A,
287  int verb) {
288  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
290  GET_DM_FIELD(dm);
291  int nb_no_changed = 0;
292  int nb_replaced = 0;
293  int nb_deleted = 0;
294  int nb_added = 0;
295  std::vector<IS>::iterator it;
296  it = dm_field->coarseningIS.begin();
297  int ii = 0;
298  for (; it != dm_field->coarseningIS.end(); it++, ii++) {
299  if (ii < nb_elems) {
300  PetscBool flg;
301  CHKERR ISEqual(*it, is_vec[ii], &flg);
302  if (!flg) {
303  CHKERR ISDestroy(&*it);
304  CHKERR MatDestroy(&dm_field->kspOperators[ii]);
305  *it = is_vec[ii];
306  CHKERR PetscObjectReference((PetscObject)is_vec[ii]);
307  if (ii < nb_elems - 1) {
308  IS is = is_vec[ii];
309  if (dm_field->aO) {
310  CHKERR ISDuplicate(is_vec[ii], &is);
311  CHKERR ISCopy(is_vec[ii], is);
312  CHKERR AOApplicationToPetscIS(dm_field->aO, is);
313  }
314  Mat subA;
315 #if PETSC_VERSION_GE(3, 8, 0)
316  CHKERR MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &subA);
317 #else
318  CHKERR MatGetSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &subA);
319 #endif
320  CHKERR PetscObjectReference((PetscObject)subA);
321  dm_field->kspOperators[ii] = subA;
322  CHKERR MatDestroy(&subA);
323  if (dm_field->aO) {
324  CHKERR ISDestroy(&is);
325  }
326  } else {
327  CHKERR PetscObjectReference((PetscObject)A);
328  dm_field->kspOperators[ii] = A;
329  }
330  nb_replaced++;
331  }
332  } else {
333  nb_no_changed++;
334  continue;
335  }
336  }
337  if (static_cast<int>(dm_field->coarseningIS.size()) < nb_elems) {
338  for (; ii < nb_elems - 1; ii++) {
339  Mat subA;
340  CHKERR DMMGViaApproxOrdersPushBackCoarseningIS(dm, is_vec[ii], A, &subA,
341  true, false);
342  CHKERR MatDestroy(&subA);
343  nb_added++;
344  }
345  CHKERR DMMGViaApproxOrdersPushBackCoarseningIS(dm, is_vec[ii], A, &A, false,
346  false);
347  nb_added++;
348  } else {
349  for (; ii < static_cast<int>(dm_field->coarseningIS.size()); ii++) {
351  nb_deleted++;
352  }
353  }
354  MPI_Comm comm;
355  CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
356  if (verb > 0) {
357  PetscPrintf(comm,
358  "DMMGViaApproxOrders nb_no_changed = %d, nb_replaced = %d, "
359  "nb_added = %d, nb_deleted = %d, size = %d\n",
360  nb_no_changed, nb_replaced, nb_added, nb_deleted,
361  dm_field->coarseningIS.size());
362  }
363  PetscInfo(dm, "Replace IS to DMMGViaApproxOrders\n");
365 }
366 
368  const DMMGViaApproxOrdersCtx **ctx) {
369  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
371  GET_DM_FIELD(dm);
372  *ctx = dm_field;
374 }
375 
378  CHKERR DMRegister(sname, DMCreate_MGViaApproxOrders);
380 }
381 
382 // MoFEMErrorCode DMDestroy_MGViaApproxOrders(DM dm) {
383 // PetscValidHeaderSpecific(dm,DM_CLASSID,1);
384 // MoFEMFunctionBeginHot;
385 // if(!((DMMGViaApproxOrdersCtx*)dm->data)->referenceNumber) {
386 // DMMGViaApproxOrdersCtx *dm_field = (DMMGViaApproxOrdersCtx*)dm->data;
387 // if(dm_field->destroyProblem) {
388 // if(dm_field->mField_ptr->check_problem(dm_field->problemName)) {
389 // dm_field->mField_ptr->delete_problem(dm_field->problemName);
390 // } // else problem has to be deleted by the user
391 // }
392 // cerr << "Destroy " << dm_field->problemName << endl;
393 // delete (DMMGViaApproxOrdersCtx*)dm->data;
394 // } else {
395 // DMMGViaApproxOrdersCtx *dm_field = (DMMGViaApproxOrdersCtx*)dm->data;
396 //
397 // cerr << "Derefrence " << dm_field->problemName << " " <<
398 // ((DMCtx*)dm->data)->referenceNumber << endl;
399 // (((DMMGViaApproxOrdersCtx*)dm->data)->referenceNumber)--;
400 // }
401 // MoFEMFunctionReturnHot(0);
402 // }
403 
404 static MoFEMErrorCode ksp_set_operators(KSP ksp, Mat A, Mat B, void *ctx) {
407 }
408 
410  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
412  if (!dm->data) {
413  dm->data = new DMMGViaApproxOrdersCtx();
414  } else {
415  ((DMCtx *)(dm->data))->referenceNumber++;
416  }
417  // cerr << "Create " << ((DMCtx*)(dm->data))->referenceNumber << endl;
419  dm->ops->creatematrix = DMCreateMatrix_MGViaApproxOrders;
420  dm->ops->createglobalvector = DMCreateGlobalVector_MGViaApproxOrders;
421  dm->ops->coarsen = DMCoarsen_MGViaApproxOrders;
422  // dm->ops->destroy = DMDestroy_MGViaApproxOrders;
423  dm->ops->createinterpolation = DMCreateInterpolation_MGViaApproxOrders;
424  CHKERR DMKSPSetComputeOperators(dm, ksp_set_operators, NULL);
425  PetscInfo1(dm, "Create DMMGViaApproxOrders reference = %d\n",
426  ((DMCtx *)(dm->data))->referenceNumber);
428 }
429 
431 
432  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
434  GET_DM_FIELD(dm);
435 
436  int leveldown = dm->leveldown;
437 
438  if (dm_field->kspOperators.empty()) {
440  } else {
441  MPI_Comm comm;
442  CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
443  if (dm_field->kspOperators.empty()) {
444  SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
445  "data inconsistency, operator can not be set");
446  }
447  if (static_cast<int>(dm_field->kspOperators.size()) < leveldown) {
448  SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
449  "data inconsistency, no IS for that level");
450  }
451  *M = dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown];
452  CHKERR PetscObjectReference((PetscObject)*M);
453  CHKERR MatSetDM(*M, dm);
454  }
455 
456  PetscInfo1(dm, "Create Matrix DMMGViaApproxOrders leveldown = %d\n",
457  leveldown);
458 
460 }
461 
462 MoFEMErrorCode DMCoarsen_MGViaApproxOrders(DM dm, MPI_Comm comm, DM *dmc) {
463  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
465  GET_DM_FIELD(dm);
466  CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
467  CHKERR DMCreate(comm, dmc);
468  (*dmc)->data = dm->data;
469  DMType type;
470  CHKERR DMGetType(dm, &type);
471  CHKERR DMSetType(*dmc, type);
472  CHKERR PetscObjectReference((PetscObject)(*dmc));
473  PetscInfo1(dm, "Coarsen DMMGViaApproxOrders leveldown = %d\n", dm->leveldown);
475 }
476 
478  Vec *vec) {
479  PetscValidHeaderSpecific(dm1, DM_CLASSID, 1);
480  PetscValidHeaderSpecific(dm2, DM_CLASSID, 1);
482 
483  MPI_Comm comm;
484  CHKERR PetscObjectGetComm((PetscObject)dm1, &comm);
485 
486  int m, n, M, N;
487 
488  DM dm_down = dm1;
489  DM dm_up = dm2;
490 
491  int dm_down_leveldown = dm_down->leveldown;
492  int dm_up_leveldown = dm_up->leveldown;
493 
494  PetscInfo2(dm1,
495  "Create interpolation DMMGViaApproxOrders dm1_leveldown = %d "
496  "dm2_leveldown = %d\n",
497  dm_down_leveldown, dm_up_leveldown);
498 
499  IS is_down, is_up;
500  {
501  // Coarser mesh
502  GET_DM_FIELD(dm_down);
503  if (static_cast<int>(dm_field->coarseningIS.size()) < dm_down_leveldown) {
504  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
505  }
506  is_down = dm_field->coarseningIS[dm_field->coarseningIS.size() - 1 -
507  dm_down_leveldown];
508  CHKERR ISGetSize(is_down, &M);
509  CHKERR ISGetLocalSize(is_down, &m);
510  }
511  {
512  // Finer mesh
513  GET_DM_FIELD(dm_up);
514  if (static_cast<int>(dm_field->coarseningIS.size()) < dm_up_leveldown) {
515  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
516  }
517  is_up =
518  dm_field
519  ->coarseningIS[dm_field->coarseningIS.size() - 1 - dm_up_leveldown];
520  CHKERR ISGetSize(is_up, &N);
521  CHKERR ISGetLocalSize(is_up, &n);
522  }
523 
524  // is_dow rows
525  // is_up columns
526 
527  CHKERR MatCreate(comm, mat);
528  CHKERR MatSetSizes(*mat, m, n, M, N);
529  CHKERR MatSetType(*mat, MATMPIAIJ);
530  CHKERR MatMPIAIJSetPreallocation(*mat, 1, PETSC_NULL, 0, PETSC_NULL);
531 
532  // get matrix layout
533  PetscLayout rmap, cmap;
534  CHKERR MatGetLayouts(*mat, &rmap, &cmap);
535  int rstart, rend, cstart, cend;
536  CHKERR PetscLayoutGetRange(rmap, &rstart, &rend);
537  CHKERR PetscLayoutGetRange(cmap, &cstart, &cend);
538 
539  // if(verb>0) {
540  // PetscSynchronizedPrintf(comm,"level %d row start %d row end
541  // %d\n",kk,rstart,rend); PetscSynchronizedPrintf(comm,"level %d col start
542  // %d col end %d\n",kk,cstart,cend);
543  // }
544 
545  const int *row_indices_ptr, *col_indices_ptr;
546  CHKERR ISGetIndices(is_down, &row_indices_ptr);
547  CHKERR ISGetIndices(is_up, &col_indices_ptr);
548 
549  map<int, int> idx_map;
550  for (int ii = 0; ii < m; ii++) {
551  idx_map[row_indices_ptr[ii]] = rstart + ii;
552  }
553 
554  CHKERR MatZeroEntries(*mat);
555  // FIXME: Use MatCreateMPIAIJWithArrays and set array directly
556  for (int jj = 0; jj < n; jj++) {
557  map<int, int>::iterator mit = idx_map.find(col_indices_ptr[jj]);
558  if (mit != idx_map.end()) {
559  CHKERR MatSetValue(*mat, mit->second, cstart + jj, 1, INSERT_VALUES);
560  }
561  }
562 
563  CHKERR ISRestoreIndices(is_down, &row_indices_ptr);
564  CHKERR ISRestoreIndices(is_up, &col_indices_ptr);
565 
566  CHKERR MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
567  CHKERR MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);
568 
569  if (vec != NULL) {
570  *vec = PETSC_NULL;
571  }
572 
574 }
575 
577 
578  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
580  int leveldown = dm->leveldown;
581  GET_DM_FIELD(dm);
582  if (dm_field->kspOperators.empty()) {
584  } else {
585 #if PETSC_VERSION_GE(3, 5, 3)
586  CHKERR MatCreateVecs(
587  dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
588  g, NULL);
589 #else
590  CHKERR MatGetVecs(
591  dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
592  g, NULL);
593 #endif
594  CHKERR VecSetDM(*g, dm);
595  }
596  PetscInfo1(dm, "Create global vector DMMGViaApproxOrders leveldown = %d\n",
597  dm->leveldown);
599 }
600 
603  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
604  "MOFEM Multi-Grid (Orders) pre-conditioner", "none");
605  CHKERRG(ierr);
606 
607  CHKERR PetscOptionsInt("-mofem_mg_levels", "nb levels of multi-grid solver",
608  "", 2, &nbLevels, PETSC_NULL);
609  CHKERR PetscOptionsInt("-mofem_mg_coarse_order",
610  "approximation order of coarse level", "", 2,
611  &coarseOrder, PETSC_NULL);
612  CHKERR PetscOptionsInt("-mofem_mg_order_at_last_level", "order at last level",
613  "", 100, &orderAtLastLevel, PETSC_NULL);
614  CHKERR PetscOptionsInt("-mofem_mg_verbose", "nb levels of multi-grid solver",
615  "", 0, &verboseLevel, PETSC_NULL);
616  PetscBool shell_sub_a = shellSubA ? PETSC_TRUE : PETSC_FALSE;
617  CHKERR PetscOptionsBool("-mofem_mg_shell_a", "use shell matrix as sub matrix",
618  "", shell_sub_a, &shell_sub_a, NULL);
619  shellSubA = (shellSubA == PETSC_TRUE);
620 
621  ierr = PetscOptionsEnd();
622  CHKERRG(ierr);
624 }
625 
627  MoFEM::Interface *m_field_ptr;
628  MoFEM::ISManager *is_manager_ptr;
630  // if is last level, take all remaining orders dofs, if any left
631  CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
632  CHKERR m_field_ptr->getInterface(is_manager_ptr);
633  const Problem *problem_ptr;
634  CHKERR DMMoFEMGetProblemPtr(dM, &problem_ptr);
635  int order_at_next_level = kk + coarseOrder;
636  if (kk == nbLevels - 1) {
637  int first = problem_ptr->getNumeredRowDofsPtr()
638  ->get<PetscLocalIdx_mi_tag>()
639  .find(0)
640  ->get()
641  ->getPetscGlobalDofIdx();
642  CHKERR ISCreateStride(PETSC_COMM_WORLD, problem_ptr->getNbLocalDofsRow(),
643  first, 1, is);
645  // order_at_next_level = orderAtLastLevel;
646  }
647  string problem_name = problem_ptr->getName();
648  CHKERR is_manager_ptr->isCreateProblemOrder(problem_name, ROW, 0,
649  order_at_next_level, is);
651 }
652 
655  CHKERR ISDestroy(is);
657 }
658 
661  int verb) {
663  verb = verb > verboseLevel ? verb : verboseLevel;
664 
665  MPI_Comm comm;
666  CHKERR PetscObjectGetComm((PetscObject)dM, &comm);
667 
668  if (verb > QUIET) {
669  PetscPrintf(comm, "set MG levels %u\n", nbLevels);
670  }
671 
672  std::vector<IS> is_vec(nbLevels + 1);
673  std::vector<int> is_glob_size(nbLevels + 1), is_loc_size(nbLevels + 1);
674 
675  for (int kk = 0; kk < nbLevels; kk++) {
676 
677  // get indices up to up to give approximation order
678  CHKERR createIsAtLevel(kk, &is_vec[kk]);
679  CHKERR ISGetSize(is_vec[kk], &is_glob_size[kk]);
680  CHKERR ISGetLocalSize(is_vec[kk], &is_loc_size[kk]);
681 
682  if (verb > QUIET) {
683  PetscSynchronizedPrintf(comm,
684  "Nb. dofs at level [ %d ] global %u local %d\n",
685  kk, is_glob_size[kk], is_loc_size[kk]);
686  }
687 
688  // if no dofs on level kk finish here
689  if (is_glob_size[kk] == 0) {
690  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "no dofs at level");
691  }
692  }
693 
694  for (int kk = 0; kk != nbLevels; kk++) {
695  Mat subA;
696  if (kk == nbLevels - 1 && use_mat_a) {
697  subA = A;
699  false, false);
700  } else {
701  if (kk > 0) {
702  // Not coarse level
704  true, shellSubA);
705  } else {
706  // Coarse lave is compressed matrix allowing for factorization when
707  // needed
709  true, false);
710  }
711  if (subA) {
712  CHKERR MatDestroy(&subA);
713  }
714  }
715  }
716 
717  for (unsigned int kk = 0; kk < is_vec.size(); kk++) {
718  CHKERR destroyIsAtLevel(kk, &is_vec[kk]);
719  }
720 
721  if (verb > QUIET) {
722  PetscSynchronizedFlush(comm, PETSC_STDOUT);
723  }
724 
726 }
727 
729  int verb) {
730 
731  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
733 
734  MPI_Comm comm;
735  CHKERR PetscObjectGetComm((PetscObject)pc, &comm);
736  if (verb > 0) {
737  PetscPrintf(comm, "Start PCMGSetUpViaApproxOrders\n");
738  }
739 
740  CHKERR ctx->getOptions();
741  CHKERR ctx->buildProlongationOperator(true, verb);
742 
743 #if PETSC_VERSION_GE(3, 8, 0)
744  CHKERR PCMGSetGalerkin(pc, PC_MG_GALERKIN_NONE);
745 #else
746  CHKERR PCMGSetGalerkin(pc, PETSC_FALSE);
747 #endif
748 
749  CHKERR PCMGSetLevels(pc, ctx->nbLevels, NULL);
750 
751  if (verb > 0) {
752  PetscPrintf(comm, "End PCMGSetUpViaApproxOrders\n");
753  }
754 
756 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
PCMGSubMatrixCtx::~PCMGSubMatrixCtx
virtual ~PCMGSubMatrixCtx()
Definition: PCMGSetUpViaApproxOrders.cpp:41
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
PCMGSetUpViaApproxOrdersCtx::getOptions
virtual MoFEMErrorCode getOptions()
get options from line command
Definition: PCMGSetUpViaApproxOrders.cpp:601
PCMGSubMatrixCtx::X
Vec X
Definition: PCMGSetUpViaApproxOrders.hpp:11
DMMGViaApproxOrdersCtx::~DMMGViaApproxOrdersCtx
virtual ~DMMGViaApproxOrdersCtx()
Definition: PCMGSetUpViaApproxOrders.cpp:138
DMMGViaApproxOrdersClearCoarseningIS
MoFEMErrorCode DMMGViaApproxOrdersClearCoarseningIS(DM dm)
Clear approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:276
DMMGViaApproxOrdersCtx::shellMatrixCtxPtr
boost::ptr_vector< PCMGSubMatrixCtx > shellMatrixCtxPtr
Shell sub-matrix context.
Definition: PCMGSetUpViaApproxOrders.hpp:36
PCMGSetUpViaApproxOrdersCtx::shellSubA
bool shellSubA
Definition: PCMGSetUpViaApproxOrders.hpp:208
PCMGSubMatrixCtx_private::~PCMGSubMatrixCtx_private
~PCMGSubMatrixCtx_private()
Definition: PCMGSetUpViaApproxOrders.cpp:54
DMCoarsen_MGViaApproxOrders
MoFEMErrorCode DMCoarsen_MGViaApproxOrders(DM dm, MPI_Comm comm, DM *dmc)
Coarsen DM.
Definition: PCMGSetUpViaApproxOrders.cpp:462
DMMGViaApproxOrdersCtx::kspOperators
std::vector< Mat > kspOperators
Get KSP operators.
Definition: PCMGSetUpViaApproxOrders.hpp:34
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
PCMGSetUpViaApproxOrdersCtx
Set data structures of MG pre-conditioner via approximation orders.
Definition: PCMGSetUpViaApproxOrders.hpp:190
DMMGViaApproxOrdersCtx
Structure for DM for multi-grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.hpp:22
MoFEM::PetscLocalIdx_mi_tag
Definition: TagMultiIndices.hpp:45
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
DMMGViaApproxOrdersPushBackCoarseningIS
MoFEMErrorCode DMMGViaApproxOrdersPushBackCoarseningIS(DM dm, IS is, Mat A, Mat *subA, bool create_sub_matrix, bool shell_sub_a)
Push back coarsening level to MG via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:206
PCMGSetUpViaApproxOrders.hpp
header of multi-grid solver for p- adaptivity
PCMGSetUpViaApproxOrdersCtx::nbLevels
int nbLevels
number of multi-grid levels
Definition: PCMGSetUpViaApproxOrders.hpp:204
GET_DM_FIELD
#define GET_DM_FIELD(DM)
Definition: PCMGSetUpViaApproxOrders.cpp:170
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:123
sub_mat_mult_generic
MoFEMErrorCode sub_mat_mult_generic(Mat a, Vec x, Vec f)
Definition: PCMGSetUpViaApproxOrders.cpp:88
DMCreate_MGViaApproxOrders
MoFEMErrorCode DMCreate_MGViaApproxOrders(DM dm)
Create DM data structure for Multi-Grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:409
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1171
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1201
DMMGViaApproxOrdersCtx::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, MoFEM::UnknownInterface **iface) const
Definition: PCMGSetUpViaApproxOrders.cpp:163
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
PCMGSetUpViaApproxOrders
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, PCMGSetUpViaApproxOrdersCtx *ctx, int verb)
Function build MG structure.
Definition: PCMGSetUpViaApproxOrders.cpp:728
DMMGViaApproxOrdersCtx::DMMGViaApproxOrdersCtx
DMMGViaApproxOrdersCtx()
Definition: PCMGSetUpViaApproxOrders.cpp:134
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
sub_mat_mult
MoFEMErrorCode sub_mat_mult(Mat a, Vec x, Vec f)
Definition: PCMGSetUpViaApproxOrders.cpp:106
DMMGViaApproxOrdersSetAO
MoFEMErrorCode DMMGViaApproxOrdersSetAO(DM dm, AO ao)
Set DM ordering.
Definition: PCMGSetUpViaApproxOrders.cpp:183
convert.type
type
Definition: convert.py:64
PCMGSubMatrixCtx::F
Vec F
Definition: PCMGSetUpViaApproxOrders.hpp:11
PCMGSetUpViaApproxOrdersCtx::dM
DM dM
Distributed mesh manager.
Definition: PCMGSetUpViaApproxOrders.hpp:195
PCMGSubMatrixCtx_private
Definition: PCMGSetUpViaApproxOrders.cpp:48
DMMGViaApproxOrdersCtx::destroyCoarseningIS
MoFEMErrorCode destroyCoarseningIS()
Definition: PCMGSetUpViaApproxOrders.cpp:143
MoFEM::Problem::getNumeredRowDofsPtr
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Definition: ProblemsMultiIndices.hpp:82
UnknownInterface.hpp
MoFEM interface.
PCMGSubMatrixCtx_private::MOFEM_EVENT_mult
PetscLogEvent MOFEM_EVENT_mult
Definition: PCMGSetUpViaApproxOrders.cpp:82
PCMGSubMatrixCtx
Definition: PCMGSetUpViaApproxOrders.hpp:9
DMMGViaApproxOrdersGetCoarseningISSize
MoFEMErrorCode DMMGViaApproxOrdersGetCoarseningISSize(DM dm, int *size)
Gets size of coarseningIS in internal data struture DMMGViaApproxOrdersCtx.
Definition: PCMGSetUpViaApproxOrders.cpp:198
PCMGSetUpViaApproxOrdersCtx::createIsAtLevel
virtual MoFEMErrorCode createIsAtLevel(int kk, IS *is)
Set IS for levels.
Definition: PCMGSetUpViaApproxOrders.cpp:626
PCMGSubMatrixCtx_private::PCMGSubMatrixCtx_private
PCMGSubMatrixCtx_private(Mat a, IS is)
Definition: PCMGSetUpViaApproxOrders.cpp:49
MoFEM::DMCtx
PETSc Discrete Manager data structure.
Definition: DMMoFEM.hpp:944
PCMGSubMatrixCtx_private::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:114
PCMGSubMatrixCtx::iS
IS iS
Definition: PCMGSetUpViaApproxOrders.hpp:12
MoFEM::Problem::getNbLocalDofsRow
DofIdx getNbLocalDofsRow() const
Definition: ProblemsMultiIndices.hpp:378
DMCreateMatrix_MGViaApproxOrders
MoFEMErrorCode DMCreateMatrix_MGViaApproxOrders(DM dm, Mat *M)
Create matrix for Multi-Grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:430
DMMGViaApproxOrdersPopBackCoarseningIS
MoFEMErrorCode DMMGViaApproxOrdersPopBackCoarseningIS(DM dm)
Pop is form MG via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:260
MoFEM::DMSetOperators_MoFEM
PetscErrorCode DMSetOperators_MoFEM(DM dm)
Set operators for MoFEM dm.
Definition: DMMoFEM.cpp:53
PCMGSubMatrixCtx_private::MOFEM_EVENT_sor
PetscLogEvent MOFEM_EVENT_sor
Definition: PCMGSetUpViaApproxOrders.cpp:83
PCMGSubMatrixCtx_private::initData
MoFEMErrorCode initData(Vec x)
Definition: PCMGSetUpViaApproxOrders.cpp:71
DMCreateGlobalVector_MGViaApproxOrders
MoFEMErrorCode DMCreateGlobalVector_MGViaApproxOrders(DM dm, Vec *g)
Create global vector for DMGViaApproxOrders.
Definition: PCMGSetUpViaApproxOrders.cpp:576
PCMGSetUpViaApproxOrdersCtx::orderAtLastLevel
int orderAtLastLevel
set maximal evaluated order
Definition: PCMGSetUpViaApproxOrders.hpp:206
convert.n
n
Definition: convert.py:82
PCMGSetUpViaApproxOrdersCtx::destroyIsAtLevel
virtual MoFEMErrorCode destroyIsAtLevel(int kk, IS *is)
Destroy IS if internally created.
Definition: PCMGSetUpViaApproxOrders.cpp:653
PCMGSubMatrixCtx::A
Mat A
Definition: PCMGSetUpViaApproxOrders.hpp:10
N
const int N
Definition: speed_test.cpp:3
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
sub_mat_mult_add
MoFEMErrorCode sub_mat_mult_add(Mat a, Vec x, Vec f)
Definition: PCMGSetUpViaApproxOrders.cpp:110
DMRegister_MGViaApproxOrders
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:376
PCMGSubMatrixCtx::sCat
VecScatter sCat
Definition: PCMGSetUpViaApproxOrders.hpp:13
DMMGViaApproxOrdersGetCtx
MoFEMErrorCode DMMGViaApproxOrdersGetCtx(DM dm, DMMGViaApproxOrdersCtx **ctx)
Definition: PCMGSetUpViaApproxOrders.cpp:175
DMMGViaApproxOrdersCtx::coarseningIS
std::vector< IS > coarseningIS
Coarsening IS.
Definition: PCMGSetUpViaApproxOrders.hpp:33
PCMGSetUpViaApproxOrdersCtx::verboseLevel
int verboseLevel
Definition: PCMGSetUpViaApproxOrders.hpp:209
MoFEM::DMoFEMGetInterfacePtr
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMoFEM.cpp:418
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
PCMGSetUpViaApproxOrdersCtx::coarseOrder
int coarseOrder
approximation order of coarse level
Definition: PCMGSetUpViaApproxOrders.hpp:205
PCMGSubMatrixCtx_private::sub_mat_mult_generic
friend MoFEMErrorCode sub_mat_mult_generic(Mat a, Vec x, Vec f)
Definition: PCMGSetUpViaApproxOrders.cpp:88
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::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:430
PCMGSubMatrixCtx::PCMGSubMatrixCtx
PCMGSubMatrixCtx(Mat a, IS is)
Definition: PCMGSetUpViaApproxOrders.cpp:32
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:114
ksp_set_operators
static MoFEMErrorCode ksp_set_operators(KSP ksp, Mat A, Mat B, void *ctx)
Definition: PCMGSetUpViaApproxOrders.cpp:404
PCMGSetUpViaApproxOrdersCtx::buildProlongationOperator
virtual MoFEMErrorCode buildProlongationOperator(bool use_mat_a, int verb=0)
Set up data structures for MG.
Definition: PCMGSetUpViaApproxOrders.cpp:660
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:477
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
DMMGViaApproxOrdersReplaceCoarseningIS
MoFEMErrorCode DMMGViaApproxOrdersReplaceCoarseningIS(DM dm, IS *is_vec, int nb_elems, Mat A, int verb)
Replace coarsening IS in DM via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:285
DMMGViaApproxOrdersCtx::aO
AO aO
Definition: PCMGSetUpViaApproxOrders.hpp:32
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
QUIET
@ QUIET
Definition: definitions.h:208
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
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:197
MoFEM::Problem::getName
auto getName() const
Definition: ProblemsMultiIndices.hpp:372
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
PCMGSetUpViaApproxOrdersCtx::A
Mat A
Matrix at fine level.
Definition: PCMGSetUpViaApproxOrders.hpp:196
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
PCMGSubMatrixCtx_private::isInitisalised
bool isInitisalised
Definition: PCMGSetUpViaApproxOrders.cpp:84