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