Actual source code: mg.c
petsc-3.9.4 2018-09-11
2: /*
3: Defines the multigrid preconditioner interface.
4: */
5: #include <petsc/private/pcmgimpl.h>
6: #include <petscdm.h>
7: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC,PetscBool*);
9: PetscErrorCode PCMGMCycle_Private(PC pc,PC_MG_Levels **mglevelsin,PCRichardsonConvergedReason *reason)
10: {
11: PC_MG *mg = (PC_MG*)pc->data;
12: PC_MG_Levels *mgc,*mglevels = *mglevelsin;
14: PetscInt cycles = (mglevels->level == 1) ? 1 : (PetscInt) mglevels->cycles;
15: PC subpc;
16: PCFailedReason pcreason;
19: if (mglevels->eventsmoothsolve) {PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);}
20: KSPSolve(mglevels->smoothd,mglevels->b,mglevels->x); /* pre-smooth */
21: KSPGetPC(mglevels->smoothd,&subpc);
22: PCGetSetUpFailedReason(subpc,&pcreason);
23: if (pcreason) {
24: pc->failedreason = PC_SUBPC_ERROR;
25: }
26: if (mglevels->eventsmoothsolve) {PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);}
27: if (mglevels->level) { /* not the coarsest grid */
28: if (mglevels->eventresidual) {PetscLogEventBegin(mglevels->eventresidual,0,0,0,0);}
29: (*mglevels->residual)(mglevels->A,mglevels->b,mglevels->x,mglevels->r);
30: if (mglevels->eventresidual) {PetscLogEventEnd(mglevels->eventresidual,0,0,0,0);}
32: /* if on finest level and have convergence criteria set */
33: if (mglevels->level == mglevels->levels-1 && mg->ttol && reason) {
34: PetscReal rnorm;
35: VecNorm(mglevels->r,NORM_2,&rnorm);
36: if (rnorm <= mg->ttol) {
37: if (rnorm < mg->abstol) {
38: *reason = PCRICHARDSON_CONVERGED_ATOL;
39: PetscInfo2(pc,"Linear solver has converged. Residual norm %g is less than absolute tolerance %g\n",(double)rnorm,(double)mg->abstol);
40: } else {
41: *reason = PCRICHARDSON_CONVERGED_RTOL;
42: PetscInfo2(pc,"Linear solver has converged. Residual norm %g is less than relative tolerance times initial residual norm %g\n",(double)rnorm,(double)mg->ttol);
43: }
44: return(0);
45: }
46: }
48: mgc = *(mglevelsin - 1);
49: if (mglevels->eventinterprestrict) {PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);}
50: MatRestrict(mglevels->restrct,mglevels->r,mgc->b);
51: if (mglevels->eventinterprestrict) {PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);}
52: VecSet(mgc->x,0.0);
53: while (cycles--) {
54: PCMGMCycle_Private(pc,mglevelsin-1,reason);
55: }
56: if (mglevels->eventinterprestrict) {PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);}
57: MatInterpolateAdd(mglevels->interpolate,mgc->x,mglevels->x,mglevels->x);
58: if (mglevels->eventinterprestrict) {PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);}
59: if (mglevels->eventsmoothsolve) {PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);}
60: KSPSolve(mglevels->smoothu,mglevels->b,mglevels->x); /* post smooth */
61: if (mglevels->eventsmoothsolve) {PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);}
62: }
63: return(0);
64: }
66: 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)
67: {
68: PC_MG *mg = (PC_MG*)pc->data;
69: PC_MG_Levels **mglevels = mg->levels;
71: PC tpc;
72: PetscBool changeu,changed;
73: PetscInt levels = mglevels[0]->levels,i;
76: /* When the DM is supplying the matrix then it will not exist until here */
77: for (i=0; i<levels; i++) {
78: if (!mglevels[i]->A) {
79: KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);
80: PetscObjectReference((PetscObject)mglevels[i]->A);
81: }
82: }
84: KSPGetPC(mglevels[levels-1]->smoothd,&tpc);
85: PCPreSolveChangeRHS(tpc,&changed);
86: KSPGetPC(mglevels[levels-1]->smoothu,&tpc);
87: PCPreSolveChangeRHS(tpc,&changeu);
88: if (!changed && !changeu) {
89: VecDestroy(&mglevels[levels-1]->b);
90: mglevels[levels-1]->b = b;
91: } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */
92: if (!mglevels[levels-1]->b) {
93: Vec *vec;
95: KSPCreateVecs(mglevels[levels-1]->smoothd,1,&vec,0,NULL);
96: mglevels[levels-1]->b = *vec;
97: }
98: VecCopy(b,mglevels[levels-1]->b);
99: }
100: mglevels[levels-1]->x = x;
102: mg->rtol = rtol;
103: mg->abstol = abstol;
104: mg->dtol = dtol;
105: if (rtol) {
106: /* compute initial residual norm for relative convergence test */
107: PetscReal rnorm;
108: if (zeroguess) {
109: VecNorm(b,NORM_2,&rnorm);
110: } else {
111: (*mglevels[levels-1]->residual)(mglevels[levels-1]->A,b,x,w);
112: VecNorm(w,NORM_2,&rnorm);
113: }
114: mg->ttol = PetscMax(rtol*rnorm,abstol);
115: } else if (abstol) mg->ttol = abstol;
116: else mg->ttol = 0.0;
118: /* since smoother is applied to full system, not just residual we need to make sure that smoothers don't
119: stop prematurely due to small residual */
120: for (i=1; i<levels; i++) {
121: KSPSetTolerances(mglevels[i]->smoothu,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
122: if (mglevels[i]->smoothu != mglevels[i]->smoothd) {
123: /* For Richardson the initial guess is nonzero since it is solving in each cycle the original system not just applying as a preconditioner */
124: KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);
125: KSPSetTolerances(mglevels[i]->smoothd,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
126: }
127: }
129: *reason = (PCRichardsonConvergedReason)0;
130: for (i=0; i<its; i++) {
131: PCMGMCycle_Private(pc,mglevels+levels-1,reason);
132: if (*reason) break;
133: }
134: if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
135: *outits = i;
136: if (!changed && !changeu) mglevels[levels-1]->b = NULL;
137: return(0);
138: }
140: PetscErrorCode PCReset_MG(PC pc)
141: {
142: PC_MG *mg = (PC_MG*)pc->data;
143: PC_MG_Levels **mglevels = mg->levels;
145: PetscInt i,n;
148: if (mglevels) {
149: n = mglevels[0]->levels;
150: for (i=0; i<n-1; i++) {
151: VecDestroy(&mglevels[i+1]->r);
152: VecDestroy(&mglevels[i]->b);
153: VecDestroy(&mglevels[i]->x);
154: MatDestroy(&mglevels[i+1]->restrct);
155: MatDestroy(&mglevels[i+1]->interpolate);
156: MatDestroy(&mglevels[i+1]->inject);
157: VecDestroy(&mglevels[i+1]->rscale);
158: }
159: /* this is not null only if the smoother on the finest level
160: changes the rhs during PreSolve */
161: VecDestroy(&mglevels[n-1]->b);
163: for (i=0; i<n; i++) {
164: MatDestroy(&mglevels[i]->A);
165: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
166: KSPReset(mglevels[i]->smoothd);
167: }
168: KSPReset(mglevels[i]->smoothu);
169: }
170: }
171: return(0);
172: }
174: /*@C
175: PCMGSetLevels - Sets the number of levels to use with MG.
176: Must be called before any other MG routine.
178: Logically Collective on PC
180: Input Parameters:
181: + pc - the preconditioner context
182: . levels - the number of levels
183: - comms - optional communicators for each level; this is to allow solving the coarser problems
184: on smaller sets of processors.
186: Level: intermediate
188: Notes:
189: If the number of levels is one then the multigrid uses the -mg_levels prefix
190: for setting the level options rather than the -mg_coarse prefix.
192: .keywords: MG, set, levels, multigrid
194: .seealso: PCMGSetType(), PCMGGetLevels()
195: @*/
196: PetscErrorCode PCMGSetLevels(PC pc,PetscInt levels,MPI_Comm *comms)
197: {
199: PC_MG *mg = (PC_MG*)pc->data;
200: MPI_Comm comm;
201: PC_MG_Levels **mglevels = mg->levels;
202: PCMGType mgtype = mg->am;
203: PetscInt mgctype = (PetscInt) PC_MG_CYCLE_V;
204: PetscInt i;
205: PetscMPIInt size;
206: const char *prefix;
207: PC ipc;
208: PetscInt n;
213: if (mg->nlevels == levels) return(0);
214: PetscObjectGetComm((PetscObject)pc,&comm);
215: if (mglevels) {
216: mgctype = mglevels[0]->cycles;
217: /* changing the number of levels so free up the previous stuff */
218: PCReset_MG(pc);
219: n = mglevels[0]->levels;
220: for (i=0; i<n; i++) {
221: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
222: KSPDestroy(&mglevels[i]->smoothd);
223: }
224: KSPDestroy(&mglevels[i]->smoothu);
225: PetscFree(mglevels[i]);
226: }
227: PetscFree(mg->levels);
228: }
230: mg->nlevels = levels;
232: PetscMalloc1(levels,&mglevels);
233: PetscLogObjectMemory((PetscObject)pc,levels*(sizeof(PC_MG*)));
235: PCGetOptionsPrefix(pc,&prefix);
237: mg->stageApply = 0;
238: for (i=0; i<levels; i++) {
239: PetscNewLog(pc,&mglevels[i]);
241: mglevels[i]->level = i;
242: mglevels[i]->levels = levels;
243: mglevels[i]->cycles = mgctype;
244: mg->default_smoothu = 2;
245: mg->default_smoothd = 2;
246: mglevels[i]->eventsmoothsetup = 0;
247: mglevels[i]->eventsmoothsolve = 0;
248: mglevels[i]->eventresidual = 0;
249: mglevels[i]->eventinterprestrict = 0;
251: if (comms) comm = comms[i];
252: KSPCreate(comm,&mglevels[i]->smoothd);
253: KSPSetErrorIfNotConverged(mglevels[i]->smoothd,pc->erroriffailure);
254: PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd,(PetscObject)pc,levels-i);
255: KSPSetOptionsPrefix(mglevels[i]->smoothd,prefix);
256: PetscObjectComposedDataSetInt((PetscObject) mglevels[i]->smoothd, PetscMGLevelId, mglevels[i]->level);
257: if (i || levels == 1) {
258: char tprefix[128];
260: KSPSetType(mglevels[i]->smoothd,KSPCHEBYSHEV);
261: KSPSetConvergenceTest(mglevels[i]->smoothd,KSPConvergedSkip,NULL,NULL);
262: KSPSetNormType(mglevels[i]->smoothd,KSP_NORM_NONE);
263: KSPGetPC(mglevels[i]->smoothd,&ipc);
264: PCSetType(ipc,PCSOR);
265: KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, mg->default_smoothd);
267: sprintf(tprefix,"mg_levels_%d_",(int)i);
268: KSPAppendOptionsPrefix(mglevels[i]->smoothd,tprefix);
269: } else {
270: KSPAppendOptionsPrefix(mglevels[0]->smoothd,"mg_coarse_");
272: /* coarse solve is (redundant) LU by default; set shifttype NONZERO to avoid annoying zero-pivot in LU preconditioner */
273: KSPSetType(mglevels[0]->smoothd,KSPPREONLY);
274: KSPGetPC(mglevels[0]->smoothd,&ipc);
275: MPI_Comm_size(comm,&size);
276: if (size > 1) {
277: PCSetType(ipc,PCREDUNDANT);
278: } else {
279: PCSetType(ipc,PCLU);
280: }
281: PCFactorSetShiftType(ipc,MAT_SHIFT_INBLOCKS);
282: }
283: PetscLogObjectParent((PetscObject)pc,(PetscObject)mglevels[i]->smoothd);
285: mglevels[i]->smoothu = mglevels[i]->smoothd;
286: mg->rtol = 0.0;
287: mg->abstol = 0.0;
288: mg->dtol = 0.0;
289: mg->ttol = 0.0;
290: mg->cyclesperpcapply = 1;
291: }
292: mg->levels = mglevels;
293: PCMGSetType(pc,mgtype);
294: return(0);
295: }
298: PetscErrorCode PCDestroy_MG(PC pc)
299: {
301: PC_MG *mg = (PC_MG*)pc->data;
302: PC_MG_Levels **mglevels = mg->levels;
303: PetscInt i,n;
306: PCReset_MG(pc);
307: if (mglevels) {
308: n = mglevels[0]->levels;
309: for (i=0; i<n; i++) {
310: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
311: KSPDestroy(&mglevels[i]->smoothd);
312: }
313: KSPDestroy(&mglevels[i]->smoothu);
314: PetscFree(mglevels[i]);
315: }
316: PetscFree(mg->levels);
317: }
318: PetscFree(pc->data);
319: return(0);
320: }
324: extern PetscErrorCode PCMGACycle_Private(PC,PC_MG_Levels**);
325: extern PetscErrorCode PCMGFCycle_Private(PC,PC_MG_Levels**);
326: extern PetscErrorCode PCMGKCycle_Private(PC,PC_MG_Levels**);
328: /*
329: PCApply_MG - Runs either an additive, multiplicative, Kaskadic
330: or full cycle of multigrid.
332: Note:
333: A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle().
334: */
335: static PetscErrorCode PCApply_MG(PC pc,Vec b,Vec x)
336: {
337: PC_MG *mg = (PC_MG*)pc->data;
338: PC_MG_Levels **mglevels = mg->levels;
340: PC tpc;
341: PetscInt levels = mglevels[0]->levels,i;
342: PetscBool changeu,changed;
345: if (mg->stageApply) {PetscLogStagePush(mg->stageApply);}
346: /* When the DM is supplying the matrix then it will not exist until here */
347: for (i=0; i<levels; i++) {
348: if (!mglevels[i]->A) {
349: KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);
350: PetscObjectReference((PetscObject)mglevels[i]->A);
351: }
352: }
354: KSPGetPC(mglevels[levels-1]->smoothd,&tpc);
355: PCPreSolveChangeRHS(tpc,&changed);
356: KSPGetPC(mglevels[levels-1]->smoothu,&tpc);
357: PCPreSolveChangeRHS(tpc,&changeu);
358: if (!changeu && !changed) {
359: VecDestroy(&mglevels[levels-1]->b);
360: mglevels[levels-1]->b = b;
361: } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */
362: if (!mglevels[levels-1]->b) {
363: Vec *vec;
365: KSPCreateVecs(mglevels[levels-1]->smoothd,1,&vec,0,NULL);
366: mglevels[levels-1]->b = *vec;
367: }
368: VecCopy(b,mglevels[levels-1]->b);
369: }
370: mglevels[levels-1]->x = x;
372: if (mg->am == PC_MG_MULTIPLICATIVE) {
373: VecSet(x,0.0);
374: for (i=0; i<mg->cyclesperpcapply; i++) {
375: PCMGMCycle_Private(pc,mglevels+levels-1,NULL);
376: }
377: } else if (mg->am == PC_MG_ADDITIVE) {
378: PCMGACycle_Private(pc,mglevels);
379: } else if (mg->am == PC_MG_KASKADE) {
380: PCMGKCycle_Private(pc,mglevels);
381: } else {
382: PCMGFCycle_Private(pc,mglevels);
383: }
384: if (mg->stageApply) {PetscLogStagePop();}
385: if (!changeu && !changed) mglevels[levels-1]->b = NULL;
386: return(0);
387: }
390: PetscErrorCode PCSetFromOptions_MG(PetscOptionItems *PetscOptionsObject,PC pc)
391: {
392: PetscErrorCode ierr;
393: PetscInt levels,cycles;
394: PetscBool flg;
395: PC_MG *mg = (PC_MG*)pc->data;
396: PC_MG_Levels **mglevels;
397: PCMGType mgtype;
398: PCMGCycleType mgctype;
399: PCMGGalerkinType gtype;
402: levels = PetscMax(mg->nlevels,1);
403: PetscOptionsHead(PetscOptionsObject,"Multigrid options");
404: PetscOptionsInt("-pc_mg_levels","Number of Levels","PCMGSetLevels",levels,&levels,&flg);
405: if (!flg && !mg->levels && pc->dm) {
406: DMGetRefineLevel(pc->dm,&levels);
407: levels++;
408: mg->usedmfornumberoflevels = PETSC_TRUE;
409: }
410: PCMGSetLevels(pc,levels,NULL);
411: mglevels = mg->levels;
413: mgctype = (PCMGCycleType) mglevels[0]->cycles;
414: PetscOptionsEnum("-pc_mg_cycle_type","V cycle or for W-cycle","PCMGSetCycleType",PCMGCycleTypes,(PetscEnum)mgctype,(PetscEnum*)&mgctype,&flg);
415: if (flg) {
416: PCMGSetCycleType(pc,mgctype);
417: }
418: gtype = mg->galerkin;
419: PetscOptionsEnum("-pc_mg_galerkin","Use Galerkin process to compute coarser operators","PCMGSetGalerkin",PCMGGalerkinTypes,(PetscEnum)gtype,(PetscEnum*)>ype,&flg);
420: if (flg) {
421: PCMGSetGalerkin(pc,gtype);
422: }
423: flg = PETSC_FALSE;
424: PetscOptionsBool("-pc_mg_distinct_smoothup","Create seperate smoothup KSP and append the prefix _up","PCMGSetDistinctSmoothUp",PETSC_FALSE,&flg,NULL);
425: if (flg) {
426: PCMGSetDistinctSmoothUp(pc);
427: }
428: mgtype = mg->am;
429: PetscOptionsEnum("-pc_mg_type","Multigrid type","PCMGSetType",PCMGTypes,(PetscEnum)mgtype,(PetscEnum*)&mgtype,&flg);
430: if (flg) {
431: PCMGSetType(pc,mgtype);
432: }
433: if (mg->am == PC_MG_MULTIPLICATIVE) {
434: PetscOptionsInt("-pc_mg_multiplicative_cycles","Number of cycles for each preconditioner step","PCMGMultiplicativeSetCycles",mg->cyclesperpcapply,&cycles,&flg);
435: if (flg) {
436: PCMGMultiplicativeSetCycles(pc,cycles);
437: }
438: }
439: flg = PETSC_FALSE;
440: PetscOptionsBool("-pc_mg_log","Log times for each multigrid level","None",flg,&flg,NULL);
441: if (flg) {
442: PetscInt i;
443: char eventname[128];
445: levels = mglevels[0]->levels;
446: for (i=0; i<levels; i++) {
447: sprintf(eventname,"MGSetup Level %d",(int)i);
448: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsetup);
449: sprintf(eventname,"MGSmooth Level %d",(int)i);
450: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsolve);
451: if (i) {
452: sprintf(eventname,"MGResid Level %d",(int)i);
453: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventresidual);
454: sprintf(eventname,"MGInterp Level %d",(int)i);
455: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventinterprestrict);
456: }
457: }
459: #if defined(PETSC_USE_LOG)
460: {
461: const char *sname = "MG Apply";
462: PetscStageLog stageLog;
463: PetscInt st;
465: PetscLogGetStageLog(&stageLog);
466: for (st = 0; st < stageLog->numStages; ++st) {
467: PetscBool same;
469: PetscStrcmp(stageLog->stageInfo[st].name, sname, &same);
470: if (same) mg->stageApply = st;
471: }
472: if (!mg->stageApply) {
473: PetscLogStageRegister(sname, &mg->stageApply);
474: }
475: }
476: #endif
477: }
478: PetscOptionsTail();
479: return(0);
480: }
482: const char *const PCMGTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","PCMGType","PC_MG",0};
483: const char *const PCMGCycleTypes[] = {"invalid","v","w","PCMGCycleType","PC_MG_CYCLE",0};
484: const char *const PCMGGalerkinTypes[] = {"both","pmat","mat","none","external","PCMGGalerkinType","PC_MG_GALERKIN",0};
486: #include <petscdraw.h>
487: PetscErrorCode PCView_MG(PC pc,PetscViewer viewer)
488: {
489: PC_MG *mg = (PC_MG*)pc->data;
490: PC_MG_Levels **mglevels = mg->levels;
492: PetscInt levels = mglevels ? mglevels[0]->levels : 0,i;
493: PetscBool iascii,isbinary,isdraw;
496: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
497: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
498: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
499: if (iascii) {
500: const char *cyclename = levels ? (mglevels[0]->cycles == PC_MG_CYCLE_V ? "v" : "w") : "unknown";
501: PetscViewerASCIIPrintf(viewer," type is %s, levels=%D cycles=%s\n", PCMGTypes[mg->am],levels,cyclename);
502: if (mg->am == PC_MG_MULTIPLICATIVE) {
503: PetscViewerASCIIPrintf(viewer," Cycles per PCApply=%d\n",mg->cyclesperpcapply);
504: }
505: if (mg->galerkin == PC_MG_GALERKIN_BOTH) {
506: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices\n");
507: } else if (mg->galerkin == PC_MG_GALERKIN_PMAT) {
508: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices for pmat\n");
509: } else if (mg->galerkin == PC_MG_GALERKIN_MAT) {
510: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices for mat\n");
511: } else if (mg->galerkin == PC_MG_GALERKIN_EXTERNAL) {
512: PetscViewerASCIIPrintf(viewer," Using externally compute Galerkin coarse grid matrices\n");
513: } else {
514: PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid matrices\n");
515: }
516: if (mg->view){
517: (*mg->view)(pc,viewer);
518: }
519: for (i=0; i<levels; i++) {
520: if (!i) {
521: PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level -------------------------------\n",i);
522: } else {
523: PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
524: }
525: PetscViewerASCIIPushTab(viewer);
526: KSPView(mglevels[i]->smoothd,viewer);
527: PetscViewerASCIIPopTab(viewer);
528: if (i && mglevels[i]->smoothd == mglevels[i]->smoothu) {
529: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
530: } else if (i) {
531: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
532: PetscViewerASCIIPushTab(viewer);
533: KSPView(mglevels[i]->smoothu,viewer);
534: PetscViewerASCIIPopTab(viewer);
535: }
536: }
537: } else if (isbinary) {
538: for (i=levels-1; i>=0; i--) {
539: KSPView(mglevels[i]->smoothd,viewer);
540: if (i && mglevels[i]->smoothd != mglevels[i]->smoothu) {
541: KSPView(mglevels[i]->smoothu,viewer);
542: }
543: }
544: } else if (isdraw) {
545: PetscDraw draw;
546: PetscReal x,w,y,bottom,th;
547: PetscViewerDrawGetDraw(viewer,0,&draw);
548: PetscDrawGetCurrentPoint(draw,&x,&y);
549: PetscDrawStringGetSize(draw,NULL,&th);
550: bottom = y - th;
551: for (i=levels-1; i>=0; i--) {
552: if (!mglevels[i]->smoothu || (mglevels[i]->smoothu == mglevels[i]->smoothd)) {
553: PetscDrawPushCurrentPoint(draw,x,bottom);
554: KSPView(mglevels[i]->smoothd,viewer);
555: PetscDrawPopCurrentPoint(draw);
556: } else {
557: w = 0.5*PetscMin(1.0-x,x);
558: PetscDrawPushCurrentPoint(draw,x+w,bottom);
559: KSPView(mglevels[i]->smoothd,viewer);
560: PetscDrawPopCurrentPoint(draw);
561: PetscDrawPushCurrentPoint(draw,x-w,bottom);
562: KSPView(mglevels[i]->smoothu,viewer);
563: PetscDrawPopCurrentPoint(draw);
564: }
565: PetscDrawGetBoundingBox(draw,NULL,&bottom,NULL,NULL);
566: bottom -= th;
567: }
568: }
569: return(0);
570: }
572: #include <petsc/private/dmimpl.h>
573: #include <petsc/private/kspimpl.h>
575: /*
576: Calls setup for the KSP on each level
577: */
578: PetscErrorCode PCSetUp_MG(PC pc)
579: {
580: PC_MG *mg = (PC_MG*)pc->data;
581: PC_MG_Levels **mglevels = mg->levels;
583: PetscInt i,n;
584: PC cpc;
585: PetscBool dump = PETSC_FALSE,opsset,use_amat,missinginterpolate = PETSC_FALSE;
586: Mat dA,dB;
587: Vec tvec;
588: DM *dms;
589: PetscViewer viewer = 0;
590: PetscBool dAeqdB = PETSC_FALSE, needRestricts = PETSC_FALSE;
593: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels with PCMGSetLevels() before setting up");
594: n = mglevels[0]->levels;
595: /* FIX: Move this to PCSetFromOptions_MG? */
596: if (mg->usedmfornumberoflevels) {
597: PetscInt levels;
598: DMGetRefineLevel(pc->dm,&levels);
599: levels++;
600: if (levels > n) { /* the problem is now being solved on a finer grid */
601: PCMGSetLevels(pc,levels,NULL);
602: n = levels;
603: PCSetFromOptions(pc); /* it is bad to call this here, but otherwise will never be called for the new hierarchy */
604: mglevels = mg->levels;
605: }
606: }
607: KSPGetPC(mglevels[0]->smoothd,&cpc);
610: /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */
611: /* so use those from global PC */
612: /* Is this what we always want? What if user wants to keep old one? */
613: KSPGetOperatorsSet(mglevels[n-1]->smoothd,NULL,&opsset);
614: if (opsset) {
615: Mat mmat;
616: KSPGetOperators(mglevels[n-1]->smoothd,NULL,&mmat);
617: if (mmat == pc->pmat) opsset = PETSC_FALSE;
618: }
620: if (!opsset) {
621: PCGetUseAmat(pc,&use_amat);
622: if(use_amat){
623: PetscInfo(pc,"Using outer operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
624: KSPSetOperators(mglevels[n-1]->smoothd,pc->mat,pc->pmat);
625: }
626: else {
627: PetscInfo(pc,"Using matrix (pmat) operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
628: KSPSetOperators(mglevels[n-1]->smoothd,pc->pmat,pc->pmat);
629: }
630: }
632: for (i=n-1; i>0; i--) {
633: if (!(mglevels[i]->interpolate || mglevels[i]->restrct)) {
634: missinginterpolate = PETSC_TRUE;
635: continue;
636: }
637: }
639: KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB);
640: if (dA == dB) dAeqdB = PETSC_TRUE;
641: if ((mg->galerkin == PC_MG_GALERKIN_NONE) || (((mg->galerkin == PC_MG_GALERKIN_PMAT) || (mg->galerkin == PC_MG_GALERKIN_MAT)) && !dAeqdB)) {
642: needRestricts = PETSC_TRUE; /* user must compute either mat, pmat, or both so must restrict x to coarser levels */
643: }
646: /*
647: Skipping if user has provided all interpolation/restriction needed (since DM might not be able to produce them (when coming from SNES/TS)
648: Skipping for galerkin==2 (externally managed hierarchy such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs?
649: */
650: if (missinginterpolate && pc->dm && mg->galerkin != PC_MG_GALERKIN_EXTERNAL && !pc->setupcalled) {
651: /* construct the interpolation from the DMs */
652: Mat p;
653: Vec rscale;
654: PetscMalloc1(n,&dms);
655: dms[n-1] = pc->dm;
656: /* Separately create them so we do not get DMKSP interference between levels */
657: for (i=n-2; i>-1; i--) {DMCoarsen(dms[i+1],MPI_COMM_NULL,&dms[i]);}
658: /*
659: Force the mat type of coarse level operator to be AIJ because usually we want to use LU for coarse level.
660: Notice that it can be overwritten by -mat_type because KSPSetUp() reads command line options.
661: But it is safe to use -dm_mat_type.
663: The mat type should not be hardcoded like this, we need to find a better way.
664: DMSetMatType(dms[0],MATAIJ);
665: */
666: for (i=n-2; i>-1; i--) {
667: DMKSP kdm;
668: PetscBool dmhasrestrict, dmhasinject;
669: KSPSetDM(mglevels[i]->smoothd,dms[i]);
670: if (!needRestricts) {KSPSetDMActive(mglevels[i]->smoothd,PETSC_FALSE);}
671: DMGetDMKSPWrite(dms[i],&kdm);
672: /* Ugly hack so that the next KSPSetUp() will use the RHS that we set. A better fix is to change dmActive to take
673: * a bitwise OR of computing the matrix, RHS, and initial iterate. */
674: kdm->ops->computerhs = NULL;
675: kdm->rhsctx = NULL;
676: if (!mglevels[i+1]->interpolate) {
677: DMCreateInterpolation(dms[i],dms[i+1],&p,&rscale);
678: PCMGSetInterpolation(pc,i+1,p);
679: if (rscale) {PCMGSetRScale(pc,i+1,rscale);}
680: VecDestroy(&rscale);
681: MatDestroy(&p);
682: }
683: DMHasCreateRestriction(dms[i],&dmhasrestrict);
684: if (dmhasrestrict && !mglevels[i+1]->restrct){
685: DMCreateRestriction(dms[i],dms[i+1],&p);
686: PCMGSetRestriction(pc,i+1,p);
687: MatDestroy(&p);
688: }
689: DMHasCreateInjection(dms[i],&dmhasinject);
690: if (dmhasinject && !mglevels[i+1]->inject){
691: DMCreateInjection(dms[i],dms[i+1],&p);
692: PCMGSetInjection(pc,i+1,p);
693: MatDestroy(&p);
694: }
695: }
697: for (i=n-2; i>-1; i--) {DMDestroy(&dms[i]);}
698: PetscFree(dms);
699: }
701: if (pc->dm && !pc->setupcalled) {
702: /* finest smoother also gets DM but it is not active, independent of whether galerkin==PC_MG_GALERKIN_EXTERNAL */
703: KSPSetDM(mglevels[n-1]->smoothd,pc->dm);
704: KSPSetDMActive(mglevels[n-1]->smoothd,PETSC_FALSE);
705: }
707: if (mg->galerkin < PC_MG_GALERKIN_NONE) {
708: Mat A,B;
709: PetscBool doA = PETSC_FALSE,doB = PETSC_FALSE;
710: MatReuse reuse = MAT_INITIAL_MATRIX;
712: if ((mg->galerkin == PC_MG_GALERKIN_PMAT) || (mg->galerkin == PC_MG_GALERKIN_BOTH)) doB = PETSC_TRUE;
713: if ((mg->galerkin == PC_MG_GALERKIN_MAT) || ((mg->galerkin == PC_MG_GALERKIN_BOTH) && (dA != dB))) doA = PETSC_TRUE;
714: if (pc->setupcalled) reuse = MAT_REUSE_MATRIX;
715: for (i=n-2; i>-1; i--) {
716: if (!mglevels[i+1]->restrct && !mglevels[i+1]->interpolate) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must provide interpolation or restriction for each MG level except level 0");
717: if (!mglevels[i+1]->interpolate) {
718: PCMGSetInterpolation(pc,i+1,mglevels[i+1]->restrct);
719: }
720: if (!mglevels[i+1]->restrct) {
721: PCMGSetRestriction(pc,i+1,mglevels[i+1]->interpolate);
722: }
723: if (reuse == MAT_REUSE_MATRIX) {
724: KSPGetOperators(mglevels[i]->smoothd,&A,&B);
725: }
726: if (doA) {
727: MatGalerkin(mglevels[i+1]->restrct,dA,mglevels[i+1]->interpolate,reuse,1.0,&A);
728: }
729: if (doB) {
730: MatGalerkin(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,reuse,1.0,&B);
731: }
732: /* the management of the PetscObjectReference() and PetscObjecDereference() below is rather delicate */
733: if (!doA && dAeqdB) {
734: if (reuse == MAT_INITIAL_MATRIX) {PetscObjectReference((PetscObject)B);}
735: A = B;
736: } else if (!doA && reuse == MAT_INITIAL_MATRIX ) {
737: KSPGetOperators(mglevels[i]->smoothd,&A,NULL);
738: PetscObjectReference((PetscObject)A);
739: }
740: if (!doB && dAeqdB) {
741: if (reuse == MAT_INITIAL_MATRIX) {PetscObjectReference((PetscObject)A);}
742: B = A;
743: } else if (!doB && reuse == MAT_INITIAL_MATRIX) {
744: KSPGetOperators(mglevels[i]->smoothd,NULL,&B);
745: PetscObjectReference((PetscObject)B);
746: }
747: if (reuse == MAT_INITIAL_MATRIX) {
748: KSPSetOperators(mglevels[i]->smoothd,A,B);
749: PetscObjectDereference((PetscObject)A);
750: PetscObjectDereference((PetscObject)B);
751: }
752: dA = A;
753: dB = B;
754: }
755: }
756: if (needRestricts && pc->dm && pc->dm->x) {
757: /* need to restrict Jacobian location to coarser meshes for evaluation */
758: for (i=n-2; i>-1; i--) {
759: Mat R;
760: Vec rscale;
761: if (!mglevels[i]->smoothd->dm->x) {
762: Vec *vecs;
763: KSPCreateVecs(mglevels[i]->smoothd,1,&vecs,0,NULL);
764: mglevels[i]->smoothd->dm->x = vecs[0];
765: PetscFree(vecs);
766: }
767: PCMGGetRestriction(pc,i+1,&R);
768: PCMGGetRScale(pc,i+1,&rscale);
769: MatRestrict(R,mglevels[i+1]->smoothd->dm->x,mglevels[i]->smoothd->dm->x);
770: VecPointwiseMult(mglevels[i]->smoothd->dm->x,mglevels[i]->smoothd->dm->x,rscale);
771: }
772: }
773: if (needRestricts && pc->dm) {
774: for (i=n-2; i>=0; i--) {
775: DM dmfine,dmcoarse;
776: Mat Restrict,Inject;
777: Vec rscale;
778: KSPGetDM(mglevels[i+1]->smoothd,&dmfine);
779: KSPGetDM(mglevels[i]->smoothd,&dmcoarse);
780: PCMGGetRestriction(pc,i+1,&Restrict);
781: PCMGGetRScale(pc,i+1,&rscale);
782: PCMGGetInjection(pc,i+1,&Inject);
783: DMRestrict(dmfine,Restrict,rscale,Inject,dmcoarse);
784: }
785: }
787: if (!pc->setupcalled) {
788: for (i=0; i<n; i++) {
789: KSPSetFromOptions(mglevels[i]->smoothd);
790: }
791: for (i=1; i<n; i++) {
792: if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) {
793: KSPSetFromOptions(mglevels[i]->smoothu);
794: }
795: }
796: /* insure that if either interpolation or restriction is set the other other one is set */
797: for (i=1; i<n; i++) {
798: PCMGGetInterpolation(pc,i,NULL);
799: PCMGGetRestriction(pc,i,NULL);
800: }
801: for (i=0; i<n-1; i++) {
802: if (!mglevels[i]->b) {
803: Vec *vec;
804: KSPCreateVecs(mglevels[i]->smoothd,1,&vec,0,NULL);
805: PCMGSetRhs(pc,i,*vec);
806: VecDestroy(vec);
807: PetscFree(vec);
808: }
809: if (!mglevels[i]->r && i) {
810: VecDuplicate(mglevels[i]->b,&tvec);
811: PCMGSetR(pc,i,tvec);
812: VecDestroy(&tvec);
813: }
814: if (!mglevels[i]->x) {
815: VecDuplicate(mglevels[i]->b,&tvec);
816: PCMGSetX(pc,i,tvec);
817: VecDestroy(&tvec);
818: }
819: }
820: if (n != 1 && !mglevels[n-1]->r) {
821: /* PCMGSetR() on the finest level if user did not supply it */
822: Vec *vec;
823: KSPCreateVecs(mglevels[n-1]->smoothd,1,&vec,0,NULL);
824: PCMGSetR(pc,n-1,*vec);
825: VecDestroy(vec);
826: PetscFree(vec);
827: }
828: }
830: if (pc->dm) {
831: /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */
832: for (i=0; i<n-1; i++) {
833: if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
834: }
835: }
837: for (i=1; i<n; i++) {
838: if (mglevels[i]->smoothu == mglevels[i]->smoothd || mg->am == PC_MG_FULL || mg->am == PC_MG_KASKADE || mg->cyclesperpcapply > 1){
839: /* if doing only down then initial guess is zero */
840: KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);
841: }
842: if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
843: KSPSetUp(mglevels[i]->smoothd);
844: if (mglevels[i]->smoothd->reason == KSP_DIVERGED_PCSETUP_FAILED) {
845: pc->failedreason = PC_SUBPC_ERROR;
846: }
847: if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
848: if (!mglevels[i]->residual) {
849: Mat mat;
850: KSPGetOperators(mglevels[i]->smoothd,&mat,NULL);
851: PCMGSetResidual(pc,i,PCMGResidualDefault,mat);
852: }
853: }
854: for (i=1; i<n; i++) {
855: if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) {
856: Mat downmat,downpmat;
858: /* check if operators have been set for up, if not use down operators to set them */
859: KSPGetOperatorsSet(mglevels[i]->smoothu,&opsset,NULL);
860: if (!opsset) {
861: KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat);
862: KSPSetOperators(mglevels[i]->smoothu,downmat,downpmat);
863: }
865: KSPSetInitialGuessNonzero(mglevels[i]->smoothu,PETSC_TRUE);
866: if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
867: KSPSetUp(mglevels[i]->smoothu);
868: if (mglevels[i]->smoothu->reason == KSP_DIVERGED_PCSETUP_FAILED) {
869: pc->failedreason = PC_SUBPC_ERROR;
870: }
871: if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
872: }
873: }
875: if (mglevels[0]->eventsmoothsetup) {PetscLogEventBegin(mglevels[0]->eventsmoothsetup,0,0,0,0);}
876: KSPSetUp(mglevels[0]->smoothd);
877: if (mglevels[0]->smoothd->reason == KSP_DIVERGED_PCSETUP_FAILED) {
878: pc->failedreason = PC_SUBPC_ERROR;
879: }
880: if (mglevels[0]->eventsmoothsetup) {PetscLogEventEnd(mglevels[0]->eventsmoothsetup,0,0,0,0);}
882: /*
883: Dump the interpolation/restriction matrices plus the
884: Jacobian/stiffness on each level. This allows MATLAB users to
885: easily check if the Galerkin condition A_c = R A_f R^T is satisfied.
887: Only support one or the other at the same time.
888: */
889: #if defined(PETSC_USE_SOCKET_VIEWER)
890: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_mg_dump_matlab",&dump,NULL);
891: if (dump) viewer = PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc));
892: dump = PETSC_FALSE;
893: #endif
894: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_mg_dump_binary",&dump,NULL);
895: if (dump) viewer = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc));
897: if (viewer) {
898: for (i=1; i<n; i++) {
899: MatView(mglevels[i]->restrct,viewer);
900: }
901: for (i=0; i<n; i++) {
902: KSPGetPC(mglevels[i]->smoothd,&pc);
903: MatView(pc->mat,viewer);
904: }
905: }
906: return(0);
907: }
909: /* -------------------------------------------------------------------------------------*/
911: /*@
912: PCMGGetLevels - Gets the number of levels to use with MG.
914: Not Collective
916: Input Parameter:
917: . pc - the preconditioner context
919: Output parameter:
920: . levels - the number of levels
922: Level: advanced
924: .keywords: MG, get, levels, multigrid
926: .seealso: PCMGSetLevels()
927: @*/
928: PetscErrorCode PCMGGetLevels(PC pc,PetscInt *levels)
929: {
930: PC_MG *mg = (PC_MG*)pc->data;
935: *levels = mg->nlevels;
936: return(0);
937: }
939: /*@
940: PCMGSetType - Determines the form of multigrid to use:
941: multiplicative, additive, full, or the Kaskade algorithm.
943: Logically Collective on PC
945: Input Parameters:
946: + pc - the preconditioner context
947: - form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,
948: PC_MG_FULL, PC_MG_KASKADE
950: Options Database Key:
951: . -pc_mg_type <form> - Sets <form>, one of multiplicative,
952: additive, full, kaskade
954: Level: advanced
956: .keywords: MG, set, method, multiplicative, additive, full, Kaskade, multigrid
958: .seealso: PCMGSetLevels()
959: @*/
960: PetscErrorCode PCMGSetType(PC pc,PCMGType form)
961: {
962: PC_MG *mg = (PC_MG*)pc->data;
967: mg->am = form;
968: if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
969: else pc->ops->applyrichardson = NULL;
970: return(0);
971: }
973: /*@
974: PCMGGetType - Determines the form of multigrid to use:
975: multiplicative, additive, full, or the Kaskade algorithm.
977: Logically Collective on PC
979: Input Parameter:
980: . pc - the preconditioner context
982: Output Parameter:
983: . type - one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,PC_MG_FULL, PC_MG_KASKADE
986: Level: advanced
988: .keywords: MG, set, method, multiplicative, additive, full, Kaskade, multigrid
990: .seealso: PCMGSetLevels()
991: @*/
992: PetscErrorCode PCMGGetType(PC pc,PCMGType *type)
993: {
994: PC_MG *mg = (PC_MG*)pc->data;
998: *type = mg->am;
999: return(0);
1000: }
1002: /*@
1003: PCMGSetCycleType - Sets the type cycles to use. Use PCMGSetCycleTypeOnLevel() for more
1004: complicated cycling.
1006: Logically Collective on PC
1008: Input Parameters:
1009: + pc - the multigrid context
1010: - n - either PC_MG_CYCLE_V or PC_MG_CYCLE_W
1012: Options Database Key:
1013: . -pc_mg_cycle_type <v,w> - provide the cycle desired
1015: Level: advanced
1017: .keywords: MG, set, cycles, V-cycle, W-cycle, multigrid
1019: .seealso: PCMGSetCycleTypeOnLevel()
1020: @*/
1021: PetscErrorCode PCMGSetCycleType(PC pc,PCMGCycleType n)
1022: {
1023: PC_MG *mg = (PC_MG*)pc->data;
1024: PC_MG_Levels **mglevels = mg->levels;
1025: PetscInt i,levels;
1030: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling");
1031: levels = mglevels[0]->levels;
1032: for (i=0; i<levels; i++) mglevels[i]->cycles = n;
1033: return(0);
1034: }
1036: /*@
1037: PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step
1038: of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used
1040: Logically Collective on PC
1042: Input Parameters:
1043: + pc - the multigrid context
1044: - n - number of cycles (default is 1)
1046: Options Database Key:
1047: . -pc_mg_multiplicative_cycles n
1049: Level: advanced
1051: Notes: This is not associated with setting a v or w cycle, that is set with PCMGSetCycleType()
1053: .keywords: MG, set, cycles, V-cycle, W-cycle, multigrid
1055: .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType()
1056: @*/
1057: PetscErrorCode PCMGMultiplicativeSetCycles(PC pc,PetscInt n)
1058: {
1059: PC_MG *mg = (PC_MG*)pc->data;
1064: mg->cyclesperpcapply = n;
1065: return(0);
1066: }
1068: PetscErrorCode PCMGSetGalerkin_MG(PC pc,PCMGGalerkinType use)
1069: {
1070: PC_MG *mg = (PC_MG*)pc->data;
1073: mg->galerkin = use;
1074: return(0);
1075: }
1077: /*@
1078: PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the
1079: finest grid via the Galerkin process: A_i-1 = r_i * A_i * p_i
1081: Logically Collective on PC
1083: Input Parameters:
1084: + pc - the multigrid context
1085: - use - one of PC_MG_GALERKIN_BOTH,PC_MG_GALERKIN_PMAT,PC_MG_GALERKIN_MAT, or PC_MG_GALERKIN_NONE
1087: Options Database Key:
1088: . -pc_mg_galerkin <both,pmat,mat,none>
1090: Level: intermediate
1092: Notes: Some codes that use PCMG such as PCGAMG use Galerkin internally while constructing the hierarchy and thus do not
1093: use the PCMG construction of the coarser grids.
1095: .keywords: MG, set, Galerkin
1097: .seealso: PCMGGetGalerkin(), PCMGGalerkinType
1099: @*/
1100: PetscErrorCode PCMGSetGalerkin(PC pc,PCMGGalerkinType use)
1101: {
1106: PetscTryMethod(pc,"PCMGSetGalerkin_C",(PC,PCMGGalerkinType),(pc,use));
1107: return(0);
1108: }
1110: /*@
1111: PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e.
1112: A_i-1 = r_i * A_i * p_i
1114: Not Collective
1116: Input Parameter:
1117: . pc - the multigrid context
1119: Output Parameter:
1120: . galerkin - one of PC_MG_GALERKIN_BOTH,PC_MG_GALERKIN_PMAT,PC_MG_GALERKIN_MAT, PC_MG_GALERKIN_NONE, or PC_MG_GALERKIN_EXTERNAL
1122: Level: intermediate
1124: .keywords: MG, set, Galerkin
1126: .seealso: PCMGSetGalerkin(), PCMGGalerkinType
1128: @*/
1129: PetscErrorCode PCMGGetGalerkin(PC pc,PCMGGalerkinType *galerkin)
1130: {
1131: PC_MG *mg = (PC_MG*)pc->data;
1135: *galerkin = mg->galerkin;
1136: return(0);
1137: }
1139: /*@
1140: PCMGSetNumberSmooth - Sets the number of pre and post-smoothing steps to use
1141: on all levels. Use PCMGDistinctSmoothUp() to create separate up and down smoothers if you want different numbers of
1142: pre- and post-smoothing steps.
1144: Logically Collective on PC
1146: Input Parameters:
1147: + mg - the multigrid context
1148: - n - the number of smoothing steps
1150: Options Database Key:
1151: + -mg_levels_ksp_max_it <n> - Sets number of pre and post-smoothing steps
1153: Level: advanced
1155: Notes: this does not set a value on the coarsest grid, since we assume that
1156: there is no separate smooth up on the coarsest grid.
1158: .keywords: MG, smooth, up, post-smoothing, steps, multigrid
1160: .seealso: PCMGSetDistinctSmoothUp()
1161: @*/
1162: PetscErrorCode PCMGSetNumberSmooth(PC pc,PetscInt n)
1163: {
1164: PC_MG *mg = (PC_MG*)pc->data;
1165: PC_MG_Levels **mglevels = mg->levels;
1167: PetscInt i,levels;
1172: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling");
1173: levels = mglevels[0]->levels;
1175: for (i=1; i<levels; i++) {
1176: KSPSetTolerances(mglevels[i]->smoothu,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
1177: KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
1178: mg->default_smoothu = n;
1179: mg->default_smoothd = n;
1180: }
1181: return(0);
1182: }
1184: /*@
1185: PCMGSetDistinctSmoothUp - sets the up (post) smoother to be a seperate KSP from the down (pre) smoother on all levels
1186: and adds the suffix _up to the options name
1188: Logically Collective on PC
1190: Input Parameters:
1191: . pc - the preconditioner context
1193: Options Database Key:
1194: . -pc_mg_distinct_smoothup
1196: Level: advanced
1198: Notes: this does not set a value on the coarsest grid, since we assume that
1199: there is no separate smooth up on the coarsest grid.
1201: .keywords: MG, smooth, up, post-smoothing, steps, multigrid
1203: .seealso: PCMGSetNumberSmooth()
1204: @*/
1205: PetscErrorCode PCMGSetDistinctSmoothUp(PC pc)
1206: {
1207: PC_MG *mg = (PC_MG*)pc->data;
1208: PC_MG_Levels **mglevels = mg->levels;
1210: PetscInt i,levels;
1211: KSP subksp;
1215: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling");
1216: levels = mglevels[0]->levels;
1218: for (i=1; i<levels; i++) {
1219: const char *prefix = NULL;
1220: /* make sure smoother up and down are different */
1221: PCMGGetSmootherUp(pc,i,&subksp);
1222: KSPGetOptionsPrefix(mglevels[i]->smoothd,&prefix);
1223: KSPSetOptionsPrefix(subksp,prefix);
1224: KSPAppendOptionsPrefix(subksp,"up_");
1225: }
1226: return(0);
1227: }
1229: /* ----------------------------------------------------------------------------------------*/
1231: /*MC
1232: PCMG - Use multigrid preconditioning. This preconditioner requires you provide additional
1233: information about the coarser grid matrices and restriction/interpolation operators.
1235: Options Database Keys:
1236: + -pc_mg_levels <nlevels> - number of levels including finest
1237: . -pc_mg_cycle_type <v,w> - provide the cycle desired
1238: . -pc_mg_type <additive,multiplicative,full,kaskade> - multiplicative is the default
1239: . -pc_mg_log - log information about time spent on each level of the solver
1240: . -pc_mg_distinct_smoothup - configure up (after interpolation) and down (before restriction) smoothers separately (with different options prefixes)
1241: . -pc_mg_galerkin <both,pmat,mat,none> - use Galerkin process to compute coarser operators, i.e. Acoarse = R A R'
1242: . -pc_mg_multiplicative_cycles - number of cycles to use as the preconditioner (defaults to 1)
1243: . -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices
1244: to the Socket viewer for reading from MATLAB.
1245: - -pc_mg_dump_binary - dumps the matrices for each level and the restriction/interpolation matrices
1246: to the binary output file called binaryoutput
1248: Notes: If one uses a Krylov method such GMRES or CG as the smoother than one must use KSPFGMRES, KSPGCG, or KSPRICHARDSON as the outer Krylov method
1250: When run with a single level the smoother options are used on that level NOT the coarse grid solver options
1252: When run with KSPRICHARDSON the convergence test changes slightly if monitor is turned on. The iteration count may change slightly. This
1253: is because without monitoring the residual norm is computed WITHIN each multigrid cycle on the finest level after the pre-smoothing
1254: (because the residual has just been computed for the multigrid algorithm and is hence available for free) while with monitoring the
1255: residual is computed at the end of each cycle.
1257: Level: intermediate
1259: Concepts: multigrid/multilevel
1261: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, PCEXOTIC, PCGAMG, PCML, PCHYPRE
1262: PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(),
1263: PCMGSetDistinctSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
1264: PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
1265: PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()
1266: M*/
1268: PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc)
1269: {
1270: PC_MG *mg;
1274: PetscNewLog(pc,&mg);
1275: pc->data = (void*)mg;
1276: mg->nlevels = -1;
1277: mg->am = PC_MG_MULTIPLICATIVE;
1278: mg->galerkin = PC_MG_GALERKIN_NONE;
1280: pc->useAmat = PETSC_TRUE;
1282: pc->ops->apply = PCApply_MG;
1283: pc->ops->setup = PCSetUp_MG;
1284: pc->ops->reset = PCReset_MG;
1285: pc->ops->destroy = PCDestroy_MG;
1286: pc->ops->setfromoptions = PCSetFromOptions_MG;
1287: pc->ops->view = PCView_MG;
1289: PetscObjectComposeFunction((PetscObject)pc,"PCMGSetGalerkin_C",PCMGSetGalerkin_MG);
1290: return(0);
1291: }