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