Actual source code: fas.c
1: /* Defines the basic SNES object */
2: #include <../src/snes/impls/fas/fasimpls.h>
4: const char *const SNESFASTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","SNESFASType","SNES_FAS",NULL};
6: static PetscErrorCode SNESReset_FAS(SNES snes)
7: {
8: SNES_FAS *fas = (SNES_FAS*)snes->data;
12: SNESDestroy(&fas->smoothu);
13: SNESDestroy(&fas->smoothd);
14: MatDestroy(&fas->inject);
15: MatDestroy(&fas->interpolate);
16: MatDestroy(&fas->restrct);
17: VecDestroy(&fas->rscale);
18: VecDestroy(&fas->Xg);
19: VecDestroy(&fas->Fg);
20: if (fas->next) {SNESReset(fas->next);}
21: return(0);
22: }
24: static PetscErrorCode SNESDestroy_FAS(SNES snes)
25: {
26: SNES_FAS *fas = (SNES_FAS*)snes->data;
30: /* recursively resets and then destroys */
31: SNESReset_FAS(snes);
32: SNESDestroy(&fas->next);
33: PetscFree(fas);
34: return(0);
35: }
37: static PetscErrorCode SNESFASSetUpLineSearch_Private(SNES snes, SNES smooth)
38: {
39: SNESLineSearch linesearch;
40: SNESLineSearch slinesearch;
41: void *lsprectx,*lspostctx;
42: PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*);
43: PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
47: if (!snes->linesearch) return(0);
48: SNESGetLineSearch(snes,&linesearch);
49: SNESGetLineSearch(smooth,&slinesearch);
50: SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
51: SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
52: SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
53: SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
54: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
55: return(0);
56: }
58: static PetscErrorCode SNESFASCycleSetUpSmoother_Private(SNES snes, SNES smooth)
59: {
60: SNES_FAS *fas = (SNES_FAS*) snes->data;
64: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)smooth);
65: SNESSetFromOptions(smooth);
66: SNESFASSetUpLineSearch_Private(snes, smooth);
68: PetscObjectReference((PetscObject)snes->vec_sol);
69: PetscObjectReference((PetscObject)snes->vec_sol_update);
70: PetscObjectReference((PetscObject)snes->vec_func);
71: smooth->vec_sol = snes->vec_sol;
72: smooth->vec_sol_update = snes->vec_sol_update;
73: smooth->vec_func = snes->vec_func;
75: if (fas->eventsmoothsetup) {PetscLogEventBegin(fas->eventsmoothsetup,smooth,0,0,0);}
76: SNESSetUp(smooth);
77: if (fas->eventsmoothsetup) {PetscLogEventEnd(fas->eventsmoothsetup,smooth,0,0,0);}
78: return(0);
79: }
81: static PetscErrorCode SNESSetUp_FAS(SNES snes)
82: {
83: SNES_FAS *fas = (SNES_FAS*) snes->data;
85: PetscInt dm_levels;
86: SNES next;
87: PetscBool isFine, hasCreateRestriction, hasCreateInjection;
90: SNESFASCycleIsFine(snes, &isFine);
91: if (fas->usedmfornumberoflevels && isFine) {
92: DMGetRefineLevel(snes->dm,&dm_levels);
93: dm_levels++;
94: if (dm_levels > fas->levels) {
95: /* reset the number of levels */
96: SNESFASSetLevels(snes,dm_levels,NULL);
97: SNESSetFromOptions(snes);
98: }
99: }
100: SNESFASCycleGetCorrection(snes, &next);
101: if (!isFine) snes->gridsequence = 0; /* no grid sequencing inside the multigrid hierarchy! */
103: SNESSetWorkVecs(snes, 2); /* work vectors used for intergrid transfers */
105: /* set up the smoothers if they haven't already been set up */
106: if (!fas->smoothd) {
107: SNESFASCycleCreateSmoother_Private(snes, &fas->smoothd);
108: }
110: if (snes->dm) {
111: /* set the smoother DMs properly */
112: if (fas->smoothu) {SNESSetDM(fas->smoothu, snes->dm);}
113: SNESSetDM(fas->smoothd, snes->dm);
114: /* construct EVERYTHING from the DM -- including the progressive set of smoothers */
115: if (next) {
116: /* for now -- assume the DM and the evaluation functions have been set externally */
117: if (!next->dm) {
118: DMCoarsen(snes->dm, PetscObjectComm((PetscObject)next), &next->dm);
119: SNESSetDM(next, next->dm);
120: }
121: /* set the interpolation and restriction from the DM */
122: if (!fas->interpolate) {
123: DMCreateInterpolation(next->dm, snes->dm, &fas->interpolate, &fas->rscale);
124: if (!fas->restrct) {
125: DMHasCreateRestriction(next->dm, &hasCreateRestriction);
126: /* DM can create restrictions, use that */
127: if (hasCreateRestriction) {
128: DMCreateRestriction(next->dm, snes->dm, &fas->restrct);
129: } else {
130: PetscObjectReference((PetscObject)fas->interpolate);
131: fas->restrct = fas->interpolate;
132: }
133: }
134: }
135: /* set the injection from the DM */
136: if (!fas->inject) {
137: DMHasCreateInjection(next->dm, &hasCreateInjection);
138: if (hasCreateInjection) {
139: DMCreateInjection(next->dm, snes->dm, &fas->inject);
140: }
141: }
142: }
143: }
145: /*pass the smoother, function, and jacobian up to the next level if it's not user set already */
146: if (fas->galerkin) {
147: if (next) {
148: SNESSetFunction(next, NULL, SNESFASGalerkinFunctionDefault, next);
149: }
150: if (fas->smoothd && fas->level != fas->levels - 1) {
151: SNESSetFunction(fas->smoothd, NULL, SNESFASGalerkinFunctionDefault, snes);
152: }
153: if (fas->smoothu && fas->level != fas->levels - 1) {
154: SNESSetFunction(fas->smoothu, NULL, SNESFASGalerkinFunctionDefault, snes);
155: }
156: }
158: /* sets the down (pre) smoother's default norm and sets it from options */
159: if (fas->smoothd) {
160: if (fas->level == 0 && fas->levels != 1) {
161: SNESSetNormSchedule(fas->smoothd, SNES_NORM_NONE);
162: } else {
163: SNESSetNormSchedule(fas->smoothd, SNES_NORM_FINAL_ONLY);
164: }
165: SNESFASCycleSetUpSmoother_Private(snes, fas->smoothd);
166: }
168: /* sets the up (post) smoother's default norm and sets it from options */
169: if (fas->smoothu) {
170: if (fas->level != fas->levels - 1) {
171: SNESSetNormSchedule(fas->smoothu, SNES_NORM_NONE);
172: } else {
173: SNESSetNormSchedule(fas->smoothu, SNES_NORM_FINAL_ONLY);
174: }
175: SNESFASCycleSetUpSmoother_Private(snes, fas->smoothu);
176: }
178: if (next) {
179: /* gotta set up the solution vector for this to work */
180: if (!next->vec_sol) {SNESFASCreateCoarseVec(snes,&next->vec_sol);}
181: if (!next->vec_rhs) {SNESFASCreateCoarseVec(snes,&next->vec_rhs);}
182: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next);
183: SNESFASSetUpLineSearch_Private(snes, next);
184: SNESSetUp(next);
185: }
187: /* setup FAS work vectors */
188: if (fas->galerkin) {
189: VecDuplicate(snes->vec_sol, &fas->Xg);
190: VecDuplicate(snes->vec_sol, &fas->Fg);
191: }
192: return(0);
193: }
195: static PetscErrorCode SNESSetFromOptions_FAS(PetscOptionItems *PetscOptionsObject,SNES snes)
196: {
197: SNES_FAS *fas = (SNES_FAS*) snes->data;
198: PetscInt levels = 1;
199: PetscBool flg = PETSC_FALSE, upflg = PETSC_FALSE, downflg = PETSC_FALSE, monflg = PETSC_FALSE, galerkinflg = PETSC_FALSE,continuationflg = PETSC_FALSE;
201: SNESFASType fastype;
202: const char *optionsprefix;
203: SNESLineSearch linesearch;
204: PetscInt m, n_up, n_down;
205: SNES next;
206: PetscBool isFine;
209: SNESFASCycleIsFine(snes, &isFine);
210: PetscOptionsHead(PetscOptionsObject,"SNESFAS Options-----------------------------------");
212: /* number of levels -- only process most options on the finest level */
213: if (isFine) {
214: PetscOptionsInt("-snes_fas_levels", "Number of Levels", "SNESFASSetLevels", levels, &levels, &flg);
215: if (!flg && snes->dm) {
216: DMGetRefineLevel(snes->dm,&levels);
217: levels++;
218: fas->usedmfornumberoflevels = PETSC_TRUE;
219: }
220: SNESFASSetLevels(snes, levels, NULL);
221: fastype = fas->fastype;
222: PetscOptionsEnum("-snes_fas_type","FAS correction type","SNESFASSetType",SNESFASTypes,(PetscEnum)fastype,(PetscEnum*)&fastype,&flg);
223: if (flg) {
224: SNESFASSetType(snes, fastype);
225: }
227: SNESGetOptionsPrefix(snes, &optionsprefix);
228: PetscOptionsInt("-snes_fas_cycles","Number of cycles","SNESFASSetCycles",fas->n_cycles,&m,&flg);
229: if (flg) {
230: SNESFASSetCycles(snes, m);
231: }
232: PetscOptionsBool("-snes_fas_continuation","Corrected grid-sequence continuation","SNESFASSetContinuation",fas->continuation,&continuationflg,&flg);
233: if (flg) {
234: SNESFASSetContinuation(snes,continuationflg);
235: }
237: PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin","SNESFASSetGalerkin",fas->galerkin,&galerkinflg,&flg);
238: if (flg) {
239: SNESFASSetGalerkin(snes, galerkinflg);
240: }
242: if (fas->fastype == SNES_FAS_FULL) {
243: PetscOptionsBool("-snes_fas_full_downsweep","Smooth on the initial down sweep for full FAS cycles","SNESFASFullSetDownSweep",fas->full_downsweep,&fas->full_downsweep,&flg);
244: if (flg) {SNESFASFullSetDownSweep(snes,fas->full_downsweep);}
245: PetscOptionsBool("-snes_fas_full_total","Use total restriction and interpolaton on the indial down and up sweeps for the full FAS cycle","SNESFASFullSetUseTotal",fas->full_total,&fas->full_total,&flg);
246: if (flg) {SNESFASFullSetTotal(snes,fas->full_total);}
247: }
249: PetscOptionsInt("-snes_fas_smoothup","Number of post-smoothing steps","SNESFASSetNumberSmoothUp",fas->max_up_it,&n_up,&upflg);
251: PetscOptionsInt("-snes_fas_smoothdown","Number of pre-smoothing steps","SNESFASSetNumberSmoothDown",fas->max_down_it,&n_down,&downflg);
253: {
254: PetscViewer viewer;
255: PetscViewerFormat format;
256: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,"-snes_fas_monitor",&viewer,&format,&monflg);
257: if (monflg) {
258: PetscViewerAndFormat *vf;
259: PetscViewerAndFormatCreate(viewer,format,&vf);
260: PetscObjectDereference((PetscObject)viewer);
261: SNESFASSetMonitor(snes,vf,PETSC_TRUE);
262: }
263: }
264: flg = PETSC_FALSE;
265: monflg = PETSC_TRUE;
266: PetscOptionsBool("-snes_fas_log","Log times for each FAS level","SNESFASSetLog",monflg,&monflg,&flg);
267: if (flg) {SNESFASSetLog(snes,monflg);}
268: }
270: PetscOptionsTail();
272: /* setup from the determined types if there is no pointwise procedure or smoother defined */
273: if (upflg) {
274: SNESFASSetNumberSmoothUp(snes,n_up);
275: }
276: if (downflg) {
277: SNESFASSetNumberSmoothDown(snes,n_down);
278: }
280: /* set up the default line search for coarse grid corrections */
281: if (fas->fastype == SNES_FAS_ADDITIVE) {
282: if (!snes->linesearch) {
283: SNESGetLineSearch(snes, &linesearch);
284: SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
285: }
286: }
288: /* recursive option setting for the smoothers */
289: SNESFASCycleGetCorrection(snes, &next);
290: if (next) {SNESSetFromOptions(next);}
291: return(0);
292: }
294: #include <petscdraw.h>
295: static PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer)
296: {
297: SNES_FAS *fas = (SNES_FAS*) snes->data;
298: PetscBool isFine,iascii,isdraw;
299: PetscInt i;
301: SNES smoothu, smoothd, levelsnes;
304: SNESFASCycleIsFine(snes, &isFine);
305: if (isFine) {
306: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
307: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
308: if (iascii) {
309: PetscViewerASCIIPrintf(viewer, " type is %s, levels=%D, cycles=%D\n", SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles);
310: if (fas->galerkin) {
311: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid function evaluation\n");
312: } else {
313: PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid function evaluation\n");
314: }
315: for (i=0; i<fas->levels; i++) {
316: SNESFASGetCycleSNES(snes, i, &levelsnes);
317: SNESFASCycleGetSmootherUp(levelsnes, &smoothu);
318: SNESFASCycleGetSmootherDown(levelsnes, &smoothd);
319: if (!i) {
320: PetscViewerASCIIPrintf(viewer," Coarse grid solver -- level %D -------------------------------\n",i);
321: } else {
322: PetscViewerASCIIPrintf(viewer," Down solver (pre-smoother) on level %D -------------------------------\n",i);
323: }
324: PetscViewerASCIIPushTab(viewer);
325: if (smoothd) {
326: SNESView(smoothd,viewer);
327: } else {
328: PetscViewerASCIIPrintf(viewer,"Not yet available\n");
329: }
330: PetscViewerASCIIPopTab(viewer);
331: if (i && (smoothd == smoothu)) {
332: PetscViewerASCIIPrintf(viewer," Up solver (post-smoother) same as down solver (pre-smoother)\n");
333: } else if (i) {
334: PetscViewerASCIIPrintf(viewer," Up solver (post-smoother) on level %D -------------------------------\n",i);
335: PetscViewerASCIIPushTab(viewer);
336: if (smoothu) {
337: SNESView(smoothu,viewer);
338: } else {
339: PetscViewerASCIIPrintf(viewer,"Not yet available\n");
340: }
341: PetscViewerASCIIPopTab(viewer);
342: }
343: }
344: } else if (isdraw) {
345: PetscDraw draw;
346: PetscReal x,w,y,bottom,th,wth;
347: SNES_FAS *curfas = fas;
348: PetscViewerDrawGetDraw(viewer,0,&draw);
349: PetscDrawGetCurrentPoint(draw,&x,&y);
350: PetscDrawStringGetSize(draw,&wth,&th);
351: bottom = y - th;
352: while (curfas) {
353: if (!curfas->smoothu) {
354: PetscDrawPushCurrentPoint(draw,x,bottom);
355: if (curfas->smoothd) {SNESView(curfas->smoothd,viewer);}
356: PetscDrawPopCurrentPoint(draw);
357: } else {
358: w = 0.5*PetscMin(1.0-x,x);
359: PetscDrawPushCurrentPoint(draw,x-w,bottom);
360: if (curfas->smoothd) {SNESView(curfas->smoothd,viewer);}
361: PetscDrawPopCurrentPoint(draw);
362: PetscDrawPushCurrentPoint(draw,x+w,bottom);
363: if (curfas->smoothu) {SNESView(curfas->smoothu,viewer);}
364: PetscDrawPopCurrentPoint(draw);
365: }
366: /* this is totally bogus but we have no way of knowing how low the previous one was draw to */
367: bottom -= 5*th;
368: if (curfas->next) curfas = (SNES_FAS*)curfas->next->data;
369: else curfas = NULL;
370: }
371: }
372: }
373: return(0);
374: }
376: /*
377: Defines the action of the downsmoother
378: */
379: static PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
380: {
381: PetscErrorCode ierr;
382: SNESConvergedReason reason;
383: Vec FPC;
384: SNES smoothd;
385: PetscBool flg;
386: SNES_FAS *fas = (SNES_FAS*) snes->data;
389: SNESFASCycleGetSmootherDown(snes, &smoothd);
390: SNESSetInitialFunction(smoothd, F);
391: if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,smoothd,B,X,0);}
392: SNESSolve(smoothd, B, X);
393: if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,smoothd,B,X,0);}
394: /* check convergence reason for the smoother */
395: SNESGetConvergedReason(smoothd,&reason);
396: if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
397: snes->reason = SNES_DIVERGED_INNER;
398: return(0);
399: }
401: SNESGetFunction(smoothd, &FPC, NULL, NULL);
402: SNESGetAlwaysComputesFinalResidual(smoothd, &flg);
403: if (!flg) {
404: SNESComputeFunction(smoothd, X, FPC);
405: }
406: VecCopy(FPC, F);
407: if (fnorm) {VecNorm(F,NORM_2,fnorm);}
408: return(0);
409: }
411: /*
412: Defines the action of the upsmoother
413: */
414: static PetscErrorCode SNESFASUpSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
415: {
416: PetscErrorCode ierr;
417: SNESConvergedReason reason;
418: Vec FPC;
419: SNES smoothu;
420: PetscBool flg;
421: SNES_FAS *fas = (SNES_FAS*) snes->data;
424: SNESFASCycleGetSmootherUp(snes, &smoothu);
425: if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,smoothu,0,0,0);}
426: SNESSolve(smoothu, B, X);
427: if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,smoothu,0,0,0);}
428: /* check convergence reason for the smoother */
429: SNESGetConvergedReason(smoothu,&reason);
430: if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
431: snes->reason = SNES_DIVERGED_INNER;
432: return(0);
433: }
434: SNESGetFunction(smoothu, &FPC, NULL, NULL);
435: SNESGetAlwaysComputesFinalResidual(smoothu, &flg);
436: if (!flg) {
437: SNESComputeFunction(smoothu, X, FPC);
438: }
439: VecCopy(FPC, F);
440: if (fnorm) {VecNorm(F,NORM_2,fnorm);}
441: return(0);
442: }
444: /*@
445: SNESFASCreateCoarseVec - create Vec corresponding to a state vector on one level coarser than current level
447: Collective
449: Input Parameter:
450: . snes - SNESFAS
452: Output Parameter:
453: . Xcoarse - vector on level one coarser than snes
455: Level: developer
457: .seealso: SNESFASSetRestriction(), SNESFASRestrict()
458: @*/
459: PetscErrorCode SNESFASCreateCoarseVec(SNES snes,Vec *Xcoarse)
460: {
462: SNES_FAS *fas;
467: fas = (SNES_FAS*)snes->data;
468: if (fas->rscale) {
469: VecDuplicate(fas->rscale,Xcoarse);
470: } else if (fas->inject) {
471: MatCreateVecs(fas->inject,Xcoarse,NULL);
472: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or injection");
473: return(0);
474: }
476: /*@
477: SNESFASRestrict - restrict a Vec to the next coarser level
479: Collective
481: Input Parameters:
482: + fine - SNES from which to restrict
483: - Xfine - vector to restrict
485: Output Parameter:
486: . Xcoarse - result of restriction
488: Level: developer
490: .seealso: SNESFASSetRestriction(), SNESFASSetInjection()
491: @*/
492: PetscErrorCode SNESFASRestrict(SNES fine,Vec Xfine,Vec Xcoarse)
493: {
495: SNES_FAS *fas;
501: fas = (SNES_FAS*)fine->data;
502: if (fas->inject) {
503: MatRestrict(fas->inject,Xfine,Xcoarse);
504: } else {
505: MatRestrict(fas->restrct,Xfine,Xcoarse);
506: VecPointwiseMult(Xcoarse,fas->rscale,Xcoarse);
507: }
508: return(0);
509: }
511: /*
513: Performs a variant of FAS using the interpolated total coarse solution
515: fine problem: F(x) = b
516: coarse problem: F^c(x^c) = Rb, Initial guess Rx
517: interpolated solution: x^f = I x^c (total solution interpolation
519: */
520: static PetscErrorCode SNESFASInterpolatedCoarseSolution(SNES snes, Vec X, Vec X_new)
521: {
522: PetscErrorCode ierr;
523: Vec X_c, B_c;
524: SNESConvergedReason reason;
525: SNES next;
526: Mat restrct, interpolate;
527: SNES_FAS *fasc;
530: SNESFASCycleGetCorrection(snes, &next);
531: if (next) {
532: fasc = (SNES_FAS*)next->data;
534: SNESFASCycleGetRestriction(snes, &restrct);
535: SNESFASCycleGetInterpolation(snes, &interpolate);
537: X_c = next->vec_sol;
539: if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,snes,0,0,0);}
540: /* restrict the total solution: Rb */
541: SNESFASRestrict(snes, X, X_c);
542: B_c = next->vec_rhs;
543: if (snes->vec_rhs) {
544: /* restrict the total rhs defect: Rb */
545: MatRestrict(restrct, snes->vec_rhs, B_c);
546: } else {
547: VecSet(B_c, 0.);
548: }
549: if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,snes,0,0,0);}
551: SNESSolve(next, B_c, X_c);
552: SNESGetConvergedReason(next,&reason);
553: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
554: snes->reason = SNES_DIVERGED_INNER;
555: return(0);
556: }
557: /* x^f <- Ix^c*/
558: DM dmc,dmf;
560: SNESGetDM(next, &dmc);
561: SNESGetDM(snes, &dmf);
562: if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,snes,0,0,0);}
563: DMInterpolateSolution(dmc, dmf, interpolate, X_c, X_new);
564: if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,snes,0,0,0);}
565: PetscObjectSetName((PetscObject) X_c, "Coarse solution");
566: VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view");
567: PetscObjectSetName((PetscObject) X_new, "Updated Fine solution");
568: VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view");
569: }
570: return(0);
571: }
573: /*
575: Performs the FAS coarse correction as:
577: fine problem: F(x) = b
578: coarse problem: F^c(x^c) = b^c
580: b^c = F^c(Rx) - R(F(x) - b)
582: */
583: PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new)
584: {
585: PetscErrorCode ierr;
586: Vec X_c, Xo_c, F_c, B_c;
587: SNESConvergedReason reason;
588: SNES next;
589: Mat restrct, interpolate;
590: SNES_FAS *fasc;
593: SNESFASCycleGetCorrection(snes, &next);
594: if (next) {
595: fasc = (SNES_FAS*)next->data;
597: SNESFASCycleGetRestriction(snes, &restrct);
598: SNESFASCycleGetInterpolation(snes, &interpolate);
600: X_c = next->vec_sol;
601: Xo_c = next->work[0];
602: F_c = next->vec_func;
603: B_c = next->vec_rhs;
605: if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,snes,0,0,0);}
606: SNESFASRestrict(snes, X, Xo_c);
607: /* restrict the defect: R(F(x) - b) */
608: MatRestrict(restrct, F, B_c);
609: if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,snes,0,0,0);}
611: if (fasc->eventresidual) {PetscLogEventBegin(fasc->eventresidual,next,0,0,0);}
612: /* F_c = F^c(Rx) - R(F(x) - b) since the second term was sitting in next->vec_rhs */
613: SNESComputeFunction(next, Xo_c, F_c);
614: if (fasc->eventresidual) {PetscLogEventEnd(fasc->eventresidual,next,0,0,0);}
616: /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */
617: VecCopy(B_c, X_c);
618: VecCopy(F_c, B_c);
619: VecCopy(X_c, F_c);
620: /* set initial guess of the coarse problem to the projected fine solution */
621: VecCopy(Xo_c, X_c);
623: /* recurse to the next level */
624: SNESSetInitialFunction(next, F_c);
625: SNESSolve(next, B_c, X_c);
626: SNESGetConvergedReason(next,&reason);
627: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
628: snes->reason = SNES_DIVERGED_INNER;
629: return(0);
630: }
631: /* correct as x <- x + I(x^c - Rx)*/
632: VecAXPY(X_c, -1.0, Xo_c);
634: if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,snes,0,0,0);}
635: MatInterpolateAdd(interpolate, X_c, X, X_new);
636: if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,snes,0,0,0);}
637: PetscObjectSetName((PetscObject) X_c, "Coarse correction");
638: VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view");
639: PetscObjectSetName((PetscObject) X_new, "Updated Fine solution");
640: VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view");
641: }
642: return(0);
643: }
645: /*
647: The additive cycle looks like:
649: xhat = x
650: xhat = dS(x, b)
651: x = coarsecorrection(xhat, b_d)
652: x = x + nu*(xhat - x);
653: (optional) x = uS(x, b)
655: With the coarse RHS (defect correction) as below.
657: */
658: static PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X)
659: {
660: Vec F, B, Xhat;
661: Vec X_c, Xo_c, F_c, B_c;
662: PetscErrorCode ierr;
663: SNESConvergedReason reason;
664: PetscReal xnorm, fnorm, ynorm;
665: SNESLineSearchReason lsresult;
666: SNES next;
667: Mat restrct, interpolate;
668: SNES_FAS *fas = (SNES_FAS*)snes->data,*fasc;
671: SNESFASCycleGetCorrection(snes, &next);
672: F = snes->vec_func;
673: B = snes->vec_rhs;
674: Xhat = snes->work[1];
675: VecCopy(X, Xhat);
676: /* recurse first */
677: if (next) {
678: fasc = (SNES_FAS*)next->data;
679: SNESFASCycleGetRestriction(snes, &restrct);
680: SNESFASCycleGetInterpolation(snes, &interpolate);
681: if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,snes,0,0,0);}
682: SNESComputeFunction(snes, Xhat, F);
683: if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,snes,0,0,0);}
684: VecNorm(F, NORM_2, &fnorm);
685: X_c = next->vec_sol;
686: Xo_c = next->work[0];
687: F_c = next->vec_func;
688: B_c = next->vec_rhs;
690: SNESFASRestrict(snes,Xhat,Xo_c);
691: /* restrict the defect */
692: MatRestrict(restrct, F, B_c);
694: /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */
695: if (fasc->eventresidual) {PetscLogEventBegin(fasc->eventresidual,next,0,0,0);}
696: SNESComputeFunction(next, Xo_c, F_c);
697: if (fasc->eventresidual) {PetscLogEventEnd(fasc->eventresidual,next,0,0,0);}
698: VecCopy(B_c, X_c);
699: VecCopy(F_c, B_c);
700: VecCopy(X_c, F_c);
701: /* set initial guess of the coarse problem to the projected fine solution */
702: VecCopy(Xo_c, X_c);
704: /* recurse */
705: SNESSetInitialFunction(next, F_c);
706: SNESSolve(next, B_c, X_c);
708: /* smooth on this level */
709: SNESFASDownSmooth_Private(snes, B, X, F, &fnorm);
711: SNESGetConvergedReason(next,&reason);
712: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
713: snes->reason = SNES_DIVERGED_INNER;
714: return(0);
715: }
717: /* correct as x <- x + I(x^c - Rx)*/
718: VecAYPX(X_c, -1.0, Xo_c);
719: MatInterpolate(interpolate, X_c, Xhat);
721: /* additive correction of the coarse direction*/
722: SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat);
723: SNESLineSearchGetReason(snes->linesearch, &lsresult);
724: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm);
725: if (lsresult) {
726: if (++snes->numFailures >= snes->maxFailures) {
727: snes->reason = SNES_DIVERGED_LINE_SEARCH;
728: return(0);
729: }
730: }
731: } else {
732: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
733: }
734: return(0);
735: }
737: /*
739: Defines the FAS cycle as:
741: fine problem: F(x) = b
742: coarse problem: F^c(x) = b^c
744: b^c = F^c(Rx) - R(F(x) - b)
746: correction:
748: x = x + I(x^c - Rx)
750: */
751: static PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X)
752: {
755: Vec F,B;
756: SNES next;
759: F = snes->vec_func;
760: B = snes->vec_rhs;
761: /* pre-smooth -- just update using the pre-smoother */
762: SNESFASCycleGetCorrection(snes, &next);
763: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
764: if (next) {
765: SNESFASCoarseCorrection(snes, X, F, X);
766: SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
767: }
768: return(0);
769: }
771: static PetscErrorCode SNESFASCycleSetupPhase_Full(SNES snes)
772: {
773: SNES next;
774: SNES_FAS *fas = (SNES_FAS*)snes->data;
775: PetscBool isFine;
779: /* pre-smooth -- just update using the pre-smoother */
780: SNESFASCycleIsFine(snes,&isFine);
781: SNESFASCycleGetCorrection(snes,&next);
782: fas->full_stage = 0;
783: if (next) {SNESFASCycleSetupPhase_Full(next);}
784: return(0);
785: }
787: static PetscErrorCode SNESFASCycle_Full(SNES snes, Vec X)
788: {
790: Vec F,B;
791: SNES_FAS *fas = (SNES_FAS*)snes->data;
792: PetscBool isFine;
793: SNES next;
796: F = snes->vec_func;
797: B = snes->vec_rhs;
798: SNESFASCycleIsFine(snes,&isFine);
799: SNESFASCycleGetCorrection(snes,&next);
801: if (isFine) {
802: SNESFASCycleSetupPhase_Full(snes);
803: }
805: if (fas->full_stage == 0) {
806: /* downsweep */
807: if (next) {
808: if (fas->level != 1) next->max_its += 1;
809: if (fas->full_downsweep) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
810: fas->full_downsweep = PETSC_TRUE;
811: if (fas->full_total) {SNESFASInterpolatedCoarseSolution(snes,X,X);}
812: else {SNESFASCoarseCorrection(snes,X,F,X);}
813: fas->full_total = PETSC_FALSE;
814: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
815: if (fas->level != 1) next->max_its -= 1;
816: } else {
817: /* The smoother on the coarse level is the coarse solver */
818: SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
819: }
820: fas->full_stage = 1;
821: } else if (fas->full_stage == 1) {
822: if (snes->iter == 0) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
823: if (next) {
824: SNESFASCoarseCorrection(snes,X,F,X);
825: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
826: }
827: }
828: /* final v-cycle */
829: if (isFine) {
830: if (next) {
831: SNESFASCoarseCorrection(snes,X,F,X);
832: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
833: }
834: }
835: return(0);
836: }
838: static PetscErrorCode SNESFASCycle_Kaskade(SNES snes, Vec X)
839: {
841: Vec F,B;
842: SNES next;
845: F = snes->vec_func;
846: B = snes->vec_rhs;
847: SNESFASCycleGetCorrection(snes,&next);
848: if (next) {
849: SNESFASCoarseCorrection(snes,X,F,X);
850: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
851: } else {
852: SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
853: }
854: return(0);
855: }
857: PetscBool SNEScite = PETSC_FALSE;
858: const char SNESCitation[] = "@techreport{pbmkbsxt2012,\n"
859: " title = {Composing Scalable Nonlinear Algebraic Solvers},\n"
860: " author = {Peter Brune and Mathew Knepley and Barry Smith and Xuemin Tu},\n"
861: " year = 2013,\n"
862: " type = Preprint,\n"
863: " number = {ANL/MCS-P2010-0112},\n"
864: " institution = {Argonne National Laboratory}\n}\n";
866: static PetscErrorCode SNESSolve_FAS(SNES snes)
867: {
869: PetscInt i;
870: Vec X, F;
871: PetscReal fnorm;
872: SNES_FAS *fas = (SNES_FAS*)snes->data,*ffas;
873: DM dm;
874: PetscBool isFine;
877: if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
879: PetscCitationsRegister(SNESCitation,&SNEScite);
880: snes->reason = SNES_CONVERGED_ITERATING;
881: X = snes->vec_sol;
882: F = snes->vec_func;
884: SNESFASCycleIsFine(snes, &isFine);
885: /* norm setup */
886: PetscObjectSAWsTakeAccess((PetscObject)snes);
887: snes->iter = 0;
888: snes->norm = 0;
889: PetscObjectSAWsGrantAccess((PetscObject)snes);
890: if (!snes->vec_func_init_set) {
891: if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,snes,0,0,0);}
892: SNESComputeFunction(snes,X,F);
893: if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,snes,0,0,0);}
894: } else snes->vec_func_init_set = PETSC_FALSE;
896: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
897: SNESCheckFunctionNorm(snes,fnorm);
898: PetscObjectSAWsTakeAccess((PetscObject)snes);
899: snes->norm = fnorm;
900: PetscObjectSAWsGrantAccess((PetscObject)snes);
901: SNESLogConvergenceHistory(snes,fnorm,0);
902: SNESMonitor(snes,snes->iter,fnorm);
904: /* test convergence */
905: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
906: if (snes->reason) return(0);
908: if (isFine) {
909: /* propagate scale-dependent data up the hierarchy */
910: SNESGetDM(snes,&dm);
911: for (ffas=fas; ffas->next; ffas=(SNES_FAS*)ffas->next->data) {
912: DM dmcoarse;
913: SNESGetDM(ffas->next,&dmcoarse);
914: DMRestrict(dm,ffas->restrct,ffas->rscale,ffas->inject,dmcoarse);
915: dm = dmcoarse;
916: }
917: }
919: for (i = 0; i < snes->max_its; i++) {
920: /* Call general purpose update function */
921: if (snes->ops->update) {
922: (*snes->ops->update)(snes, snes->iter);
923: }
925: if (fas->fastype == SNES_FAS_MULTIPLICATIVE) {
926: SNESFASCycle_Multiplicative(snes, X);
927: } else if (fas->fastype == SNES_FAS_ADDITIVE) {
928: SNESFASCycle_Additive(snes, X);
929: } else if (fas->fastype == SNES_FAS_FULL) {
930: SNESFASCycle_Full(snes, X);
931: } else if (fas->fastype == SNES_FAS_KASKADE) {
932: SNESFASCycle_Kaskade(snes, X);
933: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported FAS type");
935: /* check for FAS cycle divergence */
936: if (snes->reason != SNES_CONVERGED_ITERATING) return(0);
938: /* Monitor convergence */
939: PetscObjectSAWsTakeAccess((PetscObject)snes);
940: snes->iter = i+1;
941: PetscObjectSAWsGrantAccess((PetscObject)snes);
942: SNESLogConvergenceHistory(snes,snes->norm,0);
943: SNESMonitor(snes,snes->iter,snes->norm);
944: /* Test for convergence */
945: if (isFine) {
946: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,snes->norm,&snes->reason,snes->cnvP);
947: if (snes->reason) break;
948: }
949: }
950: if (i == snes->max_its) {
951: PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", i);
952: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
953: }
954: return(0);
955: }
957: /*MC
959: SNESFAS - Full Approximation Scheme nonlinear multigrid solver.
961: The nonlinear problem is solved by correction using coarse versions
962: of the nonlinear problem. This problem is perturbed so that a projected
963: solution of the fine problem elicits no correction from the coarse problem.
965: Options Database:
966: + -snes_fas_levels - The number of levels
967: . -snes_fas_cycles<1> - The number of cycles -- 1 for V, 2 for W
968: . -snes_fas_type<additive,multiplicative,full,kaskade> - Additive or multiplicative cycle
969: . -snes_fas_galerkin<PETSC_FALSE> - Form coarse problems by projection back upon the fine problem
970: . -snes_fas_smoothup<1> - The number of iterations of the post-smoother
971: . -snes_fas_smoothdown<1> - The number of iterations of the pre-smoother
972: . -snes_fas_monitor - Monitor progress of all of the levels
973: . -snes_fas_full_downsweep<PETSC_FALSE> - call the downsmooth on the initial downsweep of full FAS
974: . -fas_levels_snes_ - SNES options for all smoothers
975: . -fas_levels_cycle_snes_ - SNES options for all cycles
976: . -fas_levels_i_snes_ - SNES options for the smoothers on level i
977: . -fas_levels_i_cycle_snes_ - SNES options for the cycle on level i
978: - -fas_coarse_snes_ - SNES options for the coarsest smoother
980: Notes:
981: The organization of the FAS solver is slightly different from the organization of PCMG
982: As each level has smoother SNES instances(down and potentially up) and a cycle SNES instance.
983: The cycle SNES instance may be used for monitoring convergence on a particular level.
985: Level: beginner
987: References:
988: . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
989: SIAM Review, 57(4), 2015
991: .seealso: PCMG, SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types)
992: M*/
994: PETSC_EXTERN PetscErrorCode SNESCreate_FAS(SNES snes)
995: {
996: SNES_FAS *fas;
1000: snes->ops->destroy = SNESDestroy_FAS;
1001: snes->ops->setup = SNESSetUp_FAS;
1002: snes->ops->setfromoptions = SNESSetFromOptions_FAS;
1003: snes->ops->view = SNESView_FAS;
1004: snes->ops->solve = SNESSolve_FAS;
1005: snes->ops->reset = SNESReset_FAS;
1007: snes->usesksp = PETSC_FALSE;
1008: snes->usesnpc = PETSC_FALSE;
1010: if (!snes->tolerancesset) {
1011: snes->max_funcs = 30000;
1012: snes->max_its = 10000;
1013: }
1015: snes->alwayscomputesfinalresidual = PETSC_TRUE;
1017: PetscNewLog(snes,&fas);
1019: snes->data = (void*) fas;
1020: fas->level = 0;
1021: fas->levels = 1;
1022: fas->n_cycles = 1;
1023: fas->max_up_it = 1;
1024: fas->max_down_it = 1;
1025: fas->smoothu = NULL;
1026: fas->smoothd = NULL;
1027: fas->next = NULL;
1028: fas->previous = NULL;
1029: fas->fine = snes;
1030: fas->interpolate = NULL;
1031: fas->restrct = NULL;
1032: fas->inject = NULL;
1033: fas->usedmfornumberoflevels = PETSC_FALSE;
1034: fas->fastype = SNES_FAS_MULTIPLICATIVE;
1035: fas->full_downsweep = PETSC_FALSE;
1036: fas->full_total = PETSC_FALSE;
1038: fas->eventsmoothsetup = 0;
1039: fas->eventsmoothsolve = 0;
1040: fas->eventresidual = 0;
1041: fas->eventinterprestrict = 0;
1042: return(0);
1043: }