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);
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);
87template <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);
140 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
145 for (
unsigned int ii = 0; ii <
coarseningIS.size(); ii++) {
149 for (
unsigned int ii = 0; ii <
kspOperators.size(); ii++) {
170#define GET_DM_FIELD(DM) \
172 static_cast<DMCtx *>(DM->data)->getInterface<DMMGViaApproxOrdersCtx>(); \
176 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
184 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
189 CHKERR AODestroy(&dm_field->aO);
193 CHKERR PetscObjectReference((PetscObject)ao);
199 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
202 *size = dm_field->coarseningIS.size();
208 bool create_sub_matrix,
210 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
213 dm_field->coarseningIS.push_back(is);
216 CHKERR PetscObjectReference((PetscObject)is);
221 CHKERR ISDuplicate(is, &is2);
223 CHKERR AOApplicationToPetscIS(dm_field->aO, is2);
225 if (create_sub_matrix) {
231 CHKERR PetscObjectGetComm((PetscObject)
A, &comm);
233 &(dm_field->shellMatrixCtxPtr.back()), subA);
234 CHKERR MatShellSetOperation(*subA, MATOP_MULT,
236 CHKERR MatShellSetOperation(*subA, MATOP_MULT_ADD,
238 CHKERR MatShellSetOperation(*subA, MATOP_SOR,
241#if PETSC_VERSION_GE(3, 8, 0)
242 CHKERR MatCreateSubMatrix(
A, is2, is2, MAT_INITIAL_MATRIX, subA);
244 CHKERR MatGetSubMatrix(
A, is2, is2, MAT_INITIAL_MATRIX, subA);
251 dm_field->kspOperators.push_back(*subA);
252 CHKERR PetscObjectReference((PetscObject)(*subA));
256 PetscInfo(dm,
"Push back IS to DMMGViaApproxOrders\n");
261 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
264 if (dm_field->coarseningIS.back()) {
265 CHKERR ISDestroy(&dm_field->coarseningIS.back());
266 dm_field->coarseningIS.pop_back();
268 if (dm_field->kspOperators.back()) {
269 CHKERR MatDestroy(&dm_field->kspOperators.back());
271 dm_field->kspOperators.pop_back();
272 PetscInfo(dm,
"Pop back IS to DMMGViaApproxOrders\n");
277 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
280 CHKERR dm_field->destroyCoarseningIS();
281 PetscInfo(dm,
"Clear DMs data structures\n");
288 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
291 int nb_no_changed = 0;
295 std::vector<IS>::iterator it;
296 it = dm_field->coarseningIS.begin();
298 for (; it != dm_field->coarseningIS.end(); it++, ii++) {
301 CHKERR ISEqual(*it, is_vec[ii], &flg);
304 CHKERR MatDestroy(&dm_field->kspOperators[ii]);
306 CHKERR PetscObjectReference((PetscObject)is_vec[ii]);
307 if (ii < nb_elems - 1) {
310 CHKERR ISDuplicate(is_vec[ii], &is);
311 CHKERR ISCopy(is_vec[ii], is);
312 CHKERR AOApplicationToPetscIS(dm_field->aO, is);
315#if PETSC_VERSION_GE(3, 8, 0)
316 CHKERR MatCreateSubMatrix(
A, is, is, MAT_INITIAL_MATRIX, &subA);
318 CHKERR MatGetSubMatrix(
A, is, is, MAT_INITIAL_MATRIX, &subA);
320 CHKERR PetscObjectReference((PetscObject)subA);
321 dm_field->kspOperators[ii] = subA;
327 CHKERR PetscObjectReference((PetscObject)
A);
328 dm_field->kspOperators[ii] =
A;
337 if (
static_cast<int>(dm_field->coarseningIS.size()) < nb_elems) {
338 for (; ii < nb_elems - 1; ii++) {
349 for (; ii < static_cast<int>(dm_field->coarseningIS.size()); ii++) {
355 CHKERR PetscObjectGetComm((PetscObject)dm, &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());
363 PetscInfo(dm,
"Replace IS to DMMGViaApproxOrders\n");
369 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
410 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
415 ((
DMCtx *)(dm->data))->referenceNumber++;
425 PetscInfo1(dm,
"Create DMMGViaApproxOrders reference = %d\n",
426 ((
DMCtx *)(dm->data))->referenceNumber);
432 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
436 int leveldown = dm->leveldown;
438 if (dm_field->kspOperators.empty()) {
442 CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
443 if (dm_field->kspOperators.empty()) {
445 "data inconsistency, operator can not be set");
447 if (
static_cast<int>(dm_field->kspOperators.size()) < leveldown) {
449 "data inconsistency, no IS for that level");
451 *
M = dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown];
452 CHKERR PetscObjectReference((PetscObject)*
M);
456 PetscInfo1(dm,
"Create Matrix DMMGViaApproxOrders leveldown = %d\n",
463 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
466 CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
467 CHKERR DMCreate(comm, dmc);
468 (*dmc)->data = dm->data;
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);
479 PetscValidHeaderSpecific(dm1, DM_CLASSID, 1);
480 PetscValidHeaderSpecific(dm2, DM_CLASSID, 1);
484 CHKERR PetscObjectGetComm((PetscObject)dm1, &comm);
491 int dm_down_leveldown = dm_down->leveldown;
492 int dm_up_leveldown = dm_up->leveldown;
495 "Create interpolation DMMGViaApproxOrders dm1_leveldown = %d "
496 "dm2_leveldown = %d\n",
497 dm_down_leveldown, dm_up_leveldown);
503 if (
static_cast<int>(dm_field->coarseningIS.size()) < dm_down_leveldown) {
506 is_down = dm_field->coarseningIS[dm_field->coarseningIS.size() - 1 -
509 CHKERR ISGetLocalSize(is_down, &
m);
514 if (
static_cast<int>(dm_field->coarseningIS.size()) < dm_up_leveldown) {
519 ->coarseningIS[dm_field->coarseningIS.size() - 1 - dm_up_leveldown];
521 CHKERR ISGetLocalSize(is_up, &
n);
527 CHKERR MatCreate(comm, mat);
529 CHKERR MatSetType(*mat, MATMPIAIJ);
530 CHKERR MatMPIAIJSetPreallocation(*mat, 1, PETSC_NULL, 0, PETSC_NULL);
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);
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);
549 map<int, int> idx_map;
550 for (
int ii = 0; ii <
m; ii++) {
551 idx_map[row_indices_ptr[ii]] = rstart + ii;
554 CHKERR MatZeroEntries(*mat);
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);
563 CHKERR ISRestoreIndices(is_down, &row_indices_ptr);
564 CHKERR ISRestoreIndices(is_up, &col_indices_ptr);
566 CHKERR MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
567 CHKERR MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);
578 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
580 int leveldown = dm->leveldown;
582 if (dm_field->kspOperators.empty()) {
585#if PETSC_VERSION_GE(3, 5, 3)
587 dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
591 dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
596 PetscInfo1(dm,
"Create global vector DMMGViaApproxOrders leveldown = %d\n",
603 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
604 "MOFEM Multi-Grid (Orders) pre-conditioner",
"none");
607 CHKERR PetscOptionsInt(
"-mofem_mg_levels",
"nb levels of multi-grid solver",
609 CHKERR PetscOptionsInt(
"-mofem_mg_coarse_order",
610 "approximation order of coarse level",
"", 2,
612 CHKERR PetscOptionsInt(
"-mofem_mg_order_at_last_level",
"order at last level",
614 CHKERR PetscOptionsInt(
"-mofem_mg_verbose",
"nb levels of multi-grid solver",
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);
621 ierr = PetscOptionsEnd();
641 ->getPetscGlobalDofIdx();
647 string problem_name = problem_ptr->
getName();
649 order_at_next_level, is);
666 CHKERR PetscObjectGetComm((PetscObject)
dM, &comm);
669 PetscPrintf(comm,
"set MG levels %u\n",
nbLevels);
672 std::vector<IS> is_vec(
nbLevels + 1);
675 for (
int kk = 0; kk <
nbLevels; kk++) {
679 CHKERR ISGetSize(is_vec[kk], &is_glob_size[kk]);
680 CHKERR ISGetLocalSize(is_vec[kk], &is_loc_size[kk]);
683 PetscSynchronizedPrintf(comm,
684 "Nb. dofs at level [ %d ] global %u local %d\n",
685 kk, is_glob_size[kk], is_loc_size[kk]);
689 if (is_glob_size[kk] == 0) {
694 for (
int kk = 0; kk !=
nbLevels; kk++) {
696 if (kk ==
nbLevels - 1 && use_mat_a) {
717 for (
unsigned int kk = 0; kk < is_vec.size(); kk++) {
722 PetscSynchronizedFlush(comm, PETSC_STDOUT);
731 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
735 CHKERR PetscObjectGetComm((PetscObject)pc, &comm);
737 PetscPrintf(comm,
"Start PCMGSetUpViaApproxOrders\n");
743#if PETSC_VERSION_GE(3, 8, 0)
744 CHKERR PCMGSetGalerkin(pc, PC_MG_GALERKIN_NONE);
746 CHKERR PCMGSetGalerkin(pc, PETSC_FALSE);
752 PetscPrintf(comm,
"End PCMGSetUpViaApproxOrders\n");
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)
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
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
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.
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)
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.
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
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.
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)
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Structure for DM for multi-grid via approximation orders.
MoFEMErrorCode destroyCoarseningIS()
std::vector< Mat > kspOperators
Get KSP operators.
boost::ptr_vector< PCMGSubMatrixCtx > shellMatrixCtxPtr
Shell sub-matrix context.
virtual ~DMMGViaApproxOrdersCtx()
std::vector< IS > coarseningIS
Coarsening IS.
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, MoFEM::UnknownInterface **iface) const
PETSc Discrete Manager data structure.
Deprecated interface functions.
Section manager is used to create indexes and sections.
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
Mat A
Matrix at fine level.
virtual MoFEMErrorCode buildProlongationOperator(bool use_mat_a, int verb=0)
Set up data structures for MG.
MoFEMErrorCode initData(Vec x)
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)
PCMGSubMatrixCtx_private(Mat a, IS is)
PetscLogEvent MOFEM_EVENT_sor
~PCMGSubMatrixCtx_private()
PetscLogEvent MOFEM_EVENT_mult
virtual ~PCMGSubMatrixCtx()
PCMGSubMatrixCtx(Mat a, IS is)