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