Actual source code: fas.c

petsc-3.9.2 2018-05-20
Report Typos and Errors
  1: /* Defines the basic SNES object */
  2:  #include <../src/snes/impls/fas/fasimpls.h>

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

  6: extern PetscErrorCode SNESFASGalerkinDefaultFunction(SNES,Vec,Vec,void*);

  8: static PetscErrorCode SNESReset_FAS(SNES snes)
  9: {
 10:   PetscErrorCode 0;
 11:   SNES_FAS       * fas = (SNES_FAS*)snes->data;

 14:   SNESDestroy(&fas->smoothu);
 15:   SNESDestroy(&fas->smoothd);
 16:   MatDestroy(&fas->inject);
 17:   MatDestroy(&fas->interpolate);
 18:   MatDestroy(&fas->restrct);
 19:   VecDestroy(&fas->rscale);
 20:   if (fas->galerkin) {
 21:     VecDestroy(&fas->Xg);
 22:     VecDestroy(&fas->Fg);
 23:   }
 24:   if (fas->next) {SNESReset(fas->next);}
 25:   return(0);
 26: }

 28: static PetscErrorCode SNESDestroy_FAS(SNES snes)
 29: {
 30:   SNES_FAS       * fas = (SNES_FAS*)snes->data;
 31:   PetscErrorCode 0;

 34:   /* recursively resets and then destroys */
 35:   SNESReset(snes);
 36:   if (fas->next) {
 37:     SNESDestroy(&fas->next);
 38:   }
 39:   PetscFree(fas);
 40:   return(0);
 41: }

 43: static PetscErrorCode SNESSetUp_FAS(SNES snes)
 44: {
 45:   SNES_FAS       *fas = (SNES_FAS*) snes->data;
 47:   PetscInt       dm_levels;
 48:   Vec            vec_sol, vec_func, vec_sol_update, vec_rhs; /* preserve these if they're set through the reset */
 49:   SNES           next;
 50:   PetscBool      isFine, hasCreateRestriction;
 51:   SNESLineSearch linesearch;
 52:   SNESLineSearch slinesearch;
 53:   void           *lsprectx,*lspostctx;
 54:   PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*);
 55:   PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);

 58:   SNESFASCycleIsFine(snes, &isFine);
 59:   if (fas->usedmfornumberoflevels && isFine) {
 60:     DMGetRefineLevel(snes->dm,&dm_levels);
 61:     dm_levels++;
 62:     if (dm_levels > fas->levels) {
 63:       /* we don't want the solution and func vectors to be destroyed in the SNESReset when it's called in SNESFASSetLevels_FAS*/
 64:       vec_sol              = snes->vec_sol;
 65:       vec_func             = snes->vec_func;
 66:       vec_sol_update       = snes->vec_sol_update;
 67:       vec_rhs              = snes->vec_rhs;
 68:       snes->vec_sol        = NULL;
 69:       snes->vec_func       = NULL;
 70:       snes->vec_sol_update = NULL;
 71:       snes->vec_rhs        = NULL;

 73:       /* reset the number of levels */
 74:       SNESFASSetLevels(snes,dm_levels,NULL);
 75:       SNESSetFromOptions(snes);

 77:       snes->vec_sol        = vec_sol;
 78:       snes->vec_func       = vec_func;
 79:       snes->vec_rhs        = vec_rhs;
 80:       snes->vec_sol_update = vec_sol_update;
 81:     }
 82:   }
 83:   SNESFASCycleGetCorrection(snes, &next);
 84:   if (!isFine) snes->gridsequence = 0; /* no grid sequencing inside the multigrid hierarchy! */

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

 88:   /* set up the smoothers if they haven't already been set up */
 89:   if (!fas->smoothd) {
 90:     SNESFASCycleCreateSmoother_Private(snes, &fas->smoothd);
 91:   }

 93:   if (snes->dm) {
 94:     /* set the smoother DMs properly */
 95:     if (fas->smoothu) SNESSetDM(fas->smoothu, snes->dm);
 96:     SNESSetDM(fas->smoothd, snes->dm);
 97:     /* construct EVERYTHING from the DM -- including the progressive set of smoothers */
 98:     if (next) {
 99:       /* for now -- assume the DM and the evaluation functions have been set externally */
100:       if (!next->dm) {
101:         DMCoarsen(snes->dm, PetscObjectComm((PetscObject)next), &next->dm);
102:         SNESSetDM(next, next->dm);
103:       }
104:       /* set the interpolation and restriction from the DM */
105:       if (!fas->interpolate) {
106:         DMCreateInterpolation(next->dm, snes->dm, &fas->interpolate, &fas->rscale);
107:         if (!fas->restrct) {
108:           DMHasCreateRestriction(next->dm, &hasCreateRestriction);
109:           /* DM can create restrictions, use that */
110:           if (hasCreateRestriction) {
111:             DMCreateRestriction(next->dm, snes->dm, &fas->restrct);
112:           } else {
113:             PetscObjectReference((PetscObject)fas->interpolate);
114:             fas->restrct = fas->interpolate;
115:           }
116:         }
117:       }
118:       /* set the injection from the DM */
119:       if (!fas->inject) {
120:         PetscBool flg;
121:         DMHasCreateInjection(next->dm, &flg);
122:         if (flg) {
123:           DMCreateInjection(next->dm, snes->dm, &fas->inject);
124:         }
125:       }
126:     }
127:   }
128:   /*pass the smoother, function, and jacobian up to the next level if it's not user set already */
129:   if (fas->galerkin) {
130:     if (next) {
131:       SNESSetFunction(next, NULL, SNESFASGalerkinDefaultFunction, next);
132:     }
133:     if (fas->smoothd && fas->level != fas->levels - 1) {
134:       SNESSetFunction(fas->smoothd, NULL, SNESFASGalerkinDefaultFunction, snes);
135:     }
136:     if (fas->smoothu && fas->level != fas->levels - 1) {
137:       SNESSetFunction(fas->smoothu, NULL, SNESFASGalerkinDefaultFunction, snes);
138:     }
139:   }

141:   /* sets the down (pre) smoother's default norm and sets it from options */
142:   if (fas->smoothd) {
143:     if (fas->level == 0 && fas->levels != 1) {
144:       SNESSetNormSchedule(fas->smoothd, SNES_NORM_NONE);
145:     } else {
146:       SNESSetNormSchedule(fas->smoothd, SNES_NORM_FINAL_ONLY);
147:     }
148:     PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothd);
149:     SNESSetFromOptions(fas->smoothd);
150:     SNESGetLineSearch(snes,&linesearch);
151:     SNESGetLineSearch(fas->smoothd,&slinesearch);
152:     SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
153:     SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
154:     SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
155:     SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
156:     PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);

158:     fas->smoothd->vec_sol        = snes->vec_sol;
159:     PetscObjectReference((PetscObject)snes->vec_sol);
160:     fas->smoothd->vec_sol_update = snes->vec_sol_update;
161:     PetscObjectReference((PetscObject)snes->vec_sol_update);
162:     fas->smoothd->vec_func       = snes->vec_func;
163:     PetscObjectReference((PetscObject)snes->vec_func);

165:     if (fas->eventsmoothsetup) {PetscLogEventBegin(fas->eventsmoothsetup,0,0,0,0);}
166:     SNESSetUp(fas->smoothd);
167:     if (fas->eventsmoothsetup) {PetscLogEventEnd(fas->eventsmoothsetup,0,0,0,0);}
168:   }

170:   /* sets the up (post) smoother's default norm and sets it from options */
171:   if (fas->smoothu) {
172:     if (fas->level != fas->levels - 1) {
173:       SNESSetNormSchedule(fas->smoothu, SNES_NORM_NONE);
174:     } else {
175:       SNESSetNormSchedule(fas->smoothu, SNES_NORM_FINAL_ONLY);
176:     }
177:     PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothu);
178:     SNESSetFromOptions(fas->smoothu);
179:     SNESGetLineSearch(snes,&linesearch);
180:     SNESGetLineSearch(fas->smoothu,&slinesearch);
181:     SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
182:     SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
183:     SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
184:     SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
185:     PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);

187:     fas->smoothu->vec_sol        = snes->vec_sol;
188:     PetscObjectReference((PetscObject)snes->vec_sol);
189:     fas->smoothu->vec_sol_update = snes->vec_sol_update;
190:     PetscObjectReference((PetscObject)snes->vec_sol_update);
191:     fas->smoothu->vec_func       = snes->vec_func;
192:     PetscObjectReference((PetscObject)snes->vec_func);

194:     if (fas->eventsmoothsetup) {PetscLogEventBegin(fas->eventsmoothsetup,0,0,0,0);}
195:     SNESSetUp(fas->smoothu);
196:     if (fas->eventsmoothsetup) {PetscLogEventEnd(fas->eventsmoothsetup,0,0,0,0);}

198:   }

200:   if (next) {
201:     /* gotta set up the solution vector for this to work */
202:     if (!next->vec_sol) {SNESFASCreateCoarseVec(snes,&next->vec_sol);}
203:     if (!next->vec_rhs) {SNESFASCreateCoarseVec(snes,&next->vec_rhs);}
204:     PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next);
205:     SNESGetLineSearch(snes,&linesearch);
206:     SNESGetLineSearch(fas->next,&slinesearch);
207:     SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
208:     SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
209:     SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
210:     SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
211:     PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
212:     SNESSetUp(next);
213:   }
214:   /* setup FAS work vectors */
215:   if (fas->galerkin) {
216:     VecDuplicate(snes->vec_sol, &fas->Xg);
217:     VecDuplicate(snes->vec_sol, &fas->Fg);
218:   }
219:   return(0);
220: }

222: static PetscErrorCode SNESSetFromOptions_FAS(PetscOptionItems *PetscOptionsObject,SNES snes)
223: {
224:   SNES_FAS       *fas   = (SNES_FAS*) snes->data;
225:   PetscInt       levels = 1;
226:   PetscBool      flg    = PETSC_FALSE, upflg = PETSC_FALSE, downflg = PETSC_FALSE, monflg = PETSC_FALSE, galerkinflg = PETSC_FALSE,continuationflg = PETSC_FALSE;
228:   SNESFASType    fastype;
229:   const char     *optionsprefix;
230:   SNESLineSearch linesearch;
231:   PetscInt       m, n_up, n_down;
232:   SNES           next;
233:   PetscBool      isFine;

236:   SNESFASCycleIsFine(snes, &isFine);
237:   PetscOptionsHead(PetscOptionsObject,"SNESFAS Options-----------------------------------");

239:   /* number of levels -- only process most options on the finest level */
240:   if (isFine) {
241:     PetscOptionsInt("-snes_fas_levels", "Number of Levels", "SNESFASSetLevels", levels, &levels, &flg);
242:     if (!flg && snes->dm) {
243:       DMGetRefineLevel(snes->dm,&levels);
244:       levels++;
245:       fas->usedmfornumberoflevels = PETSC_TRUE;
246:     }
247:     SNESFASSetLevels(snes, levels, NULL);
248:     fastype = fas->fastype;
249:     PetscOptionsEnum("-snes_fas_type","FAS correction type","SNESFASSetType",SNESFASTypes,(PetscEnum)fastype,(PetscEnum*)&fastype,&flg);
250:     if (flg) {
251:       SNESFASSetType(snes, fastype);
252:     }

254:     SNESGetOptionsPrefix(snes, &optionsprefix);
255:     PetscOptionsInt("-snes_fas_cycles","Number of cycles","SNESFASSetCycles",fas->n_cycles,&m,&flg);
256:     if (flg) {
257:       SNESFASSetCycles(snes, m);
258:     }
259:     PetscOptionsBool("-snes_fas_continuation","Corrected grid-sequence continuation","SNESFASSetContinuation",fas->continuation,&continuationflg,&flg);
260:     if (flg) {
261:       SNESFASSetContinuation(snes,continuationflg);
262:     }

264:     PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin","SNESFASSetGalerkin",fas->galerkin,&galerkinflg,&flg);
265:     if (flg) {
266:       SNESFASSetGalerkin(snes, galerkinflg);
267:     }

269:     if (fas->fastype == SNES_FAS_FULL) {
270:       PetscOptionsBool("-snes_fas_full_downsweep","Smooth on the initial upsweep for full FAS cycles","SNESFASFullSetDownSweep",fas->full_downsweep,&fas->full_downsweep,&flg);
271:       if (flg) {SNESFASFullSetDownSweep(snes,fas->full_downsweep);}
272:     }

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

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

278:     {
279:       PetscViewer viewer;
280:       PetscViewerFormat format;
281:       PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,
282:                                    "-snes_fas_monitor",&viewer,&format,&monflg);
283:       if (monflg) {
284:         PetscViewerAndFormat *vf;
285:         PetscViewerAndFormatCreate(viewer,format,&vf);
286:         PetscObjectDereference((PetscObject)viewer);
287:         SNESFASSetMonitor(snes,vf,PETSC_TRUE);
288:       }
289:     }
290:     flg    = PETSC_FALSE;
291:     monflg = PETSC_TRUE;
292:     PetscOptionsBool("-snes_fas_log","Log times for each FAS level","SNESFASSetLog",monflg,&monflg,&flg);
293:     if (flg) {SNESFASSetLog(snes,monflg);}
294:   }

296:   PetscOptionsTail();
297:   /* setup from the determined types if there is no pointwise procedure or smoother defined */
298:   if (upflg) {
299:     SNESFASSetNumberSmoothUp(snes,n_up);
300:   }
301:   if (downflg) {
302:     SNESFASSetNumberSmoothDown(snes,n_down);
303:   }

305:   /* set up the default line search for coarse grid corrections */
306:   if (fas->fastype == SNES_FAS_ADDITIVE) {
307:     if (!snes->linesearch) {
308:       SNESGetLineSearch(snes, &linesearch);
309:       SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
310:     }
311:   }

313:   SNESFASCycleGetCorrection(snes, &next);
314:   /* recursive option setting for the smoothers */
315:   if (next) {SNESSetFromOptions(next);}
316:   return(0);
317: }

319:  #include <petscdraw.h>
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;
326:   SNES           smoothu, smoothd, levelsnes;

329:   SNESFASCycleIsFine(snes, &isFine);
330:   if (isFine) {
331:     PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
332:     PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
333:     if (iascii) {
334:       PetscViewerASCIIPrintf(viewer, "  type is %s, levels=%D, cycles=%D\n",  SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles);
335:       if (fas->galerkin) {
336:         PetscViewerASCIIPrintf(viewer,"  Using Galerkin computed coarse grid function evaluation\n");
337:       } else {
338:         PetscViewerASCIIPrintf(viewer,"  Not using Galerkin computed coarse grid function evaluation\n");
339:       }
340:       for (i=0; i<fas->levels; i++) {
341:         SNESFASGetCycleSNES(snes, i, &levelsnes);
342:         SNESFASCycleGetSmootherUp(levelsnes, &smoothu);
343:         SNESFASCycleGetSmootherDown(levelsnes, &smoothd);
344:         if (!i) {
345:           PetscViewerASCIIPrintf(viewer,"  Coarse grid solver -- level %D -------------------------------\n",i);
346:         } else {
347:           PetscViewerASCIIPrintf(viewer,"  Down solver (pre-smoother) on level %D -------------------------------\n",i);
348:         }
349:         PetscViewerASCIIPushTab(viewer);
350:         if (smoothd) {
351:           SNESView(smoothd,viewer);
352:         } else {
353:           PetscViewerASCIIPrintf(viewer,"Not yet available\n");
354:         }
355:         PetscViewerASCIIPopTab(viewer);
356:         if (i && (smoothd == smoothu)) {
357:           PetscViewerASCIIPrintf(viewer,"  Up solver (post-smoother) same as down solver (pre-smoother)\n");
358:         } else if (i) {
359:           PetscViewerASCIIPrintf(viewer,"  Up solver (post-smoother) on level %D -------------------------------\n",i);
360:           PetscViewerASCIIPushTab(viewer);
361:           if (smoothu) {
362:             SNESView(smoothu,viewer);
363:           } else {
364:             PetscViewerASCIIPrintf(viewer,"Not yet available\n");
365:           }
366:           PetscViewerASCIIPopTab(viewer);
367:         }
368:       }
369:     } else if (isdraw) {
370:       PetscDraw draw;
371:       PetscReal x,w,y,bottom,th,wth;
372:       SNES_FAS  *curfas = fas;
373:       PetscViewerDrawGetDraw(viewer,0,&draw);
374:       PetscDrawGetCurrentPoint(draw,&x,&y);
375:       PetscDrawStringGetSize(draw,&wth,&th);
376:       bottom = y - th;
377:       while (curfas) {
378:         if (!curfas->smoothu) {
379:           PetscDrawPushCurrentPoint(draw,x,bottom);
380:           if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
381:           PetscDrawPopCurrentPoint(draw);
382:         } else {
383:           w    = 0.5*PetscMin(1.0-x,x);
384:           PetscDrawPushCurrentPoint(draw,x-w,bottom);
385:           if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
386:           PetscDrawPopCurrentPoint(draw);
387:           PetscDrawPushCurrentPoint(draw,x+w,bottom);
388:           if (curfas->smoothu) SNESView(curfas->smoothu,viewer);
389:           PetscDrawPopCurrentPoint(draw);
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:   return(0);
399: }

401: /*
402: Defines the action of the downsmoother
403:  */
404: static PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
405: {
406:   PetscErrorCode      0;
407:   SNESConvergedReason reason;
408:   Vec                 FPC;
409:   SNES                smoothd;
410:   PetscBool           flg;
411:   SNES_FAS            *fas = (SNES_FAS*) snes->data;

414:   SNESFASCycleGetSmootherDown(snes, &smoothd);
415:   SNESSetInitialFunction(smoothd, F);
416:   if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,0,0,0,0);}
417:   SNESSolve(smoothd, B, X);
418:   if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,0,0,0,0);}
419:   /* check convergence reason for the smoother */
420:   SNESGetConvergedReason(smoothd,&reason);
421:   if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
422:     snes->reason = SNES_DIVERGED_INNER;
423:     return(0);
424:   }

426:   SNESGetFunction(smoothd, &FPC, NULL, NULL);
427:   SNESGetAlwaysComputesFinalResidual(smoothd, &flg);
428:   if (!flg) {
429:     SNESComputeFunction(smoothd, X, FPC);
430:   }
431:   VecCopy(FPC, F);
432:   if (fnorm) {VecNorm(F,NORM_2,fnorm);}
433:   return(0);
434: }


437: /*
438: Defines the action of the upsmoother
439:  */
440: static PetscErrorCode SNESFASUpSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
441: {
442:   PetscErrorCode      0;
443:   SNESConvergedReason reason;
444:   Vec                 FPC;
445:   SNES                smoothu;
446:   PetscBool           flg;
447:   SNES_FAS            *fas = (SNES_FAS*) snes->data;

450:   SNESFASCycleGetSmootherUp(snes, &smoothu);
451:   if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,0,0,0,0);}
452:   SNESSolve(smoothu, B, X);
453:   if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,0,0,0,0);}
454:   /* check convergence reason for the smoother */
455:   SNESGetConvergedReason(smoothu,&reason);
456:   if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
457:     snes->reason = SNES_DIVERGED_INNER;
458:     return(0);
459:   }
460:   SNESGetFunction(smoothu, &FPC, NULL, NULL);
461:   SNESGetAlwaysComputesFinalResidual(smoothu, &flg);
462:   if (!flg) {
463:     SNESComputeFunction(smoothu, X, FPC);
464:   }
465:   VecCopy(FPC, F);
466:   if (fnorm) {VecNorm(F,NORM_2,fnorm);}
467:   return(0);
468: }

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

473:    Collective

475:    Input Arguments:
476: .  snes - SNESFAS

478:    Output Arguments:
479: .  Xcoarse - vector on level one coarser than snes

481:    Level: developer

483: .seealso: SNESFASSetRestriction(), SNESFASRestrict()
484: @*/
485: PetscErrorCode SNESFASCreateCoarseVec(SNES snes,Vec *Xcoarse)
486: {
488:   SNES_FAS       *fas = (SNES_FAS*)snes->data;

491:   if (fas->rscale) {
492:     VecDuplicate(fas->rscale,Xcoarse);
493:   } else if (fas->inject) {
494:     MatCreateVecs(fas->inject,Xcoarse,NULL);
495:   } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or injection");
496:   return(0);
497: }

499: /*@
500:    SNESFASRestrict - restrict a Vec to the next coarser level

502:    Collective

504:    Input Arguments:
505: +  fine - SNES from which to restrict
506: -  Xfine - vector to restrict

508:    Output Arguments:
509: .  Xcoarse - result of restriction

511:    Level: developer

513: .seealso: SNESFASSetRestriction(), SNESFASSetInjection()
514: @*/
515: PetscErrorCode SNESFASRestrict(SNES fine,Vec Xfine,Vec Xcoarse)
516: {
518:   SNES_FAS       *fas = (SNES_FAS*)fine->data;

524:   if (fas->inject) {
525:     MatRestrict(fas->inject,Xfine,Xcoarse);
526:   } else {
527:     MatRestrict(fas->restrct,Xfine,Xcoarse);
528:     VecPointwiseMult(Xcoarse,fas->rscale,Xcoarse);
529:   }
530:   return(0);
531: }

533: /*

535: Performs the FAS coarse correction as:

537: fine problem:   F(x) = b
538: coarse problem: F^c(x^c) = b^c

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

542:  */
543: PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new)
544: {
545:   PetscErrorCode      ierr;
546:   Vec                 X_c, Xo_c, F_c, B_c;
547:   SNESConvergedReason reason;
548:   SNES                next;
549:   Mat                 restrct, interpolate;
550:   SNES_FAS            *fasc;

553:   SNESFASCycleGetCorrection(snes, &next);
554:   if (next) {
555:     fasc = (SNES_FAS*)next->data;

557:     SNESFASCycleGetRestriction(snes, &restrct);
558:     SNESFASCycleGetInterpolation(snes, &interpolate);

560:     X_c  = next->vec_sol;
561:     Xo_c = next->work[0];
562:     F_c  = next->vec_func;
563:     B_c  = next->vec_rhs;

565:     if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,0,0,0,0);}
566:     SNESFASRestrict(snes,X,Xo_c);
567:     /* restrict the defect: R(F(x) - b) */
568:     MatRestrict(restrct, F, B_c);
569:     if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,0,0,0,0);}

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

576:     /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */
577:     VecCopy(B_c, X_c);
578:     VecCopy(F_c, B_c);
579:     VecCopy(X_c, F_c);
580:     /* set initial guess of the coarse problem to the projected fine solution */
581:     VecCopy(Xo_c, X_c);

583:     /* recurse to the next level */
584:     SNESSetInitialFunction(next, F_c);
585:     SNESSolve(next, B_c, X_c);
586:     SNESGetConvergedReason(next,&reason);
587:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
588:       snes->reason = SNES_DIVERGED_INNER;
589:       return(0);
590:     }
591:     /* correct as x <- x + I(x^c - Rx)*/
592:     VecAXPY(X_c, -1.0, Xo_c);

594:     if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,0,0,0,0);}
595:     MatInterpolateAdd(interpolate, X_c, X, X_new);
596:     if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,0,0,0,0);}
597:   }
598:   return(0);
599: }

601: /*

603: The additive cycle looks like:

605: xhat = x
606: xhat = dS(x, b)
607: x = coarsecorrection(xhat, b_d)
608: x = x + nu*(xhat - x);
609: (optional) x = uS(x, b)

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

613:  */
614: static PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X)
615: {
616:   Vec                  F, B, Xhat;
617:   Vec                  X_c, Xo_c, F_c, B_c;
618:   PetscErrorCode       ierr;
619:   SNESConvergedReason  reason;
620:   PetscReal            xnorm, fnorm, ynorm;
621:   SNESLineSearchReason lsresult;
622:   SNES                 next;
623:   Mat                  restrct, interpolate;
624:   SNES_FAS             *fas = (SNES_FAS*)snes->data,*fasc;

627:   SNESFASCycleGetCorrection(snes, &next);
628:   F    = snes->vec_func;
629:   B    = snes->vec_rhs;
630:   Xhat = snes->work[1];
631:   VecCopy(X, Xhat);
632:   /* recurse first */
633:   if (next) {
634:     fasc = (SNES_FAS*)next->data;
635:     SNESFASCycleGetRestriction(snes, &restrct);
636:     SNESFASCycleGetInterpolation(snes, &interpolate);
637:     if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,0,0,0,0);}
638:     SNESComputeFunction(snes, Xhat, F);
639:     if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,0,0,0,0);}
640:     VecNorm(F, NORM_2, &fnorm);
641:     X_c  = next->vec_sol;
642:     Xo_c = next->work[0];
643:     F_c  = next->vec_func;
644:     B_c  = next->vec_rhs;

646:     SNESFASRestrict(snes,Xhat,Xo_c);
647:     /* restrict the defect */
648:     MatRestrict(restrct, F, B_c);

650:     /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */
651:     if (fasc->eventresidual) {PetscLogEventBegin(fasc->eventresidual,0,0,0,0);}
652:     SNESComputeFunction(next, Xo_c, F_c);
653:     if (fasc->eventresidual) {PetscLogEventEnd(fasc->eventresidual,0,0,0,0);}
654:     VecCopy(B_c, X_c);
655:     VecCopy(F_c, B_c);
656:     VecCopy(X_c, F_c);
657:     /* set initial guess of the coarse problem to the projected fine solution */
658:     VecCopy(Xo_c, X_c);

660:     /* recurse */
661:     SNESSetInitialFunction(next, F_c);
662:     SNESSolve(next, B_c, X_c);

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

667:     SNESGetConvergedReason(next,&reason);
668:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
669:       snes->reason = SNES_DIVERGED_INNER;
670:       return(0);
671:     }

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

677:     /* additive correction of the coarse direction*/
678:     SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat);
679:     SNESLineSearchGetReason(snes->linesearch, &lsresult);
680:     SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm);
681:     if (lsresult) {
682:       if (++snes->numFailures >= snes->maxFailures) {
683:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
684:         return(0);
685:       }
686:     }
687:   } else {
688:     SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
689:   }
690:   return(0);
691: }

693: /*

695: Defines the FAS cycle as:

697: fine problem: F(x) = b
698: coarse problem: F^c(x) = b^c

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

702: correction:

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

706:  */
707: static PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X)
708: {

711:   Vec            F,B;
712:   SNES           next;

715:   F = snes->vec_func;
716:   B = snes->vec_rhs;
717:   /* pre-smooth -- just update using the pre-smoother */
718:   SNESFASCycleGetCorrection(snes,&next);
719:   SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
720:   if (next) {
721:     SNESFASCoarseCorrection(snes, X, F, X);
722:     SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
723:   }
724:   return(0);
725: }

727: static PetscErrorCode SNESFASCycleSetupPhase_Full(SNES snes)
728: {
729:   SNES           next;
730:   SNES_FAS       *fas = (SNES_FAS*)snes->data;
731:   PetscBool      isFine;

735:   /* pre-smooth -- just update using the pre-smoother */
736:   SNESFASCycleIsFine(snes,&isFine);
737:   SNESFASCycleGetCorrection(snes,&next);
738:   fas->full_stage = 0;
739:   if (next) {SNESFASCycleSetupPhase_Full(next);}
740:   return(0);
741: }

743: static PetscErrorCode SNESFASCycle_Full(SNES snes, Vec X)
744: {
746:   Vec            F,B;
747:   SNES_FAS       *fas = (SNES_FAS*)snes->data;
748:   PetscBool      isFine;
749:   SNES           next;

752:   F = snes->vec_func;
753:   B = snes->vec_rhs;
754:   SNESFASCycleIsFine(snes,&isFine);
755:   SNESFASCycleGetCorrection(snes,&next);

757:   if (isFine) {
758:     SNESFASCycleSetupPhase_Full(snes);
759:   }

761:   if (fas->full_stage == 0) {
762:     /* downsweep */
763:     if (next) {
764:       if (fas->level != 1) next->max_its += 1;
765:       if (fas->full_downsweep||isFine) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
766:       SNESFASCoarseCorrection(snes,X,F,X);
767:       SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
768:       if (fas->level != 1) next->max_its -= 1;
769:     } else {
770:       /* The smoother on the coarse level is the coarse solver */
771:       SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
772:     }
773:     fas->full_stage = 1;
774:   } else if (fas->full_stage == 1) {
775:     if (snes->iter == 0) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
776:     if (next) {
777:       SNESFASCoarseCorrection(snes,X,F,X);
778:       SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
779:     }
780:   }
781:   /* final v-cycle */
782:   if (isFine) {
783:     if (next) {
784:       SNESFASCoarseCorrection(snes,X,F,X);
785:       SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
786:     }
787:   }
788:   return(0);
789: }

791: static PetscErrorCode SNESFASCycle_Kaskade(SNES snes, Vec X)
792: {
794:   Vec            F,B;
795:   SNES           next;

798:   F = snes->vec_func;
799:   B = snes->vec_rhs;
800:   SNESFASCycleGetCorrection(snes,&next);
801:   if (next) {
802:     SNESFASCoarseCorrection(snes,X,F,X);
803:     SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
804:   } else {
805:     SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
806:   }
807:   return(0);
808: }

810: PetscBool SNEScite = PETSC_FALSE;
811: const char SNESCitation[] = "@techreport{pbmkbsxt2012,\n"
812:                             "  title = {Composing Scalable Nonlinear Algebraic Solvers},\n"
813:                             "  author = {Peter Brune and Mathew Knepley and Barry Smith and Xuemin Tu},\n"
814:                             "  year = 2013,\n"
815:                             "  type = Preprint,\n"
816:                             "  number = {ANL/MCS-P2010-0112},\n"
817:                             "  institution = {Argonne National Laboratory}\n}\n";

819: static PetscErrorCode SNESSolve_FAS(SNES snes)
820: {
822:   PetscInt       i, maxits;
823:   Vec            X, F;
824:   PetscReal      fnorm;
825:   SNES_FAS       *fas = (SNES_FAS*)snes->data,*ffas;
826:   DM             dm;
827:   PetscBool      isFine;


831:   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);

833:   PetscCitationsRegister(SNESCitation,&SNEScite);
834:   maxits       = snes->max_its;      /* maximum number of iterations */
835:   snes->reason = SNES_CONVERGED_ITERATING;
836:   X            = snes->vec_sol;
837:   F            = snes->vec_func;

839:   SNESFASCycleIsFine(snes, &isFine);
840:   /*norm setup */
841:   PetscObjectSAWsTakeAccess((PetscObject)snes);
842:   snes->iter = 0;
843:   snes->norm = 0.;
844:   PetscObjectSAWsGrantAccess((PetscObject)snes);
845:   if (!snes->vec_func_init_set) {
846:     if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,0,0,0,0);}
847:     SNESComputeFunction(snes,X,F);
848:     if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,0,0,0,0);}
849:   } else snes->vec_func_init_set = PETSC_FALSE;

851:   VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
852:   SNESCheckFunctionNorm(snes,fnorm);
853:   PetscObjectSAWsTakeAccess((PetscObject)snes);
854:   snes->norm = fnorm;
855:   PetscObjectSAWsGrantAccess((PetscObject)snes);
856:   SNESLogConvergenceHistory(snes,fnorm,0);
857:   SNESMonitor(snes,0,fnorm);

859:   /* test convergence */
860:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
861:   if (snes->reason) return(0);


864:   if (isFine) {
865:     /* propagate scale-dependent data up the hierarchy */
866:     SNESGetDM(snes,&dm);
867:     for (ffas=fas; ffas->next; ffas=(SNES_FAS*)ffas->next->data) {
868:       DM dmcoarse;
869:       SNESGetDM(ffas->next,&dmcoarse);
870:       DMRestrict(dm,ffas->restrct,ffas->rscale,ffas->inject,dmcoarse);
871:       dm   = dmcoarse;
872:     }
873:   }

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

878:     if (snes->ops->update) {
879:       (*snes->ops->update)(snes, snes->iter);
880:     }
881:     if (fas->fastype == SNES_FAS_MULTIPLICATIVE) {
882:       SNESFASCycle_Multiplicative(snes, X);
883:     } else if (fas->fastype == SNES_FAS_ADDITIVE) {
884:       SNESFASCycle_Additive(snes, X);
885:     } else if (fas->fastype == SNES_FAS_FULL) {
886:       SNESFASCycle_Full(snes, X);
887:     } else if (fas->fastype ==SNES_FAS_KASKADE) {
888:       SNESFASCycle_Kaskade(snes, X);
889:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported FAS type");

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

894:     /* Monitor convergence */
895:     PetscObjectSAWsTakeAccess((PetscObject)snes);
896:     snes->iter = i+1;
897:     PetscObjectSAWsGrantAccess((PetscObject)snes);
898:     SNESLogConvergenceHistory(snes,snes->norm,0);
899:     SNESMonitor(snes,snes->iter,snes->norm);
900:     /* Test for convergence */
901:     if (isFine) {
902:       (*snes->ops->converged)(snes,snes->iter,0.0,0.0,snes->norm,&snes->reason,snes->cnvP);
903:       if (snes->reason) break;
904:     }
905:   }
906:   if (i == maxits) {
907:     PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);
908:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
909:   }
910:   return(0);
911: }

913: /*MC

915: SNESFAS - Full Approximation Scheme nonlinear multigrid solver.

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

921: Options Database:
922: +   -snes_fas_levels -  The number of levels
923: .   -snes_fas_cycles<1> -  The number of cycles -- 1 for V, 2 for W
924: .   -snes_fas_type<additive,multiplicative,full,kaskade>  -  Additive or multiplicative cycle
925: .   -snes_fas_galerkin<PETSC_FALSE> -  Form coarse problems by projection back upon the fine problem
926: .   -snes_fas_smoothup<1> -  The number of iterations of the post-smoother
927: .   -snes_fas_smoothdown<1> -  The number of iterations of the pre-smoother
928: .   -snes_fas_monitor -  Monitor progress of all of the levels
929: .   -snes_fas_full_downsweep<PETSC_FALSE> - call the downsmooth on the initial downsweep of full FAS
930: .   -fas_levels_snes_ -  SNES options for all smoothers
931: .   -fas_levels_cycle_snes_ -  SNES options for all cycles
932: .   -fas_levels_i_snes_ -  SNES options for the smoothers on level i
933: .   -fas_levels_i_cycle_snes_ - SNES options for the cycle on level i
934: -   -fas_coarse_snes_ -  SNES options for the coarsest smoother

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

941: Level: beginner

943:    References:
944: . 1. -  Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
945:    SIAM Review, 57(4), 2015

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

950: PETSC_EXTERN PetscErrorCode SNESCreate_FAS(SNES snes)
951: {
952:   SNES_FAS       *fas;

956:   snes->ops->destroy        = SNESDestroy_FAS;
957:   snes->ops->setup          = SNESSetUp_FAS;
958:   snes->ops->setfromoptions = SNESSetFromOptions_FAS;
959:   snes->ops->view           = SNESView_FAS;
960:   snes->ops->solve          = SNESSolve_FAS;
961:   snes->ops->reset          = SNESReset_FAS;

963:   snes->usesksp = PETSC_FALSE;
964:   snes->usesnpc = PETSC_FALSE;

966:   if (!snes->tolerancesset) {
967:     snes->max_funcs = 30000;
968:     snes->max_its   = 10000;
969:   }

971:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

973:   PetscNewLog(snes,&fas);

975:   snes->data                  = (void*) fas;
976:   fas->level                  = 0;
977:   fas->levels                 = 1;
978:   fas->n_cycles               = 1;
979:   fas->max_up_it              = 1;
980:   fas->max_down_it            = 1;
981:   fas->smoothu                = NULL;
982:   fas->smoothd                = NULL;
983:   fas->next                   = NULL;
984:   fas->previous               = NULL;
985:   fas->fine                   = snes;
986:   fas->interpolate            = NULL;
987:   fas->restrct                = NULL;
988:   fas->inject                 = NULL;
989:   fas->usedmfornumberoflevels = PETSC_FALSE;
990:   fas->fastype                = SNES_FAS_MULTIPLICATIVE;
991:   fas->full_downsweep         = PETSC_FALSE;

993:   fas->eventsmoothsetup    = 0;
994:   fas->eventsmoothsolve    = 0;
995:   fas->eventresidual       = 0;
996:   fas->eventinterprestrict = 0;
997:   return(0);
998: }