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