Actual source code: fas.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: /* Defines the basic SNES object */
  2: #include <../src/snes/impls/fas/fasimpls.h>    /*I  "petscsnes.h"  I*/

  4: const char *const SNESFASTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","SNESFASType","SNES_FAS",0};

  6: extern PetscErrorCode SNESDestroy_FAS(SNES snes);
  7: extern PetscErrorCode SNESSetUp_FAS(SNES snes);
  8: extern PetscErrorCode SNESSetFromOptions_FAS(PetscOptions *PetscOptionsObject,SNES snes);
  9: extern PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer);
 10: extern PetscErrorCode SNESSolve_FAS(SNES snes);
 11: extern PetscErrorCode SNESReset_FAS(SNES snes);
 12: extern PetscErrorCode SNESFASGalerkinDefaultFunction(SNES, Vec, Vec, void*);
 13: extern PetscErrorCode SNESFASCycleCreateSmoother_Private(SNES, SNES*);

 15: /*MC

 17: SNESFAS - Full Approximation Scheme nonlinear multigrid solver.

 19:    The nonlinear problem is solved by correction using coarse versions
 20:    of the nonlinear problem.  This problem is perturbed so that a projected
 21:    solution of the fine problem elicits no correction from the coarse problem.

 23: Options Database:
 24: +   -snes_fas_levels -  The number of levels
 25: .   -snes_fas_cycles<1> -  The number of cycles -- 1 for V, 2 for W
 26: .   -snes_fas_type<additive,multiplicative,full,kaskade>  -  Additive or multiplicative cycle
 27: .   -snes_fas_galerkin<PETSC_FALSE> -  Form coarse problems by projection back upon the fine problem
 28: .   -snes_fas_smoothup<1> -  The number of iterations of the post-smoother
 29: .   -snes_fas_smoothdown<1> -  The number of iterations of the pre-smoother
 30: .   -snes_fas_monitor -  Monitor progress of all of the levels
 31: .   -snes_fas_full_downsweepsmooth<PETSC_FALSE> - call the downsmooth on the initial downsweep of full FAS
 32: .   -fas_levels_snes_ -  SNES options for all smoothers
 33: .   -fas_levels_cycle_snes_ -  SNES options for all cycles
 34: .   -fas_levels_i_snes_ -  SNES options for the smoothers on level i
 35: .   -fas_levels_i_cycle_snes_ - SNES options for the cycle on level i
 36: -   -fas_coarse_snes_ -  SNES options for the coarsest smoother

 38: Notes:
 39:    The organization of the FAS solver is slightly different from the organization of PCMG
 40:    As each level has smoother SNES instances(down and potentially up) and a cycle SNES instance.
 41:    The cycle SNES instance may be used for monitoring convergence on a particular level.

 43: Level: beginner

 45: .seealso: PCMG, SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types)
 46: M*/

 50: PETSC_EXTERN PetscErrorCode SNESCreate_FAS(SNES snes)
 51: {
 52:   SNES_FAS       *fas;

 56:   snes->ops->destroy        = SNESDestroy_FAS;
 57:   snes->ops->setup          = SNESSetUp_FAS;
 58:   snes->ops->setfromoptions = SNESSetFromOptions_FAS;
 59:   snes->ops->view           = SNESView_FAS;
 60:   snes->ops->solve          = SNESSolve_FAS;
 61:   snes->ops->reset          = SNESReset_FAS;

 63:   snes->usesksp = PETSC_FALSE;
 64:   snes->usespc  = PETSC_FALSE;

 66:   if (!snes->tolerancesset) {
 67:     snes->max_funcs = 30000;
 68:     snes->max_its   = 10000;
 69:   }

 71:   PetscNewLog(snes,&fas);

 73:   snes->data                  = (void*) fas;
 74:   fas->level                  = 0;
 75:   fas->levels                 = 1;
 76:   fas->n_cycles               = 1;
 77:   fas->max_up_it              = 1;
 78:   fas->max_down_it            = 1;
 79:   fas->smoothu                = NULL;
 80:   fas->smoothd                = NULL;
 81:   fas->next                   = NULL;
 82:   fas->previous               = NULL;
 83:   fas->fine                   = snes;
 84:   fas->interpolate            = NULL;
 85:   fas->restrct                = NULL;
 86:   fas->inject                 = NULL;
 87:   fas->monitor                = NULL;
 88:   fas->usedmfornumberoflevels = PETSC_FALSE;
 89:   fas->fastype                = SNES_FAS_MULTIPLICATIVE;
 90:   fas->full_downsweep         = PETSC_FALSE;

 92:   fas->eventsmoothsetup    = 0;
 93:   fas->eventsmoothsolve    = 0;
 94:   fas->eventresidual       = 0;
 95:   fas->eventinterprestrict = 0;
 96:   return(0);
 97: }

101: PetscErrorCode SNESReset_FAS(SNES snes)
102: {
103:   PetscErrorCode 0;
104:   SNES_FAS       * fas = (SNES_FAS*)snes->data;

107:   SNESDestroy(&fas->smoothu);
108:   SNESDestroy(&fas->smoothd);
109:   MatDestroy(&fas->inject);
110:   MatDestroy(&fas->interpolate);
111:   MatDestroy(&fas->restrct);
112:   VecDestroy(&fas->rscale);
113:   if (fas->galerkin) {
114:     VecDestroy(&fas->Xg);
115:     VecDestroy(&fas->Fg);
116:   }
117:   if (fas->next) {SNESReset(fas->next);}
118:   return(0);
119: }

123: PetscErrorCode SNESDestroy_FAS(SNES snes)
124: {
125:   SNES_FAS       * fas = (SNES_FAS*)snes->data;
126:   PetscErrorCode 0;

129:   /* recursively resets and then destroys */
130:   SNESReset(snes);
131:   if (fas->next) {
132:     SNESDestroy(&fas->next);
133:   }
134:   PetscFree(fas);
135:   return(0);
136: }

140: PetscErrorCode SNESSetUp_FAS(SNES snes)
141: {
142:   SNES_FAS       *fas = (SNES_FAS*) snes->data;
144:   PetscInt       dm_levels;
145:   Vec            vec_sol, vec_func, vec_sol_update, vec_rhs; /* preserve these if they're set through the reset */
146:   SNES           next;
147:   PetscBool      isFine;
148:   SNESLineSearch linesearch;
149:   SNESLineSearch slinesearch;
150:   void           *lsprectx,*lspostctx;
151:   PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*);
152:   PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);

155:   SNESFASCycleIsFine(snes, &isFine);
156:   if (fas->usedmfornumberoflevels && isFine) {
157:     DMGetRefineLevel(snes->dm,&dm_levels);
158:     dm_levels++;
159:     if (dm_levels > fas->levels) {
160:       /* we don't want the solution and func vectors to be destroyed in the SNESReset when it's called in SNESFASSetLevels_FAS*/
161:       vec_sol              = snes->vec_sol;
162:       vec_func             = snes->vec_func;
163:       vec_sol_update       = snes->vec_sol_update;
164:       vec_rhs              = snes->vec_rhs;
165:       snes->vec_sol        = NULL;
166:       snes->vec_func       = NULL;
167:       snes->vec_sol_update = NULL;
168:       snes->vec_rhs        = NULL;

170:       /* reset the number of levels */
171:       SNESFASSetLevels(snes,dm_levels,NULL);
172:       SNESSetFromOptions(snes);

174:       snes->vec_sol        = vec_sol;
175:       snes->vec_func       = vec_func;
176:       snes->vec_rhs        = vec_rhs;
177:       snes->vec_sol_update = vec_sol_update;
178:     }
179:   }
180:   SNESFASCycleGetCorrection(snes, &next);
181:   if (!isFine) snes->gridsequence = 0; /* no grid sequencing inside the multigrid hierarchy! */

183:   SNESSetWorkVecs(snes, 2); /* work vectors used for intergrid transfers */

185:   /* set up the smoothers if they haven't already been set up */
186:   if (!fas->smoothd) {
187:     SNESFASCycleCreateSmoother_Private(snes, &fas->smoothd);
188:   }

190:   if (snes->dm) {
191:     /* set the smoother DMs properly */
192:     if (fas->smoothu) SNESSetDM(fas->smoothu, snes->dm);
193:     SNESSetDM(fas->smoothd, snes->dm);
194:     /* construct EVERYTHING from the DM -- including the progressive set of smoothers */
195:     if (next) {
196:       /* for now -- assume the DM and the evaluation functions have been set externally */
197:       if (!next->dm) {
198:         DMCoarsen(snes->dm, PetscObjectComm((PetscObject)next), &next->dm);
199:         SNESSetDM(next, next->dm);
200:       }
201:       /* set the interpolation and restriction from the DM */
202:       if (!fas->interpolate) {
203:         DMCreateInterpolation(next->dm, snes->dm, &fas->interpolate, &fas->rscale);
204:         if (!fas->restrct) {
205:           PetscObjectReference((PetscObject)fas->interpolate);
206:           fas->restrct = fas->interpolate;
207:         }
208:       }
209:       /* set the injection from the DM */
210:       if (!fas->inject) {
211:         DMCreateInjection(next->dm, snes->dm, &fas->inject);
212:       }
213:     }
214:   }
215:   /*pass the smoother, function, and jacobian up to the next level if it's not user set already */
216:   if (fas->galerkin) {
217:     if (next) {
218:       SNESSetFunction(next, NULL, SNESFASGalerkinDefaultFunction, next);
219:     }
220:     if (fas->smoothd && fas->level != fas->levels - 1) {
221:       SNESSetFunction(fas->smoothd, NULL, SNESFASGalerkinDefaultFunction, snes);
222:     }
223:     if (fas->smoothu && fas->level != fas->levels - 1) {
224:       SNESSetFunction(fas->smoothu, NULL, SNESFASGalerkinDefaultFunction, snes);
225:     }
226:   }

228:   /* sets the down (pre) smoother's default norm and sets it from options */
229:   if (fas->smoothd) {
230:     if (fas->level == 0 && fas->levels != 1) {
231:       SNESSetNormSchedule(fas->smoothd, SNES_NORM_NONE);
232:     } else {
233:       SNESSetNormSchedule(fas->smoothd, SNES_NORM_FINAL_ONLY);
234:     }
235:     PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothd);
236:     SNESSetFromOptions(fas->smoothd);
237:     SNESGetLineSearch(snes,&linesearch);
238:     SNESGetLineSearch(fas->smoothd,&slinesearch);
239:     SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
240:     SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
241:     SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
242:     SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
243:     PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);

245:     fas->smoothd->vec_sol        = snes->vec_sol;
246:     PetscObjectReference((PetscObject)snes->vec_sol);
247:     fas->smoothd->vec_sol_update = snes->vec_sol_update;
248:     PetscObjectReference((PetscObject)snes->vec_sol_update);
249:     fas->smoothd->vec_func       = snes->vec_func;
250:     PetscObjectReference((PetscObject)snes->vec_func);

252:     if (fas->eventsmoothsetup) {PetscLogEventBegin(fas->eventsmoothsetup,0,0,0,0);}
253:     SNESSetUp(fas->smoothd);
254:     if (fas->eventsmoothsetup) {PetscLogEventEnd(fas->eventsmoothsetup,0,0,0,0);}
255:   }

257:   /* sets the up (post) smoother's default norm and sets it from options */
258:   if (fas->smoothu) {
259:     if (fas->level != fas->levels - 1) {
260:       SNESSetNormSchedule(fas->smoothu, SNES_NORM_NONE);
261:     } else {
262:       SNESSetNormSchedule(fas->smoothu, SNES_NORM_FINAL_ONLY);
263:     }
264:     PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothu);
265:     SNESSetFromOptions(fas->smoothu);
266:     SNESGetLineSearch(snes,&linesearch);
267:     SNESGetLineSearch(fas->smoothu,&slinesearch);
268:     SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
269:     SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
270:     SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
271:     SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
272:     PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);

274:     fas->smoothu->vec_sol        = snes->vec_sol;
275:     PetscObjectReference((PetscObject)snes->vec_sol);
276:     fas->smoothu->vec_sol_update = snes->vec_sol_update;
277:     PetscObjectReference((PetscObject)snes->vec_sol_update);
278:     fas->smoothu->vec_func       = snes->vec_func;
279:     PetscObjectReference((PetscObject)snes->vec_func);

281:     if (fas->eventsmoothsetup) {PetscLogEventBegin(fas->eventsmoothsetup,0,0,0,0);}
282:     SNESSetUp(fas->smoothu);
283:     if (fas->eventsmoothsetup) {PetscLogEventEnd(fas->eventsmoothsetup,0,0,0,0);}

285:   }

287:   if (next) {
288:     /* gotta set up the solution vector for this to work */
289:     if (!next->vec_sol) {SNESFASCreateCoarseVec(snes,&next->vec_sol);}
290:     if (!next->vec_rhs) {SNESFASCreateCoarseVec(snes,&next->vec_rhs);}
291:     PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next);
292:     SNESGetLineSearch(snes,&linesearch);
293:     SNESGetLineSearch(fas->next,&slinesearch);
294:     SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
295:     SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
296:     SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
297:     SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
298:     PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
299:     SNESSetUp(next);
300:   }
301:   /* setup FAS work vectors */
302:   if (fas->galerkin) {
303:     VecDuplicate(snes->vec_sol, &fas->Xg);
304:     VecDuplicate(snes->vec_sol, &fas->Fg);
305:   }
306:   return(0);
307: }

311: PetscErrorCode SNESSetFromOptions_FAS(PetscOptions *PetscOptionsObject,SNES snes)
312: {
313:   SNES_FAS       *fas   = (SNES_FAS*) snes->data;
314:   PetscInt       levels = 1;
315:   PetscBool      flg    = PETSC_FALSE, upflg = PETSC_FALSE, downflg = PETSC_FALSE, monflg = PETSC_FALSE, galerkinflg = PETSC_FALSE,continuationflg = PETSC_FALSE;
317:   char           monfilename[PETSC_MAX_PATH_LEN];
318:   SNESFASType    fastype;
319:   const char     *optionsprefix;
320:   SNESLineSearch linesearch;
321:   PetscInt       m, n_up, n_down;
322:   SNES           next;
323:   PetscBool      isFine;

326:   SNESFASCycleIsFine(snes, &isFine);
327:   PetscOptionsHead(PetscOptionsObject,"SNESFAS Options-----------------------------------");

329:   /* number of levels -- only process most options on the finest level */
330:   if (isFine) {
331:     PetscOptionsInt("-snes_fas_levels", "Number of Levels", "SNESFASSetLevels", levels, &levels, &flg);
332:     if (!flg && snes->dm) {
333:       DMGetRefineLevel(snes->dm,&levels);
334:       levels++;
335:       fas->usedmfornumberoflevels = PETSC_TRUE;
336:     }
337:     SNESFASSetLevels(snes, levels, NULL);
338:     fastype = fas->fastype;
339:     PetscOptionsEnum("-snes_fas_type","FAS correction type","SNESFASSetType",SNESFASTypes,(PetscEnum)fastype,(PetscEnum*)&fastype,&flg);
340:     if (flg) {
341:       SNESFASSetType(snes, fastype);
342:     }

344:     SNESGetOptionsPrefix(snes, &optionsprefix);
345:     PetscOptionsInt("-snes_fas_cycles","Number of cycles","SNESFASSetCycles",fas->n_cycles,&m,&flg);
346:     if (flg) {
347:       SNESFASSetCycles(snes, m);
348:     }
349:     PetscOptionsBool("-snes_fas_continuation","Corrected grid-sequence continuation","SNESFASSetContinuation",fas->continuation,&continuationflg,&flg);
350:     if (flg) {
351:       SNESFASSetContinuation(snes,continuationflg);
352:     }

354:     PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin","SNESFASSetGalerkin",fas->galerkin,&galerkinflg,&flg);
355:     if (flg) {
356:       SNESFASSetGalerkin(snes, galerkinflg);
357:     }

359:     if (fas->fastype == SNES_FAS_FULL) {
360:       PetscOptionsBool("-snes_fas_full_downsweep","Smooth on the initial upsweep for full FAS cycles","SNESFASFullSetDownSweep",fas->full_downsweep,&fas->full_downsweep,&flg);
361:       if (flg) {SNESFASFullSetDownSweep(snes,fas->full_downsweep);}
362:     }

364:     PetscOptionsInt("-snes_fas_smoothup","Number of post-smoothing steps","SNESFASSetNumberSmoothUp",fas->max_up_it,&n_up,&upflg);

366:     PetscOptionsInt("-snes_fas_smoothdown","Number of pre-smoothing steps","SNESFASSetNumberSmoothDown",fas->max_down_it,&n_down,&downflg);

368:     PetscOptionsString("-snes_fas_monitor","Monitor FAS progress","SNESFASSetMonitor","stdout",monfilename,PETSC_MAX_PATH_LEN,&monflg);
369:     if (monflg) SNESFASSetMonitor(snes, PETSC_TRUE);

371:     flg    = PETSC_FALSE;
372:     monflg = PETSC_TRUE;
373:     PetscOptionsBool("-snes_fas_log","Log times for each FAS level","SNESFASSetLog",monflg,&monflg,&flg);
374:     if (flg) {SNESFASSetLog(snes,monflg);}
375:   }

377:   PetscOptionsTail();
378:   /* setup from the determined types if there is no pointwise procedure or smoother defined */
379:   if (upflg) {
380:     SNESFASSetNumberSmoothUp(snes,n_up);
381:   }
382:   if (downflg) {
383:     SNESFASSetNumberSmoothDown(snes,n_down);
384:   }

386:   /* set up the default line search for coarse grid corrections */
387:   if (fas->fastype == SNES_FAS_ADDITIVE) {
388:     if (!snes->linesearch) {
389:       SNESGetLineSearch(snes, &linesearch);
390:       SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
391:     }
392:   }

394:   SNESFASCycleGetCorrection(snes, &next);
395:   /* recursive option setting for the smoothers */
396:   if (next) {SNESSetFromOptions(next);}
397:   return(0);
398: }

400: #include <petscdraw.h>
403: PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer)
404: {
405:   SNES_FAS       *fas = (SNES_FAS*) snes->data;
406:   PetscBool      isFine,iascii,isdraw;
407:   PetscInt       i;
409:   SNES           smoothu, smoothd, levelsnes;

412:   SNESFASCycleIsFine(snes, &isFine);
413:   if (isFine) {
414:     PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
415:     PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
416:     if (iascii) {
417:       PetscViewerASCIIPrintf(viewer, "FAS: type is %s, levels=%D, cycles=%D\n",  SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles);
418:       if (fas->galerkin) {
419:         PetscViewerASCIIPrintf(viewer,"    Using Galerkin computed coarse grid function evaluation\n");
420:       } else {
421:         PetscViewerASCIIPrintf(viewer,"    Not using Galerkin computed coarse grid function evaluation\n");
422:       }
423:       for (i=0; i<fas->levels; i++) {
424:         SNESFASGetCycleSNES(snes, i, &levelsnes);
425:         SNESFASCycleGetSmootherUp(levelsnes, &smoothu);
426:         SNESFASCycleGetSmootherDown(levelsnes, &smoothd);
427:         if (!i) {
428:           PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level %D -------------------------------\n",i);
429:         } else {
430:           PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
431:         }
432:         PetscViewerASCIIPushTab(viewer);
433:         if (smoothd) {
434:           SNESView(smoothd,viewer);
435:         } else {
436:           PetscViewerASCIIPrintf(viewer,"Not yet available\n");
437:         }
438:         PetscViewerASCIIPopTab(viewer);
439:         if (i && (smoothd == smoothu)) {
440:           PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
441:         } else if (i) {
442:           PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
443:           PetscViewerASCIIPushTab(viewer);
444:           if (smoothu) {
445:             SNESView(smoothu,viewer);
446:           } else {
447:             PetscViewerASCIIPrintf(viewer,"Not yet available\n");
448:           }
449:           PetscViewerASCIIPopTab(viewer);
450:         }
451:       }
452:     } else if (isdraw) {
453:       PetscDraw draw;
454:       PetscReal x,w,y,bottom,th,wth;
455:       SNES_FAS  *curfas = fas;
456:       PetscViewerDrawGetDraw(viewer,0,&draw);
457:       PetscDrawGetCurrentPoint(draw,&x,&y);
458:       PetscDrawStringGetSize(draw,&wth,&th);
459:       bottom = y - th;
460:       while (curfas) {
461:         if (!curfas->smoothu) {
462:           PetscDrawPushCurrentPoint(draw,x,bottom);
463:           if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
464:           PetscDrawPopCurrentPoint(draw);
465:         } else {
466:           w    = 0.5*PetscMin(1.0-x,x);
467:           PetscDrawPushCurrentPoint(draw,x-w,bottom);
468:           if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
469:           PetscDrawPopCurrentPoint(draw);
470:           PetscDrawPushCurrentPoint(draw,x+w,bottom);
471:           if (curfas->smoothu) SNESView(curfas->smoothu,viewer);
472:           PetscDrawPopCurrentPoint(draw);
473:         }
474:         /* this is totally bogus but we have no way of knowing how low the previous one was draw to */
475:         bottom -= 5*th;
476:         if (curfas->next) curfas = (SNES_FAS*)curfas->next->data;
477:         else curfas = NULL;
478:       }
479:     }
480:   }
481:   return(0);
482: }

486: /*
487: Defines the action of the downsmoother
488:  */
489: PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
490: {
491:   PetscErrorCode      0;
492:   SNESConvergedReason reason;
493:   Vec                 FPC;
494:   SNES                smoothd;
495:   SNES_FAS            *fas = (SNES_FAS*) snes->data;

498:   SNESFASCycleGetSmootherDown(snes, &smoothd);
499:   SNESSetInitialFunction(smoothd, F);
500:   if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,0,0,0,0);}
501:   SNESSolve(smoothd, B, X);
502:   if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,0,0,0,0);}
503:   /* check convergence reason for the smoother */
504:   SNESGetConvergedReason(smoothd,&reason);
505:   if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
506:     snes->reason = SNES_DIVERGED_INNER;
507:     return(0);
508:   }
509:   SNESGetFunction(smoothd, &FPC, NULL, NULL);
510:   VecCopy(FPC, F);
511:   if (fnorm) {VecNorm(F,NORM_2,fnorm);}
512:   return(0);
513: }


518: /*
519: Defines the action of the upsmoother
520:  */
521: PetscErrorCode SNESFASUpSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
522: {
523:   PetscErrorCode      0;
524:   SNESConvergedReason reason;
525:   Vec                 FPC;
526:   SNES                smoothu;
527:   SNES_FAS            *fas = (SNES_FAS*) snes->data;

530:   SNESFASCycleGetSmootherUp(snes, &smoothu);
531:   if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,0,0,0,0);}
532:   SNESSolve(smoothu, B, X);
533:   if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,0,0,0,0);}
534:   /* check convergence reason for the smoother */
535:   SNESGetConvergedReason(smoothu,&reason);
536:   if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
537:     snes->reason = SNES_DIVERGED_INNER;
538:     return(0);
539:   }
540:   SNESGetFunction(smoothu, &FPC, NULL, NULL);
541:   VecCopy(FPC, F);
542:   if (fnorm) {VecNorm(F,NORM_2,fnorm);}
543:   return(0);
544: }

548: /*@
549:    SNESFASCreateCoarseVec - create Vec corresponding to a state vector on one level coarser than current level

551:    Collective

553:    Input Arguments:
554: .  snes - SNESFAS

556:    Output Arguments:
557: .  Xcoarse - vector on level one coarser than snes

559:    Level: developer

561: .seealso: SNESFASSetRestriction(), SNESFASRestrict()
562: @*/
563: PetscErrorCode SNESFASCreateCoarseVec(SNES snes,Vec *Xcoarse)
564: {
566:   SNES_FAS       *fas = (SNES_FAS*)snes->data;

569:   if (fas->rscale) {
570:     VecDuplicate(fas->rscale,Xcoarse);
571:   } else if (fas->inject) {
572:     MatCreateVecs(fas->inject,Xcoarse,NULL);
573:   } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or injection");
574:   return(0);
575: }

579: /*@
580:    SNESFASRestrict - restrict a Vec to the next coarser level

582:    Collective

584:    Input Arguments:
585: +  fine - SNES from which to restrict
586: -  Xfine - vector to restrict

588:    Output Arguments:
589: .  Xcoarse - result of restriction

591:    Level: developer

593: .seealso: SNESFASSetRestriction(), SNESFASSetInjection()
594: @*/
595: PetscErrorCode SNESFASRestrict(SNES fine,Vec Xfine,Vec Xcoarse)
596: {
598:   SNES_FAS       *fas = (SNES_FAS*)fine->data;

604:   if (fas->inject) {
605:     MatRestrict(fas->inject,Xfine,Xcoarse);
606:   } else {
607:     MatRestrict(fas->restrct,Xfine,Xcoarse);
608:     VecPointwiseMult(Xcoarse,fas->rscale,Xcoarse);
609:   }
610:   return(0);
611: }

615: /*

617: Performs the FAS coarse correction as:

619: fine problem:   F(x) = b
620: coarse problem: F^c(x^c) = b^c

622: b^c = F^c(Rx) - R(F(x) - b)

624:  */
625: PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new)
626: {
627:   PetscErrorCode      ierr;
628:   Vec                 X_c, Xo_c, F_c, B_c;
629:   SNESConvergedReason reason;
630:   SNES                next;
631:   Mat                 restrct, interpolate;
632:   SNES_FAS            *fasc;

635:   SNESFASCycleGetCorrection(snes, &next);
636:   if (next) {
637:     fasc = (SNES_FAS*)next->data;

639:     SNESFASCycleGetRestriction(snes, &restrct);
640:     SNESFASCycleGetInterpolation(snes, &interpolate);

642:     X_c  = next->vec_sol;
643:     Xo_c = next->work[0];
644:     F_c  = next->vec_func;
645:     B_c  = next->vec_rhs;

647:     if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,0,0,0,0);}
648:     SNESFASRestrict(snes,X,Xo_c);
649:     /* restrict the defect: R(F(x) - b) */
650:     MatRestrict(restrct, F, B_c);
651:     if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,0,0,0,0);}

653:     if (fasc->eventresidual) {PetscLogEventBegin(fasc->eventresidual,0,0,0,0);}
654:     /* F_c = F^c(Rx) - R(F(x) - b) since the second term was sitting in next->vec_rhs */
655:     SNESComputeFunction(next, Xo_c, F_c);
656:     if (fasc->eventresidual) {PetscLogEventEnd(fasc->eventresidual,0,0,0,0);}

658:     /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */
659:     VecCopy(B_c, X_c);
660:     VecCopy(F_c, B_c);
661:     VecCopy(X_c, F_c);
662:     /* set initial guess of the coarse problem to the projected fine solution */
663:     VecCopy(Xo_c, X_c);

665:     /* recurse to the next level */
666:     SNESSetInitialFunction(next, F_c);
667:     SNESSolve(next, B_c, X_c);
668:     SNESGetConvergedReason(next,&reason);
669:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
670:       snes->reason = SNES_DIVERGED_INNER;
671:       return(0);
672:     }
673:     /* correct as x <- x + I(x^c - Rx)*/
674:     VecAXPY(X_c, -1.0, Xo_c);

676:     if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,0,0,0,0);}
677:     MatInterpolateAdd(interpolate, X_c, X, X_new);
678:     if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,0,0,0,0);}
679:   }
680:   return(0);
681: }

685: /*

687: The additive cycle looks like:

689: xhat = x
690: xhat = dS(x, b)
691: x = coarsecorrection(xhat, b_d)
692: x = x + nu*(xhat - x);
693: (optional) x = uS(x, b)

695: With the coarse RHS (defect correction) as below.

697:  */
698: PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X)
699: {
700:   Vec                  F, B, Xhat;
701:   Vec                  X_c, Xo_c, F_c, B_c;
702:   PetscErrorCode       ierr;
703:   SNESConvergedReason  reason;
704:   PetscReal            xnorm, fnorm, ynorm;
705:   SNESLineSearchReason lsresult;
706:   SNES                 next;
707:   Mat                  restrct, interpolate;
708:   SNES_FAS             *fas = (SNES_FAS*)snes->data,*fasc;

711:   SNESFASCycleGetCorrection(snes, &next);
712:   F    = snes->vec_func;
713:   B    = snes->vec_rhs;
714:   Xhat = snes->work[1];
715:   VecCopy(X, Xhat);
716:   /* recurse first */
717:   if (next) {
718:     fasc = (SNES_FAS*)next->data;
719:     SNESFASCycleGetRestriction(snes, &restrct);
720:     SNESFASCycleGetInterpolation(snes, &interpolate);
721:     if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,0,0,0,0);}
722:     SNESComputeFunction(snes, Xhat, F);
723:     if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,0,0,0,0);}
724:     VecNorm(F, NORM_2, &fnorm);
725:     X_c  = next->vec_sol;
726:     Xo_c = next->work[0];
727:     F_c  = next->vec_func;
728:     B_c  = next->vec_rhs;

730:     SNESFASRestrict(snes,Xhat,Xo_c);
731:     /* restrict the defect */
732:     MatRestrict(restrct, F, B_c);

734:     /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */
735:     if (fasc->eventresidual) {PetscLogEventBegin(fasc->eventresidual,0,0,0,0);}
736:     SNESComputeFunction(next, Xo_c, F_c);
737:     if (fasc->eventresidual) {PetscLogEventEnd(fasc->eventresidual,0,0,0,0);}
738:     VecCopy(B_c, X_c);
739:     VecCopy(F_c, B_c);
740:     VecCopy(X_c, F_c);
741:     /* set initial guess of the coarse problem to the projected fine solution */
742:     VecCopy(Xo_c, X_c);

744:     /* recurse */
745:     SNESSetInitialFunction(next, F_c);
746:     SNESSolve(next, B_c, X_c);

748:     /* smooth on this level */
749:     SNESFASDownSmooth_Private(snes, B, X, F, &fnorm);

751:     SNESGetConvergedReason(next,&reason);
752:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
753:       snes->reason = SNES_DIVERGED_INNER;
754:       return(0);
755:     }

757:     /* correct as x <- x + I(x^c - Rx)*/
758:     VecAYPX(X_c, -1.0, Xo_c);
759:     MatInterpolate(interpolate, X_c, Xhat);

761:     /* additive correction of the coarse direction*/
762:     SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat);
763:     SNESLineSearchGetReason(snes->linesearch, &lsresult);
764:     SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm);
765:     if (lsresult) {
766:       if (++snes->numFailures >= snes->maxFailures) {
767:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
768:         return(0);
769:       }
770:     }
771:   } else {
772:     SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
773:   }
774:   return(0);
775: }

779: /*

781: Defines the FAS cycle as:

783: fine problem: F(x) = b
784: coarse problem: F^c(x) = b^c

786: b^c = F^c(Rx) - R(F(x) - b)

788: correction:

790: x = x + I(x^c - Rx)

792:  */
793: PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X)
794: {

797:   Vec            F,B;
798:   SNES           next;

801:   F = snes->vec_func;
802:   B = snes->vec_rhs;
803:   /* pre-smooth -- just update using the pre-smoother */
804:   SNESFASCycleGetCorrection(snes,&next);
805:   SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
806:   if (next) {
807:     SNESFASCoarseCorrection(snes, X, F, X);
808:     SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
809:   }
810:   return(0);
811: }

815: PetscErrorCode SNESFASCycleSetupPhase_Full(SNES snes)
816: {
817:   SNES           next;
818:   SNES_FAS       *fas = (SNES_FAS*)snes->data;
819:   PetscBool      isFine;

823:   /* pre-smooth -- just update using the pre-smoother */
824:   SNESFASCycleIsFine(snes,&isFine);
825:   SNESFASCycleGetCorrection(snes,&next);
826:   fas->full_stage = 0;
827:   if (next) {SNESFASCycleSetupPhase_Full(next);}
828:   return(0);
829: }

833: PetscErrorCode SNESFASCycle_Full(SNES snes, Vec X)
834: {
836:   Vec            F,B;
837:   SNES_FAS       *fas = (SNES_FAS*)snes->data;
838:   PetscBool      isFine;
839:   SNES           next;

842:   F = snes->vec_func;
843:   B = snes->vec_rhs;
844:   SNESFASCycleIsFine(snes,&isFine);
845:   SNESFASCycleGetCorrection(snes,&next);

847:   if (isFine) {
848:     SNESFASCycleSetupPhase_Full(snes);
849:   }

851:   if (fas->full_stage == 0) {
852:     /* downsweep */
853:     if (next) {
854:       if (fas->level != 1) next->max_its += 1;
855:       if (fas->full_downsweep||isFine) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
856:       SNESFASCoarseCorrection(snes,X,F,X);
857:       SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
858:       if (fas->level != 1) next->max_its -= 1;
859:     } else {
860:       SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
861:     }
862:     fas->full_stage = 1;
863:   } else if (fas->full_stage == 1) {
864:     if (snes->iter == 0) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
865:     if (next) {
866:       SNESFASCoarseCorrection(snes,X,F,X);
867:       SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
868:     }
869:   }
870:   /* final v-cycle */
871:   if (isFine) {
872:     if (next) {
873:       SNESFASCoarseCorrection(snes,X,F,X);
874:       SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
875:     }
876:   }
877:   return(0);
878: }

882: PetscErrorCode SNESFASCycle_Kaskade(SNES snes, Vec X)
883: {
885:   Vec            F,B;
886:   SNES           next;

889:   F = snes->vec_func;
890:   B = snes->vec_rhs;
891:   SNESFASCycleGetCorrection(snes,&next);
892:   if (next) {
893:     SNESFASCoarseCorrection(snes,X,F,X);
894:     SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
895:   } else {
896:     SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
897:   }
898:   return(0);
899: }

901: PetscBool SNEScite = PETSC_FALSE;
902: const char SNESCitation[] = "@techreport{pbmkbsxt2012,\n"
903:                             "  title = {Composing Scalable Nonlinear Algebraic Solvers},\n"
904:                             "  author = {Peter Brune and Mathew Knepley and Barry Smith and Xuemin Tu},\n"
905:                             "  year = 2013,\n"
906:                             "  type = Preprint,\n"
907:                             "  number = {ANL/MCS-P2010-0112},\n"
908:                             "  institution = {Argonne National Laboratory}\n}\n";

912: PetscErrorCode SNESSolve_FAS(SNES snes)
913: {
915:   PetscInt       i, maxits;
916:   Vec            X, F;
917:   PetscReal      fnorm;
918:   SNES_FAS       *fas = (SNES_FAS*)snes->data,*ffas;
919:   DM             dm;
920:   PetscBool      isFine;


924:   if (snes->xl || snes->xu || snes->ops->computevariablebounds) {
925:     SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
926:   }

928:   PetscCitationsRegister(SNESCitation,&SNEScite);
929:   maxits       = snes->max_its;      /* maximum number of iterations */
930:   snes->reason = SNES_CONVERGED_ITERATING;
931:   X            = snes->vec_sol;
932:   F            = snes->vec_func;

934:   SNESFASCycleIsFine(snes, &isFine);
935:   /*norm setup */
936:   PetscObjectSAWsTakeAccess((PetscObject)snes);
937:   snes->iter = 0;
938:   snes->norm = 0.;
939:   PetscObjectSAWsGrantAccess((PetscObject)snes);
940:   if (!snes->vec_func_init_set) {
941:     if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,0,0,0,0);}
942:     SNESComputeFunction(snes,X,F);
943:     if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,0,0,0,0);}
944:   } else snes->vec_func_init_set = PETSC_FALSE;

946:   VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
947:   SNESCheckFunctionNorm(snes,fnorm);
948:   PetscObjectSAWsTakeAccess((PetscObject)snes);
949:   snes->norm = fnorm;
950:   PetscObjectSAWsGrantAccess((PetscObject)snes);
951:   SNESLogConvergenceHistory(snes,fnorm,0);
952:   SNESMonitor(snes,0,fnorm);

954:   /* test convergence */
955:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
956:   if (snes->reason) return(0);


959:   if (isFine) {
960:     /* propagate scale-dependent data up the hierarchy */
961:     SNESGetDM(snes,&dm);
962:     for (ffas=fas; ffas->next; ffas=(SNES_FAS*)ffas->next->data) {
963:       DM dmcoarse;
964:       SNESGetDM(ffas->next,&dmcoarse);
965:       DMRestrict(dm,ffas->restrct,ffas->rscale,ffas->inject,dmcoarse);
966:       dm   = dmcoarse;
967:     }
968:   }

970:   for (i = 0; i < maxits; i++) {
971:     /* Call general purpose update function */

973:     if (snes->ops->update) {
974:       (*snes->ops->update)(snes, snes->iter);
975:     }
976:     if (fas->fastype == SNES_FAS_MULTIPLICATIVE) {
977:       SNESFASCycle_Multiplicative(snes, X);
978:     } else if (fas->fastype == SNES_FAS_ADDITIVE) {
979:       SNESFASCycle_Additive(snes, X);
980:     } else if (fas->fastype == SNES_FAS_FULL) {
981:       SNESFASCycle_Full(snes, X);
982:     } else if (fas->fastype ==SNES_FAS_KASKADE) {
983:       SNESFASCycle_Kaskade(snes, X);
984:     } else {
985:       SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported FAS type");
986:     }

988:     /* check for FAS cycle divergence */
989:     if (snes->reason != SNES_CONVERGED_ITERATING) return(0);

991:     /* Monitor convergence */
992:     PetscObjectSAWsTakeAccess((PetscObject)snes);
993:     snes->iter = i+1;
994:     PetscObjectSAWsGrantAccess((PetscObject)snes);
995:     SNESLogConvergenceHistory(snes,snes->norm,0);
996:     SNESMonitor(snes,snes->iter,snes->norm);
997:     /* Test for convergence */
998:     if (isFine) {
999:       (*snes->ops->converged)(snes,snes->iter,0.0,0.0,snes->norm,&snes->reason,snes->cnvP);
1000:       if (snes->reason) break;
1001:     }
1002:   }
1003:   if (i == maxits) {
1004:     PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);
1005:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
1006:   }
1007:   return(0);
1008: }