Actual source code: fas.c
petsc-3.6.4 2016-04-12
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(PetscOptions *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_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: PetscInt dm_levels;
145: Vec vec_sol, vec_func, vec_sol_update, vec_rhs; /* preserve these if they're set through the reset */
146: SNES next;
147: PetscBool isFine;
148: SNESLineSearch linesearch;
149: SNESLineSearch slinesearch;
150: void *lsprectx,*lspostctx;
151: PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*);
152: PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
155: SNESFASCycleIsFine(snes, &isFine);
156: if (fas->usedmfornumberoflevels && isFine) {
157: DMGetRefineLevel(snes->dm,&dm_levels);
158: dm_levels++;
159: if (dm_levels > fas->levels) {
160: /* we don't want the solution and func vectors to be destroyed in the SNESReset when it's called in SNESFASSetLevels_FAS*/
161: vec_sol = snes->vec_sol;
162: vec_func = snes->vec_func;
163: vec_sol_update = snes->vec_sol_update;
164: vec_rhs = snes->vec_rhs;
165: snes->vec_sol = NULL;
166: snes->vec_func = NULL;
167: snes->vec_sol_update = NULL;
168: snes->vec_rhs = NULL;
170: /* reset the number of levels */
171: SNESFASSetLevels(snes,dm_levels,NULL);
172: SNESSetFromOptions(snes);
174: snes->vec_sol = vec_sol;
175: snes->vec_func = vec_func;
176: snes->vec_rhs = vec_rhs;
177: snes->vec_sol_update = vec_sol_update;
178: }
179: }
180: SNESFASCycleGetCorrection(snes, &next);
181: if (!isFine) snes->gridsequence = 0; /* no grid sequencing inside the multigrid hierarchy! */
183: SNESSetWorkVecs(snes, 2); /* work vectors used for intergrid transfers */
185: /* set up the smoothers if they haven't already been set up */
186: if (!fas->smoothd) {
187: SNESFASCycleCreateSmoother_Private(snes, &fas->smoothd);
188: }
190: if (snes->dm) {
191: /* set the smoother DMs properly */
192: if (fas->smoothu) SNESSetDM(fas->smoothu, snes->dm);
193: SNESSetDM(fas->smoothd, snes->dm);
194: /* construct EVERYTHING from the DM -- including the progressive set of smoothers */
195: if (next) {
196: /* for now -- assume the DM and the evaluation functions have been set externally */
197: if (!next->dm) {
198: DMCoarsen(snes->dm, PetscObjectComm((PetscObject)next), &next->dm);
199: SNESSetDM(next, next->dm);
200: }
201: /* set the interpolation and restriction from the DM */
202: if (!fas->interpolate) {
203: DMCreateInterpolation(next->dm, snes->dm, &fas->interpolate, &fas->rscale);
204: if (!fas->restrct) {
205: PetscObjectReference((PetscObject)fas->interpolate);
206: fas->restrct = fas->interpolate;
207: }
208: }
209: /* set the injection from the DM */
210: if (!fas->inject) {
211: DMCreateInjection(next->dm, snes->dm, &fas->inject);
212: }
213: }
214: }
215: /*pass the smoother, function, and jacobian up to the next level if it's not user set already */
216: if (fas->galerkin) {
217: if (next) {
218: SNESSetFunction(next, NULL, SNESFASGalerkinDefaultFunction, next);
219: }
220: if (fas->smoothd && fas->level != fas->levels - 1) {
221: SNESSetFunction(fas->smoothd, NULL, SNESFASGalerkinDefaultFunction, snes);
222: }
223: if (fas->smoothu && fas->level != fas->levels - 1) {
224: SNESSetFunction(fas->smoothu, NULL, SNESFASGalerkinDefaultFunction, snes);
225: }
226: }
228: /* sets the down (pre) smoother's default norm and sets it from options */
229: if (fas->smoothd) {
230: if (fas->level == 0 && fas->levels != 1) {
231: SNESSetNormSchedule(fas->smoothd, SNES_NORM_NONE);
232: } else {
233: SNESSetNormSchedule(fas->smoothd, SNES_NORM_FINAL_ONLY);
234: }
235: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothd);
236: SNESSetFromOptions(fas->smoothd);
237: SNESGetLineSearch(snes,&linesearch);
238: SNESGetLineSearch(fas->smoothd,&slinesearch);
239: SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
240: SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
241: SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
242: SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
243: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
245: fas->smoothd->vec_sol = snes->vec_sol;
246: PetscObjectReference((PetscObject)snes->vec_sol);
247: fas->smoothd->vec_sol_update = snes->vec_sol_update;
248: PetscObjectReference((PetscObject)snes->vec_sol_update);
249: fas->smoothd->vec_func = snes->vec_func;
250: PetscObjectReference((PetscObject)snes->vec_func);
252: if (fas->eventsmoothsetup) {PetscLogEventBegin(fas->eventsmoothsetup,0,0,0,0);}
253: SNESSetUp(fas->smoothd);
254: if (fas->eventsmoothsetup) {PetscLogEventEnd(fas->eventsmoothsetup,0,0,0,0);}
255: }
257: /* sets the up (post) smoother's default norm and sets it from options */
258: if (fas->smoothu) {
259: if (fas->level != fas->levels - 1) {
260: SNESSetNormSchedule(fas->smoothu, SNES_NORM_NONE);
261: } else {
262: SNESSetNormSchedule(fas->smoothu, SNES_NORM_FINAL_ONLY);
263: }
264: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothu);
265: SNESSetFromOptions(fas->smoothu);
266: SNESGetLineSearch(snes,&linesearch);
267: SNESGetLineSearch(fas->smoothu,&slinesearch);
268: SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
269: SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
270: SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
271: SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
272: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
274: fas->smoothu->vec_sol = snes->vec_sol;
275: PetscObjectReference((PetscObject)snes->vec_sol);
276: fas->smoothu->vec_sol_update = snes->vec_sol_update;
277: PetscObjectReference((PetscObject)snes->vec_sol_update);
278: fas->smoothu->vec_func = snes->vec_func;
279: PetscObjectReference((PetscObject)snes->vec_func);
281: if (fas->eventsmoothsetup) {PetscLogEventBegin(fas->eventsmoothsetup,0,0,0,0);}
282: SNESSetUp(fas->smoothu);
283: if (fas->eventsmoothsetup) {PetscLogEventEnd(fas->eventsmoothsetup,0,0,0,0);}
285: }
287: if (next) {
288: /* gotta set up the solution vector for this to work */
289: if (!next->vec_sol) {SNESFASCreateCoarseVec(snes,&next->vec_sol);}
290: if (!next->vec_rhs) {SNESFASCreateCoarseVec(snes,&next->vec_rhs);}
291: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next);
292: SNESGetLineSearch(snes,&linesearch);
293: SNESGetLineSearch(fas->next,&slinesearch);
294: SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
295: SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
296: SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
297: SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
298: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
299: SNESSetUp(next);
300: }
301: /* setup FAS work vectors */
302: if (fas->galerkin) {
303: VecDuplicate(snes->vec_sol, &fas->Xg);
304: VecDuplicate(snes->vec_sol, &fas->Fg);
305: }
306: return(0);
307: }
311: PetscErrorCode SNESSetFromOptions_FAS(PetscOptions *PetscOptionsObject,SNES snes)
312: {
313: SNES_FAS *fas = (SNES_FAS*) snes->data;
314: PetscInt levels = 1;
315: PetscBool flg = PETSC_FALSE, upflg = PETSC_FALSE, downflg = PETSC_FALSE, monflg = PETSC_FALSE, galerkinflg = PETSC_FALSE,continuationflg = PETSC_FALSE;
317: char monfilename[PETSC_MAX_PATH_LEN];
318: SNESFASType fastype;
319: const char *optionsprefix;
320: SNESLineSearch linesearch;
321: PetscInt m, n_up, n_down;
322: SNES next;
323: PetscBool isFine;
326: SNESFASCycleIsFine(snes, &isFine);
327: PetscOptionsHead(PetscOptionsObject,"SNESFAS Options-----------------------------------");
329: /* number of levels -- only process most options on the finest level */
330: if (isFine) {
331: PetscOptionsInt("-snes_fas_levels", "Number of Levels", "SNESFASSetLevels", levels, &levels, &flg);
332: if (!flg && snes->dm) {
333: DMGetRefineLevel(snes->dm,&levels);
334: levels++;
335: fas->usedmfornumberoflevels = PETSC_TRUE;
336: }
337: SNESFASSetLevels(snes, levels, NULL);
338: fastype = fas->fastype;
339: PetscOptionsEnum("-snes_fas_type","FAS correction type","SNESFASSetType",SNESFASTypes,(PetscEnum)fastype,(PetscEnum*)&fastype,&flg);
340: if (flg) {
341: SNESFASSetType(snes, fastype);
342: }
344: SNESGetOptionsPrefix(snes, &optionsprefix);
345: PetscOptionsInt("-snes_fas_cycles","Number of cycles","SNESFASSetCycles",fas->n_cycles,&m,&flg);
346: if (flg) {
347: SNESFASSetCycles(snes, m);
348: }
349: PetscOptionsBool("-snes_fas_continuation","Corrected grid-sequence continuation","SNESFASSetContinuation",fas->continuation,&continuationflg,&flg);
350: if (flg) {
351: SNESFASSetContinuation(snes,continuationflg);
352: }
354: PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin","SNESFASSetGalerkin",fas->galerkin,&galerkinflg,&flg);
355: if (flg) {
356: SNESFASSetGalerkin(snes, galerkinflg);
357: }
359: if (fas->fastype == SNES_FAS_FULL) {
360: PetscOptionsBool("-snes_fas_full_downsweep","Smooth on the initial upsweep for full FAS cycles","SNESFASFullSetDownSweep",fas->full_downsweep,&fas->full_downsweep,&flg);
361: if (flg) {SNESFASFullSetDownSweep(snes,fas->full_downsweep);}
362: }
364: PetscOptionsInt("-snes_fas_smoothup","Number of post-smoothing steps","SNESFASSetNumberSmoothUp",fas->max_up_it,&n_up,&upflg);
366: PetscOptionsInt("-snes_fas_smoothdown","Number of pre-smoothing steps","SNESFASSetNumberSmoothDown",fas->max_down_it,&n_down,&downflg);
368: PetscOptionsString("-snes_fas_monitor","Monitor FAS progress","SNESFASSetMonitor","stdout",monfilename,PETSC_MAX_PATH_LEN,&monflg);
369: if (monflg) SNESFASSetMonitor(snes, PETSC_TRUE);
371: flg = PETSC_FALSE;
372: monflg = PETSC_TRUE;
373: PetscOptionsBool("-snes_fas_log","Log times for each FAS level","SNESFASSetLog",monflg,&monflg,&flg);
374: if (flg) {SNESFASSetLog(snes,monflg);}
375: }
377: PetscOptionsTail();
378: /* setup from the determined types if there is no pointwise procedure or smoother defined */
379: if (upflg) {
380: SNESFASSetNumberSmoothUp(snes,n_up);
381: }
382: if (downflg) {
383: SNESFASSetNumberSmoothDown(snes,n_down);
384: }
386: /* set up the default line search for coarse grid corrections */
387: if (fas->fastype == SNES_FAS_ADDITIVE) {
388: if (!snes->linesearch) {
389: SNESGetLineSearch(snes, &linesearch);
390: SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
391: }
392: }
394: SNESFASCycleGetCorrection(snes, &next);
395: /* recursive option setting for the smoothers */
396: if (next) {SNESSetFromOptions(next);}
397: return(0);
398: }
400: #include <petscdraw.h>
403: PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer)
404: {
405: SNES_FAS *fas = (SNES_FAS*) snes->data;
406: PetscBool isFine,iascii,isdraw;
407: PetscInt i;
409: SNES smoothu, smoothd, levelsnes;
412: SNESFASCycleIsFine(snes, &isFine);
413: if (isFine) {
414: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
415: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
416: if (iascii) {
417: PetscViewerASCIIPrintf(viewer, "FAS: type is %s, levels=%D, cycles=%D\n", SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles);
418: if (fas->galerkin) {
419: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid function evaluation\n");
420: } else {
421: PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid function evaluation\n");
422: }
423: for (i=0; i<fas->levels; i++) {
424: SNESFASGetCycleSNES(snes, i, &levelsnes);
425: SNESFASCycleGetSmootherUp(levelsnes, &smoothu);
426: SNESFASCycleGetSmootherDown(levelsnes, &smoothd);
427: if (!i) {
428: PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level %D -------------------------------\n",i);
429: } else {
430: PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
431: }
432: PetscViewerASCIIPushTab(viewer);
433: if (smoothd) {
434: SNESView(smoothd,viewer);
435: } else {
436: PetscViewerASCIIPrintf(viewer,"Not yet available\n");
437: }
438: PetscViewerASCIIPopTab(viewer);
439: if (i && (smoothd == smoothu)) {
440: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
441: } else if (i) {
442: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
443: PetscViewerASCIIPushTab(viewer);
444: if (smoothu) {
445: SNESView(smoothu,viewer);
446: } else {
447: PetscViewerASCIIPrintf(viewer,"Not yet available\n");
448: }
449: PetscViewerASCIIPopTab(viewer);
450: }
451: }
452: } else if (isdraw) {
453: PetscDraw draw;
454: PetscReal x,w,y,bottom,th,wth;
455: SNES_FAS *curfas = fas;
456: PetscViewerDrawGetDraw(viewer,0,&draw);
457: PetscDrawGetCurrentPoint(draw,&x,&y);
458: PetscDrawStringGetSize(draw,&wth,&th);
459: bottom = y - th;
460: while (curfas) {
461: if (!curfas->smoothu) {
462: PetscDrawPushCurrentPoint(draw,x,bottom);
463: if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
464: PetscDrawPopCurrentPoint(draw);
465: } else {
466: w = 0.5*PetscMin(1.0-x,x);
467: PetscDrawPushCurrentPoint(draw,x-w,bottom);
468: if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
469: PetscDrawPopCurrentPoint(draw);
470: PetscDrawPushCurrentPoint(draw,x+w,bottom);
471: if (curfas->smoothu) SNESView(curfas->smoothu,viewer);
472: PetscDrawPopCurrentPoint(draw);
473: }
474: /* this is totally bogus but we have no way of knowing how low the previous one was draw to */
475: bottom -= 5*th;
476: if (curfas->next) curfas = (SNES_FAS*)curfas->next->data;
477: else curfas = NULL;
478: }
479: }
480: }
481: return(0);
482: }
486: /*
487: Defines the action of the downsmoother
488: */
489: PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
490: {
491: PetscErrorCode 0;
492: SNESConvergedReason reason;
493: Vec FPC;
494: SNES smoothd;
495: SNES_FAS *fas = (SNES_FAS*) snes->data;
498: SNESFASCycleGetSmootherDown(snes, &smoothd);
499: SNESSetInitialFunction(smoothd, F);
500: if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,0,0,0,0);}
501: SNESSolve(smoothd, B, X);
502: if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,0,0,0,0);}
503: /* check convergence reason for the smoother */
504: SNESGetConvergedReason(smoothd,&reason);
505: if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
506: snes->reason = SNES_DIVERGED_INNER;
507: return(0);
508: }
509: SNESGetFunction(smoothd, &FPC, NULL, NULL);
510: VecCopy(FPC, F);
511: if (fnorm) {VecNorm(F,NORM_2,fnorm);}
512: return(0);
513: }
518: /*
519: Defines the action of the upsmoother
520: */
521: PetscErrorCode SNESFASUpSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
522: {
523: PetscErrorCode 0;
524: SNESConvergedReason reason;
525: Vec FPC;
526: SNES smoothu;
527: SNES_FAS *fas = (SNES_FAS*) snes->data;
530: SNESFASCycleGetSmootherUp(snes, &smoothu);
531: if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,0,0,0,0);}
532: SNESSolve(smoothu, B, X);
533: if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,0,0,0,0);}
534: /* check convergence reason for the smoother */
535: SNESGetConvergedReason(smoothu,&reason);
536: if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
537: snes->reason = SNES_DIVERGED_INNER;
538: return(0);
539: }
540: SNESGetFunction(smoothu, &FPC, NULL, NULL);
541: VecCopy(FPC, F);
542: if (fnorm) {VecNorm(F,NORM_2,fnorm);}
543: return(0);
544: }
548: /*@
549: SNESFASCreateCoarseVec - create Vec corresponding to a state vector on one level coarser than current level
551: Collective
553: Input Arguments:
554: . snes - SNESFAS
556: Output Arguments:
557: . Xcoarse - vector on level one coarser than snes
559: Level: developer
561: .seealso: SNESFASSetRestriction(), SNESFASRestrict()
562: @*/
563: PetscErrorCode SNESFASCreateCoarseVec(SNES snes,Vec *Xcoarse)
564: {
566: SNES_FAS *fas = (SNES_FAS*)snes->data;
569: if (fas->rscale) {
570: VecDuplicate(fas->rscale,Xcoarse);
571: } else if (fas->inject) {
572: MatCreateVecs(fas->inject,Xcoarse,NULL);
573: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or injection");
574: return(0);
575: }
579: /*@
580: SNESFASRestrict - restrict a Vec to the next coarser level
582: Collective
584: Input Arguments:
585: + fine - SNES from which to restrict
586: - Xfine - vector to restrict
588: Output Arguments:
589: . Xcoarse - result of restriction
591: Level: developer
593: .seealso: SNESFASSetRestriction(), SNESFASSetInjection()
594: @*/
595: PetscErrorCode SNESFASRestrict(SNES fine,Vec Xfine,Vec Xcoarse)
596: {
598: SNES_FAS *fas = (SNES_FAS*)fine->data;
604: if (fas->inject) {
605: MatRestrict(fas->inject,Xfine,Xcoarse);
606: } else {
607: MatRestrict(fas->restrct,Xfine,Xcoarse);
608: VecPointwiseMult(Xcoarse,fas->rscale,Xcoarse);
609: }
610: return(0);
611: }
615: /*
617: Performs the FAS coarse correction as:
619: fine problem: F(x) = b
620: coarse problem: F^c(x^c) = b^c
622: b^c = F^c(Rx) - R(F(x) - b)
624: */
625: PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new)
626: {
627: PetscErrorCode ierr;
628: Vec X_c, Xo_c, F_c, B_c;
629: SNESConvergedReason reason;
630: SNES next;
631: Mat restrct, interpolate;
632: SNES_FAS *fasc;
635: SNESFASCycleGetCorrection(snes, &next);
636: if (next) {
637: fasc = (SNES_FAS*)next->data;
639: SNESFASCycleGetRestriction(snes, &restrct);
640: SNESFASCycleGetInterpolation(snes, &interpolate);
642: X_c = next->vec_sol;
643: Xo_c = next->work[0];
644: F_c = next->vec_func;
645: B_c = next->vec_rhs;
647: if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,0,0,0,0);}
648: SNESFASRestrict(snes,X,Xo_c);
649: /* restrict the defect: R(F(x) - b) */
650: MatRestrict(restrct, F, B_c);
651: if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,0,0,0,0);}
653: if (fasc->eventresidual) {PetscLogEventBegin(fasc->eventresidual,0,0,0,0);}
654: /* F_c = F^c(Rx) - R(F(x) - b) since the second term was sitting in next->vec_rhs */
655: SNESComputeFunction(next, Xo_c, F_c);
656: if (fasc->eventresidual) {PetscLogEventEnd(fasc->eventresidual,0,0,0,0);}
658: /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */
659: VecCopy(B_c, X_c);
660: VecCopy(F_c, B_c);
661: VecCopy(X_c, F_c);
662: /* set initial guess of the coarse problem to the projected fine solution */
663: VecCopy(Xo_c, X_c);
665: /* recurse to the next level */
666: SNESSetInitialFunction(next, F_c);
667: SNESSolve(next, B_c, X_c);
668: SNESGetConvergedReason(next,&reason);
669: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
670: snes->reason = SNES_DIVERGED_INNER;
671: return(0);
672: }
673: /* correct as x <- x + I(x^c - Rx)*/
674: VecAXPY(X_c, -1.0, Xo_c);
676: if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,0,0,0,0);}
677: MatInterpolateAdd(interpolate, X_c, X, X_new);
678: if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,0,0,0,0);}
679: }
680: return(0);
681: }
685: /*
687: The additive cycle looks like:
689: xhat = x
690: xhat = dS(x, b)
691: x = coarsecorrection(xhat, b_d)
692: x = x + nu*(xhat - x);
693: (optional) x = uS(x, b)
695: With the coarse RHS (defect correction) as below.
697: */
698: PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X)
699: {
700: Vec F, B, Xhat;
701: Vec X_c, Xo_c, F_c, B_c;
702: PetscErrorCode ierr;
703: SNESConvergedReason reason;
704: PetscReal xnorm, fnorm, ynorm;
705: SNESLineSearchReason lsresult;
706: SNES next;
707: Mat restrct, interpolate;
708: SNES_FAS *fas = (SNES_FAS*)snes->data,*fasc;
711: SNESFASCycleGetCorrection(snes, &next);
712: F = snes->vec_func;
713: B = snes->vec_rhs;
714: Xhat = snes->work[1];
715: VecCopy(X, Xhat);
716: /* recurse first */
717: if (next) {
718: fasc = (SNES_FAS*)next->data;
719: SNESFASCycleGetRestriction(snes, &restrct);
720: SNESFASCycleGetInterpolation(snes, &interpolate);
721: if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,0,0,0,0);}
722: SNESComputeFunction(snes, Xhat, F);
723: if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,0,0,0,0);}
724: VecNorm(F, NORM_2, &fnorm);
725: X_c = next->vec_sol;
726: Xo_c = next->work[0];
727: F_c = next->vec_func;
728: B_c = next->vec_rhs;
730: SNESFASRestrict(snes,Xhat,Xo_c);
731: /* restrict the defect */
732: MatRestrict(restrct, F, B_c);
734: /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */
735: if (fasc->eventresidual) {PetscLogEventBegin(fasc->eventresidual,0,0,0,0);}
736: SNESComputeFunction(next, Xo_c, F_c);
737: if (fasc->eventresidual) {PetscLogEventEnd(fasc->eventresidual,0,0,0,0);}
738: VecCopy(B_c, X_c);
739: VecCopy(F_c, B_c);
740: VecCopy(X_c, F_c);
741: /* set initial guess of the coarse problem to the projected fine solution */
742: VecCopy(Xo_c, X_c);
744: /* recurse */
745: SNESSetInitialFunction(next, F_c);
746: SNESSolve(next, B_c, X_c);
748: /* smooth on this level */
749: SNESFASDownSmooth_Private(snes, B, X, F, &fnorm);
751: SNESGetConvergedReason(next,&reason);
752: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
753: snes->reason = SNES_DIVERGED_INNER;
754: return(0);
755: }
757: /* correct as x <- x + I(x^c - Rx)*/
758: VecAYPX(X_c, -1.0, Xo_c);
759: MatInterpolate(interpolate, X_c, Xhat);
761: /* additive correction of the coarse direction*/
762: SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat);
763: SNESLineSearchGetReason(snes->linesearch, &lsresult);
764: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm);
765: if (lsresult) {
766: if (++snes->numFailures >= snes->maxFailures) {
767: snes->reason = SNES_DIVERGED_LINE_SEARCH;
768: return(0);
769: }
770: }
771: } else {
772: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
773: }
774: return(0);
775: }
779: /*
781: Defines the FAS cycle as:
783: fine problem: F(x) = b
784: coarse problem: F^c(x) = b^c
786: b^c = F^c(Rx) - R(F(x) - b)
788: correction:
790: x = x + I(x^c - Rx)
792: */
793: PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X)
794: {
797: Vec F,B;
798: SNES next;
801: F = snes->vec_func;
802: B = snes->vec_rhs;
803: /* pre-smooth -- just update using the pre-smoother */
804: SNESFASCycleGetCorrection(snes,&next);
805: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
806: if (next) {
807: SNESFASCoarseCorrection(snes, X, F, X);
808: SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
809: }
810: return(0);
811: }
815: PetscErrorCode SNESFASCycleSetupPhase_Full(SNES snes)
816: {
817: SNES next;
818: SNES_FAS *fas = (SNES_FAS*)snes->data;
819: PetscBool isFine;
823: /* pre-smooth -- just update using the pre-smoother */
824: SNESFASCycleIsFine(snes,&isFine);
825: SNESFASCycleGetCorrection(snes,&next);
826: fas->full_stage = 0;
827: if (next) {SNESFASCycleSetupPhase_Full(next);}
828: return(0);
829: }
833: PetscErrorCode SNESFASCycle_Full(SNES snes, Vec X)
834: {
836: Vec F,B;
837: SNES_FAS *fas = (SNES_FAS*)snes->data;
838: PetscBool isFine;
839: SNES next;
842: F = snes->vec_func;
843: B = snes->vec_rhs;
844: SNESFASCycleIsFine(snes,&isFine);
845: SNESFASCycleGetCorrection(snes,&next);
847: if (isFine) {
848: SNESFASCycleSetupPhase_Full(snes);
849: }
851: if (fas->full_stage == 0) {
852: /* downsweep */
853: if (next) {
854: if (fas->level != 1) next->max_its += 1;
855: if (fas->full_downsweep||isFine) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
856: SNESFASCoarseCorrection(snes,X,F,X);
857: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
858: if (fas->level != 1) next->max_its -= 1;
859: } else {
860: SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
861: }
862: fas->full_stage = 1;
863: } else if (fas->full_stage == 1) {
864: if (snes->iter == 0) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
865: if (next) {
866: SNESFASCoarseCorrection(snes,X,F,X);
867: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
868: }
869: }
870: /* final v-cycle */
871: if (isFine) {
872: if (next) {
873: SNESFASCoarseCorrection(snes,X,F,X);
874: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
875: }
876: }
877: return(0);
878: }
882: PetscErrorCode SNESFASCycle_Kaskade(SNES snes, Vec X)
883: {
885: Vec F,B;
886: SNES next;
889: F = snes->vec_func;
890: B = snes->vec_rhs;
891: SNESFASCycleGetCorrection(snes,&next);
892: if (next) {
893: SNESFASCoarseCorrection(snes,X,F,X);
894: SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
895: } else {
896: SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
897: }
898: return(0);
899: }
901: PetscBool SNEScite = PETSC_FALSE;
902: const char SNESCitation[] = "@techreport{pbmkbsxt2012,\n"
903: " title = {Composing Scalable Nonlinear Algebraic Solvers},\n"
904: " author = {Peter Brune and Mathew Knepley and Barry Smith and Xuemin Tu},\n"
905: " year = 2013,\n"
906: " type = Preprint,\n"
907: " number = {ANL/MCS-P2010-0112},\n"
908: " institution = {Argonne National Laboratory}\n}\n";
912: PetscErrorCode SNESSolve_FAS(SNES snes)
913: {
915: PetscInt i, maxits;
916: Vec X, F;
917: PetscReal fnorm;
918: SNES_FAS *fas = (SNES_FAS*)snes->data,*ffas;
919: DM dm;
920: PetscBool isFine;
924: if (snes->xl || snes->xu || snes->ops->computevariablebounds) {
925: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
926: }
928: PetscCitationsRegister(SNESCitation,&SNEScite);
929: maxits = snes->max_its; /* maximum number of iterations */
930: snes->reason = SNES_CONVERGED_ITERATING;
931: X = snes->vec_sol;
932: F = snes->vec_func;
934: SNESFASCycleIsFine(snes, &isFine);
935: /*norm setup */
936: PetscObjectSAWsTakeAccess((PetscObject)snes);
937: snes->iter = 0;
938: snes->norm = 0.;
939: PetscObjectSAWsGrantAccess((PetscObject)snes);
940: if (!snes->vec_func_init_set) {
941: if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,0,0,0,0);}
942: SNESComputeFunction(snes,X,F);
943: if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,0,0,0,0);}
944: } else snes->vec_func_init_set = PETSC_FALSE;
946: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
947: SNESCheckFunctionNorm(snes,fnorm);
948: PetscObjectSAWsTakeAccess((PetscObject)snes);
949: snes->norm = fnorm;
950: PetscObjectSAWsGrantAccess((PetscObject)snes);
951: SNESLogConvergenceHistory(snes,fnorm,0);
952: SNESMonitor(snes,0,fnorm);
954: /* test convergence */
955: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
956: if (snes->reason) return(0);
959: if (isFine) {
960: /* propagate scale-dependent data up the hierarchy */
961: SNESGetDM(snes,&dm);
962: for (ffas=fas; ffas->next; ffas=(SNES_FAS*)ffas->next->data) {
963: DM dmcoarse;
964: SNESGetDM(ffas->next,&dmcoarse);
965: DMRestrict(dm,ffas->restrct,ffas->rscale,ffas->inject,dmcoarse);
966: dm = dmcoarse;
967: }
968: }
970: for (i = 0; i < maxits; i++) {
971: /* Call general purpose update function */
973: if (snes->ops->update) {
974: (*snes->ops->update)(snes, snes->iter);
975: }
976: if (fas->fastype == SNES_FAS_MULTIPLICATIVE) {
977: SNESFASCycle_Multiplicative(snes, X);
978: } else if (fas->fastype == SNES_FAS_ADDITIVE) {
979: SNESFASCycle_Additive(snes, X);
980: } else if (fas->fastype == SNES_FAS_FULL) {
981: SNESFASCycle_Full(snes, X);
982: } else if (fas->fastype ==SNES_FAS_KASKADE) {
983: SNESFASCycle_Kaskade(snes, X);
984: } else {
985: SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported FAS type");
986: }
988: /* check for FAS cycle divergence */
989: if (snes->reason != SNES_CONVERGED_ITERATING) return(0);
991: /* Monitor convergence */
992: PetscObjectSAWsTakeAccess((PetscObject)snes);
993: snes->iter = i+1;
994: PetscObjectSAWsGrantAccess((PetscObject)snes);
995: SNESLogConvergenceHistory(snes,snes->norm,0);
996: SNESMonitor(snes,snes->iter,snes->norm);
997: /* Test for convergence */
998: if (isFine) {
999: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,snes->norm,&snes->reason,snes->cnvP);
1000: if (snes->reason) break;
1001: }
1002: }
1003: if (i == maxits) {
1004: PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);
1005: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
1006: }
1007: return(0);
1008: }