6 #undef PETSC_VERSION_RELEASE
7 #define PETSC_VERSION_RELEASE 1
9 #if PETSC_VERSION_GE(3, 6, 0)
10 #include <petsc/private/petscimpl.h>
12 #include <petsc-private/petscimpl.h>
15 #if PETSC_VERSION_GE(3, 6, 0)
16 #include <petsc/private/dmimpl.h>
19 #include <petsc-private/dmimpl.h>
20 #include <petsc-private/vecimpl.h>
34 template <InsertMode MODE>
37 MatSORType flag, PetscReal shift,
38 PetscInt its, PetscInt lits,
Vec x);
86 boost::ptr_vector<PCMGSubMatrixCtx>
143 template <InsertMode MODE>
147 CHKERR MatShellGetContext(
a, &void_ctx);
151 CHKERR VecScatterBegin(ctx->
sCat, x, ctx->
X, INSERT_VALUES, SCATTER_REVERSE);
152 CHKERR VecScatterEnd(ctx->
sCat, x, ctx->
X, INSERT_VALUES, SCATTER_REVERSE);
154 CHKERR VecScatterBegin(ctx->
sCat, ctx->
F,
f, MODE, SCATTER_FORWARD);
155 CHKERR VecScatterEnd(ctx->
sCat, ctx->
F,
f, MODE, SCATTER_FORWARD);
161 return sub_mat_mult_generic<INSERT_VALUES>(
a, x,
f);
165 return sub_mat_mult_generic<ADD_VALUES>(
a, x,
f);
169 PetscReal shift, PetscInt its, PetscInt lits,
176 CHKERR MatShellGetContext(mat, &void_ctx);
180 CHKERR VecScatterBegin(ctx->
sCat, b, ctx->
X, INSERT_VALUES, SCATTER_REVERSE);
181 CHKERR VecScatterEnd(ctx->
sCat, b, ctx->
X, INSERT_VALUES, SCATTER_REVERSE);
182 CHKERR MatSOR(ctx->
A, ctx->
X,
omega, flag, shift, its, lits, ctx->
F);
183 CHKERR VecScatterBegin(ctx->
sCat, ctx->
F, x, INSERT_VALUES, SCATTER_FORWARD);
184 CHKERR VecScatterEnd(ctx->
sCat, ctx->
F, x, INSERT_VALUES, SCATTER_FORWARD);
193 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
214 #define GET_DM_FIELD(DM) \
216 static_cast<DMCtx *>(DM->data)->getInterface<DMMGViaApproxOrdersCtx>(); \
220 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
228 bool create_sub_matrix,
230 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
240 CHKERR AOApplicationToPetscIS(dm_field->aO, is2);
242 if (create_sub_matrix) {
248 CHKERR PetscObjectGetComm((PetscObject)
A, &comm);
251 &(dm_field->shellMatrixCtxPtr.back()),
253 CHKERR MatShellSetOperation(sub_a_raw, MATOP_MULT,
255 CHKERR MatShellSetOperation(sub_a_raw, MATOP_MULT_ADD,
257 CHKERR MatShellSetOperation(sub_a_raw, MATOP_SOR,
259 dm_field->kspOperators.emplace_back(
263 #if PETSC_VERSION_GE(3, 8, 0)
264 CHKERR MatCreateSubMatrix(
A, is2, is2, MAT_INITIAL_MATRIX, &sub_a_raw);
266 CHKERR MatGetSubMatrix(
A, is2, is2, MAT_INITIAL_MATRIX, &sub_a_raw);
268 dm_field->kspOperators.emplace_back(
277 PetscInfo(dm,
"Push back IS to DMMGViaApproxOrders\n");
282 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
285 if (dm_field->coarseningIS.back()) {
286 dm_field->coarseningIS.pop_back();
288 dm_field->kspOperators.pop_back();
289 PetscInfo(dm,
"Pop back IS to DMMGViaApproxOrders\n");
294 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
297 CHKERR dm_field->destroyCoarseningIS();
298 PetscInfo(dm,
"Clear DMs data structures\n");
314 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
319 ((
DMCtxImpl *)(dm->data))->incrementReference();
330 PetscInfo1(dm,
"Create DMMGViaApproxOrders reference = %d\n",
336 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
339 CHKERR dm_field->destroyCoarseningIS();
346 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
350 int leveldown = dm->leveldown;
352 if (dm_field->kspOperators.empty()) {
356 CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
357 if (dm_field->kspOperators.empty()) {
359 "data inconsistency, operator can not be set");
361 if (
static_cast<int>(dm_field->kspOperators.size()) < leveldown) {
363 "data inconsistency, no IS for that level");
365 *
M = dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown];
366 CHKERR PetscObjectReference((PetscObject)*
M);
370 PetscInfo1(dm,
"Create Matrix DMMGViaApproxOrders leveldown = %d\n",
377 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
380 CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
381 CHKERR DMCreate(comm, dmc);
382 (*dmc)->data = dm->data;
387 PetscInfo1(dm,
"Coarsen DMMGViaApproxOrders leveldown = %d\n", dm->leveldown);
393 PetscValidHeaderSpecific(dm1, DM_CLASSID, 1);
394 PetscValidHeaderSpecific(dm2, DM_CLASSID, 1);
398 CHKERR PetscObjectGetComm((PetscObject)dm1, &comm);
405 int dm_down_leveldown = dm_down->leveldown;
406 int dm_up_leveldown = dm_up->leveldown;
409 "Create interpolation DMMGViaApproxOrders dm1_leveldown = %d "
410 "dm2_leveldown = %d\n",
411 dm_down_leveldown, dm_up_leveldown);
417 if (
static_cast<int>(dm_field->coarseningIS.size()) < dm_down_leveldown) {
420 is_down = dm_field->coarseningIS[dm_field->coarseningIS.size() - 1 -
423 CHKERR ISGetLocalSize(is_down, &
m);
428 if (
static_cast<int>(dm_field->coarseningIS.size()) < dm_up_leveldown) {
433 ->coarseningIS[dm_field->coarseningIS.size() - 1 - dm_up_leveldown];
435 CHKERR ISGetLocalSize(is_up, &
n);
441 CHKERR MatCreate(comm, mat);
443 CHKERR MatSetType(*mat, MATMPIAIJ);
444 CHKERR MatMPIAIJSetPreallocation(*mat, 1, PETSC_NULL, 0, PETSC_NULL);
447 PetscLayout rmap, cmap;
448 CHKERR MatGetLayouts(*mat, &rmap, &cmap);
449 int rstart, rend, cstart, cend;
450 CHKERR PetscLayoutGetRange(rmap, &rstart, &rend);
451 CHKERR PetscLayoutGetRange(cmap, &cstart, &cend);
453 const int *row_indices_ptr, *col_indices_ptr;
454 CHKERR ISGetIndices(is_down, &row_indices_ptr);
455 CHKERR ISGetIndices(is_up, &col_indices_ptr);
457 map<int, int> idx_map;
458 for (
int ii = 0; ii <
m; ii++) {
459 idx_map[row_indices_ptr[ii]] = rstart + ii;
462 CHKERR MatZeroEntries(*mat);
464 for (
int jj = 0; jj <
n; jj++) {
465 map<int, int>::iterator mit = idx_map.find(col_indices_ptr[jj]);
466 if (mit != idx_map.end()) {
467 CHKERR MatSetValue(*mat, mit->second, cstart + jj, 1, INSERT_VALUES);
471 CHKERR ISRestoreIndices(is_down, &row_indices_ptr);
472 CHKERR ISRestoreIndices(is_up, &col_indices_ptr);
474 CHKERR MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
475 CHKERR MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);
485 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
487 int leveldown = dm->leveldown;
489 if (dm_field->kspOperators.empty()) {
492 #if PETSC_VERSION_GE(3, 5, 3)
494 dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
498 dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
503 PetscInfo1(dm,
"Create global vector DMMGViaApproxOrders leveldown = %d\n",
510 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
511 "MOFEM Multi-Grid (Orders) pre-conditioner",
"none");
514 CHKERR PetscOptionsInt(
"-mofem_mg_levels",
"nb levels of multi-grid solver",
516 CHKERR PetscOptionsInt(
"-mofem_mg_coarse_order",
517 "approximation order of coarse level",
"", 2,
519 CHKERR PetscOptionsInt(
"-mofem_mg_order_at_last_level",
"order at last level",
521 CHKERR PetscOptionsInt(
"-mofem_mg_verbose",
"nb levels of multi-grid solver",
523 PetscBool shell_sub_a =
shellSubA ? PETSC_TRUE : PETSC_FALSE;
524 CHKERR PetscOptionsBool(
"-mofem_mg_shell_a",
"use shell matrix as sub matrix",
525 "", shell_sub_a, &shell_sub_a, NULL);
528 ierr = PetscOptionsEnd();
548 ->getPetscGlobalDofIdx();
554 string problem_name = problem_ptr->
getName();
556 order_at_next_level, is);
573 CHKERR PetscObjectGetComm((PetscObject)
dM, &comm);
575 MOFEM_LOG_C(
"PCMGViaApproximationOrders", Sev::inform,
"set MG levels %u\n",
581 std::vector<IS> is_vec(
nbLevels + 1);
584 for (
int kk = 0; kk <
nbLevels; kk++) {
588 CHKERR ISGetSize(is_vec[kk], &is_glob_size[kk]);
589 CHKERR ISGetLocalSize(is_vec[kk], &is_loc_size[kk]);
592 "Nb. dofs at level [ %d ] global %u local %d\n", kk,
593 is_glob_size[kk], is_loc_size[kk]);
596 if (is_glob_size[kk] == 0) {
601 for (
int kk = 0; kk !=
nbLevels; kk++) {
602 if (kk ==
nbLevels - 1 && use_mat_a) {
619 for (
unsigned int kk = 0; kk < is_vec.size(); kk++) {
629 boost::shared_ptr<PCMGSetUpViaApproxOrdersCtx>
631 return boost::make_shared<PCMGSetUpViaApproxOrdersCtx>(dm,
A, use_shell_mat);
635 PC pc, boost::shared_ptr<PCMGSetUpViaApproxOrdersCtx> ctx,
int verb) {
637 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
640 auto core_log = logging::core::get();
642 "PCMGViaApproximationOrders"));
644 MOFEM_LOG_TAG(
"PCMGViaApproximationOrders",
"PCMGViaApproximationOrders");
646 MOFEM_LOG(
"PCMGViaApproximationOrders", Sev::verbose)
647 <<
"Setup PCMGSetUpViaApproxOrders";
650 CHKERR ctx->buildProlongationOperator(
true, verb);
652 #if PETSC_VERSION_GE(3, 8, 0)
653 CHKERR PCMGSetGalerkin(pc, PC_MG_GALERKIN_NONE);
655 CHKERR PCMGSetGalerkin(pc, PETSC_FALSE);
658 CHKERR PCMGSetLevels(pc, ctx->nbLevels, NULL);
660 MOFEM_LOG(
"PCMGViaApproximationOrders", Sev::noisy)
661 <<
"Setup PCMGSetUpViaApproxOrders <- end";