Actual source code: fas.c
petsc-3.3-p7 2013-05-11
1: /* Defines the basic SNES object */
2: #include <../src/snes/impls/fas/fasimpls.h> /*I "petscsnesfas.h" I*/
4: const char *SNESFASTypes[] = {"MULTIPLICATIVE","ADDITIVE","SNESFASType","SNES_FAS",0};
6: extern PetscErrorCode SNESDestroy_FAS(SNES snes);
7: extern PetscErrorCode SNESSetUp_FAS(SNES snes);
8: extern PetscErrorCode SNESSetFromOptions_FAS(SNES snes);
9: extern PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer);
10: extern PetscErrorCode SNESSolve_FAS(SNES snes);
11: extern PetscErrorCode SNESReset_FAS(SNES snes);
12: extern PetscErrorCode SNESFASGalerkinDefaultFunction(SNES, Vec, Vec, void *);
13: extern PetscErrorCode SNESFASCycleCreateSmoother_Private(SNES, SNES *);
15: /*MC
17: SNESFAS - Full Approximation Scheme nonlinear multigrid solver.
19: The nonlinear problem is solved by correction using coarse versions
20: of the nonlinear problem. This problem is perturbed so that a projected
21: solution of the fine problem elicits no correction from the coarse problem.
23: Options Database:
24: + -snes_fas_levels - The number of levels
25: . -snes_fas_cycles<1> - The number of cycles -- 1 for V, 2 for W
26: . -snes_fas_type<additive, multiplicative> - Additive or multiplicative cycle
27: . -snes_fas_galerkin<PETSC_FALSE> - Form coarse problems by projection back upon the fine problem
28: . -snes_fas_smoothup<1> - The number of iterations of the post-smoother
29: . -snes_fas_smoothdown<1> - The number of iterations of the pre-smoother
30: . -snes_fas_monitor - Monitor progress of all of the levels
31: . -fas_levels_snes_ - SNES options for all smoothers
32: . -fas_levels_cycle_snes_ - SNES options for all cycles
33: . -fas_levels_i_snes_ - SNES options for the smoothers on level i
34: . -fas_levels_i_cycle_snes_ - SNES options for the cycle on level i
35: - -fas_coarse_snes_ - SNES options for the coarsest smoother
37: Notes:
38: The organization of the FAS solver is slightly different from the organization of PCMG
39: As each level has smoother SNES instances(down and potentially up) and a cycle SNES instance.
40: The cycle SNES instance may be used for monitoring convergence on a particular level.
42: Level: beginner
44: .seealso: PCMG, SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types)
45: M*/
50: {
51: SNES_FAS * fas;
55: snes->ops->destroy = SNESDestroy_FAS;
56: snes->ops->setup = SNESSetUp_FAS;
57: snes->ops->setfromoptions = SNESSetFromOptions_FAS;
58: snes->ops->view = SNESView_FAS;
59: snes->ops->solve = SNESSolve_FAS;
60: snes->ops->reset = SNESReset_FAS;
62: snes->usesksp = PETSC_FALSE;
63: snes->usespc = PETSC_FALSE;
65: if (!snes->tolerancesset) {
66: snes->max_funcs = 30000;
67: snes->max_its = 10000;
68: }
70: PetscNewLog(snes, SNES_FAS, &fas);
71: snes->data = (void*) fas;
72: fas->level = 0;
73: fas->levels = 1;
74: fas->n_cycles = 1;
75: fas->max_up_it = 1;
76: fas->max_down_it = 1;
77: fas->smoothu = PETSC_NULL;
78: fas->smoothd = PETSC_NULL;
79: fas->next = PETSC_NULL;
80: fas->previous = PETSC_NULL;
81: fas->fine = snes;
82: fas->interpolate = PETSC_NULL;
83: fas->restrct = PETSC_NULL;
84: fas->inject = PETSC_NULL;
85: fas->monitor = PETSC_NULL;
86: fas->usedmfornumberoflevels = PETSC_FALSE;
87: fas->fastype = SNES_FAS_MULTIPLICATIVE;
88: return(0);
89: }
93: PetscErrorCode SNESReset_FAS(SNES snes)
94: {
95: PetscErrorCode 0;
96: SNES_FAS * fas = (SNES_FAS *)snes->data;
99: SNESDestroy(&fas->smoothu);
100: SNESDestroy(&fas->smoothd);
101: MatDestroy(&fas->inject);
102: MatDestroy(&fas->interpolate);
103: MatDestroy(&fas->restrct);
104: VecDestroy(&fas->rscale);
105: if (fas->next) SNESReset(fas->next);
107: return(0);
108: }
112: PetscErrorCode SNESDestroy_FAS(SNES snes)
113: {
114: SNES_FAS * fas = (SNES_FAS *)snes->data;
115: PetscErrorCode 0;
118: /* recursively resets and then destroys */
119: SNESReset(snes);
120: if (fas->next) SNESDestroy(&fas->next);
121: PetscFree(fas);
123: return(0);
124: }
128: PetscErrorCode SNESSetUp_FAS(SNES snes)
129: {
130: SNES_FAS *fas = (SNES_FAS *) snes->data;
131: PetscErrorCode ierr;
132: VecScatter injscatter;
133: PetscInt dm_levels;
134: Vec vec_sol, vec_func, vec_sol_update, vec_rhs; /* preserve these if they're set through the reset */
135: SNES next;
136: PetscBool isFine;
137: SNESLineSearch linesearch;
138: SNESLineSearch slinesearch;
139: void *lsprectx,*lspostctx;
140: SNESLineSearchPreCheckFunc lsprefunc;
141: SNESLineSearchPostCheckFunc lspostfunc;
144: SNESFASCycleIsFine(snes, &isFine);
145: if (fas->usedmfornumberoflevels && isFine) {
146: DMGetRefineLevel(snes->dm,&dm_levels);
147: dm_levels++;
148: if (dm_levels > fas->levels) {
149: /* we don't want the solution and func vectors to be destroyed in the SNESReset when it's called in SNESFASSetLevels_FAS*/
150: vec_sol = snes->vec_sol;
151: vec_func = snes->vec_func;
152: vec_sol_update = snes->vec_sol_update;
153: vec_rhs = snes->vec_rhs;
154: snes->vec_sol = PETSC_NULL;
155: snes->vec_func = PETSC_NULL;
156: snes->vec_sol_update = PETSC_NULL;
157: snes->vec_rhs = PETSC_NULL;
159: /* reset the number of levels */
160: SNESFASSetLevels(snes,dm_levels,PETSC_NULL);
161: SNESSetFromOptions(snes);
163: snes->vec_sol = vec_sol;
164: snes->vec_func = vec_func;
165: snes->vec_rhs = vec_rhs;
166: snes->vec_sol_update = vec_sol_update;
167: }
168: }
169: SNESFASCycleGetCorrection(snes, &next);
170: if (!isFine) snes->gridsequence = 0; /* no grid sequencing inside the multigrid hierarchy! */
172: if (fas->fastype == SNES_FAS_MULTIPLICATIVE) {
173: SNESDefaultGetWork(snes, 2); /* work vectors used for intergrid transfers */
174: } else {
175: SNESDefaultGetWork(snes, 2); /* work vectors used for intergrid transfers */
176: }
178: /* set up the smoothers if they haven't already been set up */
179: if (!fas->smoothd) {
180: SNESFASCycleCreateSmoother_Private(snes, &fas->smoothd);
181: }
182: if (!fas->smoothu && fas->level != fas->levels - 1) {
183: SNESFASCycleCreateSmoother_Private(snes, &fas->smoothu);
184: }
186: if (snes->dm) {
187: /* set the smoother DMs properly */
188: if (fas->smoothu) SNESSetDM(fas->smoothu, snes->dm);
189: SNESSetDM(fas->smoothd, snes->dm);
190: /* construct EVERYTHING from the DM -- including the progressive set of smoothers */
191: if (next) {
192: /* for now -- assume the DM and the evaluation functions have been set externally */
193: if (!next->dm) {
194: DMCoarsen(snes->dm, ((PetscObject)next)->comm, &next->dm);
195: SNESSetDM(next, next->dm);
196: }
197: /* set the interpolation and restriction from the DM */
198: if (!fas->interpolate) {
199: DMCreateInterpolation(next->dm, snes->dm, &fas->interpolate, &fas->rscale);
200: if (!fas->restrct) {
201: PetscObjectReference((PetscObject)fas->interpolate);
202: fas->restrct = fas->interpolate;
203: }
204: }
205: /* set the injection from the DM */
206: if (!fas->inject) {
207: DMCreateInjection(next->dm, snes->dm, &injscatter);
208: MatCreateScatter(((PetscObject)snes)->comm, injscatter, &fas->inject);
209: VecScatterDestroy(&injscatter);
210: }
211: }
212: }
213: /*pass the smoother, function, and jacobian up to the next level if it's not user set already */
214: if (fas->galerkin) {
215: if (next)
216: SNESSetFunction(next, PETSC_NULL, SNESFASGalerkinDefaultFunction, next);
217: if (fas->smoothd && fas->level != fas->levels - 1) SNESSetFunction(fas->smoothd, PETSC_NULL, SNESFASGalerkinDefaultFunction, snes);
218: if (fas->smoothu && fas->level != fas->levels - 1) SNESSetFunction(fas->smoothu, PETSC_NULL, SNESFASGalerkinDefaultFunction, snes);
219: }
221: /* sets the down (pre) smoother's default norm and sets it from options */
222: if(fas->smoothd){
223: if (fas->level == 0) {
224: SNESSetNormType(fas->smoothd, SNES_NORM_NONE);
225: } else {
226: SNESSetNormType(fas->smoothd, SNES_NORM_FINAL_ONLY);
227: }
228: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothd);
229: SNESSetFromOptions(fas->smoothd);
230: SNESGetSNESLineSearch(snes,&linesearch);
231: SNESGetSNESLineSearch(fas->smoothd,&slinesearch);
232: SNESLineSearchGetPreCheck(linesearch,&lsprefunc,&lsprectx);
233: SNESLineSearchGetPostCheck(linesearch,&lspostfunc,&lspostctx);
234: SNESLineSearchSetPreCheck(slinesearch,lsprefunc,lsprectx);
235: SNESLineSearchSetPostCheck(slinesearch,lspostfunc,lspostctx);
236: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
237: }
239: /* sets the up (post) smoother's default norm and sets it from options */
240: if(fas->smoothu){
241: if (fas->level != fas->levels - 1) {
242: SNESSetNormType(fas->smoothu, SNES_NORM_NONE);
243: } else {
244: SNESSetNormType(fas->smoothu, SNES_NORM_FINAL_ONLY);
245: }
246: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothu);
247: SNESSetFromOptions(fas->smoothu);
248: SNESGetSNESLineSearch(snes,&linesearch);
249: SNESGetSNESLineSearch(fas->smoothu,&slinesearch);
250: SNESLineSearchGetPreCheck(linesearch,&lsprefunc,&lsprectx);
251: SNESLineSearchGetPostCheck(linesearch,&lspostfunc,&lspostctx);
252: SNESLineSearchSetPreCheck(slinesearch,lsprefunc,lsprectx);
253: SNESLineSearchSetPostCheck(slinesearch,lspostfunc,lspostctx);
254: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
255: }
257: if (next) {
258: /* gotta set up the solution vector for this to work */
259: if (!next->vec_sol) {SNESFASCreateCoarseVec(snes,&next->vec_sol);}
260: if (!next->vec_rhs) {SNESFASCreateCoarseVec(snes,&next->vec_rhs);}
261: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next);
262: SNESGetSNESLineSearch(snes,&linesearch);
263: SNESGetSNESLineSearch(fas->next,&slinesearch);
264: SNESLineSearchGetPreCheck(linesearch,&lsprefunc,&lsprectx);
265: SNESLineSearchGetPostCheck(linesearch,&lspostfunc,&lspostctx);
266: SNESLineSearchSetPreCheck(slinesearch,lsprefunc,lsprectx);
267: SNESLineSearchSetPostCheck(slinesearch,lspostfunc,lspostctx);
268: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
269: SNESSetUp(next);
270: }
271: /* setup FAS work vectors */
272: if (fas->galerkin) {
273: VecDuplicate(snes->vec_sol, &fas->Xg);
274: VecDuplicate(snes->vec_sol, &fas->Fg);
275: }
276: return(0);
277: }
281: PetscErrorCode SNESSetFromOptions_FAS(SNES snes)
282: {
283: SNES_FAS *fas = (SNES_FAS *) snes->data;
284: PetscInt levels = 1;
285: PetscBool flg = PETSC_FALSE, upflg = PETSC_FALSE, downflg = PETSC_FALSE, monflg = PETSC_FALSE, galerkinflg = PETSC_FALSE;
287: char monfilename[PETSC_MAX_PATH_LEN];
288: SNESFASType fastype;
289: const char *optionsprefix;
290: SNESLineSearch linesearch;
291: PetscInt m, n_up, n_down;
292: SNES next;
293: PetscBool isFine;
296: SNESFASCycleIsFine(snes, &isFine);
297: PetscOptionsHead("SNESFAS Options-----------------------------------");
299: /* number of levels -- only process most options on the finest level */
300: if (isFine) {
301: PetscOptionsInt("-snes_fas_levels", "Number of Levels", "SNESFASSetLevels", levels, &levels, &flg);
302: if (!flg && snes->dm) {
303: DMGetRefineLevel(snes->dm,&levels);
304: levels++;
305: fas->usedmfornumberoflevels = PETSC_TRUE;
306: }
307: SNESFASSetLevels(snes, levels, PETSC_NULL);
308: fastype = fas->fastype;
309: PetscOptionsEnum("-snes_fas_type","FAS correction type","SNESFASSetType",SNESFASTypes,(PetscEnum)fastype,(PetscEnum*)&fastype,&flg);
310: if (flg) {
311: SNESFASSetType(snes, fastype);
312: }
314: SNESGetOptionsPrefix(snes, &optionsprefix);
315: PetscOptionsInt("-snes_fas_cycles","Number of cycles","SNESFASSetCycles",fas->n_cycles,&m,&flg);
316: if (flg) {
317: SNESFASSetCycles(snes, m);
318: }
320: PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin","SNESFASSetGalerkin",fas->galerkin,&galerkinflg,&flg);
321: if (flg) {
322: SNESFASSetGalerkin(snes, galerkinflg);
323: }
325: PetscOptionsInt("-snes_fas_smoothup","Number of post-smoothing steps","SNESFASSetNumberSmoothUp",fas->max_up_it,&n_up,&upflg);
327: PetscOptionsInt("-snes_fas_smoothdown","Number of pre-smoothing steps","SNESFASSetNumberSmoothDown",fas->max_down_it,&n_down,&downflg);
329: PetscOptionsString("-snes_fas_monitor","Monitor FAS progress","SNESFASSetMonitor","stdout",monfilename,PETSC_MAX_PATH_LEN,&monflg);
330: if (monflg) SNESFASSetMonitor(snes, PETSC_TRUE);
331: }
333: PetscOptionsTail();
334: /* setup from the determined types if there is no pointwise procedure or smoother defined */
335: if (upflg) {
336: SNESFASSetNumberSmoothUp(snes,n_up);
337: }
338: if (downflg) {
339: SNESFASSetNumberSmoothDown(snes,n_down);
340: }
342: /* set up the default line search for coarse grid corrections */
343: if (fas->fastype == SNES_FAS_ADDITIVE) {
344: if (!snes->linesearch) {
345: SNESGetSNESLineSearch(snes, &linesearch);
346: SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
347: }
348: }
350: SNESFASCycleGetCorrection(snes, &next);
351: /* recursive option setting for the smoothers */
352: if (next) {SNESSetFromOptions(next);}
353: return(0);
354: }
358: PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer)
359: {
360: SNES_FAS *fas = (SNES_FAS *) snes->data;
361: PetscBool isFine, iascii;
362: PetscInt i;
364: SNES smoothu, smoothd, levelsnes;
367: SNESFASCycleIsFine(snes, &isFine);
368: if (isFine) {
369: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
370: if (iascii) {
371: PetscViewerASCIIPrintf(viewer, "FAS: type is %s, levels=%D, cycles=%D\n", SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles);
372: if (fas->galerkin) {
373: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid function evaluation\n");
374: } else {
375: PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid function evaluation\n");
376: }
377: for (i=0; i<fas->levels; i++) {
378: SNESFASGetCycleSNES(snes, i, &levelsnes);
379: SNESFASCycleGetSmootherUp(levelsnes, &smoothu);
380: SNESFASCycleGetSmootherDown(levelsnes, &smoothd);
381: if (!i) {
382: PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level %D -------------------------------\n",i);
383: } else {
384: PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
385: }
386: PetscViewerASCIIPushTab(viewer);
387: SNESView(smoothd,viewer);
388: PetscViewerASCIIPopTab(viewer);
389: if (i && (smoothd == smoothu)) {
390: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
391: } else if (i) {
392: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
393: PetscViewerASCIIPushTab(viewer);
394: SNESView(smoothu,viewer);
395: PetscViewerASCIIPopTab(viewer);
396: }
397: }
398: } else {
399: SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for SNESFAS",((PetscObject)viewer)->type_name);
400: }
401: }
402: return(0);
403: }
407: /*
408: Defines the action of the downsmoother
409: */
410: PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
411: {
412: PetscErrorCode 0;
413: SNESConvergedReason reason;
414: Vec FPC;
415: SNES smoothd;
417: SNESFASCycleGetSmootherDown(snes, &smoothd);
418: SNESSetInitialFunction(smoothd, F);
419: SNESSetInitialFunctionNorm(smoothd, *fnorm);
420: SNESSolve(smoothd, B, X);
421: /* check convergence reason for the smoother */
422: SNESGetConvergedReason(smoothd,&reason);
423: if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN)) {
424: snes->reason = SNES_DIVERGED_INNER;
425: return(0);
426: }
427: SNESGetFunction(smoothd, &FPC, PETSC_NULL, PETSC_NULL);
428: VecCopy(FPC, F);
429: SNESGetFunctionNorm(smoothd, fnorm);
430: return(0);
431: }
436: /*
437: Defines the action of the upsmoother
438: */
439: PetscErrorCode SNESFASUpSmooth_Private (SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm) {
440: PetscErrorCode 0;
441: SNESConvergedReason reason;
442: Vec FPC;
443: SNES smoothu;
446: SNESFASCycleGetSmootherUp(snes, &smoothu);
447: SNESSolve(smoothu, B, X);
448: /* check convergence reason for the smoother */
449: SNESGetConvergedReason(smoothu,&reason);
450: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
451: snes->reason = SNES_DIVERGED_INNER;
452: return(0);
453: }
454: SNESGetFunction(smoothu, &FPC, PETSC_NULL, PETSC_NULL);
455: VecCopy(FPC, F);
456: SNESGetFunctionNorm(smoothu, fnorm);
457: return(0);
458: }
462: /*@
463: SNESFASCreateCoarseVec - create Vec corresponding to a state vector on one level coarser than current level
465: Collective
467: Input Arguments:
468: . snes - SNESFAS
470: Output Arguments:
471: . Xcoarse - vector on level one coarser than snes
473: Level: developer
475: .seealso: SNESFASSetRestriction(), SNESFASRestrict()
476: @*/
477: PetscErrorCode SNESFASCreateCoarseVec(SNES snes,Vec *Xcoarse)
478: {
480: SNES_FAS *fas = (SNES_FAS*)snes->data;
483: if (fas->rscale) {VecDuplicate(fas->rscale,Xcoarse);}
484: else if (fas->inject) {MatGetVecs(fas->inject,Xcoarse,PETSC_NULL);}
485: else SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or injection");
486: return(0);
487: }
491: /*@
492: SNESFASRestrict - restrict a Vec to the next coarser level
494: Collective
496: Input Arguments:
497: + fine - SNES from which to restrict
498: - Xfine - vector to restrict
500: Output Arguments:
501: . Xcoarse - result of restriction
503: Level: developer
505: .seealso: SNESFASSetRestriction(), SNESFASSetInjection()
506: @*/
507: PetscErrorCode SNESFASRestrict(SNES fine,Vec Xfine,Vec Xcoarse)
508: {
510: SNES_FAS *fas = (SNES_FAS*)fine->data;
516: if (fas->inject) {
517: MatRestrict(fas->inject,Xfine,Xcoarse);
518: } else {
519: MatRestrict(fas->restrct,Xfine,Xcoarse);
520: VecPointwiseMult(Xcoarse,fas->rscale,Xcoarse);
521: }
522: return(0);
523: }
527: /*
529: Performs the FAS coarse correction as:
531: fine problem: F(x) = 0
532: coarse problem: F^c(x) = b^c
534: b^c = F^c(I^c_fx^f - I^c_fF(x))
536: */
537: PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new) {
538: PetscErrorCode ierr;
539: Vec X_c, Xo_c, F_c, B_c;
540: SNESConvergedReason reason;
541: SNES next;
542: Mat restrct, interpolate;
544: SNESFASCycleGetCorrection(snes, &next);
545: if (next) {
546: SNESFASCycleGetRestriction(snes, &restrct);
547: SNESFASCycleGetInterpolation(snes, &interpolate);
549: X_c = next->vec_sol;
550: Xo_c = next->work[0];
551: F_c = next->vec_func;
552: B_c = next->vec_rhs;
554: SNESFASRestrict(snes,X,Xo_c);
555: /* restrict the defect */
556: MatRestrict(restrct, F, B_c);
557: /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */
558: SNESComputeFunction(next, Xo_c, F_c);
559: VecCopy(B_c, X_c);
560: VecCopy(F_c, B_c);
561: VecCopy(X_c, F_c);
562: /* set initial guess of the coarse problem to the projected fine solution */
563: VecCopy(Xo_c, X_c);
565: /* recurse to the next level */
566: SNESSetInitialFunction(next, F_c);
567: SNESSolve(next, B_c, X_c);
568: SNESGetConvergedReason(next,&reason);
569: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
570: snes->reason = SNES_DIVERGED_INNER;
571: return(0);
572: }
573: /* correct as x <- x + I(x^c - Rx)*/
574: VecAXPY(X_c, -1.0, Xo_c);
575: MatInterpolateAdd(interpolate, X_c, X, X_new);
576: }
577: return(0);
578: }
582: /*
584: The additive cycle looks like:
586: xhat = x
587: xhat = dS(x, b)
588: x = coarsecorrection(xhat, b_d)
589: x = x + nu*(xhat - x);
590: (optional) x = uS(x, b)
592: With the coarse RHS (defect correction) as below.
594: */
595: PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X) {
596: Vec F, B, Xhat;
597: Vec X_c, Xo_c, F_c, B_c;
598: PetscErrorCode ierr;
599: SNESConvergedReason reason;
600: PetscReal xnorm, fnorm, ynorm;
601: PetscBool lssuccess;
602: SNES next;
603: Mat restrct, interpolate;
605: SNESFASCycleGetCorrection(snes, &next);
606: F = snes->vec_func;
607: B = snes->vec_rhs;
608: Xhat = snes->work[1];
609: VecCopy(X, Xhat);
610: /* recurse first */
611: if (next) {
612: SNESFASCycleGetRestriction(snes, &restrct);
613: SNESFASCycleGetInterpolation(snes, &interpolate);
614: SNESComputeFunction(snes, Xhat, F);
615: VecNorm(F, NORM_2, &fnorm);
616: X_c = next->vec_sol;
617: Xo_c = next->work[0];
618: F_c = next->vec_func;
619: B_c = next->vec_rhs;
621: SNESFASRestrict(snes,Xhat,Xo_c);
622: /* restrict the defect */
623: MatRestrict(restrct, F, B_c);
625: /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */
626: SNESComputeFunction(next, Xo_c, F_c);
627: VecCopy(B_c, X_c);
628: VecCopy(F_c, B_c);
629: VecCopy(X_c, F_c);
630: /* set initial guess of the coarse problem to the projected fine solution */
631: VecCopy(Xo_c, X_c);
633: /* recurse */
634: SNESSetInitialFunction(next, F_c);
635: SNESSolve(next, B_c, X_c);
637: /* smooth on this level */
638: SNESFASDownSmooth_Private(snes, B, X, F, &fnorm);
640: SNESGetConvergedReason(next,&reason);
641: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
642: snes->reason = SNES_DIVERGED_INNER;
643: return(0);
644: }
646: /* correct as x <- x + I(x^c - Rx)*/
647: VecAYPX(X_c, -1.0, Xo_c);
648: MatInterpolate(interpolate, X_c, Xhat);
650: /* additive correction of the coarse direction*/
651: SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat);
652: SNESLineSearchGetSuccess(snes->linesearch, &lssuccess);
653: if (!lssuccess) {
654: if (++snes->numFailures >= snes->maxFailures) {
655: snes->reason = SNES_DIVERGED_LINE_SEARCH;
656: return(0);
657: }
658: }
659: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm);
660: } else {
661: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
662: }
663: return(0);
664: }
668: /*
670: Defines the FAS cycle as:
672: fine problem: F(x) = 0
673: coarse problem: F^c(x) = b^c
675: b^c = F^c(I^c_fx^f - I^c_fF(x))
677: correction:
679: x = x + I(x^c - Rx)
681: */
682: PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X) {
684: PetscErrorCode ierr;
685: Vec F,B;
686: SNES_FAS *fas = (SNES_FAS *)snes->data;
689: F = snes->vec_func;
690: B = snes->vec_rhs;
691: /* pre-smooth -- just update using the pre-smoother */
692: SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
693: if (fas->level != 0) {
694: SNESFASCoarseCorrection(snes, X, F, X);
695: SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
696: }
697: return(0);
698: }
703: PetscErrorCode SNESSolve_FAS(SNES snes)
704: {
706: PetscInt i, maxits;
707: Vec X, F;
708: PetscReal fnorm;
709: SNES_FAS *fas = (SNES_FAS *)snes->data,*ffas;
710: DM dm;
711: PetscBool isFine;
714: maxits = snes->max_its; /* maximum number of iterations */
715: snes->reason = SNES_CONVERGED_ITERATING;
716: X = snes->vec_sol;
717: F = snes->vec_func;
719: SNESFASCycleIsFine(snes, &isFine);
720: /*norm setup */
721: PetscObjectTakeAccess(snes);
722: snes->iter = 0;
723: snes->norm = 0.;
724: PetscObjectGrantAccess(snes);
725: if (!snes->vec_func_init_set) {
726: SNESComputeFunction(snes,X,F);
727: if (snes->domainerror) {
728: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
729: return(0);
730: }
731: } else {
732: snes->vec_func_init_set = PETSC_FALSE;
733: }
735: if (!snes->norm_init_set) {
736: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
737: if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
738: PetscObjectTakeAccess(snes);
739: } else {
740: fnorm = snes->norm_init;
741: snes->norm_init_set = PETSC_FALSE;
742: }
744: snes->norm = fnorm;
745: PetscObjectGrantAccess(snes);
746: SNESLogConvHistory(snes,fnorm,0);
747: SNESMonitor(snes,0,fnorm);
749: /* set parameter for default relative tolerance convergence test */
750: snes->ttol = fnorm*snes->rtol;
751: /* test convergence */
752: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
753: if (snes->reason) return(0);
756: if (isFine) {
757: /* propagate scale-dependent data up the hierarchy */
758: SNESGetDM(snes,&dm);
759: for (ffas=fas; ffas->next; ffas=(SNES_FAS*)ffas->next->data) {
760: DM dmcoarse;
761: SNESGetDM(ffas->next,&dmcoarse);
762: DMRestrict(dm,ffas->restrct,ffas->rscale,ffas->inject,dmcoarse);
763: dm = dmcoarse;
764: }
765: }
767: for (i = 0; i < maxits; i++) {
768: /* Call general purpose update function */
770: if (snes->ops->update) {
771: (*snes->ops->update)(snes, snes->iter);
772: }
773: if (fas->fastype == SNES_FAS_MULTIPLICATIVE) {
774: SNESFASCycle_Multiplicative(snes, X);
775: } else {
776: SNESFASCycle_Additive(snes, X);
777: }
779: /* check for FAS cycle divergence */
780: if (snes->reason != SNES_CONVERGED_ITERATING) {
781: return(0);
782: }
784: /* Monitor convergence */
785: PetscObjectTakeAccess(snes);
786: snes->iter = i+1;
787: PetscObjectGrantAccess(snes);
788: SNESLogConvHistory(snes,snes->norm,0);
789: SNESMonitor(snes,snes->iter,snes->norm);
790: /* Test for convergence */
791: if (isFine) {
792: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,snes->norm,&snes->reason,snes->cnvP);
793: if (snes->reason) break;
794: }
795: }
796: if (i == maxits) {
797: PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);
798: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
799: }
800: return(0);
801: }