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