Actual source code: mg.c
1: /*
2: Defines the multigrid preconditioner interface.
3: */
4: #include <petsc/private/pcmgimpl.h>
5: #include <petsc/private/kspimpl.h>
6: #include <petscdm.h>
7: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC, PetscBool *);
9: /*
10: Contains the list of registered coarse space construction routines
11: */
12: PetscFunctionList PCMGCoarseList = NULL;
14: /*
15: A down smoother that is distinct from the up smoother starts from a zero initial guess in the forward V cycle. It is given
16: a nonzero initial guess where a cycle continues from an existing solution, which KSPPREONLY cannot do
17: */
18: static PetscErrorCode PCMGCheckSmootherDownGuess_Private(PC pc, KSP smoothd, const char cycle[])
19: {
20: PC spc;
21: PetscBool ispreonly, allowed;
23: PetscFunctionBegin;
24: PetscCall(PetscObjectTypeCompareAny((PetscObject)smoothd, &ispreonly, KSPPREONLY, KSPNONE, ""));
25: if (!ispreonly) PetscFunctionReturn(PETSC_SUCCESS);
26: PetscCall(KSPGetPC(smoothd, &spc));
27: PetscCall(PetscObjectTypeCompareAny((PetscObject)spc, &allowed, PCREDISTRIBUTE, PCMPI, ""));
28: PetscCheck(allowed, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "%s with distinct up and down smoothers needs a down smoother that accepts a nonzero initial guess, which KSPPREONLY does not; use one iteration of KSPRICHARDSON instead", cycle);
29: PetscFunctionReturn(PETSC_SUCCESS);
30: }
32: PetscErrorCode PCMGMCycle_Private(PC pc, PC_MG_Levels **mglevelsin, PetscBool transpose, PetscBool matapp, PCRichardsonConvergedReason *reason)
33: {
34: PC_MG *mg = (PC_MG *)pc->data;
35: PC_MG_Levels *mgc, *mglevels = *mglevelsin;
36: PetscInt cycles = (mglevels->level == 1) ? 1 : mglevels->cycles;
38: PetscFunctionBegin;
39: if (mglevels->eventsmoothsolve) PetscCall(PetscLogEventBegin(mglevels->eventsmoothsolve, 0, 0, 0, 0));
40: if (!transpose) {
41: if (matapp) {
42: PetscCall(KSPMatSolve(mglevels->smoothd, mglevels->B, mglevels->X)); /* pre-smooth */
43: PetscCall(KSPCheckMatSolve(mglevels->smoothd, pc, mglevels->X));
44: } else {
45: PetscCall(KSPSolve(mglevels->smoothd, mglevels->b, mglevels->x)); /* pre-smooth */
46: PetscCall(KSPCheckSolve(mglevels->smoothd, pc, mglevels->x));
47: }
48: } else {
49: PetscCheck(!matapp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not supported");
50: PetscCall(KSPSolveTranspose(mglevels->smoothu, mglevels->b, mglevels->x)); /* transpose of post-smooth */
51: PetscCall(KSPCheckSolve(mglevels->smoothu, pc, mglevels->x));
52: }
53: if (mglevels->eventsmoothsolve) PetscCall(PetscLogEventEnd(mglevels->eventsmoothsolve, 0, 0, 0, 0));
54: if (mglevels->level) { /* not the coarsest grid */
55: if (mglevels->eventresidual) PetscCall(PetscLogEventBegin(mglevels->eventresidual, 0, 0, 0, 0));
56: if (matapp && !mglevels->R) PetscCall(MatDuplicate(mglevels->B, MAT_DO_NOT_COPY_VALUES, &mglevels->R));
57: if (!transpose) {
58: if (matapp) PetscCall((*mglevels->matresidual)(mglevels->A, mglevels->B, mglevels->X, mglevels->R));
59: else PetscCall((*mglevels->residual)(mglevels->A, mglevels->b, mglevels->x, mglevels->r));
60: } else {
61: if (matapp) PetscCall((*mglevels->matresidualtranspose)(mglevels->A, mglevels->B, mglevels->X, mglevels->R));
62: else PetscCall((*mglevels->residualtranspose)(mglevels->A, mglevels->b, mglevels->x, mglevels->r));
63: }
64: if (mglevels->eventresidual) PetscCall(PetscLogEventEnd(mglevels->eventresidual, 0, 0, 0, 0));
66: /* if on finest level and have convergence criteria set */
67: if (mglevels->level == mglevels->levels - 1 && mg->ttol && reason) {
68: PetscReal rnorm;
70: PetscCall(VecNorm(mglevels->r, NORM_2, &rnorm));
71: if (rnorm <= mg->ttol) {
72: if (rnorm < mg->abstol) {
73: *reason = PCRICHARDSON_CONVERGED_ATOL;
74: PetscCall(PetscInfo(pc, "Linear solver has converged. Residual norm %g is less than absolute tolerance %g\n", (double)rnorm, (double)mg->abstol));
75: } else {
76: *reason = PCRICHARDSON_CONVERGED_RTOL;
77: PetscCall(PetscInfo(pc, "Linear solver has converged. Residual norm %g is less than relative tolerance times initial residual norm %g\n", (double)rnorm, (double)mg->ttol));
78: }
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
81: }
83: mgc = *(mglevelsin - 1);
84: if (mglevels->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels->eventinterprestrict, 0, 0, 0, 0));
85: if (!transpose) {
86: if (matapp) PetscCall(MatMatRestrict(mglevels->restrct, mglevels->R, &mgc->B));
87: else PetscCall(MatRestrict(mglevels->restrct, mglevels->r, mgc->b));
88: } else {
89: if (matapp) PetscCall(MatMatRestrict(mglevels->interpolate, mglevels->R, &mgc->B));
90: else PetscCall(MatRestrict(mglevels->interpolate, mglevels->r, mgc->b));
91: }
92: if (mglevels->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels->eventinterprestrict, 0, 0, 0, 0));
93: if (matapp) {
94: if (!mgc->X) PetscCall(MatDuplicate(mgc->B, MAT_DO_NOT_COPY_VALUES, &mgc->X));
95: else {
96: PetscCall(MatZeroEntries(mgc->X));
97: }
98: } else {
99: PetscCall(VecZeroEntries(mgc->x));
100: }
101: while (cycles--) PetscCall(PCMGMCycle_Private(pc, mglevelsin - 1, transpose, matapp, reason));
102: if (mglevels->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels->eventinterprestrict, 0, 0, 0, 0));
103: if (!transpose) {
104: if (matapp) PetscCall(MatMatInterpolateAdd(mglevels->interpolate, mgc->X, mglevels->X, &mglevels->X));
105: else PetscCall(MatInterpolateAdd(mglevels->interpolate, mgc->x, mglevels->x, mglevels->x));
106: } else {
107: PetscCall(MatInterpolateAdd(mglevels->restrct, mgc->x, mglevels->x, mglevels->x));
108: }
109: if (mglevels->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels->eventinterprestrict, 0, 0, 0, 0));
110: if (mglevels->eventsmoothsolve) PetscCall(PetscLogEventBegin(mglevels->eventsmoothsolve, 0, 0, 0, 0));
111: if (!transpose) {
112: if (matapp) {
113: PetscCall(KSPMatSolve(mglevels->smoothu, mglevels->B, mglevels->X)); /* post smooth */
114: PetscCall(KSPCheckMatSolve(mglevels->smoothu, pc, mglevels->X));
115: } else {
116: PetscCall(KSPSolve(mglevels->smoothu, mglevels->b, mglevels->x)); /* post smooth */
117: PetscCall(KSPCheckSolve(mglevels->smoothu, pc, mglevels->x));
118: }
119: } else {
120: PetscBool guessnonzero;
122: PetscCheck(!matapp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not supported");
123: /* smoothd may use a zero initial guess as the forward pre-smoother, but its transpose is applied last and must start from the interpolated coarse correction */
124: PetscCall(KSPGetInitialGuessNonzero(mglevels->smoothd, &guessnonzero));
125: if (!guessnonzero) PetscCall(PCMGCheckSmootherDownGuess_Private(pc, mglevels->smoothd, "PCApplyTranspose()"));
126: PetscCall(KSPSetInitialGuessNonzero(mglevels->smoothd, PETSC_TRUE));
127: PetscCall(KSPSolveTranspose(mglevels->smoothd, mglevels->b, mglevels->x)); /* transpose of pre-smooth */
128: PetscCall(KSPCheckSolve(mglevels->smoothd, pc, mglevels->x));
129: PetscCall(KSPSetInitialGuessNonzero(mglevels->smoothd, guessnonzero));
130: }
131: if (mglevels->cr) {
132: Mat crA;
134: PetscCheck(!matapp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not supported");
135: /* TODO Turn on copy and turn off noisy if we have an exact solution
136: PetscCall(VecCopy(mglevels->x, mglevels->crx));
137: PetscCall(VecCopy(mglevels->b, mglevels->crb)); */
138: PetscCall(KSPGetOperators(mglevels->cr, &crA, NULL));
139: PetscCall(KSPSetNoisy_Private(crA, mglevels->crx));
140: PetscCall(KSPSolve(mglevels->cr, mglevels->crb, mglevels->crx)); /* compatible relaxation */
141: PetscCall(KSPCheckSolve(mglevels->cr, pc, mglevels->crx));
142: }
143: if (mglevels->eventsmoothsolve) PetscCall(PetscLogEventEnd(mglevels->eventsmoothsolve, 0, 0, 0, 0));
144: }
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: static PetscErrorCode PCApplyRichardson_MG(PC pc, Vec b, Vec x, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool zeroguess, PetscInt *outits, PCRichardsonConvergedReason *reason)
149: {
150: PC_MG *mg = (PC_MG *)pc->data;
151: PC_MG_Levels **mglevels = mg->levels;
152: PC tpc;
153: PetscBool changeu, changed;
154: PetscInt levels = mglevels[0]->levels, i;
156: PetscFunctionBegin;
157: /* When the DM is supplying the matrix then it will not exist until here */
158: for (i = 0; i < levels; i++) {
159: if (!mglevels[i]->A) {
160: PetscCall(KSPGetOperators(mglevels[i]->smoothu, &mglevels[i]->A, NULL));
161: PetscCall(PetscObjectReference((PetscObject)mglevels[i]->A));
162: }
163: }
165: PetscCall(KSPGetPC(mglevels[levels - 1]->smoothd, &tpc));
166: PetscCall(PCPreSolveChangeRHS(tpc, &changed));
167: PetscCall(KSPGetPC(mglevels[levels - 1]->smoothu, &tpc));
168: PetscCall(PCPreSolveChangeRHS(tpc, &changeu));
169: if (!changed && !changeu) {
170: PetscCall(VecDestroy(&mglevels[levels - 1]->b));
171: mglevels[levels - 1]->b = b;
172: } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */
173: if (!mglevels[levels - 1]->b) {
174: Vec *vec;
176: PetscCall(KSPCreateVecs(mglevels[levels - 1]->smoothd, 1, &vec, 0, NULL));
177: mglevels[levels - 1]->b = *vec;
178: PetscCall(PetscFree(vec));
179: }
180: PetscCall(VecCopy(b, mglevels[levels - 1]->b));
181: }
182: mglevels[levels - 1]->x = x;
184: mg->rtol = rtol;
185: mg->abstol = abstol;
186: mg->dtol = dtol;
187: if (rtol) {
188: /* compute initial residual norm for relative convergence test */
189: PetscReal rnorm;
191: if (zeroguess) {
192: PetscCall(VecNorm(b, NORM_2, &rnorm));
193: } else {
194: PetscCall((*mglevels[levels - 1]->residual)(mglevels[levels - 1]->A, b, x, w));
195: PetscCall(VecNorm(w, NORM_2, &rnorm));
196: }
197: mg->ttol = PetscMax(rtol * rnorm, abstol);
198: } else if (abstol) mg->ttol = abstol;
199: else mg->ttol = 0.0;
201: /* since smoother is applied to full system, not just residual we need to make sure that smoothers don't
202: stop prematurely due to small residual */
203: for (i = 1; i < levels; i++) {
204: PetscCall(KSPSetTolerances(mglevels[i]->smoothu, 0, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
205: if (mglevels[i]->smoothu != mglevels[i]->smoothd) {
206: /* For Richardson the initial guess is nonzero since it is solving in each cycle the original system not just applying as a preconditioner */
207: PetscCall(KSPSetInitialGuessNonzero(mglevels[i]->smoothd, PETSC_TRUE));
208: PetscCall(KSPSetTolerances(mglevels[i]->smoothd, 0, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
209: }
210: }
212: *reason = PCRICHARDSON_NOT_SET;
213: for (i = 0; i < its; i++) {
214: PetscCall(PCMGMCycle_Private(pc, mglevels + levels - 1, PETSC_FALSE, PETSC_FALSE, reason));
215: if (*reason) break;
216: }
217: if (*reason == PCRICHARDSON_NOT_SET) *reason = PCRICHARDSON_CONVERGED_ITS;
218: *outits = i;
219: if (!changed && !changeu) mglevels[levels - 1]->b = NULL;
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: PetscErrorCode PCReset_MG(PC pc)
224: {
225: PC_MG *mg = (PC_MG *)pc->data;
226: PC_MG_Levels **mglevels = mg->levels;
227: PetscInt i, n;
229: PetscFunctionBegin;
230: if (mglevels) {
231: n = mglevels[0]->levels;
232: for (i = 0; i < n - 1; i++) {
233: PetscCall(VecDestroy(&mglevels[i + 1]->r));
234: PetscCall(VecDestroy(&mglevels[i]->b));
235: PetscCall(VecDestroy(&mglevels[i]->x));
236: PetscCall(MatDestroy(&mglevels[i + 1]->R));
237: PetscCall(MatDestroy(&mglevels[i]->B));
238: PetscCall(MatDestroy(&mglevels[i]->X));
239: PetscCall(VecDestroy(&mglevels[i]->crx));
240: PetscCall(VecDestroy(&mglevels[i]->crb));
241: PetscCall(MatDestroy(&mglevels[i + 1]->restrct));
242: PetscCall(MatDestroy(&mglevels[i + 1]->interpolate));
243: PetscCall(MatDestroy(&mglevels[i + 1]->inject));
244: PetscCall(VecDestroy(&mglevels[i + 1]->rscale));
245: }
246: PetscCall(VecDestroy(&mglevels[n - 1]->crx));
247: PetscCall(VecDestroy(&mglevels[n - 1]->crb));
248: /* this is not null only if the smoother on the finest level
249: changes the rhs during PreSolve */
250: PetscCall(VecDestroy(&mglevels[n - 1]->b));
251: PetscCall(MatDestroy(&mglevels[n - 1]->B));
253: for (i = 0; i < n; i++) {
254: PetscCall(MatDestroy(&mglevels[i]->coarseSpace));
255: PetscCall(MatDestroy(&mglevels[i]->A));
256: if (mglevels[i]->smoothd != mglevels[i]->smoothu) PetscCall(KSPReset(mglevels[i]->smoothd));
257: PetscCall(KSPReset(mglevels[i]->smoothu));
258: if (mglevels[i]->cr) PetscCall(KSPReset(mglevels[i]->cr));
259: }
260: mg->Nc = 0;
261: }
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: /* Implementing CR
267: We only want to make corrections that ``do not change'' the coarse solution. What we mean by not changing is that if I prolong my coarse solution to the fine grid and then inject that fine solution back to the coarse grid, I get the same answer. Injection is what Brannick calls R. We want the complementary projector to Inj, which we will call S, after Brannick, so that Inj S = 0. Now the orthogonal projector onto the range of Inj^T is
269: Inj^T (Inj Inj^T)^{-1} Inj
271: and if Inj is a VecScatter, as it is now in PETSc, we have
273: Inj^T Inj
275: and
277: S = I - Inj^T Inj
279: since
281: Inj S = Inj - (Inj Inj^T) Inj = 0.
283: Brannick suggests
285: A \to S^T A S \qquad\mathrm{and}\qquad M \to S^T M S
287: but I do not think his :math:`S^T S = I` is correct. Our S is an orthogonal projector, so :math:`S^T S = S^2 = S`. We will use
289: M^{-1} A \to S M^{-1} A S
291: In fact, since it is somewhat hard in PETSc to do the symmetric application, we will just apply S on the left.
293: Check: || Inj P - I ||_F < tol
294: Check: In general, Inj Inj^T = I
295: */
297: typedef struct {
298: PC mg; /* The PCMG object */
299: PetscInt l; /* The multigrid level for this solver */
300: Mat Inj; /* The injection matrix */
301: Mat S; /* I - Inj^T Inj */
302: } CRContext;
304: static PetscErrorCode CRSetup_Private(PC pc)
305: {
306: CRContext *ctx;
307: Mat It;
309: PetscFunctionBeginUser;
310: PetscCall(PCShellGetContext(pc, &ctx));
311: PetscCall(PCMGGetInjection(ctx->mg, ctx->l, &It));
312: PetscCheck(It, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "CR requires that injection be defined for this PCMG");
313: PetscCall(MatCreateTranspose(It, &ctx->Inj));
314: PetscCall(MatCreateNormal(ctx->Inj, &ctx->S));
315: PetscCall(MatScale(ctx->S, -1.0));
316: PetscCall(MatShift(ctx->S, 1.0));
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: static PetscErrorCode CRApply_Private(PC pc, Vec x, Vec y)
321: {
322: CRContext *ctx;
324: PetscFunctionBeginUser;
325: PetscCall(PCShellGetContext(pc, &ctx));
326: PetscCall(MatMult(ctx->S, x, y));
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: static PetscErrorCode CRDestroy_Private(PC pc)
331: {
332: CRContext *ctx;
334: PetscFunctionBeginUser;
335: PetscCall(PCShellGetContext(pc, &ctx));
336: PetscCall(MatDestroy(&ctx->Inj));
337: PetscCall(MatDestroy(&ctx->S));
338: PetscCall(PetscFree(ctx));
339: PetscCall(PCShellSetContext(pc, NULL));
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: static PetscErrorCode CreateCR_Private(PC pc, PetscInt l, PC *cr)
344: {
345: CRContext *ctx;
347: PetscFunctionBeginUser;
348: PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), cr));
349: PetscCall(PetscObjectSetName((PetscObject)*cr, "S (complementary projector to injection)"));
350: PetscCall(PetscCalloc1(1, &ctx));
351: ctx->mg = pc;
352: ctx->l = l;
353: PetscCall(PCSetType(*cr, PCSHELL));
354: PetscCall(PCShellSetContext(*cr, ctx));
355: PetscCall(PCShellSetApply(*cr, CRApply_Private));
356: PetscCall(PCShellSetSetUp(*cr, CRSetup_Private));
357: PetscCall(PCShellSetDestroy(*cr, CRDestroy_Private));
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions, const char[], const char[], const char *[], const char *[], PetscBool *);
363: PetscErrorCode PCMGSetLevels_MG(PC pc, PetscInt levels, MPI_Comm *comms)
364: {
365: PC_MG *mg = (PC_MG *)pc->data;
366: MPI_Comm comm;
367: PC_MG_Levels **mglevels = mg->levels;
368: PCMGType mgtype = mg->am;
369: PetscInt mgctype = (PetscInt)PC_MG_CYCLE_V;
370: PetscInt i;
371: PetscMPIInt size;
372: const char *prefix;
373: PC ipc;
374: PetscInt n;
376: PetscFunctionBegin;
379: if (mg->nlevels == levels) PetscFunctionReturn(PETSC_SUCCESS);
380: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
381: if (mglevels) {
382: mgctype = mglevels[0]->cycles;
383: /* changing the number of levels so free up the previous stuff */
384: PetscCall(PCReset_MG(pc));
385: n = mglevels[0]->levels;
386: for (i = 0; i < n; i++) {
387: if (mglevels[i]->smoothd != mglevels[i]->smoothu) PetscCall(KSPDestroy(&mglevels[i]->smoothd));
388: PetscCall(KSPDestroy(&mglevels[i]->smoothu));
389: PetscCall(KSPDestroy(&mglevels[i]->cr));
390: PetscCall(PetscFree(mglevels[i]));
391: }
392: PetscCall(PetscFree(mg->levels));
393: }
395: mg->nlevels = levels;
397: PetscCall(PetscMalloc1(levels, &mglevels));
399: PetscCall(PCGetOptionsPrefix(pc, &prefix));
401: mg->stageApply = 0;
402: for (i = 0; i < levels; i++) {
403: PetscCall(PetscNew(&mglevels[i]));
405: mglevels[i]->level = i;
406: mglevels[i]->levels = levels;
407: mglevels[i]->cycles = mgctype;
408: mg->default_smoothu = 2;
409: mg->default_smoothd = 2;
410: mglevels[i]->eventsmoothsetup = 0;
411: mglevels[i]->eventsmoothsolve = 0;
412: mglevels[i]->eventresidual = 0;
413: mglevels[i]->eventinterprestrict = 0;
415: if (comms) comm = comms[i];
416: if (comm != MPI_COMM_NULL) {
417: PetscCall(KSPCreate(comm, &mglevels[i]->smoothd));
418: PetscCall(KSPSetNestLevel(mglevels[i]->smoothd, pc->kspnestlevel));
419: PetscCall(KSPSetErrorIfNotConverged(mglevels[i]->smoothd, pc->erroriffailure));
420: PetscCall(PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd, (PetscObject)pc, levels - i));
421: PetscCall(KSPSetOptionsPrefix(mglevels[i]->smoothd, prefix));
422: PetscCall(PetscObjectComposedDataSetInt((PetscObject)mglevels[i]->smoothd, PetscMGLevelId, mglevels[i]->level));
423: if (i == 0 && levels > 1) { // coarse grid
424: PetscCall(KSPAppendOptionsPrefix(mglevels[0]->smoothd, "mg_coarse_"));
426: /* coarse solve is (redundant) LU by default; set shifttype NONZERO to avoid annoying zero-pivot in LU preconditioner */
427: PetscCall(KSPSetType(mglevels[0]->smoothd, KSPPREONLY));
428: PetscCall(KSPGetPC(mglevels[0]->smoothd, &ipc));
429: PetscCallMPI(MPI_Comm_size(comm, &size));
430: if (size > 1) PetscCall(PCSetType(ipc, PCREDUNDANT));
431: else PetscCall(PCSetType(ipc, PCLU));
432: PetscCall(PCFactorSetShiftType(ipc, MAT_SHIFT_INBLOCKS));
433: } else {
434: char tprefix[128];
436: PetscCall(KSPSetType(mglevels[i]->smoothd, KSPCHEBYSHEV));
437: PetscCall(KSPSetConvergenceTest(mglevels[i]->smoothd, KSPConvergedSkip, NULL, NULL));
438: PetscCall(KSPSetNormType(mglevels[i]->smoothd, KSP_NORM_NONE));
439: PetscCall(KSPGetPC(mglevels[i]->smoothd, &ipc));
440: PetscCall(PCSetType(ipc, PCSOR));
441: PetscCall(KSPSetTolerances(mglevels[i]->smoothd, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, mg->default_smoothd));
443: if (i == levels - 1 && levels > 1) { // replace 'mg_finegrid_' with 'mg_levels_X_'
444: PetscBool set;
446: PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)mglevels[i]->smoothd)->options, ((PetscObject)mglevels[i]->smoothd)->prefix, "-mg_fine_", NULL, NULL, &set));
447: if (set) {
448: if (prefix) PetscCall(PetscSNPrintf(tprefix, 128, "%smg_fine_", prefix));
449: else PetscCall(PetscSNPrintf(tprefix, 128, "mg_fine_"));
450: PetscCall(KSPSetOptionsPrefix(mglevels[i]->smoothd, tprefix));
451: } else {
452: PetscCall(PetscSNPrintf(tprefix, 128, "mg_levels_%" PetscInt_FMT "_", i));
453: PetscCall(KSPAppendOptionsPrefix(mglevels[i]->smoothd, tprefix));
454: }
455: } else {
456: PetscCall(PetscSNPrintf(tprefix, 128, "mg_levels_%" PetscInt_FMT "_", i));
457: PetscCall(KSPAppendOptionsPrefix(mglevels[i]->smoothd, tprefix));
458: }
459: }
460: }
461: mglevels[i]->smoothu = mglevels[i]->smoothd;
462: mg->rtol = 0.0;
463: mg->abstol = 0.0;
464: mg->dtol = 0.0;
465: mg->ttol = 0.0;
466: mg->cyclesperpcapply = 1;
467: }
468: mg->levels = mglevels;
469: PetscCall(PCMGSetType(pc, mgtype));
470: PetscFunctionReturn(PETSC_SUCCESS);
471: }
473: /*@
474: PCMGSetLevels - Sets the number of levels to use with `PCMG`.
475: Must be called before any other `PCMG` routine.
477: Logically Collective
479: Input Parameters:
480: + pc - the preconditioner context
481: . levels - the number of levels
482: - comms - optional communicators for each level; this is to allow solving the coarser problems
483: on smaller sets of processes. For processes that are not included in the computation
484: you must pass `MPI_COMM_NULL`. Use comms = `NULL` to specify that all processes
485: should participate in each level of problem.
487: Options Database Key:
488: . -pc_mg_levels levels - set the number of levels to use
490: Level: intermediate
492: Notes:
493: If the number of levels is one then the multigrid uses the `-mg_levels` prefix
494: for setting the level options rather than the `-mg_coarse` or `-mg_fine` prefix.
496: You can free the information in `comms` after this routine is called.
498: The array of MPI communicators must contain `MPI_COMM_NULL` for those processes that at each level
499: are not participating in the coarser solve. For example, with 2 levels and 1 and 2 processes on
500: the two levels, rank 0 in the original communicator will pass in an array of 2 communicators
501: of size 2 and 1, while rank 1 in the original communicator will pass in array of 2 communicators
502: the first of size 2 and the second of value `MPI_COMM_NULL` since the rank 1 does not participate
503: in the coarse grid solve.
505: Since each coarser level may have a new `MPI_Comm` with fewer processes than the previous, one
506: must take special care in providing the restriction and interpolation operation. We recommend
507: providing these as two step operations; first perform a standard restriction or interpolation on
508: the full number of processes for that level and then use an MPI call to copy the resulting vector
509: array entries (after calls to `VecGetArray()`) to the smaller or larger number of processes, note in both
510: cases the MPI calls must be made on the larger of the two communicators. Traditional MPI send and
511: receives or `MPI_AlltoAllv()` could be used to do the reshuffling of the vector entries.
513: Fortran Notes:
514: Use comms = `PETSC_NULL_MPI_COMM` as the equivalent of `NULL` in the C interface. Note `PETSC_NULL_MPI_COMM`
515: is not `MPI_COMM_NULL`. It is more like `PETSC_NULL_INTEGER`, `PETSC_NULL_REAL` etc.
517: .seealso: [](ch_ksp), `PCMGSetType()`, `PCMGGetLevels()`
518: @*/
519: PetscErrorCode PCMGSetLevels(PC pc, PetscInt levels, MPI_Comm *comms)
520: {
521: PetscFunctionBegin;
523: if (comms) PetscAssertPointer(comms, 3);
524: PetscTryMethod(pc, "PCMGSetLevels_C", (PC, PetscInt, MPI_Comm *), (pc, levels, comms));
525: PetscFunctionReturn(PETSC_SUCCESS);
526: }
528: PetscErrorCode PCDestroy_MG(PC pc)
529: {
530: PC_MG *mg = (PC_MG *)pc->data;
531: PC_MG_Levels **mglevels = mg->levels;
532: PetscInt i, n;
534: PetscFunctionBegin;
535: PetscCall(PCReset_MG(pc));
536: if (mglevels) {
537: n = mglevels[0]->levels;
538: for (i = 0; i < n; i++) {
539: if (mglevels[i]->smoothd != mglevels[i]->smoothu) PetscCall(KSPDestroy(&mglevels[i]->smoothd));
540: PetscCall(KSPDestroy(&mglevels[i]->smoothu));
541: PetscCall(KSPDestroy(&mglevels[i]->cr));
542: PetscCall(PetscFree(mglevels[i]));
543: }
544: PetscCall(PetscFree(mg->levels));
545: }
546: PetscCall(PetscFree(pc->data));
547: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", NULL));
548: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", NULL));
549: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetGalerkin_C", NULL));
550: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetReusePreconditioner_C", NULL));
551: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetLevels_C", NULL));
552: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetLevels_C", NULL));
553: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", NULL));
554: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", NULL));
555: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptInterpolation_C", NULL));
556: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptInterpolation_C", NULL));
557: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptCR_C", NULL));
558: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptCR_C", NULL));
559: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptCoarseSpaceType_C", NULL));
560: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptCoarseSpaceType_C", NULL));
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }
564: /*
565: PCApply_MG - Runs either an additive, multiplicative, Kaskadic
566: or full cycle of multigrid.
568: Note:
569: A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle().
570: */
571: static PetscErrorCode PCApply_MG_Internal(PC pc, Vec b, Vec x, Mat B, Mat X, PetscBool transpose)
572: {
573: PC_MG *mg = (PC_MG *)pc->data;
574: PC_MG_Levels **mglevels = mg->levels;
575: PC tpc;
576: PetscInt levels = mglevels[0]->levels, i;
577: PetscBool changeu, changed, matapp;
579: PetscFunctionBegin;
580: matapp = (PetscBool)(B && X);
581: if (mg->stageApply) PetscCall(PetscLogStagePush(mg->stageApply));
582: /* When the DM is supplying the matrix then it will not exist until here */
583: for (i = 0; i < levels; i++) {
584: if (!mglevels[i]->A) {
585: PetscCall(KSPGetOperators(mglevels[i]->smoothu, &mglevels[i]->A, NULL));
586: PetscCall(PetscObjectReference((PetscObject)mglevels[i]->A));
587: }
588: }
590: PetscCall(KSPGetPC(mglevels[levels - 1]->smoothd, &tpc));
591: PetscCall(PCPreSolveChangeRHS(tpc, &changed));
592: PetscCall(KSPGetPC(mglevels[levels - 1]->smoothu, &tpc));
593: PetscCall(PCPreSolveChangeRHS(tpc, &changeu));
594: if (!changeu && !changed) {
595: if (matapp) {
596: PetscCall(MatDestroy(&mglevels[levels - 1]->B));
597: mglevels[levels - 1]->B = B;
598: } else {
599: PetscCall(VecDestroy(&mglevels[levels - 1]->b));
600: mglevels[levels - 1]->b = b;
601: }
602: } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */
603: if (matapp) {
604: if (mglevels[levels - 1]->B) {
605: PetscInt N1, N2;
606: PetscBool flg;
608: PetscCall(MatGetSize(mglevels[levels - 1]->B, NULL, &N1));
609: PetscCall(MatGetSize(B, NULL, &N2));
610: PetscCall(PetscObjectTypeCompare((PetscObject)mglevels[levels - 1]->B, ((PetscObject)B)->type_name, &flg));
611: if (N1 != N2 || !flg) PetscCall(MatDestroy(&mglevels[levels - 1]->B));
612: }
613: if (!mglevels[levels - 1]->B) PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &mglevels[levels - 1]->B));
614: else PetscCall(MatCopy(B, mglevels[levels - 1]->B, SAME_NONZERO_PATTERN));
615: } else {
616: if (!mglevels[levels - 1]->b) {
617: Vec *vec;
619: PetscCall(KSPCreateVecs(mglevels[levels - 1]->smoothd, 1, &vec, 0, NULL));
620: mglevels[levels - 1]->b = *vec;
621: PetscCall(PetscFree(vec));
622: }
623: PetscCall(VecCopy(b, mglevels[levels - 1]->b));
624: }
625: }
626: if (matapp) {
627: mglevels[levels - 1]->X = X;
628: } else {
629: mglevels[levels - 1]->x = x;
630: }
632: /* If coarser Xs are present, it means we have already block applied the PC at least once
633: Reset operators if sizes/type do no match */
634: if (matapp && levels > 1 && mglevels[levels - 2]->X) {
635: PetscInt Xc, Bc;
636: PetscBool flg;
638: PetscCall(MatGetSize(mglevels[levels - 2]->X, NULL, &Xc));
639: PetscCall(MatGetSize(mglevels[levels - 1]->B, NULL, &Bc));
640: PetscCall(PetscObjectTypeCompare((PetscObject)mglevels[levels - 2]->X, ((PetscObject)mglevels[levels - 1]->X)->type_name, &flg));
641: if (Xc != Bc || !flg) {
642: PetscCall(MatDestroy(&mglevels[levels - 1]->R));
643: for (i = 0; i < levels - 1; i++) {
644: PetscCall(MatDestroy(&mglevels[i]->R));
645: PetscCall(MatDestroy(&mglevels[i]->B));
646: PetscCall(MatDestroy(&mglevels[i]->X));
647: }
648: }
649: }
651: if (mg->am == PC_MG_MULTIPLICATIVE) {
652: if (matapp) PetscCall(MatZeroEntries(X));
653: else PetscCall(VecZeroEntries(x));
654: for (i = 0; i < mg->cyclesperpcapply; i++) PetscCall(PCMGMCycle_Private(pc, mglevels + levels - 1, transpose, matapp, NULL));
655: } else if (mg->am == PC_MG_ADDITIVE) {
656: PetscCall(PCMGACycle_Private(pc, mglevels, transpose, matapp));
657: } else if (mg->am == PC_MG_KASKADE) {
658: PetscCall(PCMGKCycle_Private(pc, mglevels, transpose, matapp));
659: } else {
660: PetscCall(PCMGFCycle_Private(pc, mglevels, transpose, matapp));
661: }
662: if (mg->stageApply) PetscCall(PetscLogStagePop());
663: if (!changeu && !changed) {
664: if (matapp) {
665: mglevels[levels - 1]->B = NULL;
666: } else {
667: mglevels[levels - 1]->b = NULL;
668: }
669: }
670: PetscFunctionReturn(PETSC_SUCCESS);
671: }
673: static PetscErrorCode PCApply_MG(PC pc, Vec b, Vec x)
674: {
675: PetscFunctionBegin;
676: PetscCall(PCApply_MG_Internal(pc, b, x, NULL, NULL, PETSC_FALSE));
677: PetscFunctionReturn(PETSC_SUCCESS);
678: }
680: static PetscErrorCode PCApplyTranspose_MG(PC pc, Vec b, Vec x)
681: {
682: PetscFunctionBegin;
683: PetscCall(PCApply_MG_Internal(pc, b, x, NULL, NULL, PETSC_TRUE));
684: PetscFunctionReturn(PETSC_SUCCESS);
685: }
687: static PetscErrorCode PCMatApply_MG(PC pc, Mat b, Mat x)
688: {
689: PetscFunctionBegin;
690: PetscCall(PCApply_MG_Internal(pc, NULL, NULL, b, x, PETSC_FALSE));
691: PetscFunctionReturn(PETSC_SUCCESS);
692: }
694: static PetscErrorCode PCMatApplyTranspose_MG(PC pc, Mat b, Mat x)
695: {
696: PetscFunctionBegin;
697: PetscCall(PCApply_MG_Internal(pc, NULL, NULL, b, x, PETSC_TRUE));
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: PetscErrorCode PCSetFromOptions_MG(PC pc, PetscOptionItems PetscOptionsObject)
702: {
703: PetscInt levels, cycles;
704: PetscBool flg, flg2;
705: PC_MG *mg = (PC_MG *)pc->data;
706: PC_MG_Levels **mglevels;
707: PCMGType mgtype;
708: PCMGCycleType mgctype;
709: PCMGGalerkinType gtype;
710: PCMGCoarseSpaceType coarseSpaceType;
712: PetscFunctionBegin;
713: levels = PetscMax(mg->nlevels, 1);
714: PetscOptionsHeadBegin(PetscOptionsObject, "Multigrid options");
715: PetscCall(PetscOptionsInt("-pc_mg_levels", "Number of Levels", "PCMGSetLevels", levels, &levels, &flg));
716: if (!flg && !mg->levels && pc->dm) {
717: PetscCall(DMGetRefineLevel(pc->dm, &levels));
718: levels++;
719: mg->usedmfornumberoflevels = PETSC_TRUE;
720: }
721: PetscCall(PCMGSetLevels(pc, levels, NULL));
722: mglevels = mg->levels;
724: mgctype = (PCMGCycleType)mglevels[0]->cycles;
725: PetscCall(PetscOptionsEnum("-pc_mg_cycle_type", "V cycle or for W-cycle", "PCMGSetCycleType", PCMGCycleTypes, (PetscEnum)mgctype, (PetscEnum *)&mgctype, &flg));
726: if (flg) PetscCall(PCMGSetCycleType(pc, mgctype));
727: coarseSpaceType = mg->coarseSpaceType;
728: PetscCall(PetscOptionsEnum("-pc_mg_adapt_interp_coarse_space", "Type of adaptive coarse space: none, polynomial, harmonic, eigenvector, generalized_eigenvector, gdsw", "PCMGSetAdaptCoarseSpaceType", PCMGCoarseSpaceTypes, (PetscEnum)coarseSpaceType, (PetscEnum *)&coarseSpaceType, &flg));
729: if (flg) PetscCall(PCMGSetAdaptCoarseSpaceType(pc, coarseSpaceType));
730: PetscCall(PetscOptionsInt("-pc_mg_adapt_interp_n", "Size of the coarse space for adaptive interpolation", "PCMGSetAdaptCoarseSpaceType", mg->Nc, &mg->Nc, &flg));
731: PetscCall(PetscOptionsBool("-pc_mg_mesp_monitor", "Monitor the multilevel eigensolver", "PCMGSetAdaptInterpolation", PETSC_FALSE, &mg->mespMonitor, &flg));
732: flg2 = PETSC_FALSE;
733: PetscCall(PetscOptionsBool("-pc_mg_adapt_cr", "Monitor coarse space quality using Compatible Relaxation (CR)", "PCMGSetAdaptCR", PETSC_FALSE, &flg2, &flg));
734: if (flg) PetscCall(PCMGSetAdaptCR(pc, flg2));
735: flg = PETSC_FALSE;
736: PetscCall(PetscOptionsBool("-pc_mg_distinct_smoothup", "Create separate smoothup KSP and append the prefix _up", "PCMGSetDistinctSmoothUp", PETSC_FALSE, &flg, NULL));
737: if (flg) PetscCall(PCMGSetDistinctSmoothUp(pc));
738: PetscCall(PetscOptionsEnum("-pc_mg_galerkin", "Use Galerkin process to compute coarser operators", "PCMGSetGalerkin", PCMGGalerkinTypes, (PetscEnum)mg->galerkin, (PetscEnum *)>ype, &flg));
739: if (flg) PetscCall(PCMGSetGalerkin(pc, gtype));
740: mgtype = mg->am;
741: PetscCall(PetscOptionsEnum("-pc_mg_type", "Multigrid type", "PCMGSetType", PCMGTypes, (PetscEnum)mgtype, (PetscEnum *)&mgtype, &flg));
742: if (flg) PetscCall(PCMGSetType(pc, mgtype));
743: if (mg->am == PC_MG_MULTIPLICATIVE) {
744: PetscCall(PetscOptionsInt("-pc_mg_multiplicative_cycles", "Number of cycles for each preconditioner step", "PCMGMultiplicativeSetCycles", mg->cyclesperpcapply, &cycles, &flg));
745: if (flg) PetscCall(PCMGMultiplicativeSetCycles(pc, cycles));
746: }
747: flg = PETSC_FALSE;
748: PetscCall(PetscOptionsBool("-pc_mg_log", "Log times for each multigrid level", "None", flg, &flg, NULL));
749: if (flg) {
750: PetscInt i;
751: char eventname[128];
753: levels = mglevels[0]->levels;
754: for (i = 0; i < levels; i++) {
755: PetscCall(PetscSNPrintf(eventname, PETSC_STATIC_ARRAY_LENGTH(eventname), "MGSetup Level %" PetscInt_FMT, i));
756: PetscCall(PetscLogEventRegister(eventname, ((PetscObject)pc)->classid, &mglevels[i]->eventsmoothsetup));
757: PetscCall(PetscSNPrintf(eventname, PETSC_STATIC_ARRAY_LENGTH(eventname), "MGSmooth Level %" PetscInt_FMT, i));
758: PetscCall(PetscLogEventRegister(eventname, ((PetscObject)pc)->classid, &mglevels[i]->eventsmoothsolve));
759: if (i) {
760: PetscCall(PetscSNPrintf(eventname, PETSC_STATIC_ARRAY_LENGTH(eventname), "MGResid Level %" PetscInt_FMT, i));
761: PetscCall(PetscLogEventRegister(eventname, ((PetscObject)pc)->classid, &mglevels[i]->eventresidual));
762: PetscCall(PetscSNPrintf(eventname, PETSC_STATIC_ARRAY_LENGTH(eventname), "MGInterp Level %" PetscInt_FMT, i));
763: PetscCall(PetscLogEventRegister(eventname, ((PetscObject)pc)->classid, &mglevels[i]->eventinterprestrict));
764: }
765: }
767: if (PetscDefined(USE_LOG)) {
768: const char sname[] = "MG Apply";
770: PetscCall(PetscLogStageGetId(sname, &mg->stageApply));
771: if (mg->stageApply < 0) PetscCall(PetscLogStageRegister(sname, &mg->stageApply));
772: }
773: }
774: PetscOptionsHeadEnd();
775: /* Check option consistency */
776: PetscCall(PCMGGetGalerkin(pc, >ype));
777: PetscCall(PCMGGetAdaptInterpolation(pc, &flg));
778: PetscCheck(!flg || !(gtype >= PC_MG_GALERKIN_NONE), PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use Galerkin coarse operators when adapting the interpolator");
779: PetscFunctionReturn(PETSC_SUCCESS);
780: }
782: const char *const PCMGTypes[] = {"MULTIPLICATIVE", "ADDITIVE", "FULL", "KASKADE", "PCMGType", "PC_MG", NULL};
783: const char *const PCMGCycleTypes[] = {"invalid", "v", "w", "PCMGCycleType", "PC_MG_CYCLE", NULL};
784: const char *const PCMGGalerkinTypes[] = {"both", "pmat", "mat", "none", "external", "PCMGGalerkinType", "PC_MG_GALERKIN", NULL};
785: const char *const PCMGCoarseSpaceTypes[] = {"none", "polynomial", "harmonic", "eigenvector", "generalized_eigenvector", "gdsw", "PCMGCoarseSpaceType", "PCMG_ADAPT_NONE", NULL};
787: #include <petscdraw.h>
788: PetscErrorCode PCView_MG(PC pc, PetscViewer viewer)
789: {
790: PC_MG *mg = (PC_MG *)pc->data;
791: PC_MG_Levels **mglevels = mg->levels;
792: PetscInt levels = mglevels ? mglevels[0]->levels : 0, i;
793: PetscBool isascii, isbinary, isdraw;
795: PetscFunctionBegin;
796: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
797: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
798: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
799: if (isascii) {
800: const char *cyclename = levels ? (mglevels[0]->cycles == PC_MG_CYCLE_V ? "v" : "w") : "unknown";
802: if (levels == 1) PetscCall(PetscViewerASCIIPrintf(viewer, " WARNING: Multigrid is being run with only a single level!\n"));
803: PetscCall(PetscViewerASCIIPrintf(viewer, " type is %s, levels=%" PetscInt_FMT " cycles=%s\n", PCMGTypes[mg->am], levels, cyclename));
804: if (mg->am == PC_MG_MULTIPLICATIVE) PetscCall(PetscViewerASCIIPrintf(viewer, " Cycles per PCApply=%" PetscInt_FMT "\n", mg->cyclesperpcapply));
805: if (mg->galerkin == PC_MG_GALERKIN_BOTH) {
806: PetscCall(PetscViewerASCIIPrintf(viewer, " Using Galerkin computed coarse grid matrices\n"));
807: } else if (mg->galerkin == PC_MG_GALERKIN_PMAT) {
808: PetscCall(PetscViewerASCIIPrintf(viewer, " Using Galerkin computed coarse grid matrices for pmat\n"));
809: } else if (mg->galerkin == PC_MG_GALERKIN_MAT) {
810: PetscCall(PetscViewerASCIIPrintf(viewer, " Using Galerkin computed coarse grid matrices for mat\n"));
811: } else if (mg->galerkin == PC_MG_GALERKIN_EXTERNAL) {
812: PetscCall(PetscViewerASCIIPrintf(viewer, " Using externally compute Galerkin coarse grid matrices\n"));
813: } else {
814: PetscCall(PetscViewerASCIIPrintf(viewer, " Not using Galerkin computed coarse grid matrices\n"));
815: }
816: if (mg->view) PetscCall((*mg->view)(pc, viewer));
817: for (i = 0; i < levels; i++) {
818: if (i) {
819: PetscCall(PetscViewerASCIIPrintf(viewer, "Down solver (pre-smoother) on level %" PetscInt_FMT " -------------------------------\n", i));
820: } else {
821: PetscCall(PetscViewerASCIIPrintf(viewer, "Coarse grid solver -- level %" PetscInt_FMT " -------------------------------\n", i));
822: }
823: PetscCall(PetscViewerASCIIPushTab(viewer));
824: PetscCall(KSPView(mglevels[i]->smoothd, viewer));
825: PetscCall(PetscViewerASCIIPopTab(viewer));
826: if (i && mglevels[i]->smoothd == mglevels[i]->smoothu) {
827: PetscCall(PetscViewerASCIIPrintf(viewer, "Up solver (post-smoother) same as down solver (pre-smoother)\n"));
828: } else if (i) {
829: PetscCall(PetscViewerASCIIPrintf(viewer, "Up solver (post-smoother) on level %" PetscInt_FMT " -------------------------------\n", i));
830: PetscCall(PetscViewerASCIIPushTab(viewer));
831: PetscCall(KSPView(mglevels[i]->smoothu, viewer));
832: PetscCall(PetscViewerASCIIPopTab(viewer));
833: }
834: if (i && mglevels[i]->cr) {
835: PetscCall(PetscViewerASCIIPrintf(viewer, "CR solver on level %" PetscInt_FMT " -------------------------------\n", i));
836: PetscCall(PetscViewerASCIIPushTab(viewer));
837: PetscCall(KSPView(mglevels[i]->cr, viewer));
838: PetscCall(PetscViewerASCIIPopTab(viewer));
839: }
840: }
841: } else if (isbinary) {
842: for (i = levels - 1; i >= 0; i--) {
843: PetscCall(KSPView(mglevels[i]->smoothd, viewer));
844: if (i && mglevels[i]->smoothd != mglevels[i]->smoothu) PetscCall(KSPView(mglevels[i]->smoothu, viewer));
845: }
846: } else if (isdraw) {
847: PetscDraw draw;
848: PetscReal x, w, y, bottom, th;
849: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
850: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
851: PetscCall(PetscDrawStringGetSize(draw, NULL, &th));
852: bottom = y - th;
853: for (i = levels - 1; i >= 0; i--) {
854: if (!mglevels[i]->smoothu || (mglevels[i]->smoothu == mglevels[i]->smoothd)) {
855: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
856: PetscCall(KSPView(mglevels[i]->smoothd, viewer));
857: PetscCall(PetscDrawPopCurrentPoint(draw));
858: } else {
859: w = 0.5 * PetscMin(1.0 - x, x);
860: PetscCall(PetscDrawPushCurrentPoint(draw, x + w, bottom));
861: PetscCall(KSPView(mglevels[i]->smoothd, viewer));
862: PetscCall(PetscDrawPopCurrentPoint(draw));
863: PetscCall(PetscDrawPushCurrentPoint(draw, x - w, bottom));
864: PetscCall(KSPView(mglevels[i]->smoothu, viewer));
865: PetscCall(PetscDrawPopCurrentPoint(draw));
866: }
867: PetscCall(PetscDrawGetBoundingBox(draw, NULL, &bottom, NULL, NULL));
868: bottom -= th;
869: }
870: }
871: PetscFunctionReturn(PETSC_SUCCESS);
872: }
874: #include <petsc/private/kspimpl.h>
876: /*
877: Calls setup for the KSP on each level
878: */
879: PetscErrorCode PCSetUp_MG(PC pc)
880: {
881: PC_MG *mg = (PC_MG *)pc->data;
882: PC_MG_Levels **mglevels = mg->levels;
883: PetscInt n;
884: PC cpc;
885: PetscBool dump = PETSC_FALSE, opsset, use_amat, missinginterpolate = PETSC_FALSE;
886: Mat dA, dB;
887: Vec tvec;
888: DM *dms;
889: PetscViewer viewer = NULL;
890: PetscBool dAeqdB = PETSC_FALSE, needRestricts = PETSC_FALSE, doCR = PETSC_FALSE;
891: PetscBool adaptInterpolation = mg->adaptInterpolation;
893: PetscFunctionBegin;
894: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels with PCMGSetLevels() before setting up");
895: n = mglevels[0]->levels;
896: /* FIX: Move this to PCSetFromOptions_MG? */
897: if (mg->usedmfornumberoflevels) {
898: PetscInt levels;
900: PetscCall(DMGetRefineLevel(pc->dm, &levels));
901: levels++;
902: if (levels > n) { /* the problem is now being solved on a finer grid */
903: PetscCall(PCMGSetLevels(pc, levels, NULL));
904: n = levels;
905: PetscCall(PCSetFromOptions(pc)); /* it is bad to call this here, but otherwise will never be called for the new hierarchy */
906: mglevels = mg->levels;
907: }
908: }
910: /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */
911: /* so use those from global PC */
912: /* Is this what we always want? What if user wants to keep old one? */
913: PetscCall(KSPGetOperatorsSet(mglevels[n - 1]->smoothd, NULL, &opsset));
914: if (opsset) {
915: Mat mmat;
916: PetscCall(KSPGetOperators(mglevels[n - 1]->smoothd, NULL, &mmat));
917: if (mmat == pc->pmat) opsset = PETSC_FALSE;
918: }
919: /* fine grid smoother inherits the reuse-pc flag */
920: PetscCall(KSPGetPC(mglevels[n - 1]->smoothd, &cpc));
921: cpc->reusepreconditioner = pc->reusepreconditioner;
922: PetscCall(KSPGetPC(mglevels[n - 1]->smoothu, &cpc));
923: cpc->reusepreconditioner = pc->reusepreconditioner;
925: /* Create CR solvers */
926: PetscCall(PCMGGetAdaptCR(pc, &doCR));
927: if (doCR) {
928: const char *prefix;
930: PetscCall(PCGetOptionsPrefix(pc, &prefix));
931: for (PetscInt i = 1; i < n; ++i) {
932: PC ipc, cr;
933: char crprefix[128];
935: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &mglevels[i]->cr));
936: PetscCall(KSPSetNestLevel(mglevels[i]->cr, pc->kspnestlevel));
937: PetscCall(KSPSetErrorIfNotConverged(mglevels[i]->cr, PETSC_FALSE));
938: PetscCall(PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->cr, (PetscObject)pc, n - i));
939: PetscCall(KSPSetOptionsPrefix(mglevels[i]->cr, prefix));
940: PetscCall(PetscObjectComposedDataSetInt((PetscObject)mglevels[i]->cr, PetscMGLevelId, mglevels[i]->level));
941: PetscCall(KSPSetType(mglevels[i]->cr, KSPCHEBYSHEV));
942: PetscCall(KSPSetConvergenceTest(mglevels[i]->cr, KSPConvergedSkip, NULL, NULL));
943: PetscCall(KSPSetNormType(mglevels[i]->cr, KSP_NORM_PRECONDITIONED));
944: PetscCall(KSPGetPC(mglevels[i]->cr, &ipc));
946: PetscCall(PCSetType(ipc, PCCOMPOSITE));
947: PetscCall(PCCompositeSetType(ipc, PC_COMPOSITE_MULTIPLICATIVE));
948: PetscCall(PCCompositeAddPCType(ipc, PCSOR));
949: PetscCall(CreateCR_Private(pc, i, &cr));
950: PetscCall(PCCompositeAddPC(ipc, cr));
951: PetscCall(PCDestroy(&cr));
953: PetscCall(KSPSetTolerances(mglevels[i]->cr, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, mg->default_smoothd));
954: PetscCall(KSPSetInitialGuessNonzero(mglevels[i]->cr, PETSC_TRUE));
955: PetscCall(PetscSNPrintf(crprefix, 128, "mg_levels_%" PetscInt_FMT "_cr_", i));
956: PetscCall(KSPAppendOptionsPrefix(mglevels[i]->cr, crprefix));
957: }
958: }
960: PetscCall(PCGetUseAmat(pc, &use_amat));
961: if (!opsset) {
962: if (use_amat) {
963: PetscCall(PetscInfo(pc, "Using outer operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n"));
964: PetscCall(KSPSetOperators(mglevels[n - 1]->smoothd, pc->mat, pc->pmat));
965: } else {
966: PetscCall(PetscInfo(pc, "Using matrix (pmat) operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n"));
967: PetscCall(KSPSetOperators(mglevels[n - 1]->smoothd, pc->pmat, pc->pmat));
968: }
969: }
971: for (PetscInt i = n - 1; i > 0; i--) {
972: if (!(mglevels[i]->interpolate || mglevels[i]->restrct)) {
973: missinginterpolate = PETSC_TRUE;
974: break;
975: }
976: }
978: PetscCall(KSPGetOperators(mglevels[n - 1]->smoothd, &dA, &dB));
979: if (dA == dB) dAeqdB = PETSC_TRUE;
980: if (mg->galerkin == PC_MG_GALERKIN_NONE || ((mg->galerkin == PC_MG_GALERKIN_PMAT || mg->galerkin == PC_MG_GALERKIN_MAT) && !dAeqdB)) needRestricts = PETSC_TRUE; /* user must compute either mat, pmat, or both so must restrict x to coarser levels */
982: if (pc->dm && !pc->setupcalled) {
983: /* finest smoother also gets DM but it is not active, independent of whether galerkin==PC_MG_GALERKIN_EXTERNAL */
984: PetscCall(KSPSetDM(mglevels[n - 1]->smoothd, pc->dm));
985: PetscCall(KSPSetDMActive(mglevels[n - 1]->smoothd, KSP_DMACTIVE_ALL, PETSC_FALSE));
986: if (mglevels[n - 1]->smoothd != mglevels[n - 1]->smoothu) {
987: PetscCall(KSPSetDM(mglevels[n - 1]->smoothu, pc->dm));
988: PetscCall(KSPSetDMActive(mglevels[n - 1]->smoothu, KSP_DMACTIVE_ALL, PETSC_FALSE));
989: }
990: if (mglevels[n - 1]->cr) {
991: PetscCall(KSPSetDM(mglevels[n - 1]->cr, pc->dm));
992: PetscCall(KSPSetDMActive(mglevels[n - 1]->cr, KSP_DMACTIVE_ALL, PETSC_FALSE));
993: }
994: }
996: /*
997: Skipping if user has provided all interpolation/restriction needed (since DM might not be able to produce them (when coming from SNES/TS)
998: Skipping for externally managed hierarchy (such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs?
999: */
1000: if (missinginterpolate && mg->galerkin != PC_MG_GALERKIN_EXTERNAL && !pc->setupcalled) {
1001: /* first see if we can compute a coarse space */
1002: if (mg->coarseSpaceType == PCMG_ADAPT_GDSW) {
1003: for (PetscInt i = n - 2; i > -1; i--) {
1004: if (!mglevels[i + 1]->restrct && !mglevels[i + 1]->interpolate) {
1005: PetscCall(PCMGComputeCoarseSpace_Internal(pc, i + 1, mg->coarseSpaceType, mg->Nc, NULL, &mglevels[i + 1]->coarseSpace));
1006: PetscCall(PCMGSetInterpolation(pc, i + 1, mglevels[i + 1]->coarseSpace));
1007: }
1008: }
1009: } else { /* construct the interpolation from the DMs */
1010: Mat p;
1011: Vec rscale;
1013: PetscCheck(n == 1 || pc->dm, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PC lacks a DM so cannot automatically construct a multigrid hierarchy. Number of levels requested %" PetscInt_FMT, n);
1014: PetscCall(PetscMalloc1(n, &dms));
1015: dms[n - 1] = pc->dm;
1016: /* Separately create them so we do not get DMKSP interference between levels */
1017: for (PetscInt i = n - 2; i > -1; i--) PetscCall(DMCoarsen(dms[i + 1], MPI_COMM_NULL, &dms[i]));
1018: for (PetscInt i = n - 2; i > -1; i--) {
1019: PetscBool dmhasrestrict, dmhasinject;
1021: PetscCall(KSPSetDM(mglevels[i]->smoothd, dms[i]));
1022: if (!needRestricts) PetscCall(KSPSetDMActive(mglevels[i]->smoothd, KSP_DMACTIVE_ALL, PETSC_FALSE));
1023: PetscCall(KSPSetDMActive(mglevels[i]->smoothd, KSP_DMACTIVE_RHS, PETSC_FALSE));
1024: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
1025: PetscCall(KSPSetDM(mglevels[i]->smoothu, dms[i]));
1026: if (!needRestricts) PetscCall(KSPSetDMActive(mglevels[i]->smoothu, KSP_DMACTIVE_ALL, PETSC_FALSE));
1027: PetscCall(KSPSetDMActive(mglevels[i]->smoothu, KSP_DMACTIVE_RHS, PETSC_FALSE));
1028: }
1029: if (mglevels[i]->cr) {
1030: PetscCall(KSPSetDM(mglevels[i]->cr, dms[i]));
1031: if (!needRestricts) PetscCall(KSPSetDMActive(mglevels[i]->cr, KSP_DMACTIVE_ALL, PETSC_FALSE));
1032: PetscCall(KSPSetDMActive(mglevels[i]->cr, KSP_DMACTIVE_RHS, PETSC_FALSE));
1033: }
1034: if (!mglevels[i + 1]->interpolate) {
1035: PetscCall(DMCreateInterpolation(dms[i], dms[i + 1], &p, &rscale));
1036: PetscCall(PCMGSetInterpolation(pc, i + 1, p));
1037: if (rscale) PetscCall(PCMGSetRScale(pc, i + 1, rscale));
1038: PetscCall(VecDestroy(&rscale));
1039: PetscCall(MatDestroy(&p));
1040: }
1041: PetscCall(DMHasCreateRestriction(dms[i], &dmhasrestrict));
1042: if (dmhasrestrict && !mglevels[i + 1]->restrct) {
1043: PetscCall(DMCreateRestriction(dms[i], dms[i + 1], &p));
1044: PetscCall(PCMGSetRestriction(pc, i + 1, p));
1045: PetscCall(MatDestroy(&p));
1046: }
1047: PetscCall(DMHasCreateInjection(dms[i], &dmhasinject));
1048: if (dmhasinject && !mglevels[i + 1]->inject) {
1049: PetscCall(DMCreateInjection(dms[i], dms[i + 1], &p));
1050: PetscCall(PCMGSetInjection(pc, i + 1, p));
1051: PetscCall(MatDestroy(&p));
1052: }
1053: }
1055: for (PetscInt i = n - 2; i > -1; i--) PetscCall(DMDestroy(&dms[i]));
1056: PetscCall(PetscFree(dms));
1057: }
1058: }
1060: if (mg->galerkin < PC_MG_GALERKIN_NONE) {
1061: Mat A, B;
1062: PetscBool doA = PETSC_FALSE, doB = PETSC_FALSE;
1063: MatReuse reuse = MAT_INITIAL_MATRIX;
1065: if (mg->galerkin == PC_MG_GALERKIN_PMAT || mg->galerkin == PC_MG_GALERKIN_BOTH) doB = PETSC_TRUE;
1066: if (mg->galerkin == PC_MG_GALERKIN_MAT || (mg->galerkin == PC_MG_GALERKIN_BOTH && dA != dB)) doA = PETSC_TRUE;
1067: PetscCheck(!doA || use_amat, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PC_MG_GALERKIN_MAT and PCSetUseAmat(pc, PETSC_FALSE) are incompatible options");
1068: if (pc->setupcalled) reuse = MAT_REUSE_MATRIX;
1069: for (PetscInt i = n - 2; i > -1; i--) {
1070: PetscCheck(mglevels[i + 1]->restrct || mglevels[i + 1]->interpolate, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must provide interpolation or restriction for each MG level except level 0");
1071: if (!mglevels[i + 1]->interpolate) PetscCall(PCMGSetInterpolation(pc, i + 1, mglevels[i + 1]->restrct));
1072: if (!mglevels[i + 1]->restrct) PetscCall(PCMGSetRestriction(pc, i + 1, mglevels[i + 1]->interpolate));
1073: if (reuse == MAT_REUSE_MATRIX) PetscCall(KSPGetOperators(mglevels[i]->smoothd, &A, &B));
1074: if (doA) PetscCall(MatGalerkin(mglevels[i + 1]->restrct, dA, mglevels[i + 1]->interpolate, reuse, 1.0, &A));
1075: if (doB) PetscCall(MatGalerkin(mglevels[i + 1]->restrct, dB, mglevels[i + 1]->interpolate, reuse, 1.0, &B));
1076: /* the management of the PetscObjectReference() and PetscObjecDereference() below is rather delicate */
1077: if (!doA && dAeqdB) {
1078: if (reuse == MAT_INITIAL_MATRIX) PetscCall(PetscObjectReference((PetscObject)B));
1079: A = B;
1080: } else if (!doA && reuse == MAT_INITIAL_MATRIX) {
1081: PetscCall(KSPGetOperators(mglevels[i]->smoothd, &A, NULL));
1082: PetscCall(PetscObjectReference((PetscObject)A));
1083: }
1084: if (!doB && dAeqdB) {
1085: if (reuse == MAT_INITIAL_MATRIX) PetscCall(PetscObjectReference((PetscObject)A));
1086: B = A;
1087: } else if (!doB && reuse == MAT_INITIAL_MATRIX) {
1088: PetscCall(KSPGetOperators(mglevels[i]->smoothd, NULL, &B));
1089: PetscCall(PetscObjectReference((PetscObject)B));
1090: }
1091: if (reuse == MAT_INITIAL_MATRIX) {
1092: PetscCall(KSPSetOperators(mglevels[i]->smoothd, A, B));
1093: PetscCall(PetscObjectDereference((PetscObject)A));
1094: PetscCall(PetscObjectDereference((PetscObject)B));
1095: }
1096: dA = A;
1097: dB = B;
1098: }
1099: } else { /* PC_MG_GALERKIN_NONE */
1100: if (!use_amat) { /* force KSP(P, P) at all levels */
1101: for (PetscInt i = n - 1; i > -1; i--) {
1102: Mat B;
1103: PetscBool Bopset;
1104: KSP smoothd = mglevels[i]->smoothd;
1106: PetscCall(KSPGetOperatorsSet(smoothd, NULL, &Bopset));
1107: /* This is a chicken-and-egg problem when DMKSP has a create operator callback.
1108: It would be called at KSPSetUp stage, but then it is too late to handle the amat = False case */
1109: if (!Bopset && (smoothd->dmActive & KSP_DMACTIVE_OPERATOR) && smoothd->dm) {
1110: DMKSP kdm;
1112: PetscCall(DMGetDMKSP(smoothd->dm, &kdm));
1113: if (kdm->ops->createoperators) {
1114: Mat A;
1116: A = B = NULL;
1117: PetscCallBack("KSP callback create operators", (*kdm->ops->createoperators)(smoothd, &A, &B, kdm->createoperatorsctx));
1118: PetscCheck(A, PetscObjectComm((PetscObject)smoothd), PETSC_ERR_ARG_WRONGSTATE, "Missing A operator from DMKSPSetCreateOperators() callback");
1119: if (!B) B = A;
1120: if (B == A) PetscCall(PetscObjectReference((PetscObject)B));
1121: PetscCall(KSPSetOperators(smoothd, B, B));
1122: PetscCall(MatDestroy(&A));
1123: PetscCall(MatDestroy(&B));
1124: }
1125: } else {
1126: PetscCheck(Bopset, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Missing Pmat level %" PetscInt_FMT, i);
1127: PetscCall(KSPGetOperators(smoothd, NULL, &B));
1128: PetscCall(KSPSetOperators(smoothd, B, B));
1129: }
1130: }
1131: }
1132: }
1134: /* Adapt interpolation matrices */
1135: if (adaptInterpolation) {
1136: for (PetscInt i = 0; i < n; ++i) {
1137: if (!mglevels[i]->coarseSpace) PetscCall(PCMGComputeCoarseSpace_Internal(pc, i, mg->coarseSpaceType, mg->Nc, !i ? NULL : mglevels[i - 1]->coarseSpace, &mglevels[i]->coarseSpace));
1138: if (i) PetscCall(PCMGAdaptInterpolator_Internal(pc, i, mglevels[i - 1]->smoothu, mglevels[i]->smoothu, mglevels[i - 1]->coarseSpace, mglevels[i]->coarseSpace));
1139: }
1140: for (PetscInt i = n - 2; i > -1; --i) PetscCall(PCMGRecomputeLevelOperators_Internal(pc, i));
1141: }
1143: if (needRestricts && pc->dm) {
1144: for (PetscInt i = n - 2; i >= 0; i--) {
1145: DM dmfine, dmcoarse;
1146: Mat Restrict, Inject;
1147: Vec rscale;
1149: PetscCall(KSPGetDM(mglevels[i + 1]->smoothd, &dmfine));
1150: PetscCall(KSPGetDM(mglevels[i]->smoothd, &dmcoarse));
1151: PetscCall(PCMGGetRestriction(pc, i + 1, &Restrict));
1152: PetscCall(PCMGGetRScale(pc, i + 1, &rscale));
1153: PetscCall(PCMGGetInjection(pc, i + 1, &Inject));
1154: PetscCall(DMRestrict(dmfine, Restrict, rscale, Inject, dmcoarse));
1155: }
1156: }
1158: if (!pc->setupcalled) {
1159: for (PetscInt i = 0; i < n; i++) PetscCall(KSPSetFromOptions(mglevels[i]->smoothd));
1160: for (PetscInt i = 1; i < n; i++) {
1161: if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) PetscCall(KSPSetFromOptions(mglevels[i]->smoothu));
1162: if (mglevels[i]->cr) PetscCall(KSPSetFromOptions(mglevels[i]->cr));
1163: }
1164: /* insure that if either interpolation or restriction is set the other one is set */
1165: for (PetscInt i = 1; i < n; i++) {
1166: PetscCall(PCMGGetInterpolation(pc, i, NULL));
1167: PetscCall(PCMGGetRestriction(pc, i, NULL));
1168: }
1169: for (PetscInt i = 0; i < n - 1; i++) {
1170: if (!mglevels[i]->b) {
1171: Vec *vec;
1172: PetscCall(KSPCreateVecs(mglevels[i]->smoothd, 1, &vec, 0, NULL));
1173: PetscCall(PCMGSetRhs(pc, i, *vec));
1174: PetscCall(VecDestroy(vec));
1175: PetscCall(PetscFree(vec));
1176: }
1177: if (!mglevels[i]->r && i) {
1178: PetscCall(VecDuplicate(mglevels[i]->b, &tvec));
1179: PetscCall(PCMGSetR(pc, i, tvec));
1180: PetscCall(VecDestroy(&tvec));
1181: }
1182: if (!mglevels[i]->x) {
1183: PetscCall(VecDuplicate(mglevels[i]->b, &tvec));
1184: PetscCall(PCMGSetX(pc, i, tvec));
1185: PetscCall(VecDestroy(&tvec));
1186: }
1187: if (doCR) {
1188: PetscCall(VecDuplicate(mglevels[i]->b, &mglevels[i]->crx));
1189: PetscCall(VecDuplicate(mglevels[i]->b, &mglevels[i]->crb));
1190: }
1191: }
1192: if (n != 1 && !mglevels[n - 1]->r) {
1193: /* PCMGSetR() on the finest level if user did not supply it */
1194: Vec *vec;
1196: PetscCall(KSPCreateVecs(mglevels[n - 1]->smoothd, 1, &vec, 0, NULL));
1197: PetscCall(PCMGSetR(pc, n - 1, *vec));
1198: PetscCall(VecDestroy(vec));
1199: PetscCall(PetscFree(vec));
1200: }
1201: if (doCR) {
1202: PetscCall(VecDuplicate(mglevels[n - 1]->r, &mglevels[n - 1]->crx));
1203: PetscCall(VecDuplicate(mglevels[n - 1]->r, &mglevels[n - 1]->crb));
1204: }
1205: }
1207: if (pc->dm) {
1208: /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */
1209: for (PetscInt i = 0; i < n - 1; i++) {
1210: if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
1211: }
1212: }
1213: // We got here (PCSetUp_MG) because the matrix has changed, which means the smoother needs to be set up again (e.g.,
1214: // new diagonal for Jacobi). Setting it here allows it to be logged under PCSetUp rather than deep inside a PCApply.
1215: if (mglevels[n - 1]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[n - 1]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
1217: for (PetscInt i = 1; i < n; i++) {
1218: PetscBool wcycle = (PetscBool)(mg->am == PC_MG_MULTIPLICATIVE && i < n - 1 && mglevels[i + 1]->cycles > 1);
1220: if (mglevels[i]->smoothu == mglevels[i]->smoothd || mg->am == PC_MG_FULL || mg->am == PC_MG_KASKADE || mg->cyclesperpcapply > 1 || wcycle) {
1221: /* if going only down then initial guess is zero, unless a W cycle re-enters this level with the previous sub-cycle's solution */
1222: if (wcycle && mglevels[i]->smoothu != mglevels[i]->smoothd) PetscCall(PCMGCheckSmootherDownGuess_Private(pc, mglevels[i]->smoothd, "A W cycle"));
1223: PetscCall(KSPSetInitialGuessNonzero(mglevels[i]->smoothd, PETSC_TRUE));
1224: }
1225: if (mglevels[i]->cr) PetscCall(KSPSetInitialGuessNonzero(mglevels[i]->cr, PETSC_TRUE));
1226: if (mglevels[i]->eventsmoothsetup) PetscCall(PetscLogEventBegin(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0));
1227: PetscCall(KSPSetUp(mglevels[i]->smoothd));
1228: if (mglevels[i]->smoothd->reason) pc->failedreason = PC_SUBPC_ERROR;
1229: if (mglevels[i]->eventsmoothsetup) PetscCall(PetscLogEventEnd(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0));
1230: if (!mglevels[i]->residual) {
1231: Mat mat;
1233: PetscCall(KSPGetOperators(mglevels[i]->smoothd, &mat, NULL));
1234: PetscCall(PCMGSetResidual(pc, i, PCMGResidualDefault, mat));
1235: }
1236: if (!mglevels[i]->residualtranspose) {
1237: Mat mat;
1239: PetscCall(KSPGetOperators(mglevels[i]->smoothd, &mat, NULL));
1240: PetscCall(PCMGSetResidualTranspose(pc, i, PCMGResidualTransposeDefault, mat));
1241: }
1242: }
1243: for (PetscInt i = 1; i < n; i++) {
1244: if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) {
1245: Mat downmat, downpmat;
1247: /* check if operators have been set for up, if not use down operators to set them */
1248: PetscCall(KSPGetOperatorsSet(mglevels[i]->smoothu, &opsset, NULL));
1249: if (!opsset) {
1250: PetscCall(KSPGetOperators(mglevels[i]->smoothd, &downmat, &downpmat));
1251: PetscCall(KSPSetOperators(mglevels[i]->smoothu, downmat, downpmat));
1252: }
1254: PetscCall(KSPSetInitialGuessNonzero(mglevels[i]->smoothu, PETSC_TRUE));
1255: if (mglevels[i]->eventsmoothsetup) PetscCall(PetscLogEventBegin(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0));
1256: PetscCall(KSPSetUp(mglevels[i]->smoothu));
1257: if (mglevels[i]->smoothu->reason) pc->failedreason = PC_SUBPC_ERROR;
1258: if (mglevels[i]->eventsmoothsetup) PetscCall(PetscLogEventEnd(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0));
1259: }
1260: if (mglevels[i]->cr) {
1261: Mat downmat, downpmat;
1263: /* check if operators have been set for up, if not use down operators to set them */
1264: PetscCall(KSPGetOperatorsSet(mglevels[i]->cr, &opsset, NULL));
1265: if (!opsset) {
1266: PetscCall(KSPGetOperators(mglevels[i]->smoothd, &downmat, &downpmat));
1267: PetscCall(KSPSetOperators(mglevels[i]->cr, downmat, downpmat));
1268: }
1270: PetscCall(KSPSetInitialGuessNonzero(mglevels[i]->cr, PETSC_TRUE));
1271: if (mglevels[i]->eventsmoothsetup) PetscCall(PetscLogEventBegin(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0));
1272: PetscCall(KSPSetUp(mglevels[i]->cr));
1273: if (mglevels[i]->cr->reason) pc->failedreason = PC_SUBPC_ERROR;
1274: if (mglevels[i]->eventsmoothsetup) PetscCall(PetscLogEventEnd(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0));
1275: }
1276: }
1278: if (mglevels[0]->eventsmoothsetup) PetscCall(PetscLogEventBegin(mglevels[0]->eventsmoothsetup, 0, 0, 0, 0));
1279: PetscCall(KSPSetUp(mglevels[0]->smoothd));
1280: if (mglevels[0]->smoothd->reason) pc->failedreason = PC_SUBPC_ERROR;
1281: if (mglevels[0]->eventsmoothsetup) PetscCall(PetscLogEventEnd(mglevels[0]->eventsmoothsetup, 0, 0, 0, 0));
1283: /*
1284: Dump the interpolation/restriction matrices plus the
1285: Jacobian/stiffness on each level. This allows MATLAB users to
1286: easily check if the Galerkin condition A_c = R A_f R^T is satisfied.
1288: Only support one or the other at the same time.
1289: */
1290: #if PetscDefined(USE_SOCKET_VIEWER)
1291: PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_mg_dump_matlab", &dump, NULL));
1292: if (dump) viewer = PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc));
1293: dump = PETSC_FALSE;
1294: #endif
1295: PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_mg_dump_binary", &dump, NULL));
1296: if (dump) viewer = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc));
1298: if (viewer) {
1299: for (PetscInt i = 1; i < n; i++) PetscCall(MatView(mglevels[i]->restrct, viewer));
1300: for (PetscInt i = 0; i < n; i++) {
1301: PetscCall(KSPGetPC(mglevels[i]->smoothd, &pc));
1302: PetscCall(MatView(pc->mat, viewer));
1303: }
1304: }
1305: PetscFunctionReturn(PETSC_SUCCESS);
1306: }
1308: PetscErrorCode PCMGGetLevels_MG(PC pc, PetscInt *levels)
1309: {
1310: PC_MG *mg = (PC_MG *)pc->data;
1312: PetscFunctionBegin;
1313: *levels = mg->nlevels;
1314: PetscFunctionReturn(PETSC_SUCCESS);
1315: }
1317: /*@
1318: PCMGGetLevels - Gets the number of levels to use with `PCMG`.
1320: Not Collective
1322: Input Parameter:
1323: . pc - the preconditioner context
1325: Output Parameter:
1326: . levels - the number of levels
1328: Level: advanced
1330: .seealso: [](ch_ksp), `PCMG`, `PCMGSetLevels()`
1331: @*/
1332: PetscErrorCode PCMGGetLevels(PC pc, PetscInt *levels)
1333: {
1334: PetscFunctionBegin;
1336: PetscAssertPointer(levels, 2);
1337: *levels = 0;
1338: PetscTryMethod(pc, "PCMGGetLevels_C", (PC, PetscInt *), (pc, levels));
1339: PetscFunctionReturn(PETSC_SUCCESS);
1340: }
1342: /*@
1343: PCMGGetGridComplexity - compute operator and grid complexity of the `PCMG` hierarchy
1345: Input Parameter:
1346: . pc - the preconditioner context
1348: Output Parameters:
1349: + gc - grid complexity, $\frac{\sum_i n_i}{n_0}$, where $n_0$ is the number of unknowns on the finest grid
1350: - oc - operator complexity, $\frac{\sum_i nnz_i}{nnz_0}$, where $nnz_0$ is the number of nonzeros on the finest grid
1352: Level: advanced
1354: .seealso: [](ch_ksp), `PCMG`, `PCMGGetLevels()`, `PCMGSetLevels()`
1355: @*/
1356: PetscErrorCode PCMGGetGridComplexity(PC pc, PetscReal *gc, PetscReal *oc)
1357: {
1358: PC_MG *mg = (PC_MG *)pc->data;
1359: PC_MG_Levels **mglevels = mg->levels;
1360: PetscInt N;
1361: PetscLogDouble nnz0 = 0, sgc = 0, soc = 0, n0 = 0;
1362: MatInfo info;
1364: PetscFunctionBegin;
1366: if (gc) PetscAssertPointer(gc, 2);
1367: if (oc) PetscAssertPointer(oc, 3);
1368: if (!pc->setupcalled) {
1369: if (gc) *gc = 0;
1370: if (oc) *oc = 0;
1371: PetscFunctionReturn(PETSC_SUCCESS);
1372: }
1373: PetscCheck(mg->nlevels > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MG has no levels");
1374: for (PetscInt lev = 0; lev < mg->nlevels; lev++) {
1375: Mat dB;
1376: PetscCall(KSPGetOperators(mglevels[lev]->smoothd, NULL, &dB));
1377: PetscCall(MatGetInfo(dB, MAT_GLOBAL_SUM, &info)); /* global reduction */
1378: PetscCall(MatGetSize(dB, &N, NULL));
1379: sgc += N;
1380: soc += info.nz_used;
1381: if (lev == mg->nlevels - 1) {
1382: nnz0 = info.nz_used;
1383: n0 = N;
1384: }
1385: }
1386: PetscCheck(n0 > 0 && gc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number for grid points on finest level is not available");
1387: *gc = (PetscReal)(sgc / n0);
1388: if (nnz0 > 0 && oc) *oc = (PetscReal)(soc / nnz0);
1389: PetscFunctionReturn(PETSC_SUCCESS);
1390: }
1392: /*@
1393: PCMGSetType - Determines the type of multigrid to use, either
1394: multiplicative, additive, full, or the Kaskade algorithm.
1396: Logically Collective
1398: Input Parameters:
1399: + pc - the preconditioner context
1400: - form - multigrid form, one of `PC_MG_MULTIPLICATIVE`, `PC_MG_ADDITIVE`, `PC_MG_FULL`, `PC_MG_KASKADE`
1402: Options Database Key:
1403: . -pc_mg_type (multiplicative|additive|full|kaskade) - Sets the type
1405: Level: advanced
1407: .seealso: [](ch_ksp), `PCMGType`, `PCMG`, `PCMGGetLevels()`, `PCMGSetLevels()`, `PCMGGetType()`, `PCMGCycleType`,
1408: `PC_MG_MULTIPLICATIVE`, `PC_MG_ADDITIVE`, `PC_MG_FULL`, `PC_MG_KASKADE`
1409: @*/
1410: PetscErrorCode PCMGSetType(PC pc, PCMGType form)
1411: {
1412: PC_MG *mg = (PC_MG *)pc->data;
1414: PetscFunctionBegin;
1417: mg->am = form;
1418: if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
1419: else pc->ops->applyrichardson = NULL;
1420: PetscFunctionReturn(PETSC_SUCCESS);
1421: }
1423: /*@
1424: PCMGGetType - Finds the form of multigrid the `PCMG` is using multiplicative, additive, full, or the Kaskade algorithm.
1426: Logically Collective
1428: Input Parameter:
1429: . pc - the preconditioner context
1431: Output Parameter:
1432: . type - one of `PC_MG_MULTIPLICATIVE`, `PC_MG_ADDITIVE`, `PC_MG_FULL`, `PC_MG_KASKADE`, `PCMGCycleType`
1434: Level: advanced
1436: .seealso: [](ch_ksp), `PCMGType`, `PCMG`, `PCMGGetLevels()`, `PCMGSetLevels()`, `PCMGSetType()`,
1437: `PC_MG_MULTIPLICATIVE`, `PC_MG_ADDITIVE`, `PC_MG_FULL`, `PC_MG_KASKADE`
1438: @*/
1439: PetscErrorCode PCMGGetType(PC pc, PCMGType *type)
1440: {
1441: PC_MG *mg = (PC_MG *)pc->data;
1443: PetscFunctionBegin;
1445: *type = mg->am;
1446: PetscFunctionReturn(PETSC_SUCCESS);
1447: }
1449: /*@
1450: PCMGSetCycleType - Sets the type of cycles to use. Use `PCMGSetCycleTypeOnLevel()` for more
1451: complicated cycling.
1453: Logically Collective
1455: Input Parameters:
1456: + pc - the multigrid context
1457: - n - either `PC_MG_CYCLE_V` or `PC_MG_CYCLE_W`
1459: Options Database Key:
1460: . -pc_mg_cycle_type (v|w) - provide the cycle desired
1462: Level: advanced
1464: .seealso: [](ch_ksp), `PCMG`, `PCMGSetCycleTypeOnLevel()`, `PCMGType`, `PCMGCycleType`, `PC_MG_CYCLE_V`, `PC_MG_CYCLE_W`
1465: @*/
1466: PetscErrorCode PCMGSetCycleType(PC pc, PCMGCycleType n)
1467: {
1468: PC_MG *mg = (PC_MG *)pc->data;
1469: PC_MG_Levels **mglevels = mg->levels;
1470: PetscInt i, levels;
1472: PetscFunctionBegin;
1475: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Must set MG levels with PCMGSetLevels() before calling");
1476: levels = mglevels[0]->levels;
1477: for (i = 0; i < levels; i++) mglevels[i]->cycles = n;
1478: PetscFunctionReturn(PETSC_SUCCESS);
1479: }
1481: /*@
1482: PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step
1483: of multigrid when `PCMGType` is `PC_MG_MULTIPLICATIVE`
1485: Logically Collective
1487: Input Parameters:
1488: + pc - the multigrid context
1489: - n - number of cycles (default is 1)
1491: Options Database Key:
1492: . -pc_mg_multiplicative_cycles n - set the number of cycles
1494: Level: advanced
1496: Note:
1497: This is not associated with setting a V or W cycle, that is set with `PCMGSetCycleType()`
1499: .seealso: [](ch_ksp), `PCMGSetCycleTypeOnLevel()`, `PCMGSetCycleType()`, `PCMGCycleType`, `PCMGType`, `PC_MG_MULTIPLICATIVE`
1500: @*/
1501: PetscErrorCode PCMGMultiplicativeSetCycles(PC pc, PetscInt n)
1502: {
1503: PC_MG *mg = (PC_MG *)pc->data;
1505: PetscFunctionBegin;
1508: mg->cyclesperpcapply = n;
1509: PetscFunctionReturn(PETSC_SUCCESS);
1510: }
1512: /*
1513: Since the finest level KSP shares the original matrix (of the entire system), it's preconditioner
1514: should not be updated if the whole PC is supposed to reuse the preconditioner
1515: */
1516: static PetscErrorCode PCSetReusePreconditioner_MG(PC pc, PetscBool flag)
1517: {
1518: PC_MG *mg = (PC_MG *)pc->data;
1519: PC_MG_Levels **mglevels = mg->levels;
1520: PetscInt levels;
1521: PC tpc;
1523: PetscFunctionBegin;
1526: if (mglevels) {
1527: levels = mglevels[0]->levels;
1528: PetscCall(KSPGetPC(mglevels[levels - 1]->smoothd, &tpc));
1529: tpc->reusepreconditioner = flag;
1530: PetscCall(KSPGetPC(mglevels[levels - 1]->smoothu, &tpc));
1531: tpc->reusepreconditioner = flag;
1532: }
1533: PetscFunctionReturn(PETSC_SUCCESS);
1534: }
1536: static PetscErrorCode PCMGSetGalerkin_MG(PC pc, PCMGGalerkinType use)
1537: {
1538: PC_MG *mg = (PC_MG *)pc->data;
1540: PetscFunctionBegin;
1541: mg->galerkin = use;
1542: PetscFunctionReturn(PETSC_SUCCESS);
1543: }
1545: /*@
1546: PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the
1547: finest grid via the Galerkin process: $A_{i-1} = r_i * A_i * p_i$.
1549: Logically Collective
1551: Input Parameters:
1552: + pc - the multigrid context
1553: - use - one of `PC_MG_GALERKIN_BOTH`, `PC_MG_GALERKIN_PMAT`, `PC_MG_GALERKIN_MAT`, or `PC_MG_GALERKIN_NONE`
1555: Options Database Key:
1556: . -pc_mg_galerkin (both|pmat|mat|none) - set the matrices to form via the Galerkin process
1558: Level: intermediate
1560: Notes:
1561: Some codes that use `PCMG` such as `PCGAMG` use Galerkin internally while constructing the hierarchy and thus do not
1562: use the `PCMG` construction of the coarser grids.
1564: If this is not used the coarser grid matrices are computed via re-discretization. That is by calling the function associated with the `DM` attached
1565: to the `PC`. For example, for nonlinear solves the function provided with `SNESSetJacobian()`.
1567: .seealso: [](ch_ksp), `PCMG`, `PCMGGetGalerkin()`, `PCMGGalerkinType`, `PC_MG_GALERKIN_BOTH`, `PC_MG_GALERKIN_PMAT`, `PC_MG_GALERKIN_MAT`, `PC_MG_GALERKIN_NONE`
1568: @*/
1569: PetscErrorCode PCMGSetGalerkin(PC pc, PCMGGalerkinType use)
1570: {
1571: PetscFunctionBegin;
1573: PetscTryMethod(pc, "PCMGSetGalerkin_C", (PC, PCMGGalerkinType), (pc, use));
1574: PetscFunctionReturn(PETSC_SUCCESS);
1575: }
1577: /*@
1578: PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e. $A_{i-1} = r_i * A_i * p_i$.
1580: Not Collective
1582: Input Parameter:
1583: . pc - the multigrid context
1585: Output Parameter:
1586: . galerkin - one of `PC_MG_GALERKIN_BOTH`,`PC_MG_GALERKIN_PMAT`,`PC_MG_GALERKIN_MAT`, `PC_MG_GALERKIN_NONE`, or `PC_MG_GALERKIN_EXTERNAL`
1588: Level: intermediate
1590: .seealso: [](ch_ksp), `PCMG`, `PCMGSetGalerkin()`, `PCMGGalerkinType`, `PC_MG_GALERKIN_BOTH`, `PC_MG_GALERKIN_PMAT`, `PC_MG_GALERKIN_MAT`, `PC_MG_GALERKIN_NONE`, `PC_MG_GALERKIN_EXTERNAL`
1591: @*/
1592: PetscErrorCode PCMGGetGalerkin(PC pc, PCMGGalerkinType *galerkin)
1593: {
1594: PC_MG *mg = (PC_MG *)pc->data;
1596: PetscFunctionBegin;
1598: *galerkin = mg->galerkin;
1599: PetscFunctionReturn(PETSC_SUCCESS);
1600: }
1602: static PetscErrorCode PCMGSetAdaptInterpolation_MG(PC pc, PetscBool adapt)
1603: {
1604: PC_MG *mg = (PC_MG *)pc->data;
1606: PetscFunctionBegin;
1607: mg->adaptInterpolation = adapt;
1608: PetscFunctionReturn(PETSC_SUCCESS);
1609: }
1611: static PetscErrorCode PCMGGetAdaptInterpolation_MG(PC pc, PetscBool *adapt)
1612: {
1613: PC_MG *mg = (PC_MG *)pc->data;
1615: PetscFunctionBegin;
1616: *adapt = mg->adaptInterpolation;
1617: PetscFunctionReturn(PETSC_SUCCESS);
1618: }
1620: static PetscErrorCode PCMGSetAdaptCoarseSpaceType_MG(PC pc, PCMGCoarseSpaceType ctype)
1621: {
1622: PC_MG *mg = (PC_MG *)pc->data;
1624: PetscFunctionBegin;
1625: mg->adaptInterpolation = ctype != PCMG_ADAPT_NONE ? PETSC_TRUE : PETSC_FALSE;
1626: mg->coarseSpaceType = ctype;
1627: PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_BOTH));
1628: PetscFunctionReturn(PETSC_SUCCESS);
1629: }
1631: static PetscErrorCode PCMGGetAdaptCoarseSpaceType_MG(PC pc, PCMGCoarseSpaceType *ctype)
1632: {
1633: PC_MG *mg = (PC_MG *)pc->data;
1635: PetscFunctionBegin;
1636: *ctype = mg->coarseSpaceType;
1637: PetscFunctionReturn(PETSC_SUCCESS);
1638: }
1640: static PetscErrorCode PCMGSetAdaptCR_MG(PC pc, PetscBool cr)
1641: {
1642: PC_MG *mg = (PC_MG *)pc->data;
1644: PetscFunctionBegin;
1645: mg->compatibleRelaxation = cr;
1646: PetscFunctionReturn(PETSC_SUCCESS);
1647: }
1649: static PetscErrorCode PCMGGetAdaptCR_MG(PC pc, PetscBool *cr)
1650: {
1651: PC_MG *mg = (PC_MG *)pc->data;
1653: PetscFunctionBegin;
1654: *cr = mg->compatibleRelaxation;
1655: PetscFunctionReturn(PETSC_SUCCESS);
1656: }
1658: /*@
1659: PCMGSetAdaptCoarseSpaceType - Set the type of adaptive coarse space. Adapts or creates the interpolator based upon a vector space which should be accurately
1660: captured by the next coarser mesh, and thus accurately interpolated.
1662: Logically Collective
1664: Input Parameters:
1665: + pc - the multigrid context
1666: - ctype - the type of coarse space
1668: Options Database Keys:
1669: + -pc_mg_adapt_interp_n nmodes - The number of modes to use
1670: - -pc_mg_adapt_interp_coarse_space (none|polynomial|harmonic|eigenvector|generalized_eigenvector|gdsw) - The type of coarse space to use
1672: Level: intermediate
1674: Note:
1675: Requires a `DM` with specific functionality be attached to the `PC`.
1677: Developer Notes:
1678: The options database key `-pc_mg_adapt_interp_coarse_space` should not have interp in it since this function does not have interp in it.
1680: The options database key `-pc_mg_adapt_interp_n` has no functional call equivalent.
1682: .seealso: [](ch_ksp), `PCMG`, `PCMGCoarseSpaceType`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetGalerkin()`, `PCMGSetAdaptInterpolation()`, `DM`,
1683: `PCMG_ADAPT_NONE`, `PCMG_ADAPT_POLYNOMIAL`, `PCMG_ADAPT_HARMONIC`, `PCMG_ADAPT_EIGENVECTOR`, `PCMG_ADAPT_GENERALIZED_EIGENVECTOR`,
1684: `PCMG_ADAPT_GDSW`
1685: @*/
1686: PetscErrorCode PCMGSetAdaptCoarseSpaceType(PC pc, PCMGCoarseSpaceType ctype)
1687: {
1688: PetscFunctionBegin;
1691: PetscTryMethod(pc, "PCMGSetAdaptCoarseSpaceType_C", (PC, PCMGCoarseSpaceType), (pc, ctype));
1692: PetscFunctionReturn(PETSC_SUCCESS);
1693: }
1695: /*@
1696: PCMGGetAdaptCoarseSpaceType - Get the type of adaptive coarse space.
1698: Not Collective
1700: Input Parameter:
1701: . pc - the multigrid context
1703: Output Parameter:
1704: . ctype - the type of coarse space
1706: Level: intermediate
1708: .seealso: [](ch_ksp), `PCMG`, `PCMGCoarseSpaceType`, `PCMGSetAdaptCoarseSpaceType()`, `PCMGSetGalerkin()`, `PCMGSetAdaptInterpolation()`
1709: @*/
1710: PetscErrorCode PCMGGetAdaptCoarseSpaceType(PC pc, PCMGCoarseSpaceType *ctype)
1711: {
1712: PetscFunctionBegin;
1714: PetscAssertPointer(ctype, 2);
1715: PetscUseMethod(pc, "PCMGGetAdaptCoarseSpaceType_C", (PC, PCMGCoarseSpaceType *), (pc, ctype));
1716: PetscFunctionReturn(PETSC_SUCCESS);
1717: }
1719: /*@
1720: PCMGSetAdaptInterpolation - Adapt the interpolator based upon a vector space which should be accurately captured by the next coarser mesh, and thus accurately interpolated.
1722: Logically Collective
1724: Input Parameters:
1725: + pc - the multigrid context
1726: - adapt - flag for adaptation of the interpolator
1728: Level: intermediate
1730: Note:
1731: This routine should never be used, rather call `PCMGSetAdaptCoarseSpaceType()`
1733: .seealso: [](ch_ksp), `PCMG`, `PCMGGetAdaptInterpolation()`, `PCMGSetGalerkin()`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetAdaptCoarseSpaceType()`
1734: @*/
1735: PetscErrorCode PCMGSetAdaptInterpolation(PC pc, PetscBool adapt)
1736: {
1737: PetscFunctionBegin;
1739: PetscTryMethod(pc, "PCMGSetAdaptInterpolation_C", (PC, PetscBool), (pc, adapt));
1740: PetscFunctionReturn(PETSC_SUCCESS);
1741: }
1743: /*@
1744: PCMGGetAdaptInterpolation - Get the flag to adapt the interpolator based upon a vector space which should be accurately captured by the next coarser mesh,
1745: and thus accurately interpolated.
1747: Not Collective
1749: Input Parameter:
1750: . pc - the multigrid context
1752: Output Parameter:
1753: . adapt - flag for adaptation of the interpolator
1755: Level: intermediate
1757: Note:
1758: This routine should never be used, rather call `PCMGGetAdaptCoarseSpaceType()`
1760: .seealso: [](ch_ksp), `PCMG`, `PCMGSetAdaptInterpolation()`, `PCMGSetGalerkin()`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetAdaptCoarseSpaceType()`
1761: @*/
1762: PetscErrorCode PCMGGetAdaptInterpolation(PC pc, PetscBool *adapt)
1763: {
1764: PetscFunctionBegin;
1766: PetscAssertPointer(adapt, 2);
1767: PetscUseMethod(pc, "PCMGGetAdaptInterpolation_C", (PC, PetscBool *), (pc, adapt));
1768: PetscFunctionReturn(PETSC_SUCCESS);
1769: }
1771: /*@
1772: PCMGSetAdaptCR - Monitor the coarse space quality using an auxiliary solve with compatible relaxation.
1774: Logically Collective
1776: Input Parameters:
1777: + pc - the multigrid context
1778: - cr - flag for compatible relaxation
1780: Options Database Key:
1781: . -pc_mg_adapt_cr (true|false) - Turn on compatible relaxation
1783: Level: intermediate
1785: .seealso: [](ch_ksp), `PCMG`, `PCMGGetAdaptCR()`, `PCMGSetAdaptInterpolation()`, `PCMGSetGalerkin()`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetAdaptCoarseSpaceType()`
1786: @*/
1787: PetscErrorCode PCMGSetAdaptCR(PC pc, PetscBool cr)
1788: {
1789: PetscFunctionBegin;
1791: PetscTryMethod(pc, "PCMGSetAdaptCR_C", (PC, PetscBool), (pc, cr));
1792: PetscFunctionReturn(PETSC_SUCCESS);
1793: }
1795: /*@
1796: PCMGGetAdaptCR - Get the flag to monitor coarse space quality using an auxiliary solve with compatible relaxation.
1798: Not Collective
1800: Input Parameter:
1801: . pc - the multigrid context
1803: Output Parameter:
1804: . cr - flag for compatible relaxaion
1806: Level: intermediate
1808: .seealso: [](ch_ksp), `PCMGSetAdaptCR()`, `PCMGGetAdaptInterpolation()`, `PCMGSetGalerkin()`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetAdaptCoarseSpaceType()`
1809: @*/
1810: PetscErrorCode PCMGGetAdaptCR(PC pc, PetscBool *cr)
1811: {
1812: PetscFunctionBegin;
1814: PetscAssertPointer(cr, 2);
1815: PetscUseMethod(pc, "PCMGGetAdaptCR_C", (PC, PetscBool *), (pc, cr));
1816: PetscFunctionReturn(PETSC_SUCCESS);
1817: }
1819: /*@
1820: PCMGSetNumberSmooth - Sets the number of pre and post-smoothing steps to use
1821: on all levels. Use `PCMGDistinctSmoothUp()` to create separate up and down smoothers if you want different numbers of
1822: pre- and post-smoothing steps.
1824: Logically Collective
1826: Input Parameters:
1827: + pc - the multigrid context
1828: - n - the number of smoothing steps
1830: Options Database Key:
1831: . -mg_levels_ksp_max_it n - Sets number of pre and post-smoothing steps
1833: Level: advanced
1835: Note:
1836: This does not set a value on the coarsest grid, since we assume that there is no separate smooth up on the coarsest grid.
1838: .seealso: [](ch_ksp), `PCMG`, `PCMGSetDistinctSmoothUp()`
1839: @*/
1840: PetscErrorCode PCMGSetNumberSmooth(PC pc, PetscInt n)
1841: {
1842: PC_MG *mg = (PC_MG *)pc->data;
1843: PC_MG_Levels **mglevels = mg->levels;
1844: PetscInt levels;
1846: PetscFunctionBegin;
1849: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Must set MG levels with PCMGSetLevels() before calling");
1850: levels = mglevels[0]->levels;
1852: for (PetscInt i = 1; i < levels; i++) {
1853: PetscCall(KSPSetTolerances(mglevels[i]->smoothu, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, n));
1854: PetscCall(KSPSetTolerances(mglevels[i]->smoothd, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, n));
1855: mg->default_smoothu = n;
1856: mg->default_smoothd = n;
1857: }
1858: PetscFunctionReturn(PETSC_SUCCESS);
1859: }
1861: /*@
1862: PCMGSetDistinctSmoothUp - sets the up (post) smoother to be a separate `KSP` from the down (pre) smoother on all levels
1863: and adds the suffix _up to the options name
1865: Logically Collective
1867: Input Parameter:
1868: . pc - the preconditioner context
1870: Options Database Key:
1871: . -pc_mg_distinct_smoothup (true|false) - use distinct smoothing objects
1873: Level: advanced
1875: Note:
1876: This does not set a value on the coarsest grid, since we assume that there is no separate smooth up on the coarsest grid.
1878: .seealso: [](ch_ksp), `PCMG`, `PCMGSetNumberSmooth()`
1879: @*/
1880: PetscErrorCode PCMGSetDistinctSmoothUp(PC pc)
1881: {
1882: PC_MG *mg = (PC_MG *)pc->data;
1883: PC_MG_Levels **mglevels = mg->levels;
1884: PetscInt levels;
1885: KSP subksp;
1887: PetscFunctionBegin;
1889: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Must set MG levels with PCMGSetLevels() before calling");
1890: levels = mglevels[0]->levels;
1892: for (PetscInt i = 1; i < levels; i++) {
1893: const char *prefix = NULL;
1894: /* make sure smoother up and down are different */
1895: PetscCall(PCMGGetSmootherUp(pc, i, &subksp));
1896: PetscCall(KSPGetOptionsPrefix(mglevels[i]->smoothd, &prefix));
1897: PetscCall(KSPSetOptionsPrefix(subksp, prefix));
1898: PetscCall(KSPAppendOptionsPrefix(subksp, "up_"));
1899: }
1900: PetscFunctionReturn(PETSC_SUCCESS);
1901: }
1903: /* No new matrices are created, and the coarse operator matrices are the references to the original ones */
1904: static PetscErrorCode PCGetInterpolations_MG(PC pc, PetscInt *num_levels, Mat *interpolations[])
1905: {
1906: PC_MG *mg = (PC_MG *)pc->data;
1907: PC_MG_Levels **mglevels = mg->levels;
1908: Mat *mat;
1910: PetscFunctionBegin;
1911: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
1912: PetscCall(PetscMalloc1(mg->nlevels, &mat));
1913: for (PetscInt l = 1; l < mg->nlevels; l++) {
1914: mat[l - 1] = mglevels[l]->interpolate;
1915: PetscCall(PetscObjectReference((PetscObject)mat[l - 1]));
1916: }
1917: *num_levels = mg->nlevels;
1918: *interpolations = mat;
1919: PetscFunctionReturn(PETSC_SUCCESS);
1920: }
1922: /* No new matrices are created, and the coarse operator matrices are the references to the original ones */
1923: static PetscErrorCode PCGetCoarseOperators_MG(PC pc, PetscInt *num_levels, Mat *coarseOperators[])
1924: {
1925: PC_MG *mg = (PC_MG *)pc->data;
1926: PC_MG_Levels **mglevels = mg->levels;
1927: Mat *mat;
1929: PetscFunctionBegin;
1930: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
1931: PetscCall(PetscMalloc1(mg->nlevels, &mat));
1932: for (PetscInt l = 0; l < mg->nlevels - 1; l++) {
1933: PetscCall(KSPGetOperators(mglevels[l]->smoothd, NULL, &mat[l]));
1934: PetscCall(PetscObjectReference((PetscObject)mat[l]));
1935: }
1936: *num_levels = mg->nlevels;
1937: *coarseOperators = mat;
1938: PetscFunctionReturn(PETSC_SUCCESS);
1939: }
1941: /*@
1942: PCMGRegisterCoarseSpaceConstructor - Adds a method to the `PCMG` package for coarse space construction.
1944: Not Collective, No Fortran Support
1946: Input Parameters:
1947: + name - name of the constructor
1948: - function - constructor routine, see `PCMGCoarseSpaceConstructorFn`
1950: Level: advanced
1952: Developer Notes:
1953: This is used by `PCMG_ADAPT_EIGENVECTOR` and `PCMG_ADAPT_GENERALIZED_EIGENVECTOR` to utilize the BAMG package
1955: `PCMGSetAdaptCoarseSpaceType()` and `PCMGCoarseSpaceType` should be refactored to use the standard PETSc registration of
1956: types as strings instead of incorrectly using an enum for `PCMGCoarseSpaceType` and thus
1957: requiring this ad hoc nonstandard registration process for BAMG.
1959: .seealso: [](ch_ksp), `PCMGCoarseSpaceConstructorFn`, `PCMG`, `PCMGGetCoarseSpaceConstructor()`, `PCRegister()`,
1960: `PCMGSetAdaptCoarseSpaceType()`, `PCMG_ADAPT_EIGENVECTOR`, `PCMG_ADAPT_GENERALIZED_EIGENVECTOR`
1961: @*/
1962: PetscErrorCode PCMGRegisterCoarseSpaceConstructor(const char name[], PCMGCoarseSpaceConstructorFn *function)
1963: {
1964: PetscFunctionBegin;
1965: PetscCall(PCInitializePackage());
1966: PetscCall(PetscFunctionListAdd(&PCMGCoarseList, name, function));
1967: PetscFunctionReturn(PETSC_SUCCESS);
1968: }
1970: /*@
1971: PCMGGetCoarseSpaceConstructor - Returns the given coarse space construction method.
1973: Not Collective, No Fortran Support
1975: Input Parameter:
1976: . name - name of the constructor
1978: Output Parameter:
1979: . function - constructor routine
1981: Level: advanced
1983: .seealso: [](ch_ksp), `PCMGCoarseSpaceConstructorFn`, `PCMG`, `PCMGRegisterCoarseSpaceConstructor()`, `PCRegister()`
1984: @*/
1985: PetscErrorCode PCMGGetCoarseSpaceConstructor(const char name[], PCMGCoarseSpaceConstructorFn **function)
1986: {
1987: PetscFunctionBegin;
1988: PetscCall(PetscFunctionListFind(PCMGCoarseList, name, function));
1989: PetscFunctionReturn(PETSC_SUCCESS);
1990: }
1992: /*MC
1993: PCMG - Use multigrid preconditioning. This preconditioner requires you provide additional information about the restriction/interpolation
1994: operators using `PCMGSetInterpolation()` and/or `PCMGSetRestriction()`(and possibly the coarser grid matrices) or a `DM` that can provide such information.
1996: Options Database Keys:
1997: + -pc_mg_levels nlevels - number of levels including finest
1998: . -pc_mg_cycle_type (v|w) - provide the cycle desired
1999: . -pc_mg_type (additive|multiplicative|full|kaskade) - multiplicative is the default
2000: . -pc_mg_log - log information about time spent on each level of the solver
2001: . -pc_mg_distinct_smoothup - configure up (after interpolation) and down (before restriction) smoothers separately (with different options prefixes)
2002: . -pc_mg_galerkin (both|pmat|mat|none) - use the Galerkin process to compute coarser operators, i.e., $A_{coarse} = R A_{fine} R^T$
2003: . -pc_mg_multiplicative_cycles ncycles - number of cycles to use as the preconditioner (defaults to 1)
2004: . -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices
2005: to a `PETSCVIEWERSOCKET` for reading from MATLAB.
2006: - -pc_mg_dump_binary - dumps the matrices for each level and the restriction/interpolation matrices
2007: to the binary output file called binaryoutput
2009: Level: intermediate
2011: Notes:
2012: `PCMG` provides a general framework for implementing multigrid methods. Use `PCGAMG` for PETSc's algebraic multigrid preconditioner, `PCHYPRE` for hypre's
2013: BoomerAMG algebraic multigrid, and `PCML` for Trilinos's ML preconditioner. `PCAMGX` provides access to NVIDIA's AmgX algebraic multigrid.
2015: If you use `KSPSetDM()` (or `SNESSetDM()` or `TSSetDM()`) with an appropriate `DM`, such as `DMDA`, then `PCMG` will use the geometric information
2016: from the `DM` to generate appropriate restriction and interpolation information and construct a geometric multigrid.
2018: If you do not provide an appropriate `DM` and do not provide restriction or interpolation operators with `PCMGSetInterpolation()` and/or `PCMGSetRestriction()`,
2019: then `PCMG` will run multigrid with only a single level (so not really multigrid).
2021: The Krylov solver (if any) and preconditioner (smoother) and their parameters are controlled from the options database with the standard
2022: options database keywords prefixed with `-mg_levels_` to affect all the levels but the coarsest, which is controlled with `-mg_coarse_`,
2023: and the finest where `-mg_fine_` can override `-mg_levels_`. One can set different preconditioners etc on specific levels with the prefix
2024: `-mg_levels_n_` where `n` is the level number (zero being the coarse level. For example
2025: .vb
2026: -mg_levels_ksp_type gmres -mg_levels_pc_type bjacobi -mg_coarse_pc_type svd -mg_levels_2_pc_type sor
2027: .ve
2028: These options also work for controlling the smoothers etc inside `PCGAMG`
2030: If one uses a Krylov method such `KSPGMRES` or `KSPCG` as the smoother than one must use `KSPFGMRES`, `KSPGCR`, or `KSPRICHARDSON` as the outer Krylov method
2032: When run with a single level the smoother options are used on that level NOT the coarse grid solver options
2034: When run with `KSPRICHARDSON` the convergence test changes slightly if monitor is turned on. The iteration count may change slightly. This
2035: is because without monitoring the residual norm is computed WITHIN each multigrid cycle on the finest level after the pre-smoothing
2036: (because the residual has just been computed for the multigrid algorithm and is hence available for free) while with monitoring the
2037: residual is computed at the end of each cycle.
2039: .seealso: [](sec_mg), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCMGType`, `PCEXOTIC`, `PCGAMG`, `PCML`, `PCHYPRE`, `PCAIR`,
2040: `PCMGSetLevels()`, `PCMGGetLevels()`, `PCMGSetType()`, `PCMGSetCycleType()`,
2041: `PCMGSetDistinctSmoothUp()`, `PCMGGetCoarseSolve()`, `PCMGSetResidual()`, `PCMGSetInterpolation()`,
2042: `PCMGSetRestriction()`, `PCMGGetSmoother()`, `PCMGGetSmootherUp()`, `PCMGGetSmootherDown()`,
2043: `PCMGSetCycleTypeOnLevel()`, `PCMGSetRhs()`, `PCMGSetX()`, `PCMGSetR()`,
2044: `PCMGSetAdaptCR()`, `PCMGGetAdaptInterpolation()`, `PCMGSetGalerkin()`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetAdaptCoarseSpaceType()`
2045: M*/
2047: PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc)
2048: {
2049: PC_MG *mg;
2051: PetscFunctionBegin;
2052: PetscCall(PetscNew(&mg));
2053: pc->data = mg;
2054: mg->nlevels = -1;
2055: mg->am = PC_MG_MULTIPLICATIVE;
2056: mg->galerkin = PC_MG_GALERKIN_NONE;
2057: mg->adaptInterpolation = PETSC_FALSE;
2058: mg->Nc = -1;
2059: mg->eigenvalue = -1;
2061: PetscObjectParameterSetDefault(pc, useAmat, PETSC_TRUE);
2063: pc->ops->apply = PCApply_MG;
2064: pc->ops->applytranspose = PCApplyTranspose_MG;
2065: pc->ops->matapply = PCMatApply_MG;
2066: pc->ops->matapplytranspose = PCMatApplyTranspose_MG;
2067: pc->ops->setup = PCSetUp_MG;
2068: pc->ops->reset = PCReset_MG;
2069: pc->ops->destroy = PCDestroy_MG;
2070: pc->ops->setfromoptions = PCSetFromOptions_MG;
2071: pc->ops->view = PCView_MG;
2073: PetscCall(PetscObjectComposedDataRegister(&mg->eigenvalue));
2074: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetGalerkin_C", PCMGSetGalerkin_MG));
2075: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetReusePreconditioner_C", PCSetReusePreconditioner_MG));
2076: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetLevels_C", PCMGGetLevels_MG));
2077: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetLevels_C", PCMGSetLevels_MG));
2078: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", PCGetInterpolations_MG));
2079: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", PCGetCoarseOperators_MG));
2080: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptInterpolation_C", PCMGSetAdaptInterpolation_MG));
2081: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptInterpolation_C", PCMGGetAdaptInterpolation_MG));
2082: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptCR_C", PCMGSetAdaptCR_MG));
2083: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptCR_C", PCMGGetAdaptCR_MG));
2084: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptCoarseSpaceType_C", PCMGSetAdaptCoarseSpaceType_MG));
2085: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptCoarseSpaceType_C", PCMGGetAdaptCoarseSpaceType_MG));
2086: PetscFunctionReturn(PETSC_SUCCESS);
2087: }