Actual source code: mg.c
petsc-3.10.5 2019-03-28
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: PetscFree(vec);
98: }
99: VecCopy(b,mglevels[levels-1]->b);
100: }
101: mglevels[levels-1]->x = x;
103: mg->rtol = rtol;
104: mg->abstol = abstol;
105: mg->dtol = dtol;
106: if (rtol) {
107: /* compute initial residual norm for relative convergence test */
108: PetscReal rnorm;
109: if (zeroguess) {
110: VecNorm(b,NORM_2,&rnorm);
111: } else {
112: (*mglevels[levels-1]->residual)(mglevels[levels-1]->A,b,x,w);
113: VecNorm(w,NORM_2,&rnorm);
114: }
115: mg->ttol = PetscMax(rtol*rnorm,abstol);
116: } else if (abstol) mg->ttol = abstol;
117: else mg->ttol = 0.0;
119: /* since smoother is applied to full system, not just residual we need to make sure that smoothers don't
120: stop prematurely due to small residual */
121: for (i=1; i<levels; i++) {
122: KSPSetTolerances(mglevels[i]->smoothu,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
123: if (mglevels[i]->smoothu != mglevels[i]->smoothd) {
124: /* For Richardson the initial guess is nonzero since it is solving in each cycle the original system not just applying as a preconditioner */
125: KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);
126: KSPSetTolerances(mglevels[i]->smoothd,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
127: }
128: }
130: *reason = (PCRichardsonConvergedReason)0;
131: for (i=0; i<its; i++) {
132: PCMGMCycle_Private(pc,mglevels+levels-1,reason);
133: if (*reason) break;
134: }
135: if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
136: *outits = i;
137: if (!changed && !changeu) mglevels[levels-1]->b = NULL;
138: return(0);
139: }
141: PetscErrorCode PCReset_MG(PC pc)
142: {
143: PC_MG *mg = (PC_MG*)pc->data;
144: PC_MG_Levels **mglevels = mg->levels;
146: PetscInt i,n;
149: if (mglevels) {
150: n = mglevels[0]->levels;
151: for (i=0; i<n-1; i++) {
152: VecDestroy(&mglevels[i+1]->r);
153: VecDestroy(&mglevels[i]->b);
154: VecDestroy(&mglevels[i]->x);
155: MatDestroy(&mglevels[i+1]->restrct);
156: MatDestroy(&mglevels[i+1]->interpolate);
157: MatDestroy(&mglevels[i+1]->inject);
158: VecDestroy(&mglevels[i+1]->rscale);
159: }
160: /* this is not null only if the smoother on the finest level
161: changes the rhs during PreSolve */
162: VecDestroy(&mglevels[n-1]->b);
164: for (i=0; i<n; i++) {
165: MatDestroy(&mglevels[i]->A);
166: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
167: KSPReset(mglevels[i]->smoothd);
168: }
169: KSPReset(mglevels[i]->smoothu);
170: }
171: }
172: return(0);
173: }
175: /*@C
176: PCMGSetLevels - Sets the number of levels to use with MG.
177: Must be called before any other MG routine.
179: Logically Collective on PC
181: Input Parameters:
182: + pc - the preconditioner context
183: . levels - the number of levels
184: - comms - optional communicators for each level; this is to allow solving the coarser problems
185: on smaller sets of processors.
187: Level: intermediate
189: Notes:
190: If the number of levels is one then the multigrid uses the -mg_levels prefix
191: for setting the level options rather than the -mg_coarse prefix.
193: .keywords: MG, set, levels, multigrid
195: .seealso: PCMGSetType(), PCMGGetLevels()
196: @*/
197: PetscErrorCode PCMGSetLevels(PC pc,PetscInt levels,MPI_Comm *comms)
198: {
200: PC_MG *mg = (PC_MG*)pc->data;
201: MPI_Comm comm;
202: PC_MG_Levels **mglevels = mg->levels;
203: PCMGType mgtype = mg->am;
204: PetscInt mgctype = (PetscInt) PC_MG_CYCLE_V;
205: PetscInt i;
206: PetscMPIInt size;
207: const char *prefix;
208: PC ipc;
209: PetscInt n;
214: if (mg->nlevels == levels) return(0);
215: PetscObjectGetComm((PetscObject)pc,&comm);
216: if (mglevels) {
217: mgctype = mglevels[0]->cycles;
218: /* changing the number of levels so free up the previous stuff */
219: PCReset_MG(pc);
220: n = mglevels[0]->levels;
221: for (i=0; i<n; i++) {
222: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
223: KSPDestroy(&mglevels[i]->smoothd);
224: }
225: KSPDestroy(&mglevels[i]->smoothu);
226: PetscFree(mglevels[i]);
227: }
228: PetscFree(mg->levels);
229: }
231: mg->nlevels = levels;
233: PetscMalloc1(levels,&mglevels);
234: PetscLogObjectMemory((PetscObject)pc,levels*(sizeof(PC_MG*)));
236: PCGetOptionsPrefix(pc,&prefix);
238: mg->stageApply = 0;
239: for (i=0; i<levels; i++) {
240: PetscNewLog(pc,&mglevels[i]);
242: mglevels[i]->level = i;
243: mglevels[i]->levels = levels;
244: mglevels[i]->cycles = mgctype;
245: mg->default_smoothu = 2;
246: mg->default_smoothd = 2;
247: mglevels[i]->eventsmoothsetup = 0;
248: mglevels[i]->eventsmoothsolve = 0;
249: mglevels[i]->eventresidual = 0;
250: mglevels[i]->eventinterprestrict = 0;
252: if (comms) comm = comms[i];
253: KSPCreate(comm,&mglevels[i]->smoothd);
254: KSPSetErrorIfNotConverged(mglevels[i]->smoothd,pc->erroriffailure);
255: PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd,(PetscObject)pc,levels-i);
256: KSPSetOptionsPrefix(mglevels[i]->smoothd,prefix);
257: PetscObjectComposedDataSetInt((PetscObject) mglevels[i]->smoothd, PetscMGLevelId, mglevels[i]->level);
258: if (i || levels == 1) {
259: char tprefix[128];
261: KSPSetType(mglevels[i]->smoothd,KSPCHEBYSHEV);
262: KSPSetConvergenceTest(mglevels[i]->smoothd,KSPConvergedSkip,NULL,NULL);
263: KSPSetNormType(mglevels[i]->smoothd,KSP_NORM_NONE);
264: KSPGetPC(mglevels[i]->smoothd,&ipc);
265: PCSetType(ipc,PCSOR);
266: KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, mg->default_smoothd);
268: sprintf(tprefix,"mg_levels_%d_",(int)i);
269: KSPAppendOptionsPrefix(mglevels[i]->smoothd,tprefix);
270: } else {
271: KSPAppendOptionsPrefix(mglevels[0]->smoothd,"mg_coarse_");
273: /* coarse solve is (redundant) LU by default; set shifttype NONZERO to avoid annoying zero-pivot in LU preconditioner */
274: KSPSetType(mglevels[0]->smoothd,KSPPREONLY);
275: KSPGetPC(mglevels[0]->smoothd,&ipc);
276: MPI_Comm_size(comm,&size);
277: if (size > 1) {
278: PCSetType(ipc,PCREDUNDANT);
279: } else {
280: PCSetType(ipc,PCLU);
281: }
282: PCFactorSetShiftType(ipc,MAT_SHIFT_INBLOCKS);
283: }
284: PetscLogObjectParent((PetscObject)pc,(PetscObject)mglevels[i]->smoothd);
286: mglevels[i]->smoothu = mglevels[i]->smoothd;
287: mg->rtol = 0.0;
288: mg->abstol = 0.0;
289: mg->dtol = 0.0;
290: mg->ttol = 0.0;
291: mg->cyclesperpcapply = 1;
292: }
293: mg->levels = mglevels;
294: PCMGSetType(pc,mgtype);
295: return(0);
296: }
299: PetscErrorCode PCDestroy_MG(PC pc)
300: {
302: PC_MG *mg = (PC_MG*)pc->data;
303: PC_MG_Levels **mglevels = mg->levels;
304: PetscInt i,n;
307: PCReset_MG(pc);
308: if (mglevels) {
309: n = mglevels[0]->levels;
310: for (i=0; i<n; i++) {
311: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
312: KSPDestroy(&mglevels[i]->smoothd);
313: }
314: KSPDestroy(&mglevels[i]->smoothu);
315: PetscFree(mglevels[i]);
316: }
317: PetscFree(mg->levels);
318: }
319: PetscFree(pc->data);
320: return(0);
321: }
325: extern PetscErrorCode PCMGACycle_Private(PC,PC_MG_Levels**);
326: extern PetscErrorCode PCMGFCycle_Private(PC,PC_MG_Levels**);
327: extern PetscErrorCode PCMGKCycle_Private(PC,PC_MG_Levels**);
329: /*
330: PCApply_MG - Runs either an additive, multiplicative, Kaskadic
331: or full cycle of multigrid.
333: Note:
334: A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle().
335: */
336: static PetscErrorCode PCApply_MG(PC pc,Vec b,Vec x)
337: {
338: PC_MG *mg = (PC_MG*)pc->data;
339: PC_MG_Levels **mglevels = mg->levels;
341: PC tpc;
342: PetscInt levels = mglevels[0]->levels,i;
343: PetscBool changeu,changed;
346: if (mg->stageApply) {PetscLogStagePush(mg->stageApply);}
347: /* When the DM is supplying the matrix then it will not exist until here */
348: for (i=0; i<levels; i++) {
349: if (!mglevels[i]->A) {
350: KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);
351: PetscObjectReference((PetscObject)mglevels[i]->A);
352: }
353: }
355: KSPGetPC(mglevels[levels-1]->smoothd,&tpc);
356: PCPreSolveChangeRHS(tpc,&changed);
357: KSPGetPC(mglevels[levels-1]->smoothu,&tpc);
358: PCPreSolveChangeRHS(tpc,&changeu);
359: if (!changeu && !changed) {
360: VecDestroy(&mglevels[levels-1]->b);
361: mglevels[levels-1]->b = b;
362: } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */
363: if (!mglevels[levels-1]->b) {
364: Vec *vec;
366: KSPCreateVecs(mglevels[levels-1]->smoothd,1,&vec,0,NULL);
367: mglevels[levels-1]->b = *vec;
368: PetscFree(vec);
369: }
370: VecCopy(b,mglevels[levels-1]->b);
371: }
372: mglevels[levels-1]->x = x;
374: if (mg->am == PC_MG_MULTIPLICATIVE) {
375: VecSet(x,0.0);
376: for (i=0; i<mg->cyclesperpcapply; i++) {
377: PCMGMCycle_Private(pc,mglevels+levels-1,NULL);
378: }
379: } else if (mg->am == PC_MG_ADDITIVE) {
380: PCMGACycle_Private(pc,mglevels);
381: } else if (mg->am == PC_MG_KASKADE) {
382: PCMGKCycle_Private(pc,mglevels);
383: } else {
384: PCMGFCycle_Private(pc,mglevels);
385: }
386: if (mg->stageApply) {PetscLogStagePop();}
387: if (!changeu && !changed) mglevels[levels-1]->b = NULL;
388: return(0);
389: }
392: PetscErrorCode PCSetFromOptions_MG(PetscOptionItems *PetscOptionsObject,PC pc)
393: {
394: PetscErrorCode ierr;
395: PetscInt levels,cycles;
396: PetscBool flg;
397: PC_MG *mg = (PC_MG*)pc->data;
398: PC_MG_Levels **mglevels;
399: PCMGType mgtype;
400: PCMGCycleType mgctype;
401: PCMGGalerkinType gtype;
404: levels = PetscMax(mg->nlevels,1);
405: PetscOptionsHead(PetscOptionsObject,"Multigrid options");
406: PetscOptionsInt("-pc_mg_levels","Number of Levels","PCMGSetLevels",levels,&levels,&flg);
407: if (!flg && !mg->levels && pc->dm) {
408: DMGetRefineLevel(pc->dm,&levels);
409: levels++;
410: mg->usedmfornumberoflevels = PETSC_TRUE;
411: }
412: PCMGSetLevels(pc,levels,NULL);
413: mglevels = mg->levels;
415: mgctype = (PCMGCycleType) mglevels[0]->cycles;
416: PetscOptionsEnum("-pc_mg_cycle_type","V cycle or for W-cycle","PCMGSetCycleType",PCMGCycleTypes,(PetscEnum)mgctype,(PetscEnum*)&mgctype,&flg);
417: if (flg) {
418: PCMGSetCycleType(pc,mgctype);
419: }
420: gtype = mg->galerkin;
421: PetscOptionsEnum("-pc_mg_galerkin","Use Galerkin process to compute coarser operators","PCMGSetGalerkin",PCMGGalerkinTypes,(PetscEnum)gtype,(PetscEnum*)>ype,&flg);
422: if (flg) {
423: PCMGSetGalerkin(pc,gtype);
424: }
425: flg = PETSC_FALSE;
426: PetscOptionsBool("-pc_mg_distinct_smoothup","Create seperate smoothup KSP and append the prefix _up","PCMGSetDistinctSmoothUp",PETSC_FALSE,&flg,NULL);
427: if (flg) {
428: PCMGSetDistinctSmoothUp(pc);
429: }
430: mgtype = mg->am;
431: PetscOptionsEnum("-pc_mg_type","Multigrid type","PCMGSetType",PCMGTypes,(PetscEnum)mgtype,(PetscEnum*)&mgtype,&flg);
432: if (flg) {
433: PCMGSetType(pc,mgtype);
434: }
435: if (mg->am == PC_MG_MULTIPLICATIVE) {
436: PetscOptionsInt("-pc_mg_multiplicative_cycles","Number of cycles for each preconditioner step","PCMGMultiplicativeSetCycles",mg->cyclesperpcapply,&cycles,&flg);
437: if (flg) {
438: PCMGMultiplicativeSetCycles(pc,cycles);
439: }
440: }
441: flg = PETSC_FALSE;
442: PetscOptionsBool("-pc_mg_log","Log times for each multigrid level","None",flg,&flg,NULL);
443: if (flg) {
444: PetscInt i;
445: char eventname[128];
447: levels = mglevels[0]->levels;
448: for (i=0; i<levels; i++) {
449: sprintf(eventname,"MGSetup Level %d",(int)i);
450: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsetup);
451: sprintf(eventname,"MGSmooth Level %d",(int)i);
452: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsolve);
453: if (i) {
454: sprintf(eventname,"MGResid Level %d",(int)i);
455: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventresidual);
456: sprintf(eventname,"MGInterp Level %d",(int)i);
457: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventinterprestrict);
458: }
459: }
461: #if defined(PETSC_USE_LOG)
462: {
463: const char *sname = "MG Apply";
464: PetscStageLog stageLog;
465: PetscInt st;
467: PetscLogGetStageLog(&stageLog);
468: for (st = 0; st < stageLog->numStages; ++st) {
469: PetscBool same;
471: PetscStrcmp(stageLog->stageInfo[st].name, sname, &same);
472: if (same) mg->stageApply = st;
473: }
474: if (!mg->stageApply) {
475: PetscLogStageRegister(sname, &mg->stageApply);
476: }
477: }
478: #endif
479: }
480: PetscOptionsTail();
481: return(0);
482: }
484: const char *const PCMGTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","PCMGType","PC_MG",0};
485: const char *const PCMGCycleTypes[] = {"invalid","v","w","PCMGCycleType","PC_MG_CYCLE",0};
486: const char *const PCMGGalerkinTypes[] = {"both","pmat","mat","none","external","PCMGGalerkinType","PC_MG_GALERKIN",0};
488: #include <petscdraw.h>
489: PetscErrorCode PCView_MG(PC pc,PetscViewer viewer)
490: {
491: PC_MG *mg = (PC_MG*)pc->data;
492: PC_MG_Levels **mglevels = mg->levels;
494: PetscInt levels = mglevels ? mglevels[0]->levels : 0,i;
495: PetscBool iascii,isbinary,isdraw;
498: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
499: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
500: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
501: if (iascii) {
502: const char *cyclename = levels ? (mglevels[0]->cycles == PC_MG_CYCLE_V ? "v" : "w") : "unknown";
503: PetscViewerASCIIPrintf(viewer," type is %s, levels=%D cycles=%s\n", PCMGTypes[mg->am],levels,cyclename);
504: if (mg->am == PC_MG_MULTIPLICATIVE) {
505: PetscViewerASCIIPrintf(viewer," Cycles per PCApply=%d\n",mg->cyclesperpcapply);
506: }
507: if (mg->galerkin == PC_MG_GALERKIN_BOTH) {
508: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices\n");
509: } else if (mg->galerkin == PC_MG_GALERKIN_PMAT) {
510: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices for pmat\n");
511: } else if (mg->galerkin == PC_MG_GALERKIN_MAT) {
512: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices for mat\n");
513: } else if (mg->galerkin == PC_MG_GALERKIN_EXTERNAL) {
514: PetscViewerASCIIPrintf(viewer," Using externally compute Galerkin coarse grid matrices\n");
515: } else {
516: PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid matrices\n");
517: }
518: if (mg->view){
519: (*mg->view)(pc,viewer);
520: }
521: for (i=0; i<levels; i++) {
522: if (!i) {
523: PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level -------------------------------\n",i);
524: } else {
525: PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
526: }
527: PetscViewerASCIIPushTab(viewer);
528: KSPView(mglevels[i]->smoothd,viewer);
529: PetscViewerASCIIPopTab(viewer);
530: if (i && mglevels[i]->smoothd == mglevels[i]->smoothu) {
531: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
532: } else if (i) {
533: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
534: PetscViewerASCIIPushTab(viewer);
535: KSPView(mglevels[i]->smoothu,viewer);
536: PetscViewerASCIIPopTab(viewer);
537: }
538: }
539: } else if (isbinary) {
540: for (i=levels-1; i>=0; i--) {
541: KSPView(mglevels[i]->smoothd,viewer);
542: if (i && mglevels[i]->smoothd != mglevels[i]->smoothu) {
543: KSPView(mglevels[i]->smoothu,viewer);
544: }
545: }
546: } else if (isdraw) {
547: PetscDraw draw;
548: PetscReal x,w,y,bottom,th;
549: PetscViewerDrawGetDraw(viewer,0,&draw);
550: PetscDrawGetCurrentPoint(draw,&x,&y);
551: PetscDrawStringGetSize(draw,NULL,&th);
552: bottom = y - th;
553: for (i=levels-1; i>=0; i--) {
554: if (!mglevels[i]->smoothu || (mglevels[i]->smoothu == mglevels[i]->smoothd)) {
555: PetscDrawPushCurrentPoint(draw,x,bottom);
556: KSPView(mglevels[i]->smoothd,viewer);
557: PetscDrawPopCurrentPoint(draw);
558: } else {
559: w = 0.5*PetscMin(1.0-x,x);
560: PetscDrawPushCurrentPoint(draw,x+w,bottom);
561: KSPView(mglevels[i]->smoothd,viewer);
562: PetscDrawPopCurrentPoint(draw);
563: PetscDrawPushCurrentPoint(draw,x-w,bottom);
564: KSPView(mglevels[i]->smoothu,viewer);
565: PetscDrawPopCurrentPoint(draw);
566: }
567: PetscDrawGetBoundingBox(draw,NULL,&bottom,NULL,NULL);
568: bottom -= th;
569: }
570: }
571: return(0);
572: }
574: #include <petsc/private/dmimpl.h>
575: #include <petsc/private/kspimpl.h>
577: /*
578: Calls setup for the KSP on each level
579: */
580: PetscErrorCode PCSetUp_MG(PC pc)
581: {
582: PC_MG *mg = (PC_MG*)pc->data;
583: PC_MG_Levels **mglevels = mg->levels;
585: PetscInt i,n;
586: PC cpc;
587: PetscBool dump = PETSC_FALSE,opsset,use_amat,missinginterpolate = PETSC_FALSE;
588: Mat dA,dB;
589: Vec tvec;
590: DM *dms;
591: PetscViewer viewer = 0;
592: PetscBool dAeqdB = PETSC_FALSE, needRestricts = PETSC_FALSE;
595: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels with PCMGSetLevels() before setting up");
596: n = mglevels[0]->levels;
597: /* FIX: Move this to PCSetFromOptions_MG? */
598: if (mg->usedmfornumberoflevels) {
599: PetscInt levels;
600: DMGetRefineLevel(pc->dm,&levels);
601: levels++;
602: if (levels > n) { /* the problem is now being solved on a finer grid */
603: PCMGSetLevels(pc,levels,NULL);
604: n = levels;
605: PCSetFromOptions(pc); /* it is bad to call this here, but otherwise will never be called for the new hierarchy */
606: mglevels = mg->levels;
607: }
608: }
609: KSPGetPC(mglevels[0]->smoothd,&cpc);
612: /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */
613: /* so use those from global PC */
614: /* Is this what we always want? What if user wants to keep old one? */
615: KSPGetOperatorsSet(mglevels[n-1]->smoothd,NULL,&opsset);
616: if (opsset) {
617: Mat mmat;
618: KSPGetOperators(mglevels[n-1]->smoothd,NULL,&mmat);
619: if (mmat == pc->pmat) opsset = PETSC_FALSE;
620: }
622: if (!opsset) {
623: PCGetUseAmat(pc,&use_amat);
624: if(use_amat){
625: PetscInfo(pc,"Using outer operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
626: KSPSetOperators(mglevels[n-1]->smoothd,pc->mat,pc->pmat);
627: }
628: else {
629: PetscInfo(pc,"Using matrix (pmat) operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
630: KSPSetOperators(mglevels[n-1]->smoothd,pc->pmat,pc->pmat);
631: }
632: }
634: for (i=n-1; i>0; i--) {
635: if (!(mglevels[i]->interpolate || mglevels[i]->restrct)) {
636: missinginterpolate = PETSC_TRUE;
637: continue;
638: }
639: }
641: KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB);
642: if (dA == dB) dAeqdB = PETSC_TRUE;
643: if ((mg->galerkin == PC_MG_GALERKIN_NONE) || (((mg->galerkin == PC_MG_GALERKIN_PMAT) || (mg->galerkin == PC_MG_GALERKIN_MAT)) && !dAeqdB)) {
644: needRestricts = PETSC_TRUE; /* user must compute either mat, pmat, or both so must restrict x to coarser levels */
645: }
648: /*
649: Skipping if user has provided all interpolation/restriction needed (since DM might not be able to produce them (when coming from SNES/TS)
650: Skipping for galerkin==2 (externally managed hierarchy such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs?
651: */
652: if (missinginterpolate && pc->dm && mg->galerkin != PC_MG_GALERKIN_EXTERNAL && !pc->setupcalled) {
653: /* construct the interpolation from the DMs */
654: Mat p;
655: Vec rscale;
656: PetscMalloc1(n,&dms);
657: dms[n-1] = pc->dm;
658: /* Separately create them so we do not get DMKSP interference between levels */
659: for (i=n-2; i>-1; i--) {DMCoarsen(dms[i+1],MPI_COMM_NULL,&dms[i]);}
660: /*
661: Force the mat type of coarse level operator to be AIJ because usually we want to use LU for coarse level.
662: Notice that it can be overwritten by -mat_type because KSPSetUp() reads command line options.
663: But it is safe to use -dm_mat_type.
665: The mat type should not be hardcoded like this, we need to find a better way.
666: DMSetMatType(dms[0],MATAIJ);
667: */
668: for (i=n-2; i>-1; i--) {
669: DMKSP kdm;
670: PetscBool dmhasrestrict, dmhasinject;
671: KSPSetDM(mglevels[i]->smoothd,dms[i]);
672: if (!needRestricts) {KSPSetDMActive(mglevels[i]->smoothd,PETSC_FALSE);}
673: DMGetDMKSPWrite(dms[i],&kdm);
674: /* Ugly hack so that the next KSPSetUp() will use the RHS that we set. A better fix is to change dmActive to take
675: * a bitwise OR of computing the matrix, RHS, and initial iterate. */
676: kdm->ops->computerhs = NULL;
677: kdm->rhsctx = NULL;
678: if (!mglevels[i+1]->interpolate) {
679: DMCreateInterpolation(dms[i],dms[i+1],&p,&rscale);
680: PCMGSetInterpolation(pc,i+1,p);
681: if (rscale) {PCMGSetRScale(pc,i+1,rscale);}
682: VecDestroy(&rscale);
683: MatDestroy(&p);
684: }
685: DMHasCreateRestriction(dms[i],&dmhasrestrict);
686: if (dmhasrestrict && !mglevels[i+1]->restrct){
687: DMCreateRestriction(dms[i],dms[i+1],&p);
688: PCMGSetRestriction(pc,i+1,p);
689: MatDestroy(&p);
690: }
691: DMHasCreateInjection(dms[i],&dmhasinject);
692: if (dmhasinject && !mglevels[i+1]->inject){
693: DMCreateInjection(dms[i],dms[i+1],&p);
694: PCMGSetInjection(pc,i+1,p);
695: MatDestroy(&p);
696: }
697: }
699: for (i=n-2; i>-1; i--) {DMDestroy(&dms[i]);}
700: PetscFree(dms);
701: }
703: if (pc->dm && !pc->setupcalled) {
704: /* finest smoother also gets DM but it is not active, independent of whether galerkin==PC_MG_GALERKIN_EXTERNAL */
705: KSPSetDM(mglevels[n-1]->smoothd,pc->dm);
706: KSPSetDMActive(mglevels[n-1]->smoothd,PETSC_FALSE);
707: }
709: if (mg->galerkin < PC_MG_GALERKIN_NONE) {
710: Mat A,B;
711: PetscBool doA = PETSC_FALSE,doB = PETSC_FALSE;
712: MatReuse reuse = MAT_INITIAL_MATRIX;
714: if ((mg->galerkin == PC_MG_GALERKIN_PMAT) || (mg->galerkin == PC_MG_GALERKIN_BOTH)) doB = PETSC_TRUE;
715: if ((mg->galerkin == PC_MG_GALERKIN_MAT) || ((mg->galerkin == PC_MG_GALERKIN_BOTH) && (dA != dB))) doA = PETSC_TRUE;
716: if (pc->setupcalled) reuse = MAT_REUSE_MATRIX;
717: for (i=n-2; i>-1; i--) {
718: 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");
719: if (!mglevels[i+1]->interpolate) {
720: PCMGSetInterpolation(pc,i+1,mglevels[i+1]->restrct);
721: }
722: if (!mglevels[i+1]->restrct) {
723: PCMGSetRestriction(pc,i+1,mglevels[i+1]->interpolate);
724: }
725: if (reuse == MAT_REUSE_MATRIX) {
726: KSPGetOperators(mglevels[i]->smoothd,&A,&B);
727: }
728: if (doA) {
729: MatGalerkin(mglevels[i+1]->restrct,dA,mglevels[i+1]->interpolate,reuse,1.0,&A);
730: }
731: if (doB) {
732: MatGalerkin(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,reuse,1.0,&B);
733: }
734: /* the management of the PetscObjectReference() and PetscObjecDereference() below is rather delicate */
735: if (!doA && dAeqdB) {
736: if (reuse == MAT_INITIAL_MATRIX) {PetscObjectReference((PetscObject)B);}
737: A = B;
738: } else if (!doA && reuse == MAT_INITIAL_MATRIX ) {
739: KSPGetOperators(mglevels[i]->smoothd,&A,NULL);
740: PetscObjectReference((PetscObject)A);
741: }
742: if (!doB && dAeqdB) {
743: if (reuse == MAT_INITIAL_MATRIX) {PetscObjectReference((PetscObject)A);}
744: B = A;
745: } else if (!doB && reuse == MAT_INITIAL_MATRIX) {
746: KSPGetOperators(mglevels[i]->smoothd,NULL,&B);
747: PetscObjectReference((PetscObject)B);
748: }
749: if (reuse == MAT_INITIAL_MATRIX) {
750: KSPSetOperators(mglevels[i]->smoothd,A,B);
751: PetscObjectDereference((PetscObject)A);
752: PetscObjectDereference((PetscObject)B);
753: }
754: dA = A;
755: dB = B;
756: }
757: }
758: if (needRestricts && pc->dm && pc->dm->x) {
759: /* need to restrict Jacobian location to coarser meshes for evaluation */
760: for (i=n-2; i>-1; i--) {
761: Mat R;
762: Vec rscale;
763: if (!mglevels[i]->smoothd->dm->x) {
764: Vec *vecs;
765: KSPCreateVecs(mglevels[i]->smoothd,1,&vecs,0,NULL);
766: mglevels[i]->smoothd->dm->x = vecs[0];
767: PetscFree(vecs);
768: }
769: PCMGGetRestriction(pc,i+1,&R);
770: PCMGGetRScale(pc,i+1,&rscale);
771: MatRestrict(R,mglevels[i+1]->smoothd->dm->x,mglevels[i]->smoothd->dm->x);
772: VecPointwiseMult(mglevels[i]->smoothd->dm->x,mglevels[i]->smoothd->dm->x,rscale);
773: }
774: }
775: if (needRestricts && pc->dm) {
776: for (i=n-2; i>=0; i--) {
777: DM dmfine,dmcoarse;
778: Mat Restrict,Inject;
779: Vec rscale;
780: KSPGetDM(mglevels[i+1]->smoothd,&dmfine);
781: KSPGetDM(mglevels[i]->smoothd,&dmcoarse);
782: PCMGGetRestriction(pc,i+1,&Restrict);
783: PCMGGetRScale(pc,i+1,&rscale);
784: PCMGGetInjection(pc,i+1,&Inject);
785: DMRestrict(dmfine,Restrict,rscale,Inject,dmcoarse);
786: }
787: }
789: if (!pc->setupcalled) {
790: for (i=0; i<n; i++) {
791: KSPSetFromOptions(mglevels[i]->smoothd);
792: }
793: for (i=1; i<n; i++) {
794: if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) {
795: KSPSetFromOptions(mglevels[i]->smoothu);
796: }
797: }
798: /* insure that if either interpolation or restriction is set the other other one is set */
799: for (i=1; i<n; i++) {
800: PCMGGetInterpolation(pc,i,NULL);
801: PCMGGetRestriction(pc,i,NULL);
802: }
803: for (i=0; i<n-1; i++) {
804: if (!mglevels[i]->b) {
805: Vec *vec;
806: KSPCreateVecs(mglevels[i]->smoothd,1,&vec,0,NULL);
807: PCMGSetRhs(pc,i,*vec);
808: VecDestroy(vec);
809: PetscFree(vec);
810: }
811: if (!mglevels[i]->r && i) {
812: VecDuplicate(mglevels[i]->b,&tvec);
813: PCMGSetR(pc,i,tvec);
814: VecDestroy(&tvec);
815: }
816: if (!mglevels[i]->x) {
817: VecDuplicate(mglevels[i]->b,&tvec);
818: PCMGSetX(pc,i,tvec);
819: VecDestroy(&tvec);
820: }
821: }
822: if (n != 1 && !mglevels[n-1]->r) {
823: /* PCMGSetR() on the finest level if user did not supply it */
824: Vec *vec;
825: KSPCreateVecs(mglevels[n-1]->smoothd,1,&vec,0,NULL);
826: PCMGSetR(pc,n-1,*vec);
827: VecDestroy(vec);
828: PetscFree(vec);
829: }
830: }
832: if (pc->dm) {
833: /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */
834: for (i=0; i<n-1; i++) {
835: if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
836: }
837: }
839: for (i=1; i<n; i++) {
840: if (mglevels[i]->smoothu == mglevels[i]->smoothd || mg->am == PC_MG_FULL || mg->am == PC_MG_KASKADE || mg->cyclesperpcapply > 1){
841: /* if doing only down then initial guess is zero */
842: KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);
843: }
844: if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
845: KSPSetUp(mglevels[i]->smoothd);
846: if (mglevels[i]->smoothd->reason == KSP_DIVERGED_PCSETUP_FAILED) {
847: pc->failedreason = PC_SUBPC_ERROR;
848: }
849: if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
850: if (!mglevels[i]->residual) {
851: Mat mat;
852: KSPGetOperators(mglevels[i]->smoothd,&mat,NULL);
853: PCMGSetResidual(pc,i,PCMGResidualDefault,mat);
854: }
855: }
856: for (i=1; i<n; i++) {
857: if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) {
858: Mat downmat,downpmat;
860: /* check if operators have been set for up, if not use down operators to set them */
861: KSPGetOperatorsSet(mglevels[i]->smoothu,&opsset,NULL);
862: if (!opsset) {
863: KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat);
864: KSPSetOperators(mglevels[i]->smoothu,downmat,downpmat);
865: }
867: KSPSetInitialGuessNonzero(mglevels[i]->smoothu,PETSC_TRUE);
868: if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
869: KSPSetUp(mglevels[i]->smoothu);
870: if (mglevels[i]->smoothu->reason == KSP_DIVERGED_PCSETUP_FAILED) {
871: pc->failedreason = PC_SUBPC_ERROR;
872: }
873: if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
874: }
875: }
877: if (mglevels[0]->eventsmoothsetup) {PetscLogEventBegin(mglevels[0]->eventsmoothsetup,0,0,0,0);}
878: KSPSetUp(mglevels[0]->smoothd);
879: if (mglevels[0]->smoothd->reason == KSP_DIVERGED_PCSETUP_FAILED) {
880: pc->failedreason = PC_SUBPC_ERROR;
881: }
882: if (mglevels[0]->eventsmoothsetup) {PetscLogEventEnd(mglevels[0]->eventsmoothsetup,0,0,0,0);}
884: /*
885: Dump the interpolation/restriction matrices plus the
886: Jacobian/stiffness on each level. This allows MATLAB users to
887: easily check if the Galerkin condition A_c = R A_f R^T is satisfied.
889: Only support one or the other at the same time.
890: */
891: #if defined(PETSC_USE_SOCKET_VIEWER)
892: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_mg_dump_matlab",&dump,NULL);
893: if (dump) viewer = PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc));
894: dump = PETSC_FALSE;
895: #endif
896: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_mg_dump_binary",&dump,NULL);
897: if (dump) viewer = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc));
899: if (viewer) {
900: for (i=1; i<n; i++) {
901: MatView(mglevels[i]->restrct,viewer);
902: }
903: for (i=0; i<n; i++) {
904: KSPGetPC(mglevels[i]->smoothd,&pc);
905: MatView(pc->mat,viewer);
906: }
907: }
908: return(0);
909: }
911: /* -------------------------------------------------------------------------------------*/
913: /*@
914: PCMGGetLevels - Gets the number of levels to use with MG.
916: Not Collective
918: Input Parameter:
919: . pc - the preconditioner context
921: Output parameter:
922: . levels - the number of levels
924: Level: advanced
926: .keywords: MG, get, levels, multigrid
928: .seealso: PCMGSetLevels()
929: @*/
930: PetscErrorCode PCMGGetLevels(PC pc,PetscInt *levels)
931: {
932: PC_MG *mg = (PC_MG*)pc->data;
937: *levels = mg->nlevels;
938: return(0);
939: }
941: /*@
942: PCMGSetType - Determines the form of multigrid to use:
943: multiplicative, additive, full, or the Kaskade algorithm.
945: Logically Collective on PC
947: Input Parameters:
948: + pc - the preconditioner context
949: - form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,
950: PC_MG_FULL, PC_MG_KASKADE
952: Options Database Key:
953: . -pc_mg_type <form> - Sets <form>, one of multiplicative,
954: additive, full, kaskade
956: Level: advanced
958: .keywords: MG, set, method, multiplicative, additive, full, Kaskade, multigrid
960: .seealso: PCMGSetLevels()
961: @*/
962: PetscErrorCode PCMGSetType(PC pc,PCMGType form)
963: {
964: PC_MG *mg = (PC_MG*)pc->data;
969: mg->am = form;
970: if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
971: else pc->ops->applyrichardson = NULL;
972: return(0);
973: }
975: /*@
976: PCMGGetType - Determines the form of multigrid to use:
977: multiplicative, additive, full, or the Kaskade algorithm.
979: Logically Collective on PC
981: Input Parameter:
982: . pc - the preconditioner context
984: Output Parameter:
985: . type - one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,PC_MG_FULL, PC_MG_KASKADE
988: Level: advanced
990: .keywords: MG, set, method, multiplicative, additive, full, Kaskade, multigrid
992: .seealso: PCMGSetLevels()
993: @*/
994: PetscErrorCode PCMGGetType(PC pc,PCMGType *type)
995: {
996: PC_MG *mg = (PC_MG*)pc->data;
1000: *type = mg->am;
1001: return(0);
1002: }
1004: /*@
1005: PCMGSetCycleType - Sets the type cycles to use. Use PCMGSetCycleTypeOnLevel() for more
1006: complicated cycling.
1008: Logically Collective on PC
1010: Input Parameters:
1011: + pc - the multigrid context
1012: - n - either PC_MG_CYCLE_V or PC_MG_CYCLE_W
1014: Options Database Key:
1015: . -pc_mg_cycle_type <v,w> - provide the cycle desired
1017: Level: advanced
1019: .keywords: MG, set, cycles, V-cycle, W-cycle, multigrid
1021: .seealso: PCMGSetCycleTypeOnLevel()
1022: @*/
1023: PetscErrorCode PCMGSetCycleType(PC pc,PCMGCycleType n)
1024: {
1025: PC_MG *mg = (PC_MG*)pc->data;
1026: PC_MG_Levels **mglevels = mg->levels;
1027: PetscInt i,levels;
1032: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling");
1033: levels = mglevels[0]->levels;
1034: for (i=0; i<levels; i++) mglevels[i]->cycles = n;
1035: return(0);
1036: }
1038: /*@
1039: PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step
1040: of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used
1042: Logically Collective on PC
1044: Input Parameters:
1045: + pc - the multigrid context
1046: - n - number of cycles (default is 1)
1048: Options Database Key:
1049: . -pc_mg_multiplicative_cycles n
1051: Level: advanced
1053: Notes:
1054: This is not associated with setting a v or w cycle, that is set with PCMGSetCycleType()
1056: .keywords: MG, set, cycles, V-cycle, W-cycle, multigrid
1058: .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType()
1059: @*/
1060: PetscErrorCode PCMGMultiplicativeSetCycles(PC pc,PetscInt n)
1061: {
1062: PC_MG *mg = (PC_MG*)pc->data;
1067: mg->cyclesperpcapply = n;
1068: return(0);
1069: }
1071: PetscErrorCode PCMGSetGalerkin_MG(PC pc,PCMGGalerkinType use)
1072: {
1073: PC_MG *mg = (PC_MG*)pc->data;
1076: mg->galerkin = use;
1077: return(0);
1078: }
1080: /*@
1081: PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the
1082: finest grid via the Galerkin process: A_i-1 = r_i * A_i * p_i
1084: Logically Collective on PC
1086: Input Parameters:
1087: + pc - the multigrid context
1088: - use - one of PC_MG_GALERKIN_BOTH,PC_MG_GALERKIN_PMAT,PC_MG_GALERKIN_MAT, or PC_MG_GALERKIN_NONE
1090: Options Database Key:
1091: . -pc_mg_galerkin <both,pmat,mat,none>
1093: Level: intermediate
1095: Notes:
1096: Some codes that use PCMG such as PCGAMG use Galerkin internally while constructing the hierarchy and thus do not
1097: use the PCMG construction of the coarser grids.
1099: .keywords: MG, set, Galerkin
1101: .seealso: PCMGGetGalerkin(), PCMGGalerkinType
1103: @*/
1104: PetscErrorCode PCMGSetGalerkin(PC pc,PCMGGalerkinType use)
1105: {
1110: PetscTryMethod(pc,"PCMGSetGalerkin_C",(PC,PCMGGalerkinType),(pc,use));
1111: return(0);
1112: }
1114: /*@
1115: PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e.
1116: A_i-1 = r_i * A_i * p_i
1118: Not Collective
1120: Input Parameter:
1121: . pc - the multigrid context
1123: Output Parameter:
1124: . galerkin - one of PC_MG_GALERKIN_BOTH,PC_MG_GALERKIN_PMAT,PC_MG_GALERKIN_MAT, PC_MG_GALERKIN_NONE, or PC_MG_GALERKIN_EXTERNAL
1126: Level: intermediate
1128: .keywords: MG, set, Galerkin
1130: .seealso: PCMGSetGalerkin(), PCMGGalerkinType
1132: @*/
1133: PetscErrorCode PCMGGetGalerkin(PC pc,PCMGGalerkinType *galerkin)
1134: {
1135: PC_MG *mg = (PC_MG*)pc->data;
1139: *galerkin = mg->galerkin;
1140: return(0);
1141: }
1143: /*@
1144: PCMGSetNumberSmooth - Sets the number of pre and post-smoothing steps to use
1145: on all levels. Use PCMGDistinctSmoothUp() to create separate up and down smoothers if you want different numbers of
1146: pre- and post-smoothing steps.
1148: Logically Collective on PC
1150: Input Parameters:
1151: + mg - the multigrid context
1152: - n - the number of smoothing steps
1154: Options Database Key:
1155: + -mg_levels_ksp_max_it <n> - Sets number of pre and post-smoothing steps
1157: Level: advanced
1159: Notes:
1160: this does not set a value on the coarsest grid, since we assume that
1161: there is no separate smooth up on the coarsest grid.
1163: .keywords: MG, smooth, up, post-smoothing, steps, multigrid
1165: .seealso: PCMGSetDistinctSmoothUp()
1166: @*/
1167: PetscErrorCode PCMGSetNumberSmooth(PC pc,PetscInt n)
1168: {
1169: PC_MG *mg = (PC_MG*)pc->data;
1170: PC_MG_Levels **mglevels = mg->levels;
1172: PetscInt i,levels;
1177: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling");
1178: levels = mglevels[0]->levels;
1180: for (i=1; i<levels; i++) {
1181: KSPSetTolerances(mglevels[i]->smoothu,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
1182: KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
1183: mg->default_smoothu = n;
1184: mg->default_smoothd = n;
1185: }
1186: return(0);
1187: }
1189: /*@
1190: PCMGSetDistinctSmoothUp - sets the up (post) smoother to be a seperate KSP from the down (pre) smoother on all levels
1191: and adds the suffix _up to the options name
1193: Logically Collective on PC
1195: Input Parameters:
1196: . pc - the preconditioner context
1198: Options Database Key:
1199: . -pc_mg_distinct_smoothup
1201: Level: advanced
1203: Notes:
1204: this does not set a value on the coarsest grid, since we assume that
1205: there is no separate smooth up on the coarsest grid.
1207: .keywords: MG, smooth, up, post-smoothing, steps, multigrid
1209: .seealso: PCMGSetNumberSmooth()
1210: @*/
1211: PetscErrorCode PCMGSetDistinctSmoothUp(PC pc)
1212: {
1213: PC_MG *mg = (PC_MG*)pc->data;
1214: PC_MG_Levels **mglevels = mg->levels;
1216: PetscInt i,levels;
1217: KSP subksp;
1221: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling");
1222: levels = mglevels[0]->levels;
1224: for (i=1; i<levels; i++) {
1225: const char *prefix = NULL;
1226: /* make sure smoother up and down are different */
1227: PCMGGetSmootherUp(pc,i,&subksp);
1228: KSPGetOptionsPrefix(mglevels[i]->smoothd,&prefix);
1229: KSPSetOptionsPrefix(subksp,prefix);
1230: KSPAppendOptionsPrefix(subksp,"up_");
1231: }
1232: return(0);
1233: }
1235: /* ----------------------------------------------------------------------------------------*/
1237: /*MC
1238: PCMG - Use multigrid preconditioning. This preconditioner requires you provide additional
1239: information about the coarser grid matrices and restriction/interpolation operators.
1241: Options Database Keys:
1242: + -pc_mg_levels <nlevels> - number of levels including finest
1243: . -pc_mg_cycle_type <v,w> - provide the cycle desired
1244: . -pc_mg_type <additive,multiplicative,full,kaskade> - multiplicative is the default
1245: . -pc_mg_log - log information about time spent on each level of the solver
1246: . -pc_mg_distinct_smoothup - configure up (after interpolation) and down (before restriction) smoothers separately (with different options prefixes)
1247: . -pc_mg_galerkin <both,pmat,mat,none> - use Galerkin process to compute coarser operators, i.e. Acoarse = R A R'
1248: . -pc_mg_multiplicative_cycles - number of cycles to use as the preconditioner (defaults to 1)
1249: . -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices
1250: to the Socket viewer for reading from MATLAB.
1251: - -pc_mg_dump_binary - dumps the matrices for each level and the restriction/interpolation matrices
1252: to the binary output file called binaryoutput
1254: Notes:
1255: If one uses a Krylov method such GMRES or CG as the smoother then one must use KSPFGMRES, KSPGCR, or KSPRICHARDSON as the outer Krylov method
1257: When run with a single level the smoother options are used on that level NOT the coarse grid solver options
1259: When run with KSPRICHARDSON the convergence test changes slightly if monitor is turned on. The iteration count may change slightly. This
1260: is because without monitoring the residual norm is computed WITHIN each multigrid cycle on the finest level after the pre-smoothing
1261: (because the residual has just been computed for the multigrid algorithm and is hence available for free) while with monitoring the
1262: residual is computed at the end of each cycle.
1264: Level: intermediate
1266: Concepts: multigrid/multilevel
1268: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, PCEXOTIC, PCGAMG, PCML, PCHYPRE
1269: PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(),
1270: PCMGSetDistinctSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
1271: PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
1272: PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()
1273: M*/
1275: PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc)
1276: {
1277: PC_MG *mg;
1281: PetscNewLog(pc,&mg);
1282: pc->data = (void*)mg;
1283: mg->nlevels = -1;
1284: mg->am = PC_MG_MULTIPLICATIVE;
1285: mg->galerkin = PC_MG_GALERKIN_NONE;
1287: pc->useAmat = PETSC_TRUE;
1289: pc->ops->apply = PCApply_MG;
1290: pc->ops->setup = PCSetUp_MG;
1291: pc->ops->reset = PCReset_MG;
1292: pc->ops->destroy = PCDestroy_MG;
1293: pc->ops->setfromoptions = PCSetFromOptions_MG;
1294: pc->ops->view = PCView_MG;
1296: PetscObjectComposeFunction((PetscObject)pc,"PCMGSetGalerkin_C",PCMGSetGalerkin_MG);
1297: return(0);
1298: }