75 PC_MG *mg = (PC_MG *)pc->data;
77 PC_MG_Levels **mglevels = mg->levels;
78 PCMGType mgtype = mg->am;
79 PetscInt mgctype = (PetscInt)PC_MG_CYCLE_V;
87 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
88 PetscValidLogicalCollectiveInt(pc, levels, 2);
89 CHKERR PetscObjectGetComm((PetscObject)pc, &comm);
90 vector<bool> set_ksp(levels,
false);
92 mgctype = mglevels[0]->cycles;
95 n = mglevels[0]->levels;
98 for (
i = 1;
i <
n - 1;
i++) {
102 for (
i = 0;
i <
n;
i++) {
105 if (mglevels[
i]->smoothd != mglevels[
i]->smoothu) {
106 CHKERR KSPDestroy(&mglevels[
i]->smoothd);
108 CHKERR KSPDestroy(&mglevels[
i]->smoothu);
113 mg->nlevels = levels;
115 CHKERR PetscMalloc1(levels, &mglevels);
116 CHKERR PetscLogObjectMemory((PetscObject)pc, levels * (
sizeof(PC_MG *)));
118 CHKERR PCGetOptionsPrefix(pc, &prefix);
121 for (
i = 0;
i < levels;
i++) {
122 CHKERR PetscNewLog(pc, &mglevels[
i]);
124 mglevels[
i]->levels = levels;
128 if (
i != mg->levels[
i]->level) {
131 mglevels[
i]->level = mg->levels[
i]->level;
132 mglevels[
i]->cycles = mg->levels[
i]->cycles;
133 mglevels[
i]->eventsmoothsetup = mg->levels[
i]->eventsmoothsetup;
134 mglevels[
i]->eventsmoothsolve = mg->levels[
i]->eventsmoothsolve;
135 mglevels[
i]->eventresidual = mg->levels[
i]->eventresidual;
136 mglevels[
i]->eventinterprestrict = mg->levels[
i]->eventinterprestrict;
137 mglevels[
i]->restrct = mg->levels[
i]->restrct;
138 mglevels[
i]->interpolate = mg->levels[
i]->interpolate;
139 mglevels[
i]->rscale = mg->levels[
i]->rscale;
140 mglevels[
i]->r = mg->levels[
i]->r;
141 mglevels[
i]->b = mg->levels[
i]->b;
142 mglevels[
i]->x = mg->levels[
i]->x;
144 mglevels[
i]->level =
i;
145 mglevels[
i]->cycles = mgctype;
146 mg->default_smoothu = 2;
147 mg->default_smoothd = 2;
148 mglevels[
i]->eventsmoothsetup = 0;
149 mglevels[
i]->eventsmoothsolve = 0;
150 mglevels[
i]->eventresidual = 0;
151 mglevels[
i]->eventinterprestrict = 0;
156 mglevels[
i]->smoothd = mg->levels[
i]->smoothd;
157 mglevels[
i]->smoothu = mg->levels[
i]->smoothu;
164 CHKERR KSPCreate(comm, &mglevels[
i]->smoothd);
166 KSPSetErrorIfNotConverged(mglevels[
i]->smoothd, pc->erroriffailure);
167 CHKERR KSPSetType(mglevels[
i]->smoothd, KSPCHEBYSHEV);
168 CHKERR KSPSetConvergenceTest(mglevels[
i]->smoothd, KSPConvergedSkip, NULL,
170 CHKERR KSPSetNormType(mglevels[
i]->smoothd, KSP_NORM_NONE);
171 CHKERR KSPGetPC(mglevels[
i]->smoothd, &ipc);
172 CHKERR PCSetType(ipc, PCSOR);
173 CHKERR PetscObjectIncrementTabLevel((PetscObject)mglevels[
i]->smoothd,
174 (PetscObject)pc, levels -
i);
176 KSPSetTolerances(mglevels[
i]->smoothd, PETSC_DEFAULT, PETSC_DEFAULT,
177 PETSC_DEFAULT,
i ? mg->default_smoothd : 1);
178 CHKERR KSPSetOptionsPrefix(mglevels[
i]->smoothd, prefix);
181 if (!
i && levels > 1) {
182 CHKERR KSPAppendOptionsPrefix(mglevels[0]->smoothd,
"mg_coarse_");
186 CHKERR KSPSetType(mglevels[0]->smoothd, KSPPREONLY);
187 CHKERR KSPGetPC(mglevels[0]->smoothd, &ipc);
188 CHKERR MPI_Comm_size(comm, &size);
192 CHKERR PCSetType(ipc, PCREDUNDANT);
193 CHKERR PCRedundantGetKSP(ipc, &innerksp);
194 CHKERR KSPGetPC(innerksp, &innerpc);
195 CHKERR PCFactorSetShiftType(innerpc, MAT_SHIFT_INBLOCKS);
197 CHKERR PCSetType(ipc, PCLU);
198 CHKERR PCFactorSetShiftType(ipc, MAT_SHIFT_INBLOCKS);
203 sprintf(tprefix,
"mg_levels_%d_", (
int)
i);
204 CHKERR KSPAppendOptionsPrefix(mglevels[
i]->smoothd, tprefix);
206 CHKERR PetscLogObjectParent((PetscObject)pc,
207 (PetscObject)mglevels[
i]->smoothd);
209 mglevels[
i]->smoothu = mglevels[
i]->smoothd;
214 mg->cyclesperpcapply = 1;
218 CHKERR PetscFree(mg->levels);
220 mg->levels = mglevels;
221 CHKERR PCMGSetType(pc, mgtype);
222 PetscFunctionReturn(0);
228 PC_MG *mg = (PC_MG *)pc->data;
229 PC_MG_Levels **mglevels = mg->levels;
231 PetscInt
i,
n = mglevels[0]->levels;
233 PetscBool preonly, lu, redundant, cholesky, svd,
234 dump = PETSC_FALSE, opsset, use_amat, missinginterpolate = PETSC_FALSE;
237 PetscViewer viewer = 0;
241 if (mg->usedmfornumberoflevels) {
243 CHKERR DMGetRefineLevel(pc->dm, &levels);
246 CHKERR PCMGSetLevels(pc, levels, NULL);
248 CHKERR PCSetFromOptions(pc);
250 mglevels = mg->levels;
253 CHKERR KSPGetPC(mglevels[0]->smoothd, &cpc);
259 CHKERR KSPGetOperatorsSet(mglevels[
n - 1]->smoothd, NULL, &opsset);
262 CHKERR KSPGetOperators(mglevels[
n - 1]->smoothd, NULL, &mmat);
263 if (mmat == pc->pmat)
264 opsset = PETSC_FALSE;
268 CHKERR PCGetUseAmat(pc, &use_amat);
272 "Using outer operators to define finest grid operator \n because "
273 "PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was "
275 CHKERR KSPSetOperators(mglevels[
n - 1]->smoothd, pc->mat, pc->pmat);
277 CHKERR PetscInfo(pc,
"Using matrix (pmat) operators to define finest "
278 "grid operator \n because "
279 "PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators("
280 "ksp,...); was not called.\n");
281 CHKERR KSPSetOperators(mglevels[
n - 1]->smoothd, pc->pmat, pc->pmat);
285 for (
i =
n - 1;
i > 0;
i--) {
286 if (!(mglevels[
i]->interpolate || mglevels[
i]->restrct)) {
287 missinginterpolate = PETSC_TRUE;
292 vector<bool> set_dms_d(
n,
false);
293 vector<bool> set_dms_u(
n,
false);
294 for (
int i = 0;
i <
n;
i++) {
295 if (mglevels[
i]->smoothd->dm)
297 if (mglevels[
i]->smoothu->dm)
307 if (missinginterpolate && pc->dm && mg->galerkin != 2 && !pc->setupcalled) {
315 for (
i =
n - 2;
i > -1;
i--) {
317 CHKERR DMCoarsen(dms[
i + 1], MPI_COMM_NULL, &dms[
i]);
318 CHKERR KSPSetDM(mglevels[
i]->smoothd, dms[
i]);
320 CHKERR KSPSetDMActive(mglevels[
i]->smoothd, PETSC_FALSE);
323 CHKERR DMGetDMKSPWrite(dms[
i], &kdm);
327 kdm->ops->computerhs = NULL;
330 dms[
i] = mglevels[
i]->smoothd->dm;
331 dms[
i]->leveldown = dms[
i + 1]->leveldown + 1;
332 CHKERR PetscObjectReference((PetscObject)dms[
i]);
335 for (
i =
n - 2;
i > -1;
i--) {
336 if (!mglevels[
i + 1]->interpolate) {
338 CHKERR DMCreateInterpolation(dms[
i], dms[
i + 1], &
p, &rscale);
339 CHKERR PCMGSetInterpolation(pc,
i + 1,
p);
341 CHKERR PCMGSetRScale(pc,
i + 1, rscale);
343 CHKERR VecDestroy(&rscale);
349 for (
i =
n - 2;
i > -1;
i--) {
355 if (pc->dm && !pc->setupcalled) {
358 CHKERR KSPSetDM(mglevels[
n - 1]->smoothd, pc->dm);
359 CHKERR KSPSetDMActive(mglevels[
n - 1]->smoothd, PETSC_FALSE);
362 if (mg->galerkin == 1) {
364 CHKERR PetscObjectGetComm((PetscObject)pc, &comm);
388 if (!mg->galerkin && pc->dm) {
389 for (
i =
n - 2;
i >= 0;
i--) {
390 if (!mglevels[
i + 1]->restrct) {
393 Mat Restrict, Inject;
395 CHKERR KSPGetDM(mglevels[
i + 1]->smoothd, &dmfine);
396 CHKERR KSPGetDM(mglevels[
i]->smoothd, &dmcoarse);
397 CHKERR PCMGGetRestriction(pc,
i + 1, &Restrict);
398 CHKERR PCMGGetRScale(pc,
i + 1, &rscale);
400 CHKERR DMRestrict(dmfine, Restrict, rscale, Inject, dmcoarse);
405 if (!pc->setupcalled) {
406 for (
i = 0;
i <
n;
i++) {
409 CHKERR KSPSetFromOptions(mglevels[
i]->smoothd);
411 for (
i = 1;
i <
n;
i++) {
414 if (mglevels[
i]->smoothu &&
415 (mglevels[
i]->smoothu != mglevels[
i]->smoothd)) {
416 CHKERR KSPSetFromOptions(mglevels[
i]->smoothu);
419 for (
i = 1;
i <
n;
i++) {
420 CHKERR PCMGGetInterpolation(pc,
i, &mglevels[
i]->interpolate);
421 CHKERR PCMGGetRestriction(pc,
i, &mglevels[
i]->restrct);
423 for (
i = 0;
i <
n - 1;
i++) {
424 if (!mglevels[
i]->b) {
426 CHKERR KSPCreateVecs(mglevels[
i]->smoothd, 1, &vec, 0, NULL);
427 CHKERR PCMGSetRhs(pc,
i, *vec);
431 if (!mglevels[
i]->r &&
i) {
432 CHKERR VecDuplicate(mglevels[
i]->b, &tvec);
436 if (!mglevels[
i]->x) {
437 CHKERR VecDuplicate(mglevels[
i]->b, &tvec);
442 if (
n != 1 && !mglevels[
n - 1]->r) {
445 CHKERR KSPCreateVecs(mglevels[
n - 1]->smoothd, 1, &vec, 0, NULL);
446 CHKERR PCMGSetR(pc,
n - 1, *vec);
455 for (
i = 0;
i <
n - 1;
i++) {
458 if (mglevels[
i]->smoothd->setupstage != KSP_SETUP_NEW)
459 mglevels[
i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
463 for (
i = 1;
i <
n;
i++) {
464 if (mglevels[
i]->smoothu == mglevels[
i]->smoothd || mg->am == PC_MG_FULL ||
465 mg->am == PC_MG_KASKADE || mg->cyclesperpcapply > 1) {
467 CHKERR KSPSetInitialGuessNonzero(mglevels[
i]->smoothd, PETSC_TRUE);
469 if (mglevels[
i]->eventsmoothsetup) {
470 CHKERR PetscLogEventBegin(mglevels[
i]->eventsmoothsetup, 0, 0, 0, 0);
472 CHKERR KSPSetUp(mglevels[
i]->smoothd);
473 if (mglevels[
i]->eventsmoothsetup) {
474 CHKERR PetscLogEventEnd(mglevels[
i]->eventsmoothsetup, 0, 0, 0, 0);
476 if (!mglevels[
i]->residual) {
478 CHKERR KSPGetOperators(mglevels[
i]->smoothd, NULL, &mat);
479 CHKERR PCMGSetResidual(pc,
i, PCMGResidualDefault, mat);
482 for (
i = 1;
i <
n;
i++) {
483 if (mglevels[
i]->smoothu && mglevels[
i]->smoothu != mglevels[
i]->smoothd) {
484 Mat downmat, downpmat;
489 CHKERR KSPGetOperatorsSet(mglevels[
i]->smoothu, &opsset, NULL);
491 CHKERR KSPGetOperators(mglevels[
i]->smoothd, &downmat, &downpmat);
492 CHKERR KSPSetOperators(mglevels[
i]->smoothu, downmat, downpmat);
494 CHKERR KSPSetInitialGuessNonzero(mglevels[
i]->smoothu, PETSC_TRUE);
495 if (mglevels[
i]->eventsmoothsetup) {
496 CHKERR PetscLogEventBegin(mglevels[
i]->eventsmoothsetup, 0, 0, 0, 0);
498 CHKERR KSPSetUp(mglevels[
i]->smoothu);
499 if (mglevels[
i]->eventsmoothsetup) {
500 CHKERR PetscLogEventEnd(mglevels[
i]->eventsmoothsetup, 0, 0, 0, 0);
509 CHKERR PetscObjectTypeCompare((PetscObject)mglevels[0]->smoothd, KSPPREONLY,
512 CHKERR PetscObjectTypeCompare((PetscObject)cpc, PCLU, &lu);
513 CHKERR PetscObjectTypeCompare((PetscObject)cpc, PCREDUNDANT, &redundant);
514 CHKERR PetscObjectTypeCompare((PetscObject)cpc, PCCHOLESKY, &cholesky);
515 CHKERR PetscObjectTypeCompare((PetscObject)cpc, PCSVD, &svd);
516 if (!lu && !redundant && !cholesky && !svd) {
517 CHKERR KSPSetType(mglevels[0]->smoothd, KSPGMRES);
521 if (mglevels[0]->eventsmoothsetup) {
522 CHKERR PetscLogEventBegin(mglevels[0]->eventsmoothsetup, 0, 0, 0, 0);
524 CHKERR KSPSetUp(mglevels[0]->smoothd);
525 if (mglevels[0]->eventsmoothsetup) {
526 CHKERR PetscLogEventEnd(mglevels[0]->eventsmoothsetup, 0, 0, 0, 0);
549 for (
i = 1;
i <
n;
i++) {
550 CHKERR MatView(mglevels[
i]->restrct, viewer);
552 for (
i = 0;
i <
n;
i++) {
553 CHKERR KSPGetPC(mglevels[
i]->smoothd, &pc);
554 CHKERR MatView(pc->mat, viewer);
557 PetscFunctionReturn(0);