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) SNESFASCycleCreateSmoother_Private(snes, &fas->smoothd);
98: if (snes->dm) {
99: /* set the smoother DMs properly */
100: if (fas->smoothu) SNESSetDM(fas->smoothu, snes->dm);
101: SNESSetDM(fas->smoothd, snes->dm);
102: /* construct EVERYTHING from the DM -- including the progressive set of smoothers */
103: if (next) {
104: /* for now -- assume the DM and the evaluation functions have been set externally */
105: if (!next->dm) {
106: DMCoarsen(snes->dm, PetscObjectComm((PetscObject)next), &next->dm);
107: SNESSetDM(next, next->dm);
108: }
109: /* set the interpolation and restriction from the DM */
110: if (!fas->interpolate) {
111: DMCreateInterpolation(next->dm, snes->dm, &fas->interpolate, &fas->rscale);
112: if (!fas->restrct) {
113: DMHasCreateRestriction(next->dm, &hasCreateRestriction);
114: /* DM can create restrictions, use that */
115: if (hasCreateRestriction) {
116: DMCreateRestriction(next->dm, snes->dm, &fas->restrct);
117: } else {
118: PetscObjectReference((PetscObject)fas->interpolate);
119: fas->restrct = fas->interpolate;
120: }
121: }
122: }
123: /* set the injection from the DM */
124: if (!fas->inject) {
125: DMHasCreateInjection(next->dm, &hasCreateInjection);
126: if (hasCreateInjection) DMCreateInjection(next->dm, snes->dm, &fas->inject);
127: }
128: }
129: }
131: /*pass the smoother, function, and jacobian up to the next level if it's not user set already */
132: if (fas->galerkin) {
133: if (next) SNESSetFunction(next, NULL, SNESFASGalerkinFunctionDefault, next);
134: if (fas->smoothd && fas->level != fas->levels - 1) SNESSetFunction(fas->smoothd, NULL, SNESFASGalerkinFunctionDefault, snes);
135: if (fas->smoothu && fas->level != fas->levels - 1) SNESSetFunction(fas->smoothu, NULL, SNESFASGalerkinFunctionDefault, snes);
136: }
138: /* sets the down (pre) smoother's default norm and sets it from options */
139: if (fas->smoothd) {
140: if (fas->level == 0 && fas->levels != 1) {
141: SNESSetNormSchedule(fas->smoothd, SNES_NORM_NONE);
142: } else {
143: SNESSetNormSchedule(fas->smoothd, SNES_NORM_FINAL_ONLY);
144: }
145: SNESFASCycleSetUpSmoother_Private(snes, fas->smoothd);
146: }
148: /* sets the up (post) smoother's default norm and sets it from options */
149: if (fas->smoothu) {
150: if (fas->level != fas->levels - 1) {
151: SNESSetNormSchedule(fas->smoothu, SNES_NORM_NONE);
152: } else {
153: SNESSetNormSchedule(fas->smoothu, SNES_NORM_FINAL_ONLY);
154: }
155: SNESFASCycleSetUpSmoother_Private(snes, fas->smoothu);
156: }
158: if (next) {
159: /* gotta set up the solution vector for this to work */
160: if (!next->vec_sol) SNESFASCreateCoarseVec(snes, &next->vec_sol);
161: if (!next->vec_rhs) SNESFASCreateCoarseVec(snes, &next->vec_rhs);
162: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next);
163: SNESFASSetUpLineSearch_Private(snes, next);
164: SNESSetUp(next);
165: }
167: /* setup FAS work vectors */
168: if (fas->galerkin) {
169: VecDuplicate(snes->vec_sol, &fas->Xg);
170: VecDuplicate(snes->vec_sol, &fas->Fg);
171: }
172: return 0;
173: }
175: static PetscErrorCode SNESSetFromOptions_FAS(SNES snes, PetscOptionItems *PetscOptionsObject)
176: {
177: SNES_FAS *fas = (SNES_FAS *)snes->data;
178: PetscInt levels = 1;
179: PetscBool flg = PETSC_FALSE, upflg = PETSC_FALSE, downflg = PETSC_FALSE, monflg = PETSC_FALSE, galerkinflg = PETSC_FALSE, continuationflg = PETSC_FALSE;
180: SNESFASType fastype;
181: const char *optionsprefix;
182: SNESLineSearch linesearch;
183: PetscInt m, n_up, n_down;
184: SNES next;
185: PetscBool isFine;
187: SNESFASCycleIsFine(snes, &isFine);
188: PetscOptionsHeadBegin(PetscOptionsObject, "SNESFAS Options-----------------------------------");
190: /* number of levels -- only process most options on the finest level */
191: if (isFine) {
192: PetscOptionsInt("-snes_fas_levels", "Number of Levels", "SNESFASSetLevels", levels, &levels, &flg);
193: if (!flg && snes->dm) {
194: DMGetRefineLevel(snes->dm, &levels);
195: levels++;
196: fas->usedmfornumberoflevels = PETSC_TRUE;
197: }
198: SNESFASSetLevels(snes, levels, NULL);
199: fastype = fas->fastype;
200: PetscOptionsEnum("-snes_fas_type", "FAS correction type", "SNESFASSetType", SNESFASTypes, (PetscEnum)fastype, (PetscEnum *)&fastype, &flg);
201: if (flg) SNESFASSetType(snes, fastype);
203: SNESGetOptionsPrefix(snes, &optionsprefix);
204: PetscOptionsInt("-snes_fas_cycles", "Number of cycles", "SNESFASSetCycles", fas->n_cycles, &m, &flg);
205: if (flg) SNESFASSetCycles(snes, m);
206: PetscOptionsBool("-snes_fas_continuation", "Corrected grid-sequence continuation", "SNESFASSetContinuation", fas->continuation, &continuationflg, &flg);
207: if (flg) SNESFASSetContinuation(snes, continuationflg);
209: PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin", "SNESFASSetGalerkin", fas->galerkin, &galerkinflg, &flg);
210: if (flg) SNESFASSetGalerkin(snes, galerkinflg);
212: if (fas->fastype == SNES_FAS_FULL) {
213: PetscOptionsBool("-snes_fas_full_downsweep", "Smooth on the initial down sweep for full FAS cycles", "SNESFASFullSetDownSweep", fas->full_downsweep, &fas->full_downsweep, &flg);
214: if (flg) SNESFASFullSetDownSweep(snes, fas->full_downsweep);
215: 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);
216: if (flg) SNESFASFullSetTotal(snes, fas->full_total);
217: }
219: PetscOptionsInt("-snes_fas_smoothup", "Number of post-smoothing steps", "SNESFASSetNumberSmoothUp", fas->max_up_it, &n_up, &upflg);
221: PetscOptionsInt("-snes_fas_smoothdown", "Number of pre-smoothing steps", "SNESFASSetNumberSmoothDown", fas->max_down_it, &n_down, &downflg);
223: {
224: PetscViewer viewer;
225: PetscViewerFormat format;
226: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_fas_monitor", &viewer, &format, &monflg);
227: if (monflg) {
228: PetscViewerAndFormat *vf;
229: PetscViewerAndFormatCreate(viewer, format, &vf);
230: PetscObjectDereference((PetscObject)viewer);
231: SNESFASSetMonitor(snes, vf, PETSC_TRUE);
232: }
233: }
234: flg = PETSC_FALSE;
235: monflg = PETSC_TRUE;
236: PetscOptionsBool("-snes_fas_log", "Log times for each FAS level", "SNESFASSetLog", monflg, &monflg, &flg);
237: if (flg) SNESFASSetLog(snes, monflg);
238: }
240: PetscOptionsHeadEnd();
242: /* setup from the determined types if there is no pointwise procedure or smoother defined */
243: if (upflg) SNESFASSetNumberSmoothUp(snes, n_up);
244: if (downflg) SNESFASSetNumberSmoothDown(snes, n_down);
246: /* set up the default line search for coarse grid corrections */
247: if (fas->fastype == SNES_FAS_ADDITIVE) {
248: if (!snes->linesearch) {
249: SNESGetLineSearch(snes, &linesearch);
250: SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
251: }
252: }
254: /* recursive option setting for the smoothers */
255: SNESFASCycleGetCorrection(snes, &next);
256: if (next) SNESSetFromOptions(next);
257: return 0;
258: }
260: #include <petscdraw.h>
261: static PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer)
262: {
263: SNES_FAS *fas = (SNES_FAS *)snes->data;
264: PetscBool isFine, iascii, isdraw;
265: PetscInt i;
266: SNES smoothu, smoothd, levelsnes;
268: SNESFASCycleIsFine(snes, &isFine);
269: if (isFine) {
270: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
271: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
272: if (iascii) {
273: PetscViewerASCIIPrintf(viewer, " type is %s, levels=%" PetscInt_FMT ", cycles=%" PetscInt_FMT "\n", SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles);
274: if (fas->galerkin) {
275: PetscViewerASCIIPrintf(viewer, " Using Galerkin computed coarse grid function evaluation\n");
276: } else {
277: PetscViewerASCIIPrintf(viewer, " Not using Galerkin computed coarse grid function evaluation\n");
278: }
279: for (i = 0; i < fas->levels; i++) {
280: SNESFASGetCycleSNES(snes, i, &levelsnes);
281: SNESFASCycleGetSmootherUp(levelsnes, &smoothu);
282: SNESFASCycleGetSmootherDown(levelsnes, &smoothd);
283: if (!i) {
284: PetscViewerASCIIPrintf(viewer, " Coarse grid solver -- level %" PetscInt_FMT " -------------------------------\n", i);
285: } else {
286: PetscViewerASCIIPrintf(viewer, " Down solver (pre-smoother) on level %" PetscInt_FMT " -------------------------------\n", i);
287: }
288: PetscViewerASCIIPushTab(viewer);
289: if (smoothd) {
290: SNESView(smoothd, viewer);
291: } else {
292: PetscViewerASCIIPrintf(viewer, "Not yet available\n");
293: }
294: PetscViewerASCIIPopTab(viewer);
295: if (i && (smoothd == smoothu)) {
296: PetscViewerASCIIPrintf(viewer, " Up solver (post-smoother) same as down solver (pre-smoother)\n");
297: } else if (i) {
298: PetscViewerASCIIPrintf(viewer, " Up solver (post-smoother) on level %" PetscInt_FMT " -------------------------------\n", i);
299: PetscViewerASCIIPushTab(viewer);
300: if (smoothu) {
301: SNESView(smoothu, viewer);
302: } else {
303: PetscViewerASCIIPrintf(viewer, "Not yet available\n");
304: }
305: PetscViewerASCIIPopTab(viewer);
306: }
307: }
308: } else if (isdraw) {
309: PetscDraw draw;
310: PetscReal x, w, y, bottom, th, wth;
311: SNES_FAS *curfas = fas;
312: PetscViewerDrawGetDraw(viewer, 0, &draw);
313: PetscDrawGetCurrentPoint(draw, &x, &y);
314: PetscDrawStringGetSize(draw, &wth, &th);
315: bottom = y - th;
316: while (curfas) {
317: if (!curfas->smoothu) {
318: PetscDrawPushCurrentPoint(draw, x, bottom);
319: if (curfas->smoothd) SNESView(curfas->smoothd, viewer);
320: PetscDrawPopCurrentPoint(draw);
321: } else {
322: w = 0.5 * PetscMin(1.0 - x, x);
323: PetscDrawPushCurrentPoint(draw, x - w, bottom);
324: if (curfas->smoothd) SNESView(curfas->smoothd, viewer);
325: PetscDrawPopCurrentPoint(draw);
326: PetscDrawPushCurrentPoint(draw, x + w, bottom);
327: if (curfas->smoothu) SNESView(curfas->smoothu, viewer);
328: PetscDrawPopCurrentPoint(draw);
329: }
330: /* this is totally bogus but we have no way of knowing how low the previous one was draw to */
331: bottom -= 5 * th;
332: if (curfas->next) curfas = (SNES_FAS *)curfas->next->data;
333: else curfas = NULL;
334: }
335: }
336: }
337: return 0;
338: }
340: /*
341: Defines the action of the downsmoother
342: */
343: static PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
344: {
345: SNESConvergedReason reason;
346: Vec FPC;
347: SNES smoothd;
348: PetscBool flg;
349: SNES_FAS *fas = (SNES_FAS *)snes->data;
351: SNESFASCycleGetSmootherDown(snes, &smoothd);
352: SNESSetInitialFunction(smoothd, F);
353: if (fas->eventsmoothsolve) PetscLogEventBegin(fas->eventsmoothsolve, smoothd, B, X, 0);
354: SNESSolve(smoothd, B, X);
355: if (fas->eventsmoothsolve) PetscLogEventEnd(fas->eventsmoothsolve, smoothd, B, X, 0);
356: /* check convergence reason for the smoother */
357: SNESGetConvergedReason(smoothd, &reason);
358: if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
359: snes->reason = SNES_DIVERGED_INNER;
360: return 0;
361: }
363: SNESGetFunction(smoothd, &FPC, NULL, NULL);
364: SNESGetAlwaysComputesFinalResidual(smoothd, &flg);
365: if (!flg) SNESComputeFunction(smoothd, X, FPC);
366: VecCopy(FPC, F);
367: if (fnorm) VecNorm(F, NORM_2, fnorm);
368: return 0;
369: }
371: /*
372: Defines the action of the upsmoother
373: */
374: static PetscErrorCode SNESFASUpSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
375: {
376: SNESConvergedReason reason;
377: Vec FPC;
378: SNES smoothu;
379: PetscBool flg;
380: SNES_FAS *fas = (SNES_FAS *)snes->data;
382: SNESFASCycleGetSmootherUp(snes, &smoothu);
383: if (fas->eventsmoothsolve) PetscLogEventBegin(fas->eventsmoothsolve, smoothu, 0, 0, 0);
384: SNESSolve(smoothu, B, X);
385: if (fas->eventsmoothsolve) PetscLogEventEnd(fas->eventsmoothsolve, smoothu, 0, 0, 0);
386: /* check convergence reason for the smoother */
387: SNESGetConvergedReason(smoothu, &reason);
388: if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
389: snes->reason = SNES_DIVERGED_INNER;
390: return 0;
391: }
392: SNESGetFunction(smoothu, &FPC, NULL, NULL);
393: SNESGetAlwaysComputesFinalResidual(smoothu, &flg);
394: if (!flg) SNESComputeFunction(smoothu, X, FPC);
395: VecCopy(FPC, F);
396: if (fnorm) VecNorm(F, NORM_2, fnorm);
397: return 0;
398: }
400: /*@
401: SNESFASCreateCoarseVec - create `Vec` corresponding to a state vector on one level coarser than current level
403: Collective
405: Input Parameter:
406: . snes - `SNESFAS` object
408: Output Parameter:
409: . Xcoarse - vector on level one coarser than snes
411: Level: developer
413: .seealso: `SNESFASSetRestriction()`, `SNESFASRestrict()`
414: @*/
415: PetscErrorCode SNESFASCreateCoarseVec(SNES snes, Vec *Xcoarse)
416: {
417: SNES_FAS *fas;
421: fas = (SNES_FAS *)snes->data;
422: if (fas->rscale) {
423: VecDuplicate(fas->rscale, Xcoarse);
424: } else if (fas->inject) {
425: MatCreateVecs(fas->inject, Xcoarse, NULL);
426: } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must set restriction or injection");
427: return 0;
428: }
430: /*@
431: SNESFASRestrict - restrict a `Vec` to the next coarser level
433: Collective
435: Input Parameters:
436: + fine - `SNES` from which to restrict
437: - Xfine - vector to restrict
439: Output Parameter:
440: . Xcoarse - result of restriction
442: Level: developer
444: .seealso: `SNES`, `SNESFAS`, `SNESFASSetRestriction()`, `SNESFASSetInjection()`
445: @*/
446: PetscErrorCode SNESFASRestrict(SNES fine, Vec Xfine, Vec Xcoarse)
447: {
448: SNES_FAS *fas;
453: fas = (SNES_FAS *)fine->data;
454: if (fas->inject) {
455: MatRestrict(fas->inject, Xfine, Xcoarse);
456: } else {
457: MatRestrict(fas->restrct, Xfine, Xcoarse);
458: VecPointwiseMult(Xcoarse, fas->rscale, Xcoarse);
459: }
460: return 0;
461: }
463: /*
465: Performs a variant of FAS using the interpolated total coarse solution
467: fine problem: F(x) = b
468: coarse problem: F^c(x^c) = Rb, Initial guess Rx
469: interpolated solution: x^f = I x^c (total solution interpolation
471: */
472: static PetscErrorCode SNESFASInterpolatedCoarseSolution(SNES snes, Vec X, Vec X_new)
473: {
474: Vec X_c, B_c;
475: SNESConvergedReason reason;
476: SNES next;
477: Mat restrct, interpolate;
478: SNES_FAS *fasc;
480: SNESFASCycleGetCorrection(snes, &next);
481: if (next) {
482: fasc = (SNES_FAS *)next->data;
484: SNESFASCycleGetRestriction(snes, &restrct);
485: SNESFASCycleGetInterpolation(snes, &interpolate);
487: X_c = next->vec_sol;
489: if (fasc->eventinterprestrict) PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0);
490: /* restrict the total solution: Rb */
491: SNESFASRestrict(snes, X, X_c);
492: B_c = next->vec_rhs;
493: if (snes->vec_rhs) {
494: /* restrict the total rhs defect: Rb */
495: MatRestrict(restrct, snes->vec_rhs, B_c);
496: } else {
497: VecSet(B_c, 0.);
498: }
499: if (fasc->eventinterprestrict) PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0);
501: SNESSolve(next, B_c, X_c);
502: SNESGetConvergedReason(next, &reason);
503: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
504: snes->reason = SNES_DIVERGED_INNER;
505: return 0;
506: }
507: /* x^f <- Ix^c*/
508: DM dmc, dmf;
510: SNESGetDM(next, &dmc);
511: SNESGetDM(snes, &dmf);
512: if (fasc->eventinterprestrict) PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0);
513: DMInterpolateSolution(dmc, dmf, interpolate, X_c, X_new);
514: if (fasc->eventinterprestrict) PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0);
515: PetscObjectSetName((PetscObject)X_c, "Coarse solution");
516: VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view");
517: PetscObjectSetName((PetscObject)X_new, "Updated Fine solution");
518: VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view");
519: }
520: return 0;
521: }
523: /*
525: Performs the FAS coarse correction as:
527: fine problem: F(x) = b
528: coarse problem: F^c(x^c) = b^c
530: b^c = F^c(Rx) - R(F(x) - b)
532: */
533: PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new)
534: {
535: Vec X_c, Xo_c, F_c, B_c;
536: SNESConvergedReason reason;
537: SNES next;
538: Mat restrct, interpolate;
539: SNES_FAS *fasc;
541: SNESFASCycleGetCorrection(snes, &next);
542: if (next) {
543: fasc = (SNES_FAS *)next->data;
545: SNESFASCycleGetRestriction(snes, &restrct);
546: SNESFASCycleGetInterpolation(snes, &interpolate);
548: X_c = next->vec_sol;
549: Xo_c = next->work[0];
550: F_c = next->vec_func;
551: B_c = next->vec_rhs;
553: if (fasc->eventinterprestrict) PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0);
554: SNESFASRestrict(snes, X, Xo_c);
555: /* restrict the defect: R(F(x) - b) */
556: MatRestrict(restrct, F, B_c);
557: if (fasc->eventinterprestrict) PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0);
559: if (fasc->eventresidual) PetscLogEventBegin(fasc->eventresidual, next, 0, 0, 0);
560: /* F_c = F^c(Rx) - R(F(x) - b) since the second term was sitting in next->vec_rhs */
561: SNESComputeFunction(next, Xo_c, F_c);
562: if (fasc->eventresidual) PetscLogEventEnd(fasc->eventresidual, next, 0, 0, 0);
564: /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */
565: VecCopy(B_c, X_c);
566: VecCopy(F_c, B_c);
567: VecCopy(X_c, F_c);
568: /* set initial guess of the coarse problem to the projected fine solution */
569: VecCopy(Xo_c, X_c);
571: /* recurse to the next level */
572: SNESSetInitialFunction(next, F_c);
573: SNESSolve(next, B_c, X_c);
574: SNESGetConvergedReason(next, &reason);
575: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
576: snes->reason = SNES_DIVERGED_INNER;
577: return 0;
578: }
579: /* correct as x <- x + I(x^c - Rx)*/
580: VecAXPY(X_c, -1.0, Xo_c);
582: if (fasc->eventinterprestrict) PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0);
583: MatInterpolateAdd(interpolate, X_c, X, X_new);
584: if (fasc->eventinterprestrict) PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0);
585: PetscObjectSetName((PetscObject)X_c, "Coarse correction");
586: VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view");
587: PetscObjectSetName((PetscObject)X_new, "Updated Fine solution");
588: VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view");
589: }
590: return 0;
591: }
593: /*
595: The additive cycle looks like:
597: xhat = x
598: xhat = dS(x, b)
599: x = coarsecorrection(xhat, b_d)
600: x = x + nu*(xhat - x);
601: (optional) x = uS(x, b)
603: With the coarse RHS (defect correction) as below.
605: */
606: static PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X)
607: {
608: Vec F, B, Xhat;
609: Vec X_c, Xo_c, F_c, B_c;
610: SNESConvergedReason reason;
611: PetscReal xnorm, fnorm, ynorm;
612: SNESLineSearchReason lsresult;
613: SNES next;
614: Mat restrct, interpolate;
615: SNES_FAS *fas = (SNES_FAS *)snes->data, *fasc;
617: SNESFASCycleGetCorrection(snes, &next);
618: F = snes->vec_func;
619: B = snes->vec_rhs;
620: Xhat = snes->work[1];
621: VecCopy(X, Xhat);
622: /* recurse first */
623: if (next) {
624: fasc = (SNES_FAS *)next->data;
625: SNESFASCycleGetRestriction(snes, &restrct);
626: SNESFASCycleGetInterpolation(snes, &interpolate);
627: if (fas->eventresidual) PetscLogEventBegin(fas->eventresidual, snes, 0, 0, 0);
628: SNESComputeFunction(snes, Xhat, F);
629: if (fas->eventresidual) PetscLogEventEnd(fas->eventresidual, snes, 0, 0, 0);
630: VecNorm(F, NORM_2, &fnorm);
631: X_c = next->vec_sol;
632: Xo_c = next->work[0];
633: F_c = next->vec_func;
634: B_c = next->vec_rhs;
636: SNESFASRestrict(snes, Xhat, Xo_c);
637: /* restrict the defect */
638: MatRestrict(restrct, F, B_c);
640: /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */
641: if (fasc->eventresidual) PetscLogEventBegin(fasc->eventresidual, next, 0, 0, 0);
642: SNESComputeFunction(next, Xo_c, F_c);
643: if (fasc->eventresidual) PetscLogEventEnd(fasc->eventresidual, next, 0, 0, 0);
644: VecCopy(B_c, X_c);
645: VecCopy(F_c, B_c);
646: VecCopy(X_c, F_c);
647: /* set initial guess of the coarse problem to the projected fine solution */
648: VecCopy(Xo_c, X_c);
650: /* recurse */
651: SNESSetInitialFunction(next, F_c);
652: SNESSolve(next, B_c, X_c);
654: /* smooth on this level */
655: SNESFASDownSmooth_Private(snes, B, X, F, &fnorm);
657: SNESGetConvergedReason(next, &reason);
658: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
659: snes->reason = SNES_DIVERGED_INNER;
660: return 0;
661: }
663: /* correct as x <- x + I(x^c - Rx)*/
664: VecAYPX(X_c, -1.0, Xo_c);
665: MatInterpolate(interpolate, X_c, Xhat);
667: /* additive correction of the coarse direction*/
668: SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat);
669: SNESLineSearchGetReason(snes->linesearch, &lsresult);
670: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm);
671: if (lsresult) {
672: if (++snes->numFailures >= snes->maxFailures) {
673: snes->reason = SNES_DIVERGED_LINE_SEARCH;
674: return 0;
675: }
676: }
677: } else {
678: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
679: }
680: return 0;
681: }
683: /*
685: Defines the FAS cycle as:
687: fine problem: F(x) = b
688: coarse problem: F^c(x) = b^c
690: b^c = F^c(Rx) - R(F(x) - b)
692: correction:
694: x = x + I(x^c - Rx)
696: */
697: static PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X)
698: {
699: Vec F, B;
700: SNES next;
702: F = snes->vec_func;
703: B = snes->vec_rhs;
704: /* pre-smooth -- just update using the pre-smoother */
705: SNESFASCycleGetCorrection(snes, &next);
706: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
707: if (next) {
708: SNESFASCoarseCorrection(snes, X, F, X);
709: SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
710: }
711: return 0;
712: }
714: static PetscErrorCode SNESFASCycleSetupPhase_Full(SNES snes)
715: {
716: SNES next;
717: SNES_FAS *fas = (SNES_FAS *)snes->data;
718: PetscBool isFine;
720: /* pre-smooth -- just update using the pre-smoother */
721: SNESFASCycleIsFine(snes, &isFine);
722: SNESFASCycleGetCorrection(snes, &next);
723: fas->full_stage = 0;
724: if (next) SNESFASCycleSetupPhase_Full(next);
725: return 0;
726: }
728: static PetscErrorCode SNESFASCycle_Full(SNES snes, Vec X)
729: {
730: Vec F, B;
731: SNES_FAS *fas = (SNES_FAS *)snes->data;
732: PetscBool isFine;
733: SNES next;
735: F = snes->vec_func;
736: B = snes->vec_rhs;
737: SNESFASCycleIsFine(snes, &isFine);
738: SNESFASCycleGetCorrection(snes, &next);
740: if (isFine) SNESFASCycleSetupPhase_Full(snes);
742: if (fas->full_stage == 0) {
743: /* downsweep */
744: if (next) {
745: if (fas->level != 1) next->max_its += 1;
746: if (fas->full_downsweep) SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
747: fas->full_downsweep = PETSC_TRUE;
748: if (fas->full_total) SNESFASInterpolatedCoarseSolution(snes, X, X);
749: else SNESFASCoarseCorrection(snes, X, F, X);
750: fas->full_total = PETSC_FALSE;
751: SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
752: if (fas->level != 1) next->max_its -= 1;
753: } else {
754: /* The smoother on the coarse level is the coarse solver */
755: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
756: }
757: fas->full_stage = 1;
758: } else if (fas->full_stage == 1) {
759: if (snes->iter == 0) SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
760: if (next) {
761: SNESFASCoarseCorrection(snes, X, F, X);
762: SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
763: }
764: }
765: /* final v-cycle */
766: if (isFine) {
767: if (next) {
768: SNESFASCoarseCorrection(snes, X, F, X);
769: SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
770: }
771: }
772: return 0;
773: }
775: static PetscErrorCode SNESFASCycle_Kaskade(SNES snes, Vec X)
776: {
777: Vec F, B;
778: SNES next;
780: F = snes->vec_func;
781: B = snes->vec_rhs;
782: SNESFASCycleGetCorrection(snes, &next);
783: if (next) {
784: SNESFASCoarseCorrection(snes, X, F, X);
785: SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
786: } else {
787: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
788: }
789: return 0;
790: }
792: PetscBool SNEScite = PETSC_FALSE;
793: const char SNESCitation[] = "@techreport{pbmkbsxt2012,\n"
794: " title = {Composing Scalable Nonlinear Algebraic Solvers},\n"
795: " author = {Peter Brune and Mathew Knepley and Barry Smith and Xuemin Tu},\n"
796: " year = 2013,\n"
797: " type = Preprint,\n"
798: " number = {ANL/MCS-P2010-0112},\n"
799: " institution = {Argonne National Laboratory}\n}\n";
801: static PetscErrorCode SNESSolve_FAS(SNES snes)
802: {
803: PetscInt i;
804: Vec X, F;
805: PetscReal fnorm;
806: SNES_FAS *fas = (SNES_FAS *)snes->data, *ffas;
807: DM dm;
808: PetscBool isFine;
812: PetscCitationsRegister(SNESCitation, &SNEScite);
813: snes->reason = SNES_CONVERGED_ITERATING;
814: X = snes->vec_sol;
815: F = snes->vec_func;
817: SNESFASCycleIsFine(snes, &isFine);
818: /* norm setup */
819: PetscObjectSAWsTakeAccess((PetscObject)snes);
820: snes->iter = 0;
821: snes->norm = 0;
822: PetscObjectSAWsGrantAccess((PetscObject)snes);
823: if (!snes->vec_func_init_set) {
824: if (fas->eventresidual) PetscLogEventBegin(fas->eventresidual, snes, 0, 0, 0);
825: SNESComputeFunction(snes, X, F);
826: if (fas->eventresidual) PetscLogEventEnd(fas->eventresidual, snes, 0, 0, 0);
827: } else snes->vec_func_init_set = PETSC_FALSE;
829: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
830: SNESCheckFunctionNorm(snes, fnorm);
831: PetscObjectSAWsTakeAccess((PetscObject)snes);
832: snes->norm = fnorm;
833: PetscObjectSAWsGrantAccess((PetscObject)snes);
834: SNESLogConvergenceHistory(snes, fnorm, 0);
835: SNESMonitor(snes, snes->iter, fnorm);
837: /* test convergence */
838: PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
839: if (snes->reason) return 0;
841: if (isFine) {
842: /* propagate scale-dependent data up the hierarchy */
843: SNESGetDM(snes, &dm);
844: for (ffas = fas; ffas->next; ffas = (SNES_FAS *)ffas->next->data) {
845: DM dmcoarse;
846: SNESGetDM(ffas->next, &dmcoarse);
847: DMRestrict(dm, ffas->restrct, ffas->rscale, ffas->inject, dmcoarse);
848: dm = dmcoarse;
849: }
850: }
852: for (i = 0; i < snes->max_its; i++) {
853: /* Call general purpose update function */
854: PetscTryTypeMethod(snes, update, snes->iter);
856: if (fas->fastype == SNES_FAS_MULTIPLICATIVE) {
857: SNESFASCycle_Multiplicative(snes, X);
858: } else if (fas->fastype == SNES_FAS_ADDITIVE) {
859: SNESFASCycle_Additive(snes, X);
860: } else if (fas->fastype == SNES_FAS_FULL) {
861: SNESFASCycle_Full(snes, X);
862: } else if (fas->fastype == SNES_FAS_KASKADE) {
863: SNESFASCycle_Kaskade(snes, X);
864: } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Unsupported FAS type");
866: /* check for FAS cycle divergence */
867: if (snes->reason != SNES_CONVERGED_ITERATING) return 0;
869: /* Monitor convergence */
870: PetscObjectSAWsTakeAccess((PetscObject)snes);
871: snes->iter = i + 1;
872: PetscObjectSAWsGrantAccess((PetscObject)snes);
873: SNESLogConvergenceHistory(snes, snes->norm, 0);
874: SNESMonitor(snes, snes->iter, snes->norm);
875: /* Test for convergence */
876: if (isFine) {
877: PetscUseTypeMethod(snes, converged, snes->iter, 0.0, 0.0, snes->norm, &snes->reason, snes->cnvP);
878: if (snes->reason) break;
879: }
880: }
881: if (i == snes->max_its) {
882: PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", i);
883: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
884: }
885: return 0;
886: }
888: /*MC
890: SNESFAS - Full Approximation Scheme nonlinear multigrid solver.
892: The nonlinear problem is solved by correction using coarse versions
893: of the nonlinear problem. This problem is perturbed so that a projected
894: solution of the fine problem elicits no correction from the coarse problem.
896: Options Database Keys and Prefixes:
897: + -snes_fas_levels - The number of levels
898: . -snes_fas_cycles<1> - The number of cycles -- 1 for V, 2 for W
899: . -snes_fas_type<additive,multiplicative,full,kaskade> - Additive or multiplicative cycle
900: . -snes_fas_galerkin<`PETSC_FALSE`> - Form coarse problems by projection back upon the fine problem
901: . -snes_fas_smoothup<1> - The number of iterations of the post-smoother
902: . -snes_fas_smoothdown<1> - The number of iterations of the pre-smoother
903: . -snes_fas_monitor - Monitor progress of all of the levels
904: . -snes_fas_full_downsweep<`PETSC_FALSE`> - call the downsmooth on the initial downsweep of full FAS
905: . -fas_levels_snes_ - `SNES` options for all smoothers
906: . -fas_levels_cycle_snes_ - `SNES` options for all cycles
907: . -fas_levels_i_snes_ - `SNES` options for the smoothers on level i
908: . -fas_levels_i_cycle_snes_ - `SNES` options for the cycle on level i
909: - -fas_coarse_snes_ - `SNES` options for the coarsest smoother
911: Note:
912: The organization of the FAS solver is slightly different from the organization of `PCMG`
913: As each level has smoother `SNES` instances(down and potentially up) and a cycle `SNES` instance.
914: The cycle `SNES` instance may be used for monitoring convergence on a particular level.
916: Level: beginner
918: References:
919: . * - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
920: SIAM Review, 57(4), 2015
922: .seealso: `PCMG`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESType`, `SNESFASSetRestriction()`, `SNESFASSetInjection()`,
923: `SNESFASFullGetTotal()`
924: M*/
926: PETSC_EXTERN PetscErrorCode SNESCreate_FAS(SNES snes)
927: {
928: SNES_FAS *fas;
930: snes->ops->destroy = SNESDestroy_FAS;
931: snes->ops->setup = SNESSetUp_FAS;
932: snes->ops->setfromoptions = SNESSetFromOptions_FAS;
933: snes->ops->view = SNESView_FAS;
934: snes->ops->solve = SNESSolve_FAS;
935: snes->ops->reset = SNESReset_FAS;
937: snes->usesksp = PETSC_FALSE;
938: snes->usesnpc = PETSC_FALSE;
940: if (!snes->tolerancesset) {
941: snes->max_funcs = 30000;
942: snes->max_its = 10000;
943: }
945: snes->alwayscomputesfinalresidual = PETSC_TRUE;
947: PetscNew(&fas);
949: snes->data = (void *)fas;
950: fas->level = 0;
951: fas->levels = 1;
952: fas->n_cycles = 1;
953: fas->max_up_it = 1;
954: fas->max_down_it = 1;
955: fas->smoothu = NULL;
956: fas->smoothd = NULL;
957: fas->next = NULL;
958: fas->previous = NULL;
959: fas->fine = snes;
960: fas->interpolate = NULL;
961: fas->restrct = NULL;
962: fas->inject = NULL;
963: fas->usedmfornumberoflevels = PETSC_FALSE;
964: fas->fastype = SNES_FAS_MULTIPLICATIVE;
965: fas->full_downsweep = PETSC_FALSE;
966: fas->full_total = PETSC_FALSE;
968: fas->eventsmoothsetup = 0;
969: fas->eventsmoothsolve = 0;
970: fas->eventresidual = 0;
971: fas->eventinterprestrict = 0;
972: return 0;
973: }