v0.9.1
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>
99 MoFEMErrorCode sub_mat_mult_generic(Mat a, Vec x, Vec f) {
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 
117 MoFEMErrorCode sub_mat_mult(Mat a, Vec x, Vec f) {
118  return sub_mat_mult_generic<INSERT_VALUES>(a, x, f);
119 }
120 
121 MoFEMErrorCode sub_mat_mult_add(Mat a, Vec x, Vec f) {
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 
175  MoFEM::UnknownInterface **iface) const {
177  *iface = NULL;
178  if (uuid == IDD_DMMGVIAAPPROXORDERSCTX) {
179  *iface = static_cast<DMMGViaApproxOrdersCtx *>(
180  const_cast<DMMGViaApproxOrdersCtx *>(this));
182  } else {
183  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "wrong interference");
184  }
185 
186  ierr = DMCtx::query_interface(uuid, iface);
187  CHKERRG(ierr);
189 }
190 
191 #define GET_DM_FIELD(DM) \
192  MoFEM::UnknownInterface *iface; \
193  CHKERR((DMCtx *)DM->data) \
194  ->query_interface(IDD_DMMGVIAAPPROXORDERSCTX, &iface); \
195  DMMGViaApproxOrdersCtx *dm_field = \
196  static_cast<DMMGViaApproxOrdersCtx *>(iface); \
197  NOT_USED(dm_field)
198 
200  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
202  GET_DM_FIELD(dm);
203  *ctx = dm_field;
205 }
206 
208  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
210  GET_DM_FIELD(dm);
211  if (dm_field->aO) {
212  // std::cerr << dm_field->aO << std::endl;
213  CHKERR AODestroy(&dm_field->aO);
214  // std::cerr << "destroy ao when adding\n";
215  }
216  dm_field->aO = ao;
217  CHKERR PetscObjectReference((PetscObject)ao);
218  // std::cerr << "add ao\n";
220 }
221 
223  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
225  GET_DM_FIELD(dm);
226  *size = dm_field->coarseningIS.size();
228 }
229 
231  Mat *subA,
232  bool create_sub_matrix,
233  bool shell_sub_a) {
234  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
236  GET_DM_FIELD(dm);
237  dm_field->coarseningIS.push_back(is);
238  dm_field->shellMatrixCtxPtr.push_back(new PCMGSubMatrixCtx_private(A, is));
239  if (is) {
240  CHKERR PetscObjectReference((PetscObject)is);
241  }
242  if (is) {
243  IS is2 = is;
244  if (dm_field->aO) {
245  CHKERR ISDuplicate(is, &is2);
246  CHKERR ISCopy(is, is2);
247  CHKERR AOApplicationToPetscIS(dm_field->aO, is2);
248  }
249  if (create_sub_matrix) {
250  if (shell_sub_a) {
251  int n, N;
252  CHKERR ISGetSize(is, &N);
253  CHKERR ISGetLocalSize(is, &n);
254  MPI_Comm comm;
255  CHKERR PetscObjectGetComm((PetscObject)A, &comm);
256  CHKERR MatCreateShell(comm, n, n, N, N,
257  &(dm_field->shellMatrixCtxPtr.back()), subA);
258  CHKERR MatShellSetOperation(*subA, MATOP_MULT,
259  (void (*)(void))sub_mat_mult);
260  CHKERR MatShellSetOperation(*subA, MATOP_MULT_ADD,
261  (void (*)(void))sub_mat_mult_add);
262  CHKERR MatShellSetOperation(*subA, MATOP_SOR,
263  (void (*)(void))sub_mat_sor);
264  } else {
265  #if PETSC_VERSION_GE(3, 8, 0)
266  CHKERR MatCreateSubMatrix(A, is2, is2, MAT_INITIAL_MATRIX, subA);
267  #else
268  CHKERR MatGetSubMatrix(A, is2, is2, MAT_INITIAL_MATRIX, subA);
269  #endif
270  }
271  }
272  if (dm_field->aO) {
273  CHKERR ISDestroy(&is2);
274  }
275  dm_field->kspOperators.push_back(*subA);
276  CHKERR PetscObjectReference((PetscObject)(*subA));
277  } else {
278  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
279  }
280  PetscInfo(dm, "Push back IS to DMMGViaApproxOrders\n");
282 }
283 
285  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
287  GET_DM_FIELD(dm);
288  if (dm_field->coarseningIS.back()) {
289  CHKERR ISDestroy(&dm_field->coarseningIS.back());
290  dm_field->coarseningIS.pop_back();
291  }
292  if (dm_field->kspOperators.back()) {
293  CHKERR MatDestroy(&dm_field->kspOperators.back());
294  }
295  dm_field->kspOperators.pop_back();
296  PetscInfo(dm, "Pop back IS to DMMGViaApproxOrders\n");
298 }
299 
301  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
303  GET_DM_FIELD(dm);
304  CHKERR dm_field->destroyCoarseningIS();
305  PetscInfo(dm, "Clear DMs data structures\n");
307 }
308 
310  int nb_elems, Mat A,
311  int verb) {
312  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
314  GET_DM_FIELD(dm);
315  int nb_no_changed = 0;
316  int nb_replaced = 0;
317  int nb_deleted = 0;
318  int nb_added = 0;
319  std::vector<IS>::iterator it;
320  it = dm_field->coarseningIS.begin();
321  int ii = 0;
322  for (; it != dm_field->coarseningIS.end(); it++, ii++) {
323  if (ii < nb_elems) {
324  PetscBool flg;
325  CHKERR ISEqual(*it, is_vec[ii], &flg);
326  if (!flg) {
327  CHKERR ISDestroy(&*it);
328  CHKERR MatDestroy(&dm_field->kspOperators[ii]);
329  *it = is_vec[ii];
330  CHKERR PetscObjectReference((PetscObject)is_vec[ii]);
331  if (ii < nb_elems - 1) {
332  IS is = is_vec[ii];
333  if (dm_field->aO) {
334  CHKERR ISDuplicate(is_vec[ii], &is);
335  CHKERR ISCopy(is_vec[ii], is);
336  CHKERR AOApplicationToPetscIS(dm_field->aO, is);
337  }
338  Mat subA;
339  #if PETSC_VERSION_GE(3, 8, 0)
340  CHKERR MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &subA);
341  #else
342  CHKERR MatGetSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &subA);
343  #endif
344  CHKERR PetscObjectReference((PetscObject)subA);
345  dm_field->kspOperators[ii] = subA;
346  CHKERR MatDestroy(&subA);
347  if (dm_field->aO) {
348  CHKERR ISDestroy(&is);
349  }
350  } else {
351  CHKERR PetscObjectReference((PetscObject)A);
352  dm_field->kspOperators[ii] = A;
353  }
354  nb_replaced++;
355  }
356  } else {
357  nb_no_changed++;
358  continue;
359  }
360  }
361  if (static_cast<int>(dm_field->coarseningIS.size()) < nb_elems) {
362  for (; ii < nb_elems - 1; ii++) {
363  Mat subA;
364  CHKERR DMMGViaApproxOrdersPushBackCoarseningIS(dm, is_vec[ii], A, &subA,
365  true, false);
366  CHKERR MatDestroy(&subA);
367  nb_added++;
368  }
369  CHKERR DMMGViaApproxOrdersPushBackCoarseningIS(dm, is_vec[ii], A, &A, false,
370  false);
371  nb_added++;
372  } else {
373  for (; ii < static_cast<int>(dm_field->coarseningIS.size()); ii++) {
375  nb_deleted++;
376  }
377  }
378  MPI_Comm comm;
379  CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
380  if (verb > 0) {
381  PetscPrintf(comm,
382  "DMMGViaApproxOrders nb_no_changed = %d, nb_replaced = %d, "
383  "nb_added = %d, nb_deleted = %d, size = %d\n",
384  nb_no_changed, nb_replaced, nb_added, nb_deleted,
385  dm_field->coarseningIS.size());
386  }
387  PetscInfo(dm, "Replace IS to DMMGViaApproxOrders\n");
389 }
390 
392  const DMMGViaApproxOrdersCtx **ctx) {
393  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
395  GET_DM_FIELD(dm);
396  *ctx = dm_field;
398 }
399 
402  CHKERR DMRegister(sname, DMCreate_MGViaApproxOrders);
404 }
405 
406 // MoFEMErrorCode DMDestroy_MGViaApproxOrders(DM dm) {
407 // PetscValidHeaderSpecific(dm,DM_CLASSID,1);
408 // MoFEMFunctionBeginHot;
409 // if(!((DMMGViaApproxOrdersCtx*)dm->data)->referenceNumber) {
410 // DMMGViaApproxOrdersCtx *dm_field = (DMMGViaApproxOrdersCtx*)dm->data;
411 // if(dm_field->destroyProblem) {
412 // if(dm_field->mField_ptr->check_problem(dm_field->problemName)) {
413 // dm_field->mField_ptr->delete_problem(dm_field->problemName);
414 // } // else problem has to be deleted by the user
415 // }
416 // cerr << "Destroy " << dm_field->problemName << endl;
417 // delete (DMMGViaApproxOrdersCtx*)dm->data;
418 // } else {
419 // DMMGViaApproxOrdersCtx *dm_field = (DMMGViaApproxOrdersCtx*)dm->data;
420 //
421 // cerr << "Derefrence " << dm_field->problemName << " " <<
422 // ((DMCtx*)dm->data)->referenceNumber << endl;
423 // (((DMMGViaApproxOrdersCtx*)dm->data)->referenceNumber)--;
424 // }
425 // MoFEMFunctionReturnHot(0);
426 // }
427 
428 static MoFEMErrorCode ksp_set_operators(KSP ksp, Mat A, Mat B, void *ctx) {
431 }
432 
434  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
436  if (!dm->data) {
437  dm->data = new DMMGViaApproxOrdersCtx();
438  } else {
439  ((DMCtx *)(dm->data))->referenceNumber++;
440  }
441  // cerr << "Create " << ((DMCtx*)(dm->data))->referenceNumber << endl;
443  dm->ops->creatematrix = DMCreateMatrix_MGViaApproxOrders;
444  dm->ops->createglobalvector = DMCreateGlobalVector_MGViaApproxOrders;
445  dm->ops->coarsen = DMCoarsen_MGViaApproxOrders;
446  // dm->ops->destroy = DMDestroy_MGViaApproxOrders;
447  dm->ops->createinterpolation = DMCreateInterpolation_MGViaApproxOrders;
448  CHKERR DMKSPSetComputeOperators(dm, ksp_set_operators, NULL);
449  PetscInfo1(dm, "Create DMMGViaApproxOrders reference = %d\n",
450  ((DMCtx *)(dm->data))->referenceNumber);
452 }
453 
455 
456  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
458  GET_DM_FIELD(dm);
459 
460  int leveldown = dm->leveldown;
461 
462  if (dm_field->kspOperators.empty()) {
464  } else {
465  MPI_Comm comm;
466  CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
467  if (dm_field->kspOperators.empty()) {
468  SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
469  "data inconsistency, operator can not be set");
470  }
471  if (static_cast<int>(dm_field->kspOperators.size()) < leveldown) {
472  SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
473  "data inconsistency, no IS for that level");
474  }
475  *M = dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown];
476  CHKERR PetscObjectReference((PetscObject)*M);
477  }
478 
479  PetscInfo1(dm, "Create Matrix DMMGViaApproxOrders leveldown = %d\n",
480  leveldown);
481 
483 }
484 
485 MoFEMErrorCode DMCoarsen_MGViaApproxOrders(DM dm, MPI_Comm comm, DM *dmc) {
486  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
488  GET_DM_FIELD(dm);
489  CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
490  CHKERR DMCreate(comm, dmc);
491  (*dmc)->data = dm->data;
492  DMType type;
493  CHKERR DMGetType(dm, &type);
494  CHKERR DMSetType(*dmc, type);
495  CHKERR PetscObjectReference((PetscObject)(*dmc));
496  PetscInfo1(dm, "Coarsen DMMGViaApproxOrders leveldown = %d\n", dm->leveldown);
498 }
499 
501  Vec *vec) {
502  PetscValidHeaderSpecific(dm1, DM_CLASSID, 1);
503  PetscValidHeaderSpecific(dm2, DM_CLASSID, 1);
505 
506  MPI_Comm comm;
507  CHKERR PetscObjectGetComm((PetscObject)dm1, &comm);
508 
509  int m, n, M, N;
510 
511  DM dm_down = dm1;
512  DM dm_up = dm2;
513 
514  int dm_down_leveldown = dm_down->leveldown;
515  int dm_up_leveldown = dm_up->leveldown;
516 
517  PetscInfo2(dm1,
518  "Create interpolation DMMGViaApproxOrders dm1_leveldown = %d "
519  "dm2_leveldown = %d\n",
520  dm_down_leveldown, dm_up_leveldown);
521 
522  IS is_down, is_up;
523  {
524  // Coarser mesh
525  GET_DM_FIELD(dm_down);
526  if (static_cast<int>(dm_field->coarseningIS.size()) < dm_down_leveldown) {
527  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
528  }
529  is_down = dm_field->coarseningIS[dm_field->coarseningIS.size() - 1 -
530  dm_down_leveldown];
531  CHKERR ISGetSize(is_down, &M);
532  CHKERR ISGetLocalSize(is_down, &m);
533  }
534  {
535  // Finer mesh
536  GET_DM_FIELD(dm_up);
537  if (static_cast<int>(dm_field->coarseningIS.size()) < dm_up_leveldown) {
538  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
539  }
540  is_up =
541  dm_field
542  ->coarseningIS[dm_field->coarseningIS.size() - 1 - dm_up_leveldown];
543  CHKERR ISGetSize(is_up, &N);
544  CHKERR ISGetLocalSize(is_up, &n);
545  }
546 
547  // is_dow rows
548  // is_up columns
549 
550  CHKERR MatCreate(comm, mat);
551  CHKERR MatSetSizes(*mat, m, n, M, N);
552  CHKERR MatSetType(*mat, MATMPIAIJ);
553  CHKERR MatMPIAIJSetPreallocation(*mat, 1, PETSC_NULL, 0, PETSC_NULL);
554 
555  // get matrix layout
556  PetscLayout rmap, cmap;
557  CHKERR MatGetLayouts(*mat, &rmap, &cmap);
558  int rstart, rend, cstart, cend;
559  CHKERR PetscLayoutGetRange(rmap, &rstart, &rend);
560  CHKERR PetscLayoutGetRange(cmap, &cstart, &cend);
561 
562  // if(verb>0) {
563  // PetscSynchronizedPrintf(comm,"level %d row start %d row end
564  // %d\n",kk,rstart,rend); PetscSynchronizedPrintf(comm,"level %d col start
565  // %d col end %d\n",kk,cstart,cend);
566  // }
567 
568  const int *row_indices_ptr, *col_indices_ptr;
569  CHKERR ISGetIndices(is_down, &row_indices_ptr);
570  CHKERR ISGetIndices(is_up, &col_indices_ptr);
571 
572  map<int, int> idx_map;
573  for (int ii = 0; ii < m; ii++) {
574  idx_map[row_indices_ptr[ii]] = rstart + ii;
575  }
576 
577  CHKERR MatZeroEntries(*mat);
578  // FIXME: Use MatCreateMPIAIJWithArrays and set array directly
579  for (int jj = 0; jj < n; jj++) {
580  map<int, int>::iterator mit = idx_map.find(col_indices_ptr[jj]);
581  if (mit != idx_map.end()) {
582  CHKERR MatSetValue(*mat, mit->second, cstart + jj, 1, INSERT_VALUES);
583  }
584  }
585 
586  CHKERR ISRestoreIndices(is_down, &row_indices_ptr);
587  CHKERR ISRestoreIndices(is_up, &col_indices_ptr);
588 
589  CHKERR MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
590  CHKERR MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);
591 
592  if (vec != NULL) {
593  *vec = PETSC_NULL;
594  }
595 
597 }
598 
600 
601  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
603  int leveldown = dm->leveldown;
604  GET_DM_FIELD(dm);
605  if (dm_field->kspOperators.empty()) {
607  } else {
608 #if PETSC_VERSION_GE(3, 5, 3)
609  CHKERR MatCreateVecs(
610  dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
611  g, NULL);
612 #else
613  CHKERR MatGetVecs(
614  dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
615  g, NULL);
616 #endif
617  }
618  PetscInfo1(dm, "Create global vector DMMGViaApproxOrders leveldown = %d\n",
619  dm->leveldown);
621 }
622 
625  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
626  "MOFEM Multi-Grid (Orders) pre-conditioner", "none");
627  CHKERRG(ierr);
628 
629  CHKERR PetscOptionsInt("-mofem_mg_levels", "nb levels of multi-grid solver",
630  "", 2, &nbLevels, PETSC_NULL);
631  CHKERR PetscOptionsInt("-mofem_mg_coarse_order",
632  "approximation order of coarse level", "", 2,
633  &coarseOrder, PETSC_NULL);
634  CHKERR PetscOptionsInt("-mofem_mg_order_at_last_level", "order at last level",
635  "", 100, &orderAtLastLevel, PETSC_NULL);
636  CHKERR PetscOptionsInt("-mofem_mg_verbose", "nb levels of multi-grid solver",
637  "", 0, &verboseLevel, PETSC_NULL);
638  PetscBool shell_sub_a = shellSubA ? PETSC_TRUE : PETSC_FALSE;
639  CHKERR PetscOptionsBool("-mofem_mg_shell_a", "use shell matrix as sub matrix",
640  "", shell_sub_a, &shell_sub_a, NULL);
641  shellSubA = (shellSubA == PETSC_TRUE);
642 
643  ierr = PetscOptionsEnd();
644  CHKERRG(ierr);
646 }
647 
649  MoFEM::Interface *m_field_ptr;
650  MoFEM::ISManager *is_manager_ptr;
652  // if is last level, take all remaining orders dofs, if any left
653  CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
654  CHKERR m_field_ptr->getInterface(is_manager_ptr);
655  const Problem *problem_ptr;
656  CHKERR DMMoFEMGetProblemPtr(dM, &problem_ptr);
657  int order_at_next_level = kk + coarseOrder;
658  if (kk == nbLevels - 1) {
659  int first = problem_ptr->getNumeredDofsRows()
660  ->get<PetscLocalIdx_mi_tag>()
661  .find(0)
662  ->get()
663  ->getPetscGlobalDofIdx();
664  CHKERR ISCreateStride(PETSC_COMM_WORLD, problem_ptr->getNbLocalDofsRow(),
665  first, 1, is);
667  // order_at_next_level = orderAtLastLevel;
668  }
669  string problem_name = problem_ptr->getName();
670  CHKERR is_manager_ptr->isCreateProblemOrder(problem_name, ROW, 0,
671  order_at_next_level, is);
673 }
674 
677  CHKERR ISDestroy(is);
679 }
680 
683  int verb) {
685  verb = verb > verboseLevel ? verb : verboseLevel;
686 
687  MPI_Comm comm;
688  CHKERR PetscObjectGetComm((PetscObject)dM, &comm);
689 
690  if (verb > QUIET) {
691  PetscPrintf(comm, "set MG levels %u\n", nbLevels);
692  }
693 
694  std::vector<IS> is_vec(nbLevels + 1);
695  std::vector<int> is_glob_size(nbLevels + 1), is_loc_size(nbLevels + 1);
696 
697  for (int kk = 0; kk < nbLevels; kk++) {
698 
699  // get indices up to up to give approximation order
700  CHKERR createIsAtLevel(kk, &is_vec[kk]);
701  CHKERR ISGetSize(is_vec[kk], &is_glob_size[kk]);
702  CHKERR ISGetLocalSize(is_vec[kk], &is_loc_size[kk]);
703 
704  if (verb > QUIET) {
705  PetscSynchronizedPrintf(comm,
706  "Nb. dofs at level [ %d ] global %u local %d\n",
707  kk, is_glob_size[kk], is_loc_size[kk]);
708  }
709 
710  // if no dofs on level kk finish here
711  if (is_glob_size[kk] == 0) {
712  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "no dofs at level");
713  }
714  }
715 
716  for (int kk = 0; kk != nbLevels; kk++) {
717  Mat subA;
718  if (kk == nbLevels - 1 && use_mat_a) {
719  subA = A;
721  false, false);
722  } else {
723  if (kk > 0) {
724  // Not coarse level
726  true, shellSubA);
727  } else {
728  // Coarse lave is compressed matrix allowing for factorization when
729  // needed
731  true, false);
732  }
733  if (subA) {
734  CHKERR MatDestroy(&subA);
735  }
736  }
737  }
738 
739  for (unsigned int kk = 0; kk < is_vec.size(); kk++) {
740  CHKERR destroyIsAtLevel(kk, &is_vec[kk]);
741  }
742 
743  if (verb > QUIET) {
744  PetscSynchronizedFlush(comm, PETSC_STDOUT);
745  }
746 
748 }
749 
751  int verb) {
752 
753  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
755 
756  MPI_Comm comm;
757  CHKERR PetscObjectGetComm((PetscObject)pc, &comm);
758  if (verb > 0) {
759  PetscPrintf(comm, "Start PCMGSetUpViaApproxOrders\n");
760  }
761 
762  CHKERR ctx->getOptions();
763  CHKERR ctx->buildProlongationOperator(true, verb);
764 
765 #if PETSC_VERSION_GE(3, 8, 0)
766  CHKERR PCMGSetGalerkin(pc, PC_MG_GALERKIN_NONE);
767 #else
768  CHKERR PCMGSetGalerkin(pc, PETSC_FALSE);
769 #endif
770 
771  CHKERR PCMGSetLevels(pc, ctx->nbLevels, NULL);
772 
773  if (verb > 0) {
774  PetscPrintf(comm, "End PCMGSetUpViaApproxOrders\n");
775  }
776 
778 }
MoFEMErrorCode sub_mat_mult_add(Mat a, Vec x, Vec f)
Deprecated interface functions.
MoFEM interface unique ID.
MoFEMErrorCode DMCoarsen_MGViaApproxOrders(DM dm, MPI_Comm comm, DM *dmc)
Coarsen DM.
MoFEMErrorCode DMMGViaApproxOrdersReplaceCoarseningIS(DM dm, IS *is_vec, int nb_elems, Mat A, int verb)
Replace coarsening IS in DM via approximation orders.
DM dM
Distributed mesh manager.
MoFEMErrorCode DMMGViaApproxOrdersClearCoarseningIS(DM dm)
Clear approximation orders.
Section manager is used to create indexes and sectionsFIXME: ISManager is not properly testsed by ato...
Definition: ISManager.hpp:36
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
static MoFEMErrorCode ksp_set_operators(KSP ksp, Mat A, Mat B, void *ctx)
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.
base class for all interface classes
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
static const MOFEMuuid IDD_DMMGVIAAPPROXORDERSCTX
boost::ptr_vector< PCMGSubMatrixCtx > shellMatrixCtxPtr
Shell sub-matrix context.
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:550
DofIdx getNbLocalDofsRow() const
MoFEMErrorCode DMCreateInterpolation_MGViaApproxOrders(DM dm1, DM dm2, Mat *mat, Vec *vec)
Create interpolation matrix between data managers dm1 and dm2.
MoFEMErrorCode DMCreateGlobalVector_MGViaApproxOrders(DM dm, Vec *g)
Create global vector for DMGViaApproxOrders.
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, MoFEM::UnknownInterface **iface) const
int nbLevels
number of multi-grid levels
MoFEM interface.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
MoFEMErrorCode sub_mat_sor(Mat mat, Vec b, PetscReal omega, MatSORType flag, PetscReal shift, PetscInt its, PetscInt lits, Vec x)
keeps basic data about problemThis is low level structure with information about problem,...
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMMoFEM.cpp:1053
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVectorsets the routine to create a global vector associated with the shell DM.
Definition: DMMMoFEM.cpp:1026
MoFEMErrorCode DMMGViaApproxOrdersGetCoarseningISSize(DM dm, int *size)
Gets size of coarseningIS in internal data struture DMMGViaApproxOrdersCtx.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
#define GET_DM_FIELD(DM)
Structure for DM for multi-grid via approximation orders.
MoFEMErrorCode DMCreate_MGViaApproxOrders(DM dm)
Create DM data structure for Multi-Grid via approximation orders.
MoFEMErrorCode DMMGViaApproxOrdersSetAO(DM dm, AO ao)
Set DM ordering.
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
MoFEMErrorCode DMCreateMatrix_MGViaApproxOrders(DM dm, Mat *M)
Create matrix for Multi-Grid via approximation orders.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
header of multi-grid solver for p- adaptivity
int orderAtLastLevel
set maximal evaluated order
Set data structures of MG pre-conditioner via approximation orders.
virtual MoFEMErrorCode buildProlongationOperator(bool use_mat_a, int verb=0)
Set up data structures for MG.
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
MoFEMErrorCode DMMGViaApproxOrdersPopBackCoarseningIS(DM dm)
Pop is form MG via approximation orders.
PetscErrorCode DMSetOperators_MoFEM(DM dm)
Set operators for MoFEM dm.
Definition: DMMMoFEM.cpp:54
MoFEMErrorCode sub_mat_mult_generic(Mat a, Vec x, Vec f)
#define CHKERR
Inline error check.
Definition: definitions.h:602
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:72
int coarseOrder
approximation order of coarse level
friend MoFEMErrorCode sub_mat_sor(Mat mat, Vec b, PetscReal omega, MatSORType flag, PetscReal shift, PetscInt its, PetscInt lits, Vec x)
virtual MoFEMErrorCode createIsAtLevel(int kk, IS *is)
Set IS for levels.
MoFEMErrorCode DMMGViaApproxOrdersGetCtx(DM dm, DMMGViaApproxOrdersCtx **ctx)
friend MoFEMErrorCode sub_mat_mult_generic(Mat a, Vec x, Vec f)
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMMoFEM.cpp:336
PETSc Discrete Manager data structure.
Definition: DMMoFEM.hpp:877
virtual MoFEMErrorCode destroyIsAtLevel(int kk, IS *is)
Destroy IS if internally created.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:348
MoFEMErrorCode isCreateProblemOrder(const std::string &problem, RowColData rc, int min_order, int max_order, IS *is) const
create IS for given order range (collective)
Definition: ISManager.cpp:176
std::string getName() const
std::vector< Mat > kspOperators
Get KSP operators.
MoFEMErrorCode sub_mat_mult(Mat a, Vec x, Vec f)
const int N
Definition: speed_test.cpp:3
std::vector< IS > coarseningIS
Coarsening IS.
virtual MoFEMErrorCode getOptions()
get options from line command
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, PCMGSetUpViaApproxOrdersCtx *ctx, int verb)
Function build MG structure.
const boost::shared_ptr< NumeredDofEntity_multiIndex > & getNumeredDofsRows() const
get access to numeredDofsRows storing DOFs on rows