Actual source code: fas.c
petsc-3.5.4 2015-05-23
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(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_downsweepsmooth<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: .seealso: PCMG, SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types)
46: M*/
50: PETSC_EXTERN PetscErrorCode SNESCreate_FAS(SNES snes)
51: {
52: SNES_FAS *fas;
56: snes->ops->destroy = SNESDestroy_FAS;
57: snes->ops->setup = SNESSetUp_FAS;
58: snes->ops->setfromoptions = SNESSetFromOptions_FAS;
59: snes->ops->view = SNESView_FAS;
60: snes->ops->solve = SNESSolve_FAS;
61: snes->ops->reset = SNESReset_FAS;
63: snes->usesksp = PETSC_FALSE;
64: snes->usespc = PETSC_FALSE;
66: if (!snes->tolerancesset) {
67: snes->max_funcs = 30000;
68: snes->max_its = 10000;
69: }
71: PetscNewLog(snes,&fas);
73: snes->data = (void*) fas;
74: fas->level = 0;
75: fas->levels = 1;
76: fas->n_cycles = 1;
77: fas->max_up_it = 1;
78: fas->max_down_it = 1;
79: fas->smoothu = NULL;
80: fas->smoothd = NULL;
81: fas->next = NULL;
82: fas->previous = NULL;
83: fas->fine = snes;
84: fas->interpolate = NULL;
85: fas->restrct = NULL;
86: fas->inject = NULL;
87: fas->monitor = NULL;
88: fas->usedmfornumberoflevels = PETSC_FALSE;
89: fas->fastype = SNES_FAS_MULTIPLICATIVE;
90: fas->full_downsweep = PETSC_FALSE;
92: fas->eventsmoothsetup = 0;
93: fas->eventsmoothsolve = 0;
94: fas->eventresidual = 0;
95: fas->eventinterprestrict = 0;
96: return(0);
97: }
101: PetscErrorCode SNESReset_FAS(SNES snes)
102: {
103: PetscErrorCode 0;
104: SNES_FAS * fas = (SNES_FAS*)snes->data;
107: SNESDestroy(&fas->smoothu);
108: SNESDestroy(&fas->smoothd);
109: MatDestroy(&fas->inject);
110: MatDestroy(&fas->interpolate);
111: MatDestroy(&fas->restrct);
112: VecDestroy(&fas->rscale);
113: if (fas->galerkin) {
114: VecDestroy(&fas->Xg);
115: VecDestroy(&fas->Fg);
116: }
117: if (fas->next) {SNESReset(fas->next);}
118: return(0);
119: }
123: PetscErrorCode SNESDestroy_FAS(SNES snes)
124: {
125: SNES_FAS * fas = (SNES_FAS*)snes->data;
126: PetscErrorCode 0;
129: /* recursively resets and then destroys */
130: SNESReset(snes);
131: if (fas->next) {
132: SNESDestroy(&fas->next);
133: }
134: PetscFree(fas);
135: return(0);
136: }
140: PetscErrorCode SNESSetUp_FAS(SNES snes)
141: {
142: SNES_FAS *fas = (SNES_FAS*) snes->data;
144: VecScatter injscatter;
145: PetscInt dm_levels;
146: Vec vec_sol, vec_func, vec_sol_update, vec_rhs; /* preserve these if they're set through the reset */
147: SNES next;
148: PetscBool isFine;
149: SNESLineSearch linesearch;
150: SNESLineSearch slinesearch;
151: void *lsprectx,*lspostctx;
152: PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*);
153: PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
156: SNESFASCycleIsFine(snes, &isFine);
157: if (fas->usedmfornumberoflevels && isFine) {
158: DMGetRefineLevel(snes->dm,&dm_levels);
159: dm_levels++;
160: if (dm_levels > fas->levels) {
161: /* we don't want the solution and func vectors to be destroyed in the SNESReset when it's called in SNESFASSetLevels_FAS*/
162: vec_sol = snes->vec_sol;
163: vec_func = snes->vec_func;
164: vec_sol_update = snes->vec_sol_update;
165: vec_rhs = snes->vec_rhs;
166: snes->vec_sol = NULL;
167: snes->vec_func = NULL;
168: snes->vec_sol_update = NULL;
169: snes->vec_rhs = NULL;
171: /* reset the number of levels */
172: SNESFASSetLevels(snes,dm_levels,NULL);
173: SNESSetFromOptions(snes);
175: snes->vec_sol = vec_sol;
176: snes->vec_func = vec_func;
177: snes->vec_rhs = vec_rhs;
178: snes->vec_sol_update = vec_sol_update;
179: }
180: }
181: SNESFASCycleGetCorrection(snes, &next);
182: if (!isFine) snes->gridsequence = 0; /* no grid sequencing inside the multigrid hierarchy! */
184: SNESSetWorkVecs(snes, 2); /* work vectors used for intergrid transfers */
186: /* set up the smoothers if they haven't already been set up */
187: if (!fas->smoothd) {
188: SNESFASCycleCreateSmoother_Private(snes, &fas->smoothd);
189: }
191: if (snes->dm) {
192: /* set the smoother DMs properly */
193: if (fas->smoothu) SNESSetDM(fas->smoothu, snes->dm);
194: SNESSetDM(fas->smoothd, snes->dm);
195: /* construct EVERYTHING from the DM -- including the progressive set of smoothers */
196: if (next) {
197: /* for now -- assume the DM and the evaluation functions have been set externally */
198: if (!next->dm) {
199: DMCoarsen(snes->dm, PetscObjectComm((PetscObject)next), &next->dm);
200: SNESSetDM(next, next->dm);
201: }
202: /* set the interpolation and restriction from the DM */
203: if (!fas->interpolate) {
204: DMCreateInterpolation(next->dm, snes->dm, &fas->interpolate, &fas->rscale);
205: if (!fas->restrct) {
206: PetscObjectReference((PetscObject)fas->interpolate);
207: fas->restrct = fas->interpolate;
208: }
209: }
210: /* set the injection from the DM */
211: if (!fas->inject) {
212: DMCreateInjection(next->dm, snes->dm, &injscatter);
213: MatCreateScatter(PetscObjectComm((PetscObject)snes), injscatter, &fas->inject);
214: VecScatterDestroy(&injscatter);
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(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: char monfilename[PETSC_MAX_PATH_LEN];
321: SNESFASType fastype;
322: const char *optionsprefix;
323: SNESLineSearch linesearch;
324: PetscInt m, n_up, n_down;
325: SNES next;
326: PetscBool isFine;
329: SNESFASCycleIsFine(snes, &isFine);
330: PetscOptionsHead("SNESFAS Options-----------------------------------");
332: /* number of levels -- only process most options on the finest level */
333: if (isFine) {
334: PetscOptionsInt("-snes_fas_levels", "Number of Levels", "SNESFASSetLevels", levels, &levels, &flg);
335: if (!flg && snes->dm) {
336: DMGetRefineLevel(snes->dm,&levels);
337: levels++;
338: fas->usedmfornumberoflevels = PETSC_TRUE;
339: }
340: SNESFASSetLevels(snes, levels, NULL);
341: fastype = fas->fastype;
342: PetscOptionsEnum("-snes_fas_type","FAS correction type","SNESFASSetType",SNESFASTypes,(PetscEnum)fastype,(PetscEnum*)&fastype,&flg);
343: if (flg) {
344: SNESFASSetType(snes, fastype);
345: }
347: SNESGetOptionsPrefix(snes, &optionsprefix);
348: PetscOptionsInt("-snes_fas_cycles","Number of cycles","SNESFASSetCycles",fas->n_cycles,&m,&flg);
349: if (flg) {
350: SNESFASSetCycles(snes, m);
351: }
352: PetscOptionsBool("-snes_fas_continuation","Corrected grid-sequence continuation","SNESFASSetContinuation",fas->continuation,&continuationflg,&flg);
353: if (flg) {
354: SNESFASSetContinuation(snes,continuationflg);
355: }
357: PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin","SNESFASSetGalerkin",fas->galerkin,&galerkinflg,&flg);
358: if (flg) {
359: SNESFASSetGalerkin(snes, galerkinflg);
360: }
362: if (fas->fastype == SNES_FAS_FULL) {
363: PetscOptionsBool("-snes_fas_full_downsweep","Smooth on the initial upsweep for full FAS cycles","SNESFASFullSetDownSweep",fas->full_downsweep,&fas->full_downsweep,&flg);
364: if (flg) {SNESFASFullSetDownSweep(snes,fas->full_downsweep);}
365: }
367: PetscOptionsInt("-snes_fas_smoothup","Number of post-smoothing steps","SNESFASSetNumberSmoothUp",fas->max_up_it,&n_up,&upflg);
369: PetscOptionsInt("-snes_fas_smoothdown","Number of pre-smoothing steps","SNESFASSetNumberSmoothDown",fas->max_down_it,&n_down,&downflg);
371: PetscOptionsString("-snes_fas_monitor","Monitor FAS progress","SNESFASSetMonitor","stdout",monfilename,PETSC_MAX_PATH_LEN,&monflg);
372: if (monflg) SNESFASSetMonitor(snes, PETSC_TRUE);
374: flg = PETSC_FALSE;
375: monflg = PETSC_TRUE;
376: PetscOptionsBool("-snes_fas_log","Log times for each FAS level","SNESFASSetLog",monflg,&monflg,&flg);
377: if (flg) {SNESFASSetLog(snes,monflg);}
378: }
380: PetscOptionsTail();
381: /* setup from the determined types if there is no pointwise procedure or smoother defined */
382: if (upflg) {
383: SNESFASSetNumberSmoothUp(snes,n_up);
384: }
385: if (downflg) {
386: SNESFASSetNumberSmoothDown(snes,n_down);
387: }
389: /* set up the default line search for coarse grid corrections */
390: if (fas->fastype == SNES_FAS_ADDITIVE) {
391: if (!snes->linesearch) {
392: SNESGetLineSearch(snes, &linesearch);
393: SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
394: }
395: }
397: SNESFASCycleGetCorrection(snes, &next);
398: /* recursive option setting for the smoothers */
399: if (next) {SNESSetFromOptions(next);}
400: return(0);
401: }
403: #include <petscdraw.h>
406: PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer)
407: {
408: SNES_FAS *fas = (SNES_FAS*) snes->data;
409: PetscBool isFine,iascii,isdraw;
410: PetscInt i;
412: SNES smoothu, smoothd, levelsnes;
415: SNESFASCycleIsFine(snes, &isFine);
416: if (isFine) {
417: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
418: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
419: if (iascii) {
420: PetscViewerASCIIPrintf(viewer, "FAS: type is %s, levels=%D, cycles=%D\n", SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles);
421: if (fas->galerkin) {
422: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid function evaluation\n");
423: } else {
424: PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid function evaluation\n");
425: }
426: for (i=0; i<fas->levels; i++) {
427: SNESFASGetCycleSNES(snes, i, &levelsnes);
428: SNESFASCycleGetSmootherUp(levelsnes, &smoothu);
429: SNESFASCycleGetSmootherDown(levelsnes, &smoothd);
430: if (!i) {
431: PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level %D -------------------------------\n",i);
432: } else {
433: PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
434: }
435: PetscViewerASCIIPushTab(viewer);
436: if (smoothd) {
437: SNESView(smoothd,viewer);
438: } else {
439: PetscViewerASCIIPrintf(viewer,"Not yet available\n");
440: }
441: PetscViewerASCIIPopTab(viewer);
442: if (i && (smoothd == smoothu)) {
443: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
444: } else if (i) {
445: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
446: PetscViewerASCIIPushTab(viewer);
447: if (smoothu) {
448: SNESView(smoothu,viewer);
449: } else {
450: PetscViewerASCIIPrintf(viewer,"Not yet available\n");
451: }
452: PetscViewerASCIIPopTab(viewer);
453: }
454: }
455: } else if (isdraw) {
456: PetscDraw draw;
457: PetscReal x,w,y,bottom,th,wth;
458: SNES_FAS *curfas = fas;
459: PetscViewerDrawGetDraw(viewer,0,&draw);
460: PetscDrawGetCurrentPoint(draw,&x,&y);
461: PetscDrawStringGetSize(draw,&wth,&th);
462: bottom = y - th;
463: while (curfas) {
464: if (!curfas->smoothu) {
465: PetscDrawPushCurrentPoint(draw,x,bottom);
466: if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
467: PetscDrawPopCurrentPoint(draw);
468: } else {
469: w = 0.5*PetscMin(1.0-x,x);
470: PetscDrawPushCurrentPoint(draw,x-w,bottom);
471: if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
472: PetscDrawPopCurrentPoint(draw);
473: PetscDrawPushCurrentPoint(draw,x+w,bottom);
474: if (curfas->smoothu) SNESView(curfas->smoothu,viewer);
475: PetscDrawPopCurrentPoint(draw);
476: }
477: /* this is totally bogus but we have no way of knowing how low the previous one was draw to */
478: bottom -= 5*th;
479: if (curfas->next) curfas = (SNES_FAS*)curfas->next->data;
480: else curfas = NULL;
481: }
482: }
483: }
484: return(0);
485: }
489: /*
490: Defines the action of the downsmoother
491: */
492: PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
493: {
494: PetscErrorCode 0;
495: SNESConvergedReason reason;
496: Vec FPC;
497: SNES smoothd;
498: SNES_FAS *fas = (SNES_FAS*) snes->data;
501: SNESFASCycleGetSmootherDown(snes, &smoothd);
502: SNESSetInitialFunction(smoothd, F);
503: if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,0,0,0,0);}
504: SNESSolve(smoothd, B, X);
505: if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,0,0,0,0);}
506: /* check convergence reason for the smoother */
507: SNESGetConvergedReason(smoothd,&reason);
508: if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
509: snes->reason = SNES_DIVERGED_INNER;
510: return(0);
511: }
512: SNESGetFunction(smoothd, &FPC, NULL, NULL);
513: VecCopy(FPC, F);
514: if (fnorm) {VecNorm(F,NORM_2,fnorm);}
515: return(0);
516: }
521: /*
522: Defines the action of the upsmoother
523: */
524: PetscErrorCode SNESFASUpSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
525: {
526: PetscErrorCode 0;
527: SNESConvergedReason reason;
528: Vec FPC;
529: SNES smoothu;
530: SNES_FAS *fas = (SNES_FAS*) snes->data;
533: SNESFASCycleGetSmootherUp(snes, &smoothu);
534: if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,0,0,0,0);}
535: SNESSolve(smoothu, B, X);
536: if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,0,0,0,0);}
537: /* check convergence reason for the smoother */
538: SNESGetConvergedReason(smoothu,&reason);
539: if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
540: snes->reason = SNES_DIVERGED_INNER;
541: return(0);
542: }
543: SNESGetFunction(smoothu, &FPC, NULL, NULL);
544: VecCopy(FPC, F);
545: if (fnorm) {VecNorm(F,NORM_2,fnorm);}
546: return(0);
547: }
551: /*@
552: SNESFASCreateCoarseVec - create Vec corresponding to a state vector on one level coarser than current level
554: Collective
556: Input Arguments:
557: . snes - SNESFAS
559: Output Arguments:
560: . Xcoarse - vector on level one coarser than snes
562: Level: developer
564: .seealso: SNESFASSetRestriction(), SNESFASRestrict()
565: @*/
566: PetscErrorCode SNESFASCreateCoarseVec(SNES snes,Vec *Xcoarse)
567: {
569: SNES_FAS *fas = (SNES_FAS*)snes->data;
572: if (fas->rscale) {
573: VecDuplicate(fas->rscale,Xcoarse);
574: } else if (fas->inject) {
575: MatGetVecs(fas->inject,Xcoarse,NULL);
576: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or injection");
577: return(0);
578: }
582: /*@
583: SNESFASRestrict - restrict a Vec to the next coarser level
585: Collective
587: Input Arguments:
588: + fine - SNES from which to restrict
589: - Xfine - vector to restrict
591: Output Arguments:
592: . Xcoarse - result of restriction
594: Level: developer
596: .seealso: SNESFASSetRestriction(), SNESFASSetInjection()
597: @*/
598: PetscErrorCode SNESFASRestrict(SNES fine,Vec Xfine,Vec Xcoarse)
599: {
601: SNES_FAS *fas = (SNES_FAS*)fine->data;
607: if (fas->inject) {
608: MatRestrict(fas->inject,Xfine,Xcoarse);
609: } else {
610: MatRestrict(fas->restrct,Xfine,Xcoarse);
611: VecPointwiseMult(Xcoarse,fas->rscale,Xcoarse);
612: }
613: return(0);
614: }
618: /*
620: Performs the FAS coarse correction as:
622: fine problem: F(x) = b
623: coarse problem: F^c(x^c) = b^c
625: b^c = F^c(Rx) - R(F(x) - b)
627: */
628: PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new)
629: {
630: PetscErrorCode ierr;
631: Vec X_c, Xo_c, F_c, B_c;
632: SNESConvergedReason reason;
633: SNES next;
634: Mat restrct, interpolate;
635: SNES_FAS *fasc;
638: SNESFASCycleGetCorrection(snes, &next);
639: if (next) {
640: fasc = (SNES_FAS*)next->data;
642: SNESFASCycleGetRestriction(snes, &restrct);
643: SNESFASCycleGetInterpolation(snes, &interpolate);
645: X_c = next->vec_sol;
646: Xo_c = next->work[0];
647: F_c = next->vec_func;
648: B_c = next->vec_rhs;
650: if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,0,0,0,0);}
651: SNESFASRestrict(snes,X,Xo_c);
652: /* restrict the defect */
653: MatRestrict(restrct, F, B_c);
654: if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,0,0,0,0);}
656: if (fasc->eventresidual) {PetscLogEventBegin(fasc->eventresidual,0,0,0,0);}
657: SNESComputeFunction(next, Xo_c, F_c);
658: if (fasc->eventresidual) {PetscLogEventEnd(fasc->eventresidual,0,0,0,0);}
660: /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */
661: VecCopy(B_c, X_c);
662: VecCopy(F_c, B_c);
663: VecCopy(X_c, F_c);
664: /* set initial guess of the coarse problem to the projected fine solution */
665: VecCopy(Xo_c, X_c);
667: /* recurse to the next level */
668: SNESSetInitialFunction(next, F_c);
669: SNESSolve(next, B_c, X_c);
670: SNESGetConvergedReason(next,&reason);
671: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
672: snes->reason = SNES_DIVERGED_INNER;
673: return(0);
674: }
675: /* correct as x <- x + I(x^c - Rx)*/
676: VecAXPY(X_c, -1.0, Xo_c);
678: if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,0,0,0,0);}
679: MatInterpolateAdd(interpolate, X_c, X, X_new);
680: if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,0,0,0,0);}
681: }
682: return(0);
683: }
687: /*
689: The additive cycle looks like:
691: xhat = x
692: xhat = dS(x, b)
693: x = coarsecorrection(xhat, b_d)
694: x = x + nu*(xhat - x);
695: (optional) x = uS(x, b)
697: With the coarse RHS (defect correction) as below.
699: */
700: PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X)
701: {
702: Vec F, B, Xhat;
703: Vec X_c, Xo_c, F_c, B_c;
704: PetscErrorCode ierr;
705: SNESConvergedReason reason;
706: PetscReal xnorm, fnorm, ynorm;
707: PetscBool lssuccess;
708: SNES next;
709: Mat restrct, interpolate;
710: SNES_FAS *fas = (SNES_FAS*)snes->data,*fasc;
713: SNESFASCycleGetCorrection(snes, &next);
714: F = snes->vec_func;
715: B = snes->vec_rhs;
716: Xhat = snes->work[1];
717: VecCopy(X, Xhat);
718: /* recurse first */
719: if (next) {
720: fasc = (SNES_FAS*)next->data;
721: SNESFASCycleGetRestriction(snes, &restrct);
722: SNESFASCycleGetInterpolation(snes, &interpolate);
723: if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,0,0,0,0);}
724: SNESComputeFunction(snes, Xhat, F);
725: if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,0,0,0,0);}
726: VecNorm(F, NORM_2, &fnorm);
727: X_c = next->vec_sol;
728: Xo_c = next->work[0];
729: F_c = next->vec_func;
730: B_c = next->vec_rhs;
732: SNESFASRestrict(snes,Xhat,Xo_c);
733: /* restrict the defect */
734: MatRestrict(restrct, F, B_c);
736: /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */
737: if (fasc->eventresidual) {PetscLogEventBegin(fasc->eventresidual,0,0,0,0);}
738: SNESComputeFunction(next, Xo_c, F_c);
739: if (fasc->eventresidual) {PetscLogEventEnd(fasc->eventresidual,0,0,0,0);}
740: VecCopy(B_c, X_c);
741: VecCopy(F_c, B_c);
742: VecCopy(X_c, F_c);
743: /* set initial guess of the coarse problem to the projected fine solution */
744: VecCopy(Xo_c, X_c);
746: /* recurse */
747: SNESSetInitialFunction(next, F_c);
748: SNESSolve(next, B_c, X_c);
750: /* smooth on this level */
751: SNESFASDownSmooth_Private(snes, B, X, F, &fnorm);
753: SNESGetConvergedReason(next,&reason);
754: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
755: snes->reason = SNES_DIVERGED_INNER;
756: return(0);
757: }
759: /* correct as x <- x + I(x^c - Rx)*/
760: VecAYPX(X_c, -1.0, Xo_c);
761: MatInterpolate(interpolate, X_c, Xhat);
763: /* additive correction of the coarse direction*/
764: SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat);
765: SNESLineSearchGetSuccess(snes->linesearch, &lssuccess);
766: if (!lssuccess) {
767: if (++snes->numFailures >= snes->maxFailures) {
768: snes->reason = SNES_DIVERGED_LINE_SEARCH;
769: return(0);
770: }
771: }
772: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm);
773: } else {
774: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
775: }
776: return(0);
777: }
781: /*
783: Defines the FAS cycle as:
785: fine problem: F(x) = b
786: coarse problem: F^c(x) = b^c
788: b^c = F^c(Rx) - R(F(x) - b)
790: correction:
792: x = x + I(x^c - Rx)
794: */
795: PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X)
796: {
799: Vec F,B;
800: SNES next;
803: F = snes->vec_func;
804: B = snes->vec_rhs;
805: /* pre-smooth -- just update using the pre-smoother */
806: SNESFASCycleGetCorrection(snes,&next);
807: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
808: if (next) {
809: SNESFASCoarseCorrection(snes, X, F, X);
810: SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
811: }
812: return(0);
813: }
817: PetscErrorCode SNESFASCycleSetupPhase_Full(SNES snes)
818: {
819: SNES next;
820: SNES_FAS *fas = (SNES_FAS*)snes->data;
821: PetscBool isFine;
825: /* pre-smooth -- just update using the pre-smoother */
826: SNESFASCycleIsFine(snes,&isFine);
827: SNESFASCycleGetCorrection(snes,&next);
828: fas->full_stage = 0;
829: if (next) {SNESFASCycleSetupPhase_Full(next);}
830: return(0);
831: }
835: PetscErrorCode SNESFASCycle_Full(SNES snes, Vec X)
836: {
838: Vec F,B;
839: SNES_FAS *fas = (SNES_FAS*)snes->data;
840: PetscBool isFine;
841: SNES next;
844: F = snes->vec_func;
845: B = snes->vec_rhs;
846: SNESFASCycleIsFine(snes,&isFine);
847: SNESFASCycleGetCorrection(snes,&next);
849: if (isFine) {
850: SNESFASCycleSetupPhase_Full(snes);
851: }
853: if (fas->full_stage == 0) {
854: /* downsweep */
855: if (next) {
856: if (fas->level != 1) next->max_its += 1;
857: if (fas->full_downsweep||isFine) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
858: SNESFASCoarseCorrection(snes,X,F,X);
859: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
860: if (fas->level != 1) next->max_its -= 1;
861: } else {
862: SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
863: }
864: fas->full_stage = 1;
865: } else if (fas->full_stage == 1) {
866: if (snes->iter == 0) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
867: if (next) {
868: SNESFASCoarseCorrection(snes,X,F,X);
869: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
870: }
871: }
872: /* final v-cycle */
873: if (isFine) {
874: if (next) {
875: SNESFASCoarseCorrection(snes,X,F,X);
876: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
877: }
878: }
879: return(0);
880: }
884: PetscErrorCode SNESFASCycle_Kaskade(SNES snes, Vec X)
885: {
887: Vec F,B;
888: SNES next;
891: F = snes->vec_func;
892: B = snes->vec_rhs;
893: SNESFASCycleGetCorrection(snes,&next);
894: if (next) {
895: SNESFASCoarseCorrection(snes,X,F,X);
896: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
897: } else {
898: SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
899: }
900: return(0);
901: }
903: PetscBool SNEScite = PETSC_FALSE;
904: const char SNESCitation[] = "@techreport{pbmkbsxt2012,\n"
905: " title = {Composing Scalable Nonlinear Algebraic Solvers},\n"
906: " author = {Peter Brune and Mathew Knepley and Barry Smith and Xuemin Tu},\n"
907: " year = 2013,\n"
908: " type = Preprint,\n"
909: " number = {ANL/MCS-P2010-0112},\n"
910: " institution = {Argonne National Laboratory}\n}\n";
914: PetscErrorCode SNESSolve_FAS(SNES snes)
915: {
917: PetscInt i, maxits;
918: Vec X, F;
919: PetscReal fnorm;
920: SNES_FAS *fas = (SNES_FAS*)snes->data,*ffas;
921: DM dm;
922: PetscBool isFine;
925: PetscCitationsRegister(SNESCitation,&SNEScite);
926: maxits = snes->max_its; /* maximum number of iterations */
927: snes->reason = SNES_CONVERGED_ITERATING;
928: X = snes->vec_sol;
929: F = snes->vec_func;
931: SNESFASCycleIsFine(snes, &isFine);
932: /*norm setup */
933: PetscObjectSAWsTakeAccess((PetscObject)snes);
934: snes->iter = 0;
935: snes->norm = 0.;
936: PetscObjectSAWsGrantAccess((PetscObject)snes);
937: if (!snes->vec_func_init_set) {
938: if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,0,0,0,0);}
939: SNESComputeFunction(snes,X,F);
940: if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,0,0,0,0);}
941: if (snes->domainerror) {
942: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
943: return(0);
944: }
945: } else snes->vec_func_init_set = PETSC_FALSE;
947: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
948: if (PetscIsInfOrNanReal(fnorm)) {
949: snes->reason = SNES_DIVERGED_FNORM_NAN;
950: return(0);
951: }
953: PetscObjectSAWsTakeAccess((PetscObject)snes);
954: snes->norm = fnorm;
955: PetscObjectSAWsGrantAccess((PetscObject)snes);
956: SNESLogConvergenceHistory(snes,fnorm,0);
957: SNESMonitor(snes,0,fnorm);
959: /* test convergence */
960: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
961: if (snes->reason) return(0);
964: if (isFine) {
965: /* propagate scale-dependent data up the hierarchy */
966: SNESGetDM(snes,&dm);
967: for (ffas=fas; ffas->next; ffas=(SNES_FAS*)ffas->next->data) {
968: DM dmcoarse;
969: SNESGetDM(ffas->next,&dmcoarse);
970: DMRestrict(dm,ffas->restrct,ffas->rscale,ffas->inject,dmcoarse);
971: dm = dmcoarse;
972: }
973: }
975: for (i = 0; i < maxits; i++) {
976: /* Call general purpose update function */
978: if (snes->ops->update) {
979: (*snes->ops->update)(snes, snes->iter);
980: }
981: if (fas->fastype == SNES_FAS_MULTIPLICATIVE) {
982: SNESFASCycle_Multiplicative(snes, X);
983: } else if (fas->fastype == SNES_FAS_ADDITIVE) {
984: SNESFASCycle_Additive(snes, X);
985: } else if (fas->fastype == SNES_FAS_FULL) {
986: SNESFASCycle_Full(snes, X);
987: } else if (fas->fastype ==SNES_FAS_KASKADE) {
988: SNESFASCycle_Kaskade(snes, X);
989: } else {
990: SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported FAS type");
991: }
993: /* check for FAS cycle divergence */
994: if (snes->reason != SNES_CONVERGED_ITERATING) return(0);
996: /* Monitor convergence */
997: PetscObjectSAWsTakeAccess((PetscObject)snes);
998: snes->iter = i+1;
999: PetscObjectSAWsGrantAccess((PetscObject)snes);
1000: SNESLogConvergenceHistory(snes,snes->norm,0);
1001: SNESMonitor(snes,snes->iter,snes->norm);
1002: /* Test for convergence */
1003: if (isFine) {
1004: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,snes->norm,&snes->reason,snes->cnvP);
1005: if (snes->reason) break;
1006: }
1007: }
1008: if (i == maxits) {
1009: PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);
1010: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
1011: }
1012: return(0);
1013: }