v0.8.16
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  CHKERR MatGetSubMatrix(A, is2, is2, MAT_INITIAL_MATRIX, subA);
266  }
267  }
268  if (dm_field->aO) {
269  CHKERR ISDestroy(&is2);
270  }
271  dm_field->kspOperators.push_back(*subA);
272  CHKERR PetscObjectReference((PetscObject)(*subA));
273  } else {
274  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
275  }
276  PetscInfo(dm, "Push back IS to DMMGViaApproxOrders\n");
278 }
279 
281  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
283  GET_DM_FIELD(dm);
284  if (dm_field->coarseningIS.back()) {
285  CHKERR ISDestroy(&dm_field->coarseningIS.back());
286  dm_field->coarseningIS.pop_back();
287  }
288  if (dm_field->kspOperators.back()) {
289  CHKERR MatDestroy(&dm_field->kspOperators.back());
290  }
291  dm_field->kspOperators.pop_back();
292  PetscInfo(dm, "Pop back IS to DMMGViaApproxOrders\n");
294 }
295 
297  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
299  GET_DM_FIELD(dm);
300  CHKERR dm_field->destroyCoarseningIS();
301  PetscInfo(dm, "Clear DMs data structures\n");
303 }
304 
306  int nb_elems, Mat A,
307  int verb) {
308  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
310  GET_DM_FIELD(dm);
311  int nb_no_changed = 0;
312  int nb_replaced = 0;
313  int nb_deleted = 0;
314  int nb_added = 0;
315  std::vector<IS>::iterator it;
316  it = dm_field->coarseningIS.begin();
317  int ii = 0;
318  for (; it != dm_field->coarseningIS.end(); it++, ii++) {
319  if (ii < nb_elems) {
320  PetscBool flg;
321  CHKERR ISEqual(*it, is_vec[ii], &flg);
322  if (!flg) {
323  CHKERR ISDestroy(&*it);
324  CHKERR MatDestroy(&dm_field->kspOperators[ii]);
325  *it = is_vec[ii];
326  CHKERR PetscObjectReference((PetscObject)is_vec[ii]);
327  if (ii < nb_elems - 1) {
328  IS is = is_vec[ii];
329  if (dm_field->aO) {
330  CHKERR ISDuplicate(is_vec[ii], &is);
331  CHKERR ISCopy(is_vec[ii], is);
332  CHKERR AOApplicationToPetscIS(dm_field->aO, is);
333  }
334  Mat subA;
335  CHKERR MatGetSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &subA);
336  CHKERR PetscObjectReference((PetscObject)subA);
337  dm_field->kspOperators[ii] = subA;
338  CHKERR MatDestroy(&subA);
339  if (dm_field->aO) {
340  CHKERR ISDestroy(&is);
341  }
342  } else {
343  CHKERR PetscObjectReference((PetscObject)A);
344  dm_field->kspOperators[ii] = A;
345  }
346  nb_replaced++;
347  }
348  } else {
349  nb_no_changed++;
350  continue;
351  }
352  }
353  if (static_cast<int>(dm_field->coarseningIS.size()) < nb_elems) {
354  for (; ii < nb_elems - 1; ii++) {
355  Mat subA;
356  CHKERR DMMGViaApproxOrdersPushBackCoarseningIS(dm, is_vec[ii], A, &subA,
357  true, false);
358  CHKERR MatDestroy(&subA);
359  nb_added++;
360  }
361  CHKERR DMMGViaApproxOrdersPushBackCoarseningIS(dm, is_vec[ii], A, &A, false,
362  false);
363  nb_added++;
364  } else {
365  for (; ii < static_cast<int>(dm_field->coarseningIS.size()); ii++) {
367  nb_deleted++;
368  }
369  }
370  MPI_Comm comm;
371  CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
372  if (verb > 0) {
373  PetscPrintf(comm,
374  "DMMGViaApproxOrders nb_no_changed = %d, nb_replaced = %d, "
375  "nb_added = %d, nb_deleted = %d, size = %d\n",
376  nb_no_changed, nb_replaced, nb_added, nb_deleted,
377  dm_field->coarseningIS.size());
378  }
379  PetscInfo(dm, "Replace IS to DMMGViaApproxOrders\n");
381 }
382 
384  const DMMGViaApproxOrdersCtx **ctx) {
385  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
387  GET_DM_FIELD(dm);
388  *ctx = dm_field;
390 }
391 
394  CHKERR DMRegister(sname, DMCreate_MGViaApproxOrders);
396 }
397 
398 // MoFEMErrorCode DMDestroy_MGViaApproxOrders(DM dm) {
399 // PetscValidHeaderSpecific(dm,DM_CLASSID,1);
400 // MoFEMFunctionBeginHot;
401 // if(!((DMMGViaApproxOrdersCtx*)dm->data)->referenceNumber) {
402 // DMMGViaApproxOrdersCtx *dm_field = (DMMGViaApproxOrdersCtx*)dm->data;
403 // if(dm_field->destroyProblem) {
404 // if(dm_field->mField_ptr->check_problem(dm_field->problemName)) {
405 // dm_field->mField_ptr->delete_problem(dm_field->problemName);
406 // } // else problem has to be deleted by the user
407 // }
408 // cerr << "Destroy " << dm_field->problemName << endl;
409 // delete (DMMGViaApproxOrdersCtx*)dm->data;
410 // } else {
411 // DMMGViaApproxOrdersCtx *dm_field = (DMMGViaApproxOrdersCtx*)dm->data;
412 //
413 // cerr << "Derefrence " << dm_field->problemName << " " <<
414 // ((DMCtx*)dm->data)->referenceNumber << endl;
415 // (((DMMGViaApproxOrdersCtx*)dm->data)->referenceNumber)--;
416 // }
417 // MoFEMFunctionReturnHot(0);
418 // }
419 
420 static MoFEMErrorCode ksp_set_operators(KSP ksp, Mat A, Mat B, void *ctx) {
423 }
424 
426  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
428  if (!dm->data) {
429  dm->data = new DMMGViaApproxOrdersCtx();
430  } else {
431  ((DMCtx *)(dm->data))->referenceNumber++;
432  }
433  // cerr << "Create " << ((DMCtx*)(dm->data))->referenceNumber << endl;
435  dm->ops->creatematrix = DMCreateMatrix_MGViaApproxOrders;
436  dm->ops->createglobalvector = DMCreateGlobalVector_MGViaApproxOrders;
437  dm->ops->coarsen = DMCoarsen_MGViaApproxOrders;
438  // dm->ops->destroy = DMDestroy_MGViaApproxOrders;
439  dm->ops->createinterpolation = DMCreateInterpolation_MGViaApproxOrders;
440  CHKERR DMKSPSetComputeOperators(dm, ksp_set_operators, NULL);
441  PetscInfo1(dm, "Create DMMGViaApproxOrders reference = %d\n",
442  ((DMCtx *)(dm->data))->referenceNumber);
444 }
445 
447 
448  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
450  GET_DM_FIELD(dm);
451 
452  int leveldown = dm->leveldown;
453 
454  if (dm_field->kspOperators.empty()) {
456  } else {
457  MPI_Comm comm;
458  CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
459  if (dm_field->kspOperators.empty()) {
460  SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
461  "data inconsistency, operator can not be set");
462  }
463  if (static_cast<int>(dm_field->kspOperators.size()) < leveldown) {
464  SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
465  "data inconsistency, no IS for that level");
466  }
467  *M = dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown];
468  CHKERR PetscObjectReference((PetscObject)*M);
469  }
470 
471  PetscInfo1(dm, "Create Matrix DMMGViaApproxOrders leveldown = %d\n",
472  leveldown);
473 
475 }
476 
477 MoFEMErrorCode DMCoarsen_MGViaApproxOrders(DM dm, MPI_Comm comm, DM *dmc) {
478  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
480  GET_DM_FIELD(dm);
481  CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
482  CHKERR DMCreate(comm, dmc);
483  (*dmc)->data = dm->data;
484  DMType type;
485  CHKERR DMGetType(dm, &type);
486  CHKERR DMSetType(*dmc, type);
487  CHKERR PetscObjectReference((PetscObject)(*dmc));
488  PetscInfo1(dm, "Coarsen DMMGViaApproxOrders leveldown = %d\n", dm->leveldown);
490 }
491 
493  Vec *vec) {
494  PetscValidHeaderSpecific(dm1, DM_CLASSID, 1);
495  PetscValidHeaderSpecific(dm2, DM_CLASSID, 1);
497 
498  MPI_Comm comm;
499  CHKERR PetscObjectGetComm((PetscObject)dm1, &comm);
500 
501  int m, n, M, N;
502 
503  DM dm_down = dm1;
504  DM dm_up = dm2;
505 
506  int dm_down_leveldown = dm_down->leveldown;
507  int dm_up_leveldown = dm_up->leveldown;
508 
509  PetscInfo2(dm1,
510  "Create interpolation DMMGViaApproxOrders dm1_leveldown = %d "
511  "dm2_leveldown = %d\n",
512  dm_down_leveldown, dm_up_leveldown);
513 
514  IS is_down, is_up;
515  {
516  // Coarser mesh
517  GET_DM_FIELD(dm_down);
518  if (static_cast<int>(dm_field->coarseningIS.size()) < dm_down_leveldown) {
519  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
520  }
521  is_down = dm_field->coarseningIS[dm_field->coarseningIS.size() - 1 -
522  dm_down_leveldown];
523  CHKERR ISGetSize(is_down, &M);
524  CHKERR ISGetLocalSize(is_down, &m);
525  }
526  {
527  // Finer mesh
528  GET_DM_FIELD(dm_up);
529  if (static_cast<int>(dm_field->coarseningIS.size()) < dm_up_leveldown) {
530  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
531  }
532  is_up =
533  dm_field
534  ->coarseningIS[dm_field->coarseningIS.size() - 1 - dm_up_leveldown];
535  CHKERR ISGetSize(is_up, &N);
536  CHKERR ISGetLocalSize(is_up, &n);
537  }
538 
539  // is_dow rows
540  // is_up columns
541 
542  CHKERR MatCreate(comm, mat);
543  CHKERR MatSetSizes(*mat, m, n, M, N);
544  CHKERR MatSetType(*mat, MATMPIAIJ);
545  CHKERR MatMPIAIJSetPreallocation(*mat, 1, PETSC_NULL, 0, PETSC_NULL);
546 
547  // get matrix layout
548  PetscLayout rmap, cmap;
549  CHKERR MatGetLayouts(*mat, &rmap, &cmap);
550  int rstart, rend, cstart, cend;
551  CHKERR PetscLayoutGetRange(rmap, &rstart, &rend);
552  CHKERR PetscLayoutGetRange(cmap, &cstart, &cend);
553 
554  // if(verb>0) {
555  // PetscSynchronizedPrintf(comm,"level %d row start %d row end
556  // %d\n",kk,rstart,rend); PetscSynchronizedPrintf(comm,"level %d col start
557  // %d col end %d\n",kk,cstart,cend);
558  // }
559 
560  const int *row_indices_ptr, *col_indices_ptr;
561  CHKERR ISGetIndices(is_down, &row_indices_ptr);
562  CHKERR ISGetIndices(is_up, &col_indices_ptr);
563 
564  map<int, int> idx_map;
565  for (int ii = 0; ii < m; ii++) {
566  idx_map[row_indices_ptr[ii]] = rstart + ii;
567  }
568 
569  CHKERR MatZeroEntries(*mat);
570  // FIXME: Use MatCreateMPIAIJWithArrays and set array directly
571  for (int jj = 0; jj < n; jj++) {
572  map<int, int>::iterator mit = idx_map.find(col_indices_ptr[jj]);
573  if (mit != idx_map.end()) {
574  CHKERR MatSetValue(*mat, mit->second, cstart + jj, 1, INSERT_VALUES);
575  }
576  }
577 
578  CHKERR ISRestoreIndices(is_down, &row_indices_ptr);
579  CHKERR ISRestoreIndices(is_up, &col_indices_ptr);
580 
581  CHKERR MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
582  CHKERR MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);
583 
584  if (vec != NULL) {
585  *vec = PETSC_NULL;
586  }
587 
589 }
590 
592 
593  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
595  int leveldown = dm->leveldown;
596  GET_DM_FIELD(dm);
597  if (dm_field->kspOperators.empty()) {
599  } else {
600 #if PETSC_VERSION_GE(3, 5, 3)
601  CHKERR MatCreateVecs(
602  dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
603  g, NULL);
604 #else
605  CHKERR MatGetVecs(
606  dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
607  g, NULL);
608 #endif
609  }
610  PetscInfo1(dm, "Create global vector DMMGViaApproxOrders leveldown = %d\n",
611  dm->leveldown);
613 }
614 
617  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
618  "MOFEM Multi-Grid (Orders) pre-conditioner", "none");
619  CHKERRG(ierr);
620 
621  CHKERR PetscOptionsInt("-mofem_mg_levels", "nb levels of multi-grid solver",
622  "", 2, &nbLevels, PETSC_NULL);
623  CHKERR PetscOptionsInt("-mofem_mg_coarse_order",
624  "approximation order of coarse level", "", 2,
625  &coarseOrder, PETSC_NULL);
626  CHKERR PetscOptionsInt("-mofem_mg_order_at_last_level", "order at last level",
627  "", 100, &orderAtLastLevel, PETSC_NULL);
628  CHKERR PetscOptionsInt("-mofem_mg_verbose", "nb levels of multi-grid solver",
629  "", 0, &verboseLevel, PETSC_NULL);
630  PetscBool shell_sub_a = shellSubA ? PETSC_TRUE : PETSC_FALSE;
631  CHKERR PetscOptionsBool("-mofem_mg_shell_a", "use shell matrix as sub matrix",
632  "", shell_sub_a, &shell_sub_a, NULL);
633  shellSubA = (shellSubA == PETSC_TRUE);
634 
635  ierr = PetscOptionsEnd();
636  CHKERRG(ierr);
638 }
639 
641  MoFEM::Interface *m_field_ptr;
642  MoFEM::ISManager *is_manager_ptr;
644  // if is last level, take all remaining orders dofs, if any left
645  CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
646  CHKERR m_field_ptr->getInterface(is_manager_ptr);
647  const Problem *problem_ptr;
648  CHKERR DMMoFEMGetProblemPtr(dM, &problem_ptr);
649  int order_at_next_level = kk + coarseOrder;
650  if (kk == nbLevels - 1) {
651  int first = problem_ptr->getNumeredDofsRows()
652  ->get<PetscLocalIdx_mi_tag>()
653  .find(0)
654  ->get()
655  ->getPetscGlobalDofIdx();
656  CHKERR ISCreateStride(PETSC_COMM_WORLD, problem_ptr->getNbLocalDofsRow(),
657  first, 1, is);
659  // order_at_next_level = orderAtLastLevel;
660  }
661  string problem_name = problem_ptr->getName();
662  CHKERR is_manager_ptr->isCreateProblemOrder(problem_name, ROW, 0,
663  order_at_next_level, is);
665 }
666 
669  CHKERR ISDestroy(is);
671 }
672 
675  int verb) {
677  verb = verb > verboseLevel ? verb : verboseLevel;
678 
679  MPI_Comm comm;
680  CHKERR PetscObjectGetComm((PetscObject)dM, &comm);
681 
682  if (verb > QUIET) {
683  PetscPrintf(comm, "set MG levels %u\n", nbLevels);
684  }
685 
686  std::vector<IS> is_vec(nbLevels + 1);
687  std::vector<int> is_glob_size(nbLevels + 1), is_loc_size(nbLevels + 1);
688 
689  for (int kk = 0; kk < nbLevels; kk++) {
690 
691  // get indices up to up to give approximation order
692  CHKERR createIsAtLevel(kk, &is_vec[kk]);
693  CHKERR ISGetSize(is_vec[kk], &is_glob_size[kk]);
694  CHKERR ISGetLocalSize(is_vec[kk], &is_loc_size[kk]);
695 
696  if (verb > QUIET) {
697  PetscSynchronizedPrintf(comm,
698  "Nb. dofs at level [ %d ] global %u local %d\n",
699  kk, is_glob_size[kk], is_loc_size[kk]);
700  }
701 
702  // if no dofs on level kk finish here
703  if (is_glob_size[kk] == 0) {
704  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "no dofs at level");
705  }
706  }
707 
708  for (int kk = 0; kk != nbLevels; kk++) {
709  Mat subA;
710  if (kk == nbLevels - 1 && use_mat_a) {
711  subA = A;
713  false, false);
714  } else {
715  if (kk > 0) {
716  // Not coarse level
718  true, shellSubA);
719  } else {
720  // Coarse lave is compressed matrix allowing for factorization when
721  // needed
723  true, false);
724  }
725  if (subA) {
726  CHKERR MatDestroy(&subA);
727  }
728  }
729  }
730 
731  for (unsigned int kk = 0; kk < is_vec.size(); kk++) {
732  CHKERR destroyIsAtLevel(kk, &is_vec[kk]);
733  }
734 
735  if (verb > QUIET) {
736  PetscSynchronizedFlush(comm, PETSC_STDOUT);
737  }
738 
740 }
741 
743  int verb) {
744 
745  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
747 
748  MPI_Comm comm;
749  CHKERR PetscObjectGetComm((PetscObject)pc, &comm);
750  if (verb > 0) {
751  PetscPrintf(comm, "Start PCMGSetUpViaApproxOrders\n");
752  }
753 
754  CHKERR ctx->getOptions();
755  CHKERR ctx->buildProlongationOperator(true, verb);
756 
757 #if PETSC_VERSION_GE(3, 8, 0)
758  CHKERR PCMGSetGalerkin(pc, PC_MG_GALERKIN_NONE);
759 #else
760  CHKERR PCMGSetGalerkin(pc, PETSC_FALSE);
761 #endif
762 
763  CHKERR PCMGSetLevels(pc, ctx->nbLevels, NULL);
764 
765  if (verb > 0) {
766  PetscPrintf(comm, "End PCMGSetUpViaApproxOrders\n");
767  }
768 
770 }
MoFEMErrorCode sub_mat_mult_add(Mat a, Vec x, Vec f)
MoFEM interface unique ID.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Common.hpp:60
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 sections.
Definition: ISManager.hpp:33
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
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
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
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:526
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:490
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, what elements compose problem and what DOFs are on rows and columns.
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMMoFEM.cpp:907
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:889
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.
MoFEMErrorCode DMCreateMatrix_MGViaApproxOrders(DM dm, Mat *M)
Create matrix for Multi-Grid via approximation orders.
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:97
MoFEMErrorCode sub_mat_mult_generic(Mat a, Vec x, Vec f)
#define CHKERR
Inline error check.
Definition: definitions.h:578
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:347
PETSc Discrete Manager data structure.
Definition: DMMoFEM.hpp:737
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:403
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:359
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:175
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