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