v0.13.1
Loading...
Searching...
No Matches
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
10
11using 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
32PCMGSubMatrixCtx::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
70public:
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
87template <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
106MoFEMErrorCode sub_mat_mult(Mat a, Vec x, Vec f) {
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
114MoFEMErrorCode 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
163DMMGViaApproxOrdersCtx::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
404static 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
462MoFEMErrorCode 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}
static Index< 'M', 3 > M
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
@ QUIET
Definition: definitions.h:208
@ ROW
Definition: definitions.h:123
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:373
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:1144
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:361
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMMoFEM.cpp:1114
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:53
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:149
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
constexpr double omega
constexpr AssemblyType A
Definition: plastic.cpp:35
constexpr double g
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:892
Deprecated interface functions.
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
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)