Actual source code: fas.c
petsc-3.7.3 2016-08-01
1: /* Defines the basic SNES object */
2: #include <../src/snes/impls/fas/fasimpls.h> /*I "petscsnes.h" I*/
4: const char *const SNESFASTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","SNESFASType","SNES_FAS",0};
6: extern PetscErrorCode SNESDestroy_FAS(SNES snes);
7: extern PetscErrorCode SNESSetUp_FAS(SNES snes);
8: extern PetscErrorCode SNESSetFromOptions_FAS(PetscOptionItems *PetscOptionsObject,SNES snes);
9: extern PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer);
10: extern PetscErrorCode SNESSolve_FAS(SNES snes);
11: extern PetscErrorCode SNESReset_FAS(SNES snes);
12: extern PetscErrorCode SNESFASGalerkinDefaultFunction(SNES, Vec, Vec, void*);
13: extern PetscErrorCode SNESFASCycleCreateSmoother_Private(SNES, SNES*);
15: /*MC
17: SNESFAS - Full Approximation Scheme nonlinear multigrid solver.
19: The nonlinear problem is solved by correction using coarse versions
20: of the nonlinear problem. This problem is perturbed so that a projected
21: solution of the fine problem elicits no correction from the coarse problem.
23: Options Database:
24: + -snes_fas_levels - The number of levels
25: . -snes_fas_cycles<1> - The number of cycles -- 1 for V, 2 for W
26: . -snes_fas_type<additive,multiplicative,full,kaskade> - Additive or multiplicative cycle
27: . -snes_fas_galerkin<PETSC_FALSE> - Form coarse problems by projection back upon the fine problem
28: . -snes_fas_smoothup<1> - The number of iterations of the post-smoother
29: . -snes_fas_smoothdown<1> - The number of iterations of the pre-smoother
30: . -snes_fas_monitor - Monitor progress of all of the levels
31: . -snes_fas_full_downsweep<PETSC_FALSE> - call the downsmooth on the initial downsweep of full FAS
32: . -fas_levels_snes_ - SNES options for all smoothers
33: . -fas_levels_cycle_snes_ - SNES options for all cycles
34: . -fas_levels_i_snes_ - SNES options for the smoothers on level i
35: . -fas_levels_i_cycle_snes_ - SNES options for the cycle on level i
36: - -fas_coarse_snes_ - SNES options for the coarsest smoother
38: Notes:
39: The organization of the FAS solver is slightly different from the organization of PCMG
40: As each level has smoother SNES instances(down and potentially up) and a cycle SNES instance.
41: The cycle SNES instance may be used for monitoring convergence on a particular level.
43: Level: beginner
45: References:
46: . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
47: SIAM Review, 57(4), 2015
49: .seealso: PCMG, SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types)
50: M*/
54: PETSC_EXTERN PetscErrorCode SNESCreate_FAS(SNES snes)
55: {
56: SNES_FAS *fas;
60: snes->ops->destroy = SNESDestroy_FAS;
61: snes->ops->setup = SNESSetUp_FAS;
62: snes->ops->setfromoptions = SNESSetFromOptions_FAS;
63: snes->ops->view = SNESView_FAS;
64: snes->ops->solve = SNESSolve_FAS;
65: snes->ops->reset = SNESReset_FAS;
67: snes->usesksp = PETSC_FALSE;
68: snes->usespc = PETSC_FALSE;
70: if (!snes->tolerancesset) {
71: snes->max_funcs = 30000;
72: snes->max_its = 10000;
73: }
75: PetscNewLog(snes,&fas);
77: snes->data = (void*) fas;
78: fas->level = 0;
79: fas->levels = 1;
80: fas->n_cycles = 1;
81: fas->max_up_it = 1;
82: fas->max_down_it = 1;
83: fas->smoothu = NULL;
84: fas->smoothd = NULL;
85: fas->next = NULL;
86: fas->previous = NULL;
87: fas->fine = snes;
88: fas->interpolate = NULL;
89: fas->restrct = NULL;
90: fas->inject = NULL;
91: fas->usedmfornumberoflevels = PETSC_FALSE;
92: fas->fastype = SNES_FAS_MULTIPLICATIVE;
93: fas->full_downsweep = PETSC_FALSE;
95: fas->eventsmoothsetup = 0;
96: fas->eventsmoothsolve = 0;
97: fas->eventresidual = 0;
98: fas->eventinterprestrict = 0;
99: return(0);
100: }
104: PetscErrorCode SNESReset_FAS(SNES snes)
105: {
106: PetscErrorCode 0;
107: SNES_FAS * fas = (SNES_FAS*)snes->data;
110: SNESDestroy(&fas->smoothu);
111: SNESDestroy(&fas->smoothd);
112: MatDestroy(&fas->inject);
113: MatDestroy(&fas->interpolate);
114: MatDestroy(&fas->restrct);
115: VecDestroy(&fas->rscale);
116: if (fas->galerkin) {
117: VecDestroy(&fas->Xg);
118: VecDestroy(&fas->Fg);
119: }
120: if (fas->next) {SNESReset(fas->next);}
121: return(0);
122: }
126: PetscErrorCode SNESDestroy_FAS(SNES snes)
127: {
128: SNES_FAS * fas = (SNES_FAS*)snes->data;
129: PetscErrorCode 0;
132: /* recursively resets and then destroys */
133: SNESReset(snes);
134: if (fas->next) {
135: SNESDestroy(&fas->next);
136: }
137: PetscFree(fas);
138: return(0);
139: }
143: PetscErrorCode SNESSetUp_FAS(SNES snes)
144: {
145: SNES_FAS *fas = (SNES_FAS*) snes->data;
147: PetscInt dm_levels;
148: Vec vec_sol, vec_func, vec_sol_update, vec_rhs; /* preserve these if they're set through the reset */
149: SNES next;
150: PetscBool isFine;
151: SNESLineSearch linesearch;
152: SNESLineSearch slinesearch;
153: void *lsprectx,*lspostctx;
154: PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*);
155: PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
158: SNESFASCycleIsFine(snes, &isFine);
159: if (fas->usedmfornumberoflevels && isFine) {
160: DMGetRefineLevel(snes->dm,&dm_levels);
161: dm_levels++;
162: if (dm_levels > fas->levels) {
163: /* we don't want the solution and func vectors to be destroyed in the SNESReset when it's called in SNESFASSetLevels_FAS*/
164: vec_sol = snes->vec_sol;
165: vec_func = snes->vec_func;
166: vec_sol_update = snes->vec_sol_update;
167: vec_rhs = snes->vec_rhs;
168: snes->vec_sol = NULL;
169: snes->vec_func = NULL;
170: snes->vec_sol_update = NULL;
171: snes->vec_rhs = NULL;
173: /* reset the number of levels */
174: SNESFASSetLevels(snes,dm_levels,NULL);
175: SNESSetFromOptions(snes);
177: snes->vec_sol = vec_sol;
178: snes->vec_func = vec_func;
179: snes->vec_rhs = vec_rhs;
180: snes->vec_sol_update = vec_sol_update;
181: }
182: }
183: SNESFASCycleGetCorrection(snes, &next);
184: if (!isFine) snes->gridsequence = 0; /* no grid sequencing inside the multigrid hierarchy! */
186: SNESSetWorkVecs(snes, 2); /* work vectors used for intergrid transfers */
188: /* set up the smoothers if they haven't already been set up */
189: if (!fas->smoothd) {
190: SNESFASCycleCreateSmoother_Private(snes, &fas->smoothd);
191: }
193: if (snes->dm) {
194: /* set the smoother DMs properly */
195: if (fas->smoothu) SNESSetDM(fas->smoothu, snes->dm);
196: SNESSetDM(fas->smoothd, snes->dm);
197: /* construct EVERYTHING from the DM -- including the progressive set of smoothers */
198: if (next) {
199: /* for now -- assume the DM and the evaluation functions have been set externally */
200: if (!next->dm) {
201: DMCoarsen(snes->dm, PetscObjectComm((PetscObject)next), &next->dm);
202: SNESSetDM(next, next->dm);
203: }
204: /* set the interpolation and restriction from the DM */
205: if (!fas->interpolate) {
206: DMCreateInterpolation(next->dm, snes->dm, &fas->interpolate, &fas->rscale);
207: if (!fas->restrct) {
208: PetscObjectReference((PetscObject)fas->interpolate);
209: fas->restrct = fas->interpolate;
210: }
211: }
212: /* set the injection from the DM */
213: if (!fas->inject) {
214: DMCreateInjection(next->dm, snes->dm, &fas->inject);
215: }
216: }
217: }
218: /*pass the smoother, function, and jacobian up to the next level if it's not user set already */
219: if (fas->galerkin) {
220: if (next) {
221: SNESSetFunction(next, NULL, SNESFASGalerkinDefaultFunction, next);
222: }
223: if (fas->smoothd && fas->level != fas->levels - 1) {
224: SNESSetFunction(fas->smoothd, NULL, SNESFASGalerkinDefaultFunction, snes);
225: }
226: if (fas->smoothu && fas->level != fas->levels - 1) {
227: SNESSetFunction(fas->smoothu, NULL, SNESFASGalerkinDefaultFunction, snes);
228: }
229: }
231: /* sets the down (pre) smoother's default norm and sets it from options */
232: if (fas->smoothd) {
233: if (fas->level == 0 && fas->levels != 1) {
234: SNESSetNormSchedule(fas->smoothd, SNES_NORM_NONE);
235: } else {
236: SNESSetNormSchedule(fas->smoothd, SNES_NORM_FINAL_ONLY);
237: }
238: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothd);
239: SNESSetFromOptions(fas->smoothd);
240: SNESGetLineSearch(snes,&linesearch);
241: SNESGetLineSearch(fas->smoothd,&slinesearch);
242: SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
243: SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
244: SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
245: SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
246: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
248: fas->smoothd->vec_sol = snes->vec_sol;
249: PetscObjectReference((PetscObject)snes->vec_sol);
250: fas->smoothd->vec_sol_update = snes->vec_sol_update;
251: PetscObjectReference((PetscObject)snes->vec_sol_update);
252: fas->smoothd->vec_func = snes->vec_func;
253: PetscObjectReference((PetscObject)snes->vec_func);
255: if (fas->eventsmoothsetup) {PetscLogEventBegin(fas->eventsmoothsetup,0,0,0,0);}
256: SNESSetUp(fas->smoothd);
257: if (fas->eventsmoothsetup) {PetscLogEventEnd(fas->eventsmoothsetup,0,0,0,0);}
258: }
260: /* sets the up (post) smoother's default norm and sets it from options */
261: if (fas->smoothu) {
262: if (fas->level != fas->levels - 1) {
263: SNESSetNormSchedule(fas->smoothu, SNES_NORM_NONE);
264: } else {
265: SNESSetNormSchedule(fas->smoothu, SNES_NORM_FINAL_ONLY);
266: }
267: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothu);
268: SNESSetFromOptions(fas->smoothu);
269: SNESGetLineSearch(snes,&linesearch);
270: SNESGetLineSearch(fas->smoothu,&slinesearch);
271: SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
272: SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
273: SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
274: SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
275: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
277: fas->smoothu->vec_sol = snes->vec_sol;
278: PetscObjectReference((PetscObject)snes->vec_sol);
279: fas->smoothu->vec_sol_update = snes->vec_sol_update;
280: PetscObjectReference((PetscObject)snes->vec_sol_update);
281: fas->smoothu->vec_func = snes->vec_func;
282: PetscObjectReference((PetscObject)snes->vec_func);
284: if (fas->eventsmoothsetup) {PetscLogEventBegin(fas->eventsmoothsetup,0,0,0,0);}
285: SNESSetUp(fas->smoothu);
286: if (fas->eventsmoothsetup) {PetscLogEventEnd(fas->eventsmoothsetup,0,0,0,0);}
288: }
290: if (next) {
291: /* gotta set up the solution vector for this to work */
292: if (!next->vec_sol) {SNESFASCreateCoarseVec(snes,&next->vec_sol);}
293: if (!next->vec_rhs) {SNESFASCreateCoarseVec(snes,&next->vec_rhs);}
294: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next);
295: SNESGetLineSearch(snes,&linesearch);
296: SNESGetLineSearch(fas->next,&slinesearch);
297: SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
298: SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
299: SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
300: SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
301: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
302: SNESSetUp(next);
303: }
304: /* setup FAS work vectors */
305: if (fas->galerkin) {
306: VecDuplicate(snes->vec_sol, &fas->Xg);
307: VecDuplicate(snes->vec_sol, &fas->Fg);
308: }
309: return(0);
310: }
314: PetscErrorCode SNESSetFromOptions_FAS(PetscOptionItems *PetscOptionsObject,SNES snes)
315: {
316: SNES_FAS *fas = (SNES_FAS*) snes->data;
317: PetscInt levels = 1;
318: PetscBool flg = PETSC_FALSE, upflg = PETSC_FALSE, downflg = PETSC_FALSE, monflg = PETSC_FALSE, galerkinflg = PETSC_FALSE,continuationflg = PETSC_FALSE;
320: SNESFASType fastype;
321: const char *optionsprefix;
322: SNESLineSearch linesearch;
323: PetscInt m, n_up, n_down;
324: SNES next;
325: PetscBool isFine;
328: SNESFASCycleIsFine(snes, &isFine);
329: PetscOptionsHead(PetscOptionsObject,"SNESFAS Options-----------------------------------");
331: /* number of levels -- only process most options on the finest level */
332: if (isFine) {
333: PetscOptionsInt("-snes_fas_levels", "Number of Levels", "SNESFASSetLevels", levels, &levels, &flg);
334: if (!flg && snes->dm) {
335: DMGetRefineLevel(snes->dm,&levels);
336: levels++;
337: fas->usedmfornumberoflevels = PETSC_TRUE;
338: }
339: SNESFASSetLevels(snes, levels, NULL);
340: fastype = fas->fastype;
341: PetscOptionsEnum("-snes_fas_type","FAS correction type","SNESFASSetType",SNESFASTypes,(PetscEnum)fastype,(PetscEnum*)&fastype,&flg);
342: if (flg) {
343: SNESFASSetType(snes, fastype);
344: }
346: SNESGetOptionsPrefix(snes, &optionsprefix);
347: PetscOptionsInt("-snes_fas_cycles","Number of cycles","SNESFASSetCycles",fas->n_cycles,&m,&flg);
348: if (flg) {
349: SNESFASSetCycles(snes, m);
350: }
351: PetscOptionsBool("-snes_fas_continuation","Corrected grid-sequence continuation","SNESFASSetContinuation",fas->continuation,&continuationflg,&flg);
352: if (flg) {
353: SNESFASSetContinuation(snes,continuationflg);
354: }
356: PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin","SNESFASSetGalerkin",fas->galerkin,&galerkinflg,&flg);
357: if (flg) {
358: SNESFASSetGalerkin(snes, galerkinflg);
359: }
361: if (fas->fastype == SNES_FAS_FULL) {
362: PetscOptionsBool("-snes_fas_full_downsweep","Smooth on the initial upsweep for full FAS cycles","SNESFASFullSetDownSweep",fas->full_downsweep,&fas->full_downsweep,&flg);
363: if (flg) {SNESFASFullSetDownSweep(snes,fas->full_downsweep);}
364: }
366: PetscOptionsInt("-snes_fas_smoothup","Number of post-smoothing steps","SNESFASSetNumberSmoothUp",fas->max_up_it,&n_up,&upflg);
368: PetscOptionsInt("-snes_fas_smoothdown","Number of pre-smoothing steps","SNESFASSetNumberSmoothDown",fas->max_down_it,&n_down,&downflg);
370: {
371: PetscViewer viewer;
372: PetscViewerFormat format;
373: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,
374: "-snes_fas_monitor",&viewer,&format,&monflg);
375: if (monflg) {
376: PetscViewerAndFormat *vf;
377: PetscViewerAndFormatCreate(viewer,format,&vf);
378: PetscObjectDereference((PetscObject)viewer);
379: SNESFASSetMonitor(snes,vf,PETSC_TRUE);
380: }
381: }
382: flg = PETSC_FALSE;
383: monflg = PETSC_TRUE;
384: PetscOptionsBool("-snes_fas_log","Log times for each FAS level","SNESFASSetLog",monflg,&monflg,&flg);
385: if (flg) {SNESFASSetLog(snes,monflg);}
386: }
388: PetscOptionsTail();
389: /* setup from the determined types if there is no pointwise procedure or smoother defined */
390: if (upflg) {
391: SNESFASSetNumberSmoothUp(snes,n_up);
392: }
393: if (downflg) {
394: SNESFASSetNumberSmoothDown(snes,n_down);
395: }
397: /* set up the default line search for coarse grid corrections */
398: if (fas->fastype == SNES_FAS_ADDITIVE) {
399: if (!snes->linesearch) {
400: SNESGetLineSearch(snes, &linesearch);
401: SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
402: }
403: }
405: SNESFASCycleGetCorrection(snes, &next);
406: /* recursive option setting for the smoothers */
407: if (next) {SNESSetFromOptions(next);}
408: return(0);
409: }
411: #include <petscdraw.h>
414: PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer)
415: {
416: SNES_FAS *fas = (SNES_FAS*) snes->data;
417: PetscBool isFine,iascii,isdraw;
418: PetscInt i;
420: SNES smoothu, smoothd, levelsnes;
423: SNESFASCycleIsFine(snes, &isFine);
424: if (isFine) {
425: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
426: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
427: if (iascii) {
428: PetscViewerASCIIPrintf(viewer, "FAS: type is %s, levels=%D, cycles=%D\n", SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles);
429: if (fas->galerkin) {
430: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid function evaluation\n");
431: } else {
432: PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid function evaluation\n");
433: }
434: for (i=0; i<fas->levels; i++) {
435: SNESFASGetCycleSNES(snes, i, &levelsnes);
436: SNESFASCycleGetSmootherUp(levelsnes, &smoothu);
437: SNESFASCycleGetSmootherDown(levelsnes, &smoothd);
438: if (!i) {
439: PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level %D -------------------------------\n",i);
440: } else {
441: PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
442: }
443: PetscViewerASCIIPushTab(viewer);
444: if (smoothd) {
445: SNESView(smoothd,viewer);
446: } else {
447: PetscViewerASCIIPrintf(viewer,"Not yet available\n");
448: }
449: PetscViewerASCIIPopTab(viewer);
450: if (i && (smoothd == smoothu)) {
451: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
452: } else if (i) {
453: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
454: PetscViewerASCIIPushTab(viewer);
455: if (smoothu) {
456: SNESView(smoothu,viewer);
457: } else {
458: PetscViewerASCIIPrintf(viewer,"Not yet available\n");
459: }
460: PetscViewerASCIIPopTab(viewer);
461: }
462: }
463: } else if (isdraw) {
464: PetscDraw draw;
465: PetscReal x,w,y,bottom,th,wth;
466: SNES_FAS *curfas = fas;
467: PetscViewerDrawGetDraw(viewer,0,&draw);
468: PetscDrawGetCurrentPoint(draw,&x,&y);
469: PetscDrawStringGetSize(draw,&wth,&th);
470: bottom = y - th;
471: while (curfas) {
472: if (!curfas->smoothu) {
473: PetscDrawPushCurrentPoint(draw,x,bottom);
474: if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
475: PetscDrawPopCurrentPoint(draw);
476: } else {
477: w = 0.5*PetscMin(1.0-x,x);
478: PetscDrawPushCurrentPoint(draw,x-w,bottom);
479: if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
480: PetscDrawPopCurrentPoint(draw);
481: PetscDrawPushCurrentPoint(draw,x+w,bottom);
482: if (curfas->smoothu) SNESView(curfas->smoothu,viewer);
483: PetscDrawPopCurrentPoint(draw);
484: }
485: /* this is totally bogus but we have no way of knowing how low the previous one was draw to */
486: bottom -= 5*th;
487: if (curfas->next) curfas = (SNES_FAS*)curfas->next->data;
488: else curfas = NULL;
489: }
490: }
491: }
492: return(0);
493: }
497: /*
498: Defines the action of the downsmoother
499: */
500: PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
501: {
502: PetscErrorCode 0;
503: SNESConvergedReason reason;
504: Vec FPC;
505: SNES smoothd;
506: SNES_FAS *fas = (SNES_FAS*) snes->data;
509: SNESFASCycleGetSmootherDown(snes, &smoothd);
510: SNESSetInitialFunction(smoothd, F);
511: if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,0,0,0,0);}
512: SNESSolve(smoothd, B, X);
513: if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,0,0,0,0);}
514: /* check convergence reason for the smoother */
515: SNESGetConvergedReason(smoothd,&reason);
516: if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
517: snes->reason = SNES_DIVERGED_INNER;
518: return(0);
519: }
520: SNESGetFunction(smoothd, &FPC, NULL, NULL);
521: VecCopy(FPC, F);
522: if (fnorm) {VecNorm(F,NORM_2,fnorm);}
523: return(0);
524: }
529: /*
530: Defines the action of the upsmoother
531: */
532: PetscErrorCode SNESFASUpSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
533: {
534: PetscErrorCode 0;
535: SNESConvergedReason reason;
536: Vec FPC;
537: SNES smoothu;
538: SNES_FAS *fas = (SNES_FAS*) snes->data;
541: SNESFASCycleGetSmootherUp(snes, &smoothu);
542: if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,0,0,0,0);}
543: SNESSolve(smoothu, B, X);
544: if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,0,0,0,0);}
545: /* check convergence reason for the smoother */
546: SNESGetConvergedReason(smoothu,&reason);
547: if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
548: snes->reason = SNES_DIVERGED_INNER;
549: return(0);
550: }
551: SNESGetFunction(smoothu, &FPC, NULL, NULL);
552: VecCopy(FPC, F);
553: if (fnorm) {VecNorm(F,NORM_2,fnorm);}
554: return(0);
555: }
559: /*@
560: SNESFASCreateCoarseVec - create Vec corresponding to a state vector on one level coarser than current level
562: Collective
564: Input Arguments:
565: . snes - SNESFAS
567: Output Arguments:
568: . Xcoarse - vector on level one coarser than snes
570: Level: developer
572: .seealso: SNESFASSetRestriction(), SNESFASRestrict()
573: @*/
574: PetscErrorCode SNESFASCreateCoarseVec(SNES snes,Vec *Xcoarse)
575: {
577: SNES_FAS *fas = (SNES_FAS*)snes->data;
580: if (fas->rscale) {
581: VecDuplicate(fas->rscale,Xcoarse);
582: } else if (fas->inject) {
583: MatCreateVecs(fas->inject,Xcoarse,NULL);
584: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or injection");
585: return(0);
586: }
590: /*@
591: SNESFASRestrict - restrict a Vec to the next coarser level
593: Collective
595: Input Arguments:
596: + fine - SNES from which to restrict
597: - Xfine - vector to restrict
599: Output Arguments:
600: . Xcoarse - result of restriction
602: Level: developer
604: .seealso: SNESFASSetRestriction(), SNESFASSetInjection()
605: @*/
606: PetscErrorCode SNESFASRestrict(SNES fine,Vec Xfine,Vec Xcoarse)
607: {
609: SNES_FAS *fas = (SNES_FAS*)fine->data;
615: if (fas->inject) {
616: MatRestrict(fas->inject,Xfine,Xcoarse);
617: } else {
618: MatRestrict(fas->restrct,Xfine,Xcoarse);
619: VecPointwiseMult(Xcoarse,fas->rscale,Xcoarse);
620: }
621: return(0);
622: }
626: /*
628: Performs the FAS coarse correction as:
630: fine problem: F(x) = b
631: coarse problem: F^c(x^c) = b^c
633: b^c = F^c(Rx) - R(F(x) - b)
635: */
636: PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new)
637: {
638: PetscErrorCode ierr;
639: Vec X_c, Xo_c, F_c, B_c;
640: SNESConvergedReason reason;
641: SNES next;
642: Mat restrct, interpolate;
643: SNES_FAS *fasc;
646: SNESFASCycleGetCorrection(snes, &next);
647: if (next) {
648: fasc = (SNES_FAS*)next->data;
650: SNESFASCycleGetRestriction(snes, &restrct);
651: SNESFASCycleGetInterpolation(snes, &interpolate);
653: X_c = next->vec_sol;
654: Xo_c = next->work[0];
655: F_c = next->vec_func;
656: B_c = next->vec_rhs;
658: if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,0,0,0,0);}
659: SNESFASRestrict(snes,X,Xo_c);
660: /* restrict the defect: R(F(x) - b) */
661: MatRestrict(restrct, F, B_c);
662: if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,0,0,0,0);}
664: if (fasc->eventresidual) {PetscLogEventBegin(fasc->eventresidual,0,0,0,0);}
665: /* F_c = F^c(Rx) - R(F(x) - b) since the second term was sitting in next->vec_rhs */
666: SNESComputeFunction(next, Xo_c, F_c);
667: if (fasc->eventresidual) {PetscLogEventEnd(fasc->eventresidual,0,0,0,0);}
669: /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */
670: VecCopy(B_c, X_c);
671: VecCopy(F_c, B_c);
672: VecCopy(X_c, F_c);
673: /* set initial guess of the coarse problem to the projected fine solution */
674: VecCopy(Xo_c, X_c);
676: /* recurse to the next level */
677: SNESSetInitialFunction(next, F_c);
678: SNESSolve(next, B_c, X_c);
679: SNESGetConvergedReason(next,&reason);
680: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
681: snes->reason = SNES_DIVERGED_INNER;
682: return(0);
683: }
684: /* correct as x <- x + I(x^c - Rx)*/
685: VecAXPY(X_c, -1.0, Xo_c);
687: if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,0,0,0,0);}
688: MatInterpolateAdd(interpolate, X_c, X, X_new);
689: if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,0,0,0,0);}
690: }
691: return(0);
692: }
696: /*
698: The additive cycle looks like:
700: xhat = x
701: xhat = dS(x, b)
702: x = coarsecorrection(xhat, b_d)
703: x = x + nu*(xhat - x);
704: (optional) x = uS(x, b)
706: With the coarse RHS (defect correction) as below.
708: */
709: PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X)
710: {
711: Vec F, B, Xhat;
712: Vec X_c, Xo_c, F_c, B_c;
713: PetscErrorCode ierr;
714: SNESConvergedReason reason;
715: PetscReal xnorm, fnorm, ynorm;
716: SNESLineSearchReason lsresult;
717: SNES next;
718: Mat restrct, interpolate;
719: SNES_FAS *fas = (SNES_FAS*)snes->data,*fasc;
722: SNESFASCycleGetCorrection(snes, &next);
723: F = snes->vec_func;
724: B = snes->vec_rhs;
725: Xhat = snes->work[1];
726: VecCopy(X, Xhat);
727: /* recurse first */
728: if (next) {
729: fasc = (SNES_FAS*)next->data;
730: SNESFASCycleGetRestriction(snes, &restrct);
731: SNESFASCycleGetInterpolation(snes, &interpolate);
732: if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,0,0,0,0);}
733: SNESComputeFunction(snes, Xhat, F);
734: if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,0,0,0,0);}
735: VecNorm(F, NORM_2, &fnorm);
736: X_c = next->vec_sol;
737: Xo_c = next->work[0];
738: F_c = next->vec_func;
739: B_c = next->vec_rhs;
741: SNESFASRestrict(snes,Xhat,Xo_c);
742: /* restrict the defect */
743: MatRestrict(restrct, F, B_c);
745: /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */
746: if (fasc->eventresidual) {PetscLogEventBegin(fasc->eventresidual,0,0,0,0);}
747: SNESComputeFunction(next, Xo_c, F_c);
748: if (fasc->eventresidual) {PetscLogEventEnd(fasc->eventresidual,0,0,0,0);}
749: VecCopy(B_c, X_c);
750: VecCopy(F_c, B_c);
751: VecCopy(X_c, F_c);
752: /* set initial guess of the coarse problem to the projected fine solution */
753: VecCopy(Xo_c, X_c);
755: /* recurse */
756: SNESSetInitialFunction(next, F_c);
757: SNESSolve(next, B_c, X_c);
759: /* smooth on this level */
760: SNESFASDownSmooth_Private(snes, B, X, F, &fnorm);
762: SNESGetConvergedReason(next,&reason);
763: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
764: snes->reason = SNES_DIVERGED_INNER;
765: return(0);
766: }
768: /* correct as x <- x + I(x^c - Rx)*/
769: VecAYPX(X_c, -1.0, Xo_c);
770: MatInterpolate(interpolate, X_c, Xhat);
772: /* additive correction of the coarse direction*/
773: SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat);
774: SNESLineSearchGetReason(snes->linesearch, &lsresult);
775: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm);
776: if (lsresult) {
777: if (++snes->numFailures >= snes->maxFailures) {
778: snes->reason = SNES_DIVERGED_LINE_SEARCH;
779: return(0);
780: }
781: }
782: } else {
783: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
784: }
785: return(0);
786: }
790: /*
792: Defines the FAS cycle as:
794: fine problem: F(x) = b
795: coarse problem: F^c(x) = b^c
797: b^c = F^c(Rx) - R(F(x) - b)
799: correction:
801: x = x + I(x^c - Rx)
803: */
804: PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X)
805: {
808: Vec F,B;
809: SNES next;
812: F = snes->vec_func;
813: B = snes->vec_rhs;
814: /* pre-smooth -- just update using the pre-smoother */
815: SNESFASCycleGetCorrection(snes,&next);
816: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
817: if (next) {
818: SNESFASCoarseCorrection(snes, X, F, X);
819: SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
820: }
821: return(0);
822: }
826: PetscErrorCode SNESFASCycleSetupPhase_Full(SNES snes)
827: {
828: SNES next;
829: SNES_FAS *fas = (SNES_FAS*)snes->data;
830: PetscBool isFine;
834: /* pre-smooth -- just update using the pre-smoother */
835: SNESFASCycleIsFine(snes,&isFine);
836: SNESFASCycleGetCorrection(snes,&next);
837: fas->full_stage = 0;
838: if (next) {SNESFASCycleSetupPhase_Full(next);}
839: return(0);
840: }
844: PetscErrorCode SNESFASCycle_Full(SNES snes, Vec X)
845: {
847: Vec F,B;
848: SNES_FAS *fas = (SNES_FAS*)snes->data;
849: PetscBool isFine;
850: SNES next;
853: F = snes->vec_func;
854: B = snes->vec_rhs;
855: SNESFASCycleIsFine(snes,&isFine);
856: SNESFASCycleGetCorrection(snes,&next);
858: if (isFine) {
859: SNESFASCycleSetupPhase_Full(snes);
860: }
862: if (fas->full_stage == 0) {
863: /* downsweep */
864: if (next) {
865: if (fas->level != 1) next->max_its += 1;
866: if (fas->full_downsweep||isFine) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
867: SNESFASCoarseCorrection(snes,X,F,X);
868: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
869: if (fas->level != 1) next->max_its -= 1;
870: } else {
871: /* The smoother on the coarse level is the coarse solver */
872: SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
873: }
874: fas->full_stage = 1;
875: } else if (fas->full_stage == 1) {
876: if (snes->iter == 0) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
877: if (next) {
878: SNESFASCoarseCorrection(snes,X,F,X);
879: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
880: }
881: }
882: /* final v-cycle */
883: if (isFine) {
884: if (next) {
885: SNESFASCoarseCorrection(snes,X,F,X);
886: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
887: }
888: }
889: return(0);
890: }
894: PetscErrorCode SNESFASCycle_Kaskade(SNES snes, Vec X)
895: {
897: Vec F,B;
898: SNES next;
901: F = snes->vec_func;
902: B = snes->vec_rhs;
903: SNESFASCycleGetCorrection(snes,&next);
904: if (next) {
905: SNESFASCoarseCorrection(snes,X,F,X);
906: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
907: } else {
908: SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
909: }
910: return(0);
911: }
913: PetscBool SNEScite = PETSC_FALSE;
914: const char SNESCitation[] = "@techreport{pbmkbsxt2012,\n"
915: " title = {Composing Scalable Nonlinear Algebraic Solvers},\n"
916: " author = {Peter Brune and Mathew Knepley and Barry Smith and Xuemin Tu},\n"
917: " year = 2013,\n"
918: " type = Preprint,\n"
919: " number = {ANL/MCS-P2010-0112},\n"
920: " institution = {Argonne National Laboratory}\n}\n";
924: PetscErrorCode SNESSolve_FAS(SNES snes)
925: {
927: PetscInt i, maxits;
928: Vec X, F;
929: PetscReal fnorm;
930: SNES_FAS *fas = (SNES_FAS*)snes->data,*ffas;
931: DM dm;
932: PetscBool isFine;
936: 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);
938: PetscCitationsRegister(SNESCitation,&SNEScite);
939: maxits = snes->max_its; /* maximum number of iterations */
940: snes->reason = SNES_CONVERGED_ITERATING;
941: X = snes->vec_sol;
942: F = snes->vec_func;
944: SNESFASCycleIsFine(snes, &isFine);
945: /*norm setup */
946: PetscObjectSAWsTakeAccess((PetscObject)snes);
947: snes->iter = 0;
948: snes->norm = 0.;
949: PetscObjectSAWsGrantAccess((PetscObject)snes);
950: if (!snes->vec_func_init_set) {
951: if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,0,0,0,0);}
952: SNESComputeFunction(snes,X,F);
953: if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,0,0,0,0);}
954: } else snes->vec_func_init_set = PETSC_FALSE;
956: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
957: SNESCheckFunctionNorm(snes,fnorm);
958: PetscObjectSAWsTakeAccess((PetscObject)snes);
959: snes->norm = fnorm;
960: PetscObjectSAWsGrantAccess((PetscObject)snes);
961: SNESLogConvergenceHistory(snes,fnorm,0);
962: SNESMonitor(snes,0,fnorm);
964: /* test convergence */
965: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
966: if (snes->reason) return(0);
969: if (isFine) {
970: /* propagate scale-dependent data up the hierarchy */
971: SNESGetDM(snes,&dm);
972: for (ffas=fas; ffas->next; ffas=(SNES_FAS*)ffas->next->data) {
973: DM dmcoarse;
974: SNESGetDM(ffas->next,&dmcoarse);
975: DMRestrict(dm,ffas->restrct,ffas->rscale,ffas->inject,dmcoarse);
976: dm = dmcoarse;
977: }
978: }
980: for (i = 0; i < maxits; i++) {
981: /* Call general purpose update function */
983: if (snes->ops->update) {
984: (*snes->ops->update)(snes, snes->iter);
985: }
986: if (fas->fastype == SNES_FAS_MULTIPLICATIVE) {
987: SNESFASCycle_Multiplicative(snes, X);
988: } else if (fas->fastype == SNES_FAS_ADDITIVE) {
989: SNESFASCycle_Additive(snes, X);
990: } else if (fas->fastype == SNES_FAS_FULL) {
991: SNESFASCycle_Full(snes, X);
992: } else if (fas->fastype ==SNES_FAS_KASKADE) {
993: SNESFASCycle_Kaskade(snes, X);
994: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported FAS type");
996: /* check for FAS cycle divergence */
997: if (snes->reason != SNES_CONVERGED_ITERATING) return(0);
999: /* Monitor convergence */
1000: PetscObjectSAWsTakeAccess((PetscObject)snes);
1001: snes->iter = i+1;
1002: PetscObjectSAWsGrantAccess((PetscObject)snes);
1003: SNESLogConvergenceHistory(snes,snes->norm,0);
1004: SNESMonitor(snes,snes->iter,snes->norm);
1005: /* Test for convergence */
1006: if (isFine) {
1007: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,snes->norm,&snes->reason,snes->cnvP);
1008: if (snes->reason) break;
1009: }
1010: }
1011: if (i == maxits) {
1012: PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);
1013: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
1014: }
1015: return(0);
1016: }