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