11 using namespace MoFEM;
13 #include <PCMGSetUpViaApproxOrders.hpp>
15 #undef PETSC_VERSION_RELEASE
16 #define PETSC_VERSION_RELEASE 1
18 #if PETSC_VERSION_GE(3, 6, 0)
19 #include <petsc/private/petscimpl.h>
21 #include <petsc-private/petscimpl.h>
24 #if PETSC_VERSION_GE(3, 6, 0)
25 #include <petsc/private/dmimpl.h>
28 #include <petsc-private/dmimpl.h>
29 #include <petsc-private/vecimpl.h>
35 ierr = PetscObjectReference((PetscObject)
A);
36 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
37 ierr = PetscObjectReference((PetscObject)iS);
38 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
42 ierr = MatDestroy(&
A);
43 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
45 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
51 PetscLogEventRegister(
"PCMGSubMatrixCtx_mult", 0, &MOFEM_EVENT_mult);
52 PetscLogEventRegister(
"PCMGSubMatrixCtx_sor", 0, &MOFEM_EVENT_sor);
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);
64 template <InsertMode MODE>
67 MatSORType flag, PetscReal shift,
68 PetscInt its, PetscInt lits,
Vec x);
73 if (!isInitisalised) {
75 CHKERR VecScatterCreate(X, iS, x, PETSC_NULL, &sCat);
78 isInitisalised =
true;
87 template <InsertMode MODE>
91 CHKERR MatShellGetContext(
a, &void_ctx);
97 CHKERR VecScatterBegin(ctx->
sCat, x, ctx->
X, INSERT_VALUES, SCATTER_REVERSE);
98 CHKERR VecScatterEnd(ctx->
sCat, x, ctx->
X, INSERT_VALUES, SCATTER_REVERSE);
100 CHKERR VecScatterBegin(ctx->
sCat, ctx->
F,
f, MODE, SCATTER_FORWARD);
101 CHKERR VecScatterEnd(ctx->
sCat, ctx->
F,
f, MODE, SCATTER_FORWARD);
107 return sub_mat_mult_generic<INSERT_VALUES>(
a, x,
f);
111 return sub_mat_mult_generic<ADD_VALUES>(
a, x,
f);
115 PetscReal shift, PetscInt its, PetscInt lits,
119 CHKERR MatShellGetContext(mat, &void_ctx);
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);
139 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
144 for (
unsigned int ii = 0; ii <
coarseningIS.size(); ii++) {
148 for (
unsigned int ii = 0; ii <
kspOperators.size(); ii++) {
169 #define GET_DM_FIELD(DM) \
171 static_cast<DMCtx *>(DM->data)->getInterface<DMMGViaApproxOrdersCtx>(); \
175 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
183 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
188 CHKERR AODestroy(&dm_field->aO);
192 CHKERR PetscObjectReference((PetscObject)ao);
198 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
201 *size = dm_field->coarseningIS.size();
207 bool create_sub_matrix,
209 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
212 dm_field->coarseningIS.push_back(is);
215 CHKERR PetscObjectReference((PetscObject)is);
220 CHKERR ISDuplicate(is, &is2);
222 CHKERR AOApplicationToPetscIS(dm_field->aO, is2);
224 if (create_sub_matrix) {
230 CHKERR PetscObjectGetComm((PetscObject)
A, &comm);
232 &(dm_field->shellMatrixCtxPtr.back()), subA);
233 CHKERR MatShellSetOperation(*subA, MATOP_MULT,
235 CHKERR MatShellSetOperation(*subA, MATOP_MULT_ADD,
237 CHKERR MatShellSetOperation(*subA, MATOP_SOR,
240 #if PETSC_VERSION_GE(3, 8, 0)
241 CHKERR MatCreateSubMatrix(
A, is2, is2, MAT_INITIAL_MATRIX, subA);
243 CHKERR MatGetSubMatrix(
A, is2, is2, MAT_INITIAL_MATRIX, subA);
250 dm_field->kspOperators.push_back(*subA);
251 CHKERR PetscObjectReference((PetscObject)(*subA));
255 PetscInfo(dm,
"Push back IS to DMMGViaApproxOrders\n");
260 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
263 if (dm_field->coarseningIS.back()) {
264 CHKERR ISDestroy(&dm_field->coarseningIS.back());
265 dm_field->coarseningIS.pop_back();
267 if (dm_field->kspOperators.back()) {
268 CHKERR MatDestroy(&dm_field->kspOperators.back());
270 dm_field->kspOperators.pop_back();
271 PetscInfo(dm,
"Pop back IS to DMMGViaApproxOrders\n");
276 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
279 CHKERR dm_field->destroyCoarseningIS();
280 PetscInfo(dm,
"Clear DMs data structures\n");
287 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
290 int nb_no_changed = 0;
294 std::vector<IS>::iterator it;
295 it = dm_field->coarseningIS.begin();
297 for (; it != dm_field->coarseningIS.end(); it++, ii++) {
300 CHKERR ISEqual(*it, is_vec[ii], &flg);
303 CHKERR MatDestroy(&dm_field->kspOperators[ii]);
305 CHKERR PetscObjectReference((PetscObject)is_vec[ii]);
306 if (ii < nb_elems - 1) {
309 CHKERR ISDuplicate(is_vec[ii], &is);
310 CHKERR ISCopy(is_vec[ii], is);
311 CHKERR AOApplicationToPetscIS(dm_field->aO, is);
314 #if PETSC_VERSION_GE(3, 8, 0)
315 CHKERR MatCreateSubMatrix(
A, is, is, MAT_INITIAL_MATRIX, &subA);
317 CHKERR MatGetSubMatrix(
A, is, is, MAT_INITIAL_MATRIX, &subA);
319 CHKERR PetscObjectReference((PetscObject)subA);
320 dm_field->kspOperators[ii] = subA;
326 CHKERR PetscObjectReference((PetscObject)
A);
327 dm_field->kspOperators[ii] =
A;
336 if (
static_cast<int>(dm_field->coarseningIS.size()) < nb_elems) {
337 for (; ii < nb_elems - 1; ii++) {
348 for (; ii < static_cast<int>(dm_field->coarseningIS.size()); ii++) {
354 CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
357 "DMMGViaApproxOrders nb_no_changed = %d, nb_replaced = %d, "
358 "nb_added = %d, nb_deleted = %d, size = %d\n",
359 nb_no_changed, nb_replaced, nb_added, nb_deleted,
360 dm_field->coarseningIS.size());
362 PetscInfo(dm,
"Replace IS to DMMGViaApproxOrders\n");
368 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
409 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
414 ((
DMCtxImpl *)(dm->data))->incrementReference();
424 PetscInfo1(dm,
"Create DMMGViaApproxOrders reference = %d\n",
431 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
435 int leveldown = dm->leveldown;
437 if (dm_field->kspOperators.empty()) {
441 CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
442 if (dm_field->kspOperators.empty()) {
444 "data inconsistency, operator can not be set");
446 if (
static_cast<int>(dm_field->kspOperators.size()) < leveldown) {
448 "data inconsistency, no IS for that level");
450 *
M = dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown];
451 CHKERR PetscObjectReference((PetscObject)*
M);
455 PetscInfo1(dm,
"Create Matrix DMMGViaApproxOrders leveldown = %d\n",
462 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
465 CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
466 CHKERR DMCreate(comm, dmc);
467 (*dmc)->data = dm->data;
471 CHKERR PetscObjectReference((PetscObject)(*dmc));
472 PetscInfo1(dm,
"Coarsen DMMGViaApproxOrders leveldown = %d\n", dm->leveldown);
478 PetscValidHeaderSpecific(dm1, DM_CLASSID, 1);
479 PetscValidHeaderSpecific(dm2, DM_CLASSID, 1);
483 CHKERR PetscObjectGetComm((PetscObject)dm1, &comm);
490 int dm_down_leveldown = dm_down->leveldown;
491 int dm_up_leveldown = dm_up->leveldown;
494 "Create interpolation DMMGViaApproxOrders dm1_leveldown = %d "
495 "dm2_leveldown = %d\n",
496 dm_down_leveldown, dm_up_leveldown);
502 if (
static_cast<int>(dm_field->coarseningIS.size()) < dm_down_leveldown) {
505 is_down = dm_field->coarseningIS[dm_field->coarseningIS.size() - 1 -
508 CHKERR ISGetLocalSize(is_down, &
m);
513 if (
static_cast<int>(dm_field->coarseningIS.size()) < dm_up_leveldown) {
518 ->coarseningIS[dm_field->coarseningIS.size() - 1 - dm_up_leveldown];
520 CHKERR ISGetLocalSize(is_up, &
n);
526 CHKERR MatCreate(comm, mat);
528 CHKERR MatSetType(*mat, MATMPIAIJ);
529 CHKERR MatMPIAIJSetPreallocation(*mat, 1, PETSC_NULL, 0, PETSC_NULL);
532 PetscLayout rmap, cmap;
533 CHKERR MatGetLayouts(*mat, &rmap, &cmap);
534 int rstart, rend, cstart, cend;
535 CHKERR PetscLayoutGetRange(rmap, &rstart, &rend);
536 CHKERR PetscLayoutGetRange(cmap, &cstart, &cend);
544 const int *row_indices_ptr, *col_indices_ptr;
545 CHKERR ISGetIndices(is_down, &row_indices_ptr);
546 CHKERR ISGetIndices(is_up, &col_indices_ptr);
548 map<int, int> idx_map;
549 for (
int ii = 0; ii <
m; ii++) {
550 idx_map[row_indices_ptr[ii]] = rstart + ii;
553 CHKERR MatZeroEntries(*mat);
555 for (
int jj = 0; jj <
n; jj++) {
556 map<int, int>::iterator mit = idx_map.find(col_indices_ptr[jj]);
557 if (mit != idx_map.end()) {
558 CHKERR MatSetValue(*mat, mit->second, cstart + jj, 1, INSERT_VALUES);
562 CHKERR ISRestoreIndices(is_down, &row_indices_ptr);
563 CHKERR ISRestoreIndices(is_up, &col_indices_ptr);
565 CHKERR MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
566 CHKERR MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);
577 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
579 int leveldown = dm->leveldown;
581 if (dm_field->kspOperators.empty()) {
584 #if PETSC_VERSION_GE(3, 5, 3)
586 dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
590 dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
595 PetscInfo1(dm,
"Create global vector DMMGViaApproxOrders leveldown = %d\n",
602 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
603 "MOFEM Multi-Grid (Orders) pre-conditioner",
"none");
606 CHKERR PetscOptionsInt(
"-mofem_mg_levels",
"nb levels of multi-grid solver",
608 CHKERR PetscOptionsInt(
"-mofem_mg_coarse_order",
609 "approximation order of coarse level",
"", 2,
611 CHKERR PetscOptionsInt(
"-mofem_mg_order_at_last_level",
"order at last level",
613 CHKERR PetscOptionsInt(
"-mofem_mg_verbose",
"nb levels of multi-grid solver",
615 PetscBool shell_sub_a =
shellSubA ? PETSC_TRUE : PETSC_FALSE;
616 CHKERR PetscOptionsBool(
"-mofem_mg_shell_a",
"use shell matrix as sub matrix",
617 "", shell_sub_a, &shell_sub_a, NULL);
620 ierr = PetscOptionsEnd();
640 ->getPetscGlobalDofIdx();
646 string problem_name = problem_ptr->
getName();
648 order_at_next_level, is);
665 CHKERR PetscObjectGetComm((PetscObject)
dM, &comm);
668 PetscPrintf(comm,
"set MG levels %u\n",
nbLevels);
671 std::vector<IS> is_vec(
nbLevels + 1);
674 for (
int kk = 0; kk <
nbLevels; kk++) {
678 CHKERR ISGetSize(is_vec[kk], &is_glob_size[kk]);
679 CHKERR ISGetLocalSize(is_vec[kk], &is_loc_size[kk]);
682 PetscSynchronizedPrintf(comm,
683 "Nb. dofs at level [ %d ] global %u local %d\n",
684 kk, is_glob_size[kk], is_loc_size[kk]);
688 if (is_glob_size[kk] == 0) {
693 for (
int kk = 0; kk !=
nbLevels; kk++) {
695 if (kk ==
nbLevels - 1 && use_mat_a) {
716 for (
unsigned int kk = 0; kk < is_vec.size(); kk++) {
721 PetscSynchronizedFlush(comm, PETSC_STDOUT);
730 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
734 CHKERR PetscObjectGetComm((PetscObject)pc, &comm);
736 PetscPrintf(comm,
"Start PCMGSetUpViaApproxOrders\n");
742 #if PETSC_VERSION_GE(3, 8, 0)
743 CHKERR PCMGSetGalerkin(pc, PC_MG_GALERKIN_NONE);
745 CHKERR PCMGSetGalerkin(pc, PETSC_FALSE);
751 PetscPrintf(comm,
"End PCMGSetUpViaApproxOrders\n");