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