Actual source code: fas.c
petsc-3.4.5 2014-06-29
1: /* Defines the basic SNES object */
2: #include <../src/snes/impls/fas/fasimpls.h> /*I "petscsnesfas.h" I*/
4: const char *const SNESFASTypes[] = {"MULTIPLICATIVE","ADDITIVE","SNESFASType","SNES_FAS",0};
6: extern PetscErrorCode SNESDestroy_FAS(SNES snes);
7: extern PetscErrorCode SNESSetUp_FAS(SNES snes);
8: extern PetscErrorCode SNESSetFromOptions_FAS(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> - 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: . -fas_levels_snes_ - SNES options for all smoothers
32: . -fas_levels_cycle_snes_ - SNES options for all cycles
33: . -fas_levels_i_snes_ - SNES options for the smoothers on level i
34: . -fas_levels_i_cycle_snes_ - SNES options for the cycle on level i
35: - -fas_coarse_snes_ - SNES options for the coarsest smoother
37: Notes:
38: The organization of the FAS solver is slightly different from the organization of PCMG
39: As each level has smoother SNES instances(down and potentially up) and a cycle SNES instance.
40: The cycle SNES instance may be used for monitoring convergence on a particular level.
42: Level: beginner
44: .seealso: PCMG, SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types)
45: M*/
49: PETSC_EXTERN PetscErrorCode SNESCreate_FAS(SNES snes)
50: {
51: SNES_FAS *fas;
55: snes->ops->destroy = SNESDestroy_FAS;
56: snes->ops->setup = SNESSetUp_FAS;
57: snes->ops->setfromoptions = SNESSetFromOptions_FAS;
58: snes->ops->view = SNESView_FAS;
59: snes->ops->solve = SNESSolve_FAS;
60: snes->ops->reset = SNESReset_FAS;
62: snes->usesksp = PETSC_FALSE;
63: snes->usespc = PETSC_FALSE;
65: if (!snes->tolerancesset) {
66: snes->max_funcs = 30000;
67: snes->max_its = 10000;
68: }
70: PetscNewLog(snes, SNES_FAS, &fas);
72: snes->data = (void*) fas;
73: fas->level = 0;
74: fas->levels = 1;
75: fas->n_cycles = 1;
76: fas->max_up_it = 1;
77: fas->max_down_it = 1;
78: fas->smoothu = NULL;
79: fas->smoothd = NULL;
80: fas->next = NULL;
81: fas->previous = NULL;
82: fas->fine = snes;
83: fas->interpolate = NULL;
84: fas->restrct = NULL;
85: fas->inject = NULL;
86: fas->monitor = NULL;
87: fas->usedmfornumberoflevels = PETSC_FALSE;
88: fas->fastype = SNES_FAS_MULTIPLICATIVE;
90: fas->eventsmoothsetup = 0;
91: fas->eventsmoothsolve = 0;
92: fas->eventresidual = 0;
93: fas->eventinterprestrict = 0;
94: return(0);
95: }
99: PetscErrorCode SNESReset_FAS(SNES snes)
100: {
101: PetscErrorCode 0;
102: SNES_FAS * fas = (SNES_FAS*)snes->data;
105: SNESDestroy(&fas->smoothu);
106: SNESDestroy(&fas->smoothd);
107: MatDestroy(&fas->inject);
108: MatDestroy(&fas->interpolate);
109: MatDestroy(&fas->restrct);
110: VecDestroy(&fas->rscale);
111: if (fas->next) SNESReset(fas->next);
112: return(0);
113: }
117: PetscErrorCode SNESDestroy_FAS(SNES snes)
118: {
119: SNES_FAS * fas = (SNES_FAS*)snes->data;
120: PetscErrorCode 0;
123: /* recursively resets and then destroys */
124: SNESReset(snes);
125: if (fas->next) {
126: SNESDestroy(&fas->next);
127: }
128: PetscFree(fas);
129: return(0);
130: }
134: PetscErrorCode SNESSetUp_FAS(SNES snes)
135: {
136: SNES_FAS *fas = (SNES_FAS*) snes->data;
138: VecScatter injscatter;
139: PetscInt dm_levels;
140: Vec vec_sol, vec_func, vec_sol_update, vec_rhs; /* preserve these if they're set through the reset */
141: SNES next;
142: PetscBool isFine;
143: SNESLineSearch linesearch;
144: SNESLineSearch slinesearch;
145: void *lsprectx,*lspostctx;
146: PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*);
147: PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
150: SNESFASCycleIsFine(snes, &isFine);
151: if (fas->usedmfornumberoflevels && isFine) {
152: DMGetRefineLevel(snes->dm,&dm_levels);
153: dm_levels++;
154: if (dm_levels > fas->levels) {
155: /* we don't want the solution and func vectors to be destroyed in the SNESReset when it's called in SNESFASSetLevels_FAS*/
156: vec_sol = snes->vec_sol;
157: vec_func = snes->vec_func;
158: vec_sol_update = snes->vec_sol_update;
159: vec_rhs = snes->vec_rhs;
160: snes->vec_sol = NULL;
161: snes->vec_func = NULL;
162: snes->vec_sol_update = NULL;
163: snes->vec_rhs = NULL;
165: /* reset the number of levels */
166: SNESFASSetLevels(snes,dm_levels,NULL);
167: SNESSetFromOptions(snes);
169: snes->vec_sol = vec_sol;
170: snes->vec_func = vec_func;
171: snes->vec_rhs = vec_rhs;
172: snes->vec_sol_update = vec_sol_update;
173: }
174: }
175: SNESFASCycleGetCorrection(snes, &next);
176: if (!isFine) snes->gridsequence = 0; /* no grid sequencing inside the multigrid hierarchy! */
178: SNESSetWorkVecs(snes, 2); /* work vectors used for intergrid transfers */
180: /* set up the smoothers if they haven't already been set up */
181: if (!fas->smoothd) {
182: SNESFASCycleCreateSmoother_Private(snes, &fas->smoothd);
183: }
185: if (snes->dm) {
186: /* set the smoother DMs properly */
187: if (fas->smoothu) SNESSetDM(fas->smoothu, snes->dm);
188: SNESSetDM(fas->smoothd, snes->dm);
189: /* construct EVERYTHING from the DM -- including the progressive set of smoothers */
190: if (next) {
191: /* for now -- assume the DM and the evaluation functions have been set externally */
192: if (!next->dm) {
193: DMCoarsen(snes->dm, PetscObjectComm((PetscObject)next), &next->dm);
194: SNESSetDM(next, next->dm);
195: }
196: /* set the interpolation and restriction from the DM */
197: if (!fas->interpolate) {
198: DMCreateInterpolation(next->dm, snes->dm, &fas->interpolate, &fas->rscale);
199: if (!fas->restrct) {
200: PetscObjectReference((PetscObject)fas->interpolate);
201: fas->restrct = fas->interpolate;
202: }
203: }
204: /* set the injection from the DM */
205: if (!fas->inject) {
206: DMCreateInjection(next->dm, snes->dm, &injscatter);
207: MatCreateScatter(PetscObjectComm((PetscObject)snes), injscatter, &fas->inject);
208: VecScatterDestroy(&injscatter);
209: }
210: }
211: }
212: /*pass the smoother, function, and jacobian up to the next level if it's not user set already */
213: if (fas->galerkin) {
214: if (next) {
215: SNESSetFunction(next, NULL, SNESFASGalerkinDefaultFunction, next);
216: }
217: if (fas->smoothd && fas->level != fas->levels - 1) {
218: SNESSetFunction(fas->smoothd, NULL, SNESFASGalerkinDefaultFunction, snes);
219: }
220: if (fas->smoothu && fas->level != fas->levels - 1) {
221: SNESSetFunction(fas->smoothu, NULL, SNESFASGalerkinDefaultFunction, snes);
222: }
223: }
225: /* sets the down (pre) smoother's default norm and sets it from options */
226: if (fas->smoothd) {
227: if (fas->level == 0 && fas->levels != 1) {
228: SNESSetNormType(fas->smoothd, SNES_NORM_NONE);
229: } else {
230: SNESSetNormType(fas->smoothd, SNES_NORM_FINAL_ONLY);
231: }
232: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothd);
233: SNESSetFromOptions(fas->smoothd);
234: SNESGetLineSearch(snes,&linesearch);
235: SNESGetLineSearch(fas->smoothd,&slinesearch);
236: SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
237: SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
238: SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
239: SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
240: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
242: fas->smoothd->vec_sol = snes->vec_sol;
243: PetscObjectReference((PetscObject)snes->vec_sol);
244: fas->smoothd->vec_sol_update = snes->vec_sol_update;
245: PetscObjectReference((PetscObject)snes->vec_sol_update);
246: fas->smoothd->vec_func = snes->vec_func;
247: PetscObjectReference((PetscObject)snes->vec_func);
249: if (fas->eventsmoothsetup) {PetscLogEventBegin(fas->eventsmoothsetup,0,0,0,0);}
250: SNESSetUp(fas->smoothd);
251: if (fas->eventsmoothsetup) {PetscLogEventEnd(fas->eventsmoothsetup,0,0,0,0);}
252: }
254: /* sets the up (post) smoother's default norm and sets it from options */
255: if (fas->smoothu) {
256: if (fas->level != fas->levels - 1) {
257: SNESSetNormType(fas->smoothu, SNES_NORM_NONE);
258: } else {
259: SNESSetNormType(fas->smoothu, SNES_NORM_FINAL_ONLY);
260: }
261: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothu);
262: SNESSetFromOptions(fas->smoothu);
263: SNESGetLineSearch(snes,&linesearch);
264: SNESGetLineSearch(fas->smoothu,&slinesearch);
265: SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
266: SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
267: SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
268: SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
269: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
271: fas->smoothu->vec_sol = snes->vec_sol;
272: PetscObjectReference((PetscObject)snes->vec_sol);
273: fas->smoothu->vec_sol_update = snes->vec_sol_update;
274: PetscObjectReference((PetscObject)snes->vec_sol_update);
275: fas->smoothu->vec_func = snes->vec_func;
276: PetscObjectReference((PetscObject)snes->vec_func);
278: if (fas->eventsmoothsetup) {PetscLogEventBegin(fas->eventsmoothsetup,0,0,0,0);}
279: SNESSetUp(fas->smoothu);
280: if (fas->eventsmoothsetup) {PetscLogEventEnd(fas->eventsmoothsetup,0,0,0,0);}
282: }
284: if (next) {
285: /* gotta set up the solution vector for this to work */
286: if (!next->vec_sol) {SNESFASCreateCoarseVec(snes,&next->vec_sol);}
287: if (!next->vec_rhs) {SNESFASCreateCoarseVec(snes,&next->vec_rhs);}
288: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next);
289: SNESGetLineSearch(snes,&linesearch);
290: SNESGetLineSearch(fas->next,&slinesearch);
291: SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
292: SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
293: SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
294: SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
295: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
296: SNESSetUp(next);
297: }
298: /* setup FAS work vectors */
299: if (fas->galerkin) {
300: VecDuplicate(snes->vec_sol, &fas->Xg);
301: VecDuplicate(snes->vec_sol, &fas->Fg);
302: }
303: return(0);
304: }
308: PetscErrorCode SNESSetFromOptions_FAS(SNES snes)
309: {
310: SNES_FAS *fas = (SNES_FAS*) snes->data;
311: PetscInt levels = 1;
312: PetscBool flg = PETSC_FALSE, upflg = PETSC_FALSE, downflg = PETSC_FALSE, monflg = PETSC_FALSE, galerkinflg = PETSC_FALSE;
314: char monfilename[PETSC_MAX_PATH_LEN];
315: SNESFASType fastype;
316: const char *optionsprefix;
317: SNESLineSearch linesearch;
318: PetscInt m, n_up, n_down;
319: SNES next;
320: PetscBool isFine;
323: SNESFASCycleIsFine(snes, &isFine);
324: PetscOptionsHead("SNESFAS Options-----------------------------------");
326: /* number of levels -- only process most options on the finest level */
327: if (isFine) {
328: PetscOptionsInt("-snes_fas_levels", "Number of Levels", "SNESFASSetLevels", levels, &levels, &flg);
329: if (!flg && snes->dm) {
330: DMGetRefineLevel(snes->dm,&levels);
331: levels++;
332: fas->usedmfornumberoflevels = PETSC_TRUE;
333: }
334: SNESFASSetLevels(snes, levels, NULL);
335: fastype = fas->fastype;
336: PetscOptionsEnum("-snes_fas_type","FAS correction type","SNESFASSetType",SNESFASTypes,(PetscEnum)fastype,(PetscEnum*)&fastype,&flg);
337: if (flg) {
338: SNESFASSetType(snes, fastype);
339: }
341: SNESGetOptionsPrefix(snes, &optionsprefix);
342: PetscOptionsInt("-snes_fas_cycles","Number of cycles","SNESFASSetCycles",fas->n_cycles,&m,&flg);
343: if (flg) {
344: SNESFASSetCycles(snes, m);
345: }
347: PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin","SNESFASSetGalerkin",fas->galerkin,&galerkinflg,&flg);
348: if (flg) {
349: SNESFASSetGalerkin(snes, galerkinflg);
350: }
352: PetscOptionsInt("-snes_fas_smoothup","Number of post-smoothing steps","SNESFASSetNumberSmoothUp",fas->max_up_it,&n_up,&upflg);
354: PetscOptionsInt("-snes_fas_smoothdown","Number of pre-smoothing steps","SNESFASSetNumberSmoothDown",fas->max_down_it,&n_down,&downflg);
356: PetscOptionsString("-snes_fas_monitor","Monitor FAS progress","SNESFASSetMonitor","stdout",monfilename,PETSC_MAX_PATH_LEN,&monflg);
357: if (monflg) SNESFASSetMonitor(snes, PETSC_TRUE);
359: flg = PETSC_FALSE;
360: monflg = PETSC_TRUE;
361: PetscOptionsBool("-snes_fas_log","Log times for each FAS level","SNESFASSetLog",monflg,&monflg,&flg);
362: if (flg) {SNESFASSetLog(snes,monflg);}
363: }
365: PetscOptionsTail();
366: /* setup from the determined types if there is no pointwise procedure or smoother defined */
367: if (upflg) {
368: SNESFASSetNumberSmoothUp(snes,n_up);
369: }
370: if (downflg) {
371: SNESFASSetNumberSmoothDown(snes,n_down);
372: }
374: /* set up the default line search for coarse grid corrections */
375: if (fas->fastype == SNES_FAS_ADDITIVE) {
376: if (!snes->linesearch) {
377: SNESGetLineSearch(snes, &linesearch);
378: SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
379: }
380: }
382: SNESFASCycleGetCorrection(snes, &next);
383: /* recursive option setting for the smoothers */
384: if (next) {SNESSetFromOptions(next);}
385: return(0);
386: }
388: #include <petscdraw.h>
391: PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer)
392: {
393: SNES_FAS *fas = (SNES_FAS*) snes->data;
394: PetscBool isFine,iascii,isdraw;
395: PetscInt i;
397: SNES smoothu, smoothd, levelsnes;
400: SNESFASCycleIsFine(snes, &isFine);
401: if (isFine) {
402: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
403: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
404: if (iascii) {
405: PetscViewerASCIIPrintf(viewer, "FAS: type is %s, levels=%D, cycles=%D\n", SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles);
406: if (fas->galerkin) {
407: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid function evaluation\n");
408: } else {
409: PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid function evaluation\n");
410: }
411: for (i=0; i<fas->levels; i++) {
412: SNESFASGetCycleSNES(snes, i, &levelsnes);
413: SNESFASCycleGetSmootherUp(levelsnes, &smoothu);
414: SNESFASCycleGetSmootherDown(levelsnes, &smoothd);
415: if (!i) {
416: PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level %D -------------------------------\n",i);
417: } else {
418: PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
419: }
420: PetscViewerASCIIPushTab(viewer);
421: if (smoothd) {
422: SNESView(smoothd,viewer);
423: } else {
424: PetscViewerASCIIPrintf(viewer,"Not yet available\n");
425: }
426: PetscViewerASCIIPopTab(viewer);
427: if (i && (smoothd == smoothu)) {
428: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
429: } else if (i) {
430: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
431: PetscViewerASCIIPushTab(viewer);
432: if (smoothu) {
433: SNESView(smoothu,viewer);
434: } else {
435: PetscViewerASCIIPrintf(viewer,"Not yet available\n");
436: }
437: PetscViewerASCIIPopTab(viewer);
438: }
439: }
440: } else if (isdraw) {
441: PetscDraw draw;
442: PetscReal x,w,y,bottom,th,wth;
443: SNES_FAS *curfas = fas;
444: PetscViewerDrawGetDraw(viewer,0,&draw);
445: PetscDrawGetCurrentPoint(draw,&x,&y);
446: PetscDrawStringGetSize(draw,&wth,&th);
447: bottom = y - th;
448: while (curfas) {
449: if (!curfas->smoothu) {
450: PetscDrawPushCurrentPoint(draw,x,bottom);
451: if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
452: PetscDrawPopCurrentPoint(draw);
453: } else {
454: w = 0.5*PetscMin(1.0-x,x);
455: PetscDrawPushCurrentPoint(draw,x-w,bottom);
456: if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
457: PetscDrawPopCurrentPoint(draw);
458: PetscDrawPushCurrentPoint(draw,x+w,bottom);
459: if (curfas->smoothu) SNESView(curfas->smoothu,viewer);
460: PetscDrawPopCurrentPoint(draw);
461: }
462: /* this is totally bogus but we have no way of knowing how low the previous one was draw to */
463: bottom -= 5*th;
464: if (curfas->next) curfas = (SNES_FAS*)curfas->next->data;
465: else curfas = NULL;
466: }
467: }
468: }
469: return(0);
470: }
474: /*
475: Defines the action of the downsmoother
476: */
477: PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
478: {
479: PetscErrorCode 0;
480: SNESConvergedReason reason;
481: Vec FPC;
482: SNES smoothd;
483: SNES_FAS *fas = (SNES_FAS*) snes->data;
486: SNESFASCycleGetSmootherDown(snes, &smoothd);
487: SNESSetInitialFunction(smoothd, F);
488: SNESSetInitialFunctionNorm(smoothd, *fnorm);
489: if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,0,0,0,0);}
490: SNESSolve(smoothd, B, X);
491: if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,0,0,0,0);}
492: /* check convergence reason for the smoother */
493: SNESGetConvergedReason(smoothd,&reason);
494: if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN)) {
495: snes->reason = SNES_DIVERGED_INNER;
496: return(0);
497: }
498: SNESGetFunction(smoothd, &FPC, NULL, NULL);
499: VecCopy(FPC, F);
500: SNESGetFunctionNorm(smoothd, fnorm);
501: return(0);
502: }
507: /*
508: Defines the action of the upsmoother
509: */
510: PetscErrorCode SNESFASUpSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
511: {
512: PetscErrorCode 0;
513: SNESConvergedReason reason;
514: Vec FPC;
515: SNES smoothu;
516: SNES_FAS *fas = (SNES_FAS*) snes->data;
519: SNESFASCycleGetSmootherUp(snes, &smoothu);
520: if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,0,0,0,0);}
521: SNESSolve(smoothu, B, X);
522: if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,0,0,0,0);}
523: /* check convergence reason for the smoother */
524: SNESGetConvergedReason(smoothu,&reason);
525: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
526: snes->reason = SNES_DIVERGED_INNER;
527: return(0);
528: }
529: SNESGetFunction(smoothu, &FPC, NULL, NULL);
530: VecCopy(FPC, F);
531: SNESGetFunctionNorm(smoothu, fnorm);
532: return(0);
533: }
537: /*@
538: SNESFASCreateCoarseVec - create Vec corresponding to a state vector on one level coarser than current level
540: Collective
542: Input Arguments:
543: . snes - SNESFAS
545: Output Arguments:
546: . Xcoarse - vector on level one coarser than snes
548: Level: developer
550: .seealso: SNESFASSetRestriction(), SNESFASRestrict()
551: @*/
552: PetscErrorCode SNESFASCreateCoarseVec(SNES snes,Vec *Xcoarse)
553: {
555: SNES_FAS *fas = (SNES_FAS*)snes->data;
558: if (fas->rscale) {
559: VecDuplicate(fas->rscale,Xcoarse);
560: } else if (fas->inject) {
561: MatGetVecs(fas->inject,Xcoarse,NULL);
562: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or injection");
563: return(0);
564: }
568: /*@
569: SNESFASRestrict - restrict a Vec to the next coarser level
571: Collective
573: Input Arguments:
574: + fine - SNES from which to restrict
575: - Xfine - vector to restrict
577: Output Arguments:
578: . Xcoarse - result of restriction
580: Level: developer
582: .seealso: SNESFASSetRestriction(), SNESFASSetInjection()
583: @*/
584: PetscErrorCode SNESFASRestrict(SNES fine,Vec Xfine,Vec Xcoarse)
585: {
587: SNES_FAS *fas = (SNES_FAS*)fine->data;
593: if (fas->inject) {
594: MatRestrict(fas->inject,Xfine,Xcoarse);
595: } else {
596: MatRestrict(fas->restrct,Xfine,Xcoarse);
597: VecPointwiseMult(Xcoarse,fas->rscale,Xcoarse);
598: }
599: return(0);
600: }
604: /*
606: Performs the FAS coarse correction as:
608: fine problem: F(x) = 0
609: coarse problem: F^c(x) = b^c
611: b^c = F^c(I^c_fx^f - I^c_fF(x))
613: */
614: PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new)
615: {
616: PetscErrorCode ierr;
617: Vec X_c, Xo_c, F_c, B_c;
618: SNESConvergedReason reason;
619: SNES next;
620: Mat restrct, interpolate;
621: SNES_FAS *fasc;
624: SNESFASCycleGetCorrection(snes, &next);
625: if (next) {
626: fasc = (SNES_FAS*)next->data;
628: SNESFASCycleGetRestriction(snes, &restrct);
629: SNESFASCycleGetInterpolation(snes, &interpolate);
631: X_c = next->vec_sol;
632: Xo_c = next->work[0];
633: F_c = next->vec_func;
634: B_c = next->vec_rhs;
636: if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,0,0,0,0);}
637: SNESFASRestrict(snes,X,Xo_c);
638: /* restrict the defect */
639: MatRestrict(restrct, F, B_c);
640: if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,0,0,0,0);}
642: if (fasc->eventresidual) {PetscLogEventBegin(fasc->eventresidual,0,0,0,0);}
643: SNESComputeFunction(next, Xo_c, F_c);
644: if (fasc->eventresidual) {PetscLogEventEnd(fasc->eventresidual,0,0,0,0);}
646: /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */
647: VecCopy(B_c, X_c);
648: VecCopy(F_c, B_c);
649: VecCopy(X_c, F_c);
650: /* set initial guess of the coarse problem to the projected fine solution */
651: VecCopy(Xo_c, X_c);
653: /* recurse to the next level */
654: SNESSetInitialFunction(next, F_c);
655: SNESSolve(next, B_c, X_c);
656: SNESGetConvergedReason(next,&reason);
657: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
658: snes->reason = SNES_DIVERGED_INNER;
659: return(0);
660: }
661: /* correct as x <- x + I(x^c - Rx)*/
662: VecAXPY(X_c, -1.0, Xo_c);
664: if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,0,0,0,0);}
665: MatInterpolateAdd(interpolate, X_c, X, X_new);
666: if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,0,0,0,0);}
667: }
668: return(0);
669: }
673: /*
675: The additive cycle looks like:
677: xhat = x
678: xhat = dS(x, b)
679: x = coarsecorrection(xhat, b_d)
680: x = x + nu*(xhat - x);
681: (optional) x = uS(x, b)
683: With the coarse RHS (defect correction) as below.
685: */
686: PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X)
687: {
688: Vec F, B, Xhat;
689: Vec X_c, Xo_c, F_c, B_c;
690: PetscErrorCode ierr;
691: SNESConvergedReason reason;
692: PetscReal xnorm, fnorm, ynorm;
693: PetscBool lssuccess;
694: SNES next;
695: Mat restrct, interpolate;
696: SNES_FAS *fas = (SNES_FAS*)snes->data,*fasc;
699: SNESFASCycleGetCorrection(snes, &next);
700: F = snes->vec_func;
701: B = snes->vec_rhs;
702: Xhat = snes->work[1];
703: VecCopy(X, Xhat);
704: /* recurse first */
705: if (next) {
706: fasc = (SNES_FAS*)next->data;
707: SNESFASCycleGetRestriction(snes, &restrct);
708: SNESFASCycleGetInterpolation(snes, &interpolate);
709: if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,0,0,0,0);}
710: SNESComputeFunction(snes, Xhat, F);
711: if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,0,0,0,0);}
712: VecNorm(F, NORM_2, &fnorm);
713: X_c = next->vec_sol;
714: Xo_c = next->work[0];
715: F_c = next->vec_func;
716: B_c = next->vec_rhs;
718: SNESFASRestrict(snes,Xhat,Xo_c);
719: /* restrict the defect */
720: MatRestrict(restrct, F, B_c);
722: /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */
723: if (fasc->eventresidual) {PetscLogEventBegin(fasc->eventresidual,0,0,0,0);}
724: SNESComputeFunction(next, Xo_c, F_c);
725: if (fasc->eventresidual) {PetscLogEventEnd(fasc->eventresidual,0,0,0,0);}
726: VecCopy(B_c, X_c);
727: VecCopy(F_c, B_c);
728: VecCopy(X_c, F_c);
729: /* set initial guess of the coarse problem to the projected fine solution */
730: VecCopy(Xo_c, X_c);
732: /* recurse */
733: SNESSetInitialFunction(next, F_c);
734: SNESSolve(next, B_c, X_c);
736: /* smooth on this level */
737: SNESFASDownSmooth_Private(snes, B, X, F, &fnorm);
739: SNESGetConvergedReason(next,&reason);
740: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
741: snes->reason = SNES_DIVERGED_INNER;
742: return(0);
743: }
745: /* correct as x <- x + I(x^c - Rx)*/
746: VecAYPX(X_c, -1.0, Xo_c);
747: MatInterpolate(interpolate, X_c, Xhat);
749: /* additive correction of the coarse direction*/
750: SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat);
751: SNESLineSearchGetSuccess(snes->linesearch, &lssuccess);
752: if (!lssuccess) {
753: if (++snes->numFailures >= snes->maxFailures) {
754: snes->reason = SNES_DIVERGED_LINE_SEARCH;
755: return(0);
756: }
757: }
758: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm);
759: } else {
760: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
761: }
762: return(0);
763: }
767: /*
769: Defines the FAS cycle as:
771: fine problem: F(x) = 0
772: coarse problem: F^c(x) = b^c
774: b^c = F^c(I^c_fx^f - I^c_fF(x))
776: correction:
778: x = x + I(x^c - Rx)
780: */
781: PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X)
782: {
785: Vec F,B;
786: SNES_FAS *fas = (SNES_FAS*)snes->data;
789: F = snes->vec_func;
790: B = snes->vec_rhs;
791: /* pre-smooth -- just update using the pre-smoother */
792: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
793: if (fas->level != 0) {
794: SNESFASCoarseCorrection(snes, X, F, X);
795: SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
796: }
797: return(0);
798: }
803: PetscErrorCode SNESSolve_FAS(SNES snes)
804: {
806: PetscInt i, maxits;
807: Vec X, F;
808: PetscReal fnorm;
809: SNES_FAS *fas = (SNES_FAS*)snes->data,*ffas;
810: DM dm;
811: PetscBool isFine;
814: maxits = snes->max_its; /* maximum number of iterations */
815: snes->reason = SNES_CONVERGED_ITERATING;
816: X = snes->vec_sol;
817: F = snes->vec_func;
819: SNESFASCycleIsFine(snes, &isFine);
820: /*norm setup */
821: PetscObjectAMSTakeAccess((PetscObject)snes);
822: snes->iter = 0;
823: snes->norm = 0.;
824: PetscObjectAMSGrantAccess((PetscObject)snes);
825: if (!snes->vec_func_init_set) {
826: if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,0,0,0,0);}
827: SNESComputeFunction(snes,X,F);
828: if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,0,0,0,0);}
829: if (snes->domainerror) {
830: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
831: return(0);
832: }
833: } else snes->vec_func_init_set = PETSC_FALSE;
835: if (!snes->norm_init_set) {
836: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
837: if (PetscIsInfOrNanReal(fnorm)) {
838: snes->reason = SNES_DIVERGED_FNORM_NAN;
839: return(0);
840: }
841: } else {
842: fnorm = snes->norm_init;
843: snes->norm_init_set = PETSC_FALSE;
844: }
846: PetscObjectAMSTakeAccess((PetscObject)snes);
847: snes->norm = fnorm;
848: PetscObjectAMSGrantAccess((PetscObject)snes);
849: SNESLogConvergenceHistory(snes,fnorm,0);
850: SNESMonitor(snes,0,fnorm);
852: /* set parameter for default relative tolerance convergence test */
853: snes->ttol = fnorm*snes->rtol;
854: /* test convergence */
855: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
856: if (snes->reason) return(0);
859: if (isFine) {
860: /* propagate scale-dependent data up the hierarchy */
861: SNESGetDM(snes,&dm);
862: for (ffas=fas; ffas->next; ffas=(SNES_FAS*)ffas->next->data) {
863: DM dmcoarse;
864: SNESGetDM(ffas->next,&dmcoarse);
865: DMRestrict(dm,ffas->restrct,ffas->rscale,ffas->inject,dmcoarse);
866: dm = dmcoarse;
867: }
868: }
870: for (i = 0; i < maxits; i++) {
871: /* Call general purpose update function */
873: if (snes->ops->update) {
874: (*snes->ops->update)(snes, snes->iter);
875: }
876: if (fas->fastype == SNES_FAS_MULTIPLICATIVE) {
877: SNESFASCycle_Multiplicative(snes, X);
878: } else {
879: SNESFASCycle_Additive(snes, X);
880: }
882: /* check for FAS cycle divergence */
883: if (snes->reason != SNES_CONVERGED_ITERATING) return(0);
885: /* Monitor convergence */
886: PetscObjectAMSTakeAccess((PetscObject)snes);
887: snes->iter = i+1;
888: PetscObjectAMSGrantAccess((PetscObject)snes);
889: SNESLogConvergenceHistory(snes,snes->norm,0);
890: SNESMonitor(snes,snes->iter,snes->norm);
891: /* Test for convergence */
892: if (isFine) {
893: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,snes->norm,&snes->reason,snes->cnvP);
894: if (snes->reason) break;
895: }
896: }
897: if (i == maxits) {
898: PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);
899: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
900: }
901: return(0);
902: }