Actual source code: fas.c

petsc-master 2019-07-16
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: static PetscErrorCode SNESReset_FAS(SNES snes)
  7: {
  8:   PetscErrorCode 0;
  9:   SNES_FAS       * fas = (SNES_FAS*)snes->data;

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

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

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

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

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

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

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

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

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

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

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

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,fas->smoothd,0,0,0);}
166:     SNESSetUp(fas->smoothd);
167:     if (fas->eventsmoothsetup) {PetscLogEventEnd(fas->eventsmoothsetup,fas->smoothd,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:     if (snes->linesearch) {
180:       SNESGetLineSearch(snes,&linesearch);
181:       SNESGetLineSearch(fas->smoothu,&slinesearch);
182:       SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
183:       SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
184:       SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
185:       SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
186:       PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
187:     }

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

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

200:   }

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

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

240:   SNESFASCycleIsFine(snes, &isFine);
241:   PetscOptionsHead(PetscOptionsObject,"SNESFAS Options-----------------------------------");

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

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

268:     PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin","SNESFASSetGalerkin",fas->galerkin,&galerkinflg,&flg);
269:     if (flg) {
270:       SNESFASSetGalerkin(snes, galerkinflg);
271:     }

273:     if (fas->fastype == SNES_FAS_FULL) {
274:       PetscOptionsBool("-snes_fas_full_downsweep","Smooth on the initial down sweep for full FAS cycles","SNESFASFullSetDownSweep",fas->full_downsweep,&fas->full_downsweep,&flg);
275:       if (flg) {SNESFASFullSetDownSweep(snes,fas->full_downsweep);}
276:     }

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

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

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

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

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

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

322:  #include <petscdraw.h>
323: static PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer)
324: {
325:   SNES_FAS       *fas = (SNES_FAS*) snes->data;
326:   PetscBool      isFine,iascii,isdraw;
327:   PetscInt       i;
329:   SNES           smoothu, smoothd, levelsnes;

332:   SNESFASCycleIsFine(snes, &isFine);
333:   if (isFine) {
334:     PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
335:     PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
336:     if (iascii) {
337:       PetscViewerASCIIPrintf(viewer, "  type is %s, levels=%D, cycles=%D\n",  SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles);
338:       if (fas->galerkin) {
339:         PetscViewerASCIIPrintf(viewer,"  Using Galerkin computed coarse grid function evaluation\n");
340:       } else {
341:         PetscViewerASCIIPrintf(viewer,"  Not using Galerkin computed coarse grid function evaluation\n");
342:       }
343:       for (i=0; i<fas->levels; i++) {
344:         SNESFASGetCycleSNES(snes, i, &levelsnes);
345:         SNESFASCycleGetSmootherUp(levelsnes, &smoothu);
346:         SNESFASCycleGetSmootherDown(levelsnes, &smoothd);
347:         if (!i) {
348:           PetscViewerASCIIPrintf(viewer,"  Coarse grid solver -- level %D -------------------------------\n",i);
349:         } else {
350:           PetscViewerASCIIPrintf(viewer,"  Down solver (pre-smoother) on level %D -------------------------------\n",i);
351:         }
352:         PetscViewerASCIIPushTab(viewer);
353:         if (smoothd) {
354:           SNESView(smoothd,viewer);
355:         } else {
356:           PetscViewerASCIIPrintf(viewer,"Not yet available\n");
357:         }
358:         PetscViewerASCIIPopTab(viewer);
359:         if (i && (smoothd == smoothu)) {
360:           PetscViewerASCIIPrintf(viewer,"  Up solver (post-smoother) same as down solver (pre-smoother)\n");
361:         } else if (i) {
362:           PetscViewerASCIIPrintf(viewer,"  Up solver (post-smoother) on level %D -------------------------------\n",i);
363:           PetscViewerASCIIPushTab(viewer);
364:           if (smoothu) {
365:             SNESView(smoothu,viewer);
366:           } else {
367:             PetscViewerASCIIPrintf(viewer,"Not yet available\n");
368:           }
369:           PetscViewerASCIIPopTab(viewer);
370:         }
371:       }
372:     } else if (isdraw) {
373:       PetscDraw draw;
374:       PetscReal x,w,y,bottom,th,wth;
375:       SNES_FAS  *curfas = fas;
376:       PetscViewerDrawGetDraw(viewer,0,&draw);
377:       PetscDrawGetCurrentPoint(draw,&x,&y);
378:       PetscDrawStringGetSize(draw,&wth,&th);
379:       bottom = y - th;
380:       while (curfas) {
381:         if (!curfas->smoothu) {
382:           PetscDrawPushCurrentPoint(draw,x,bottom);
383:           if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
384:           PetscDrawPopCurrentPoint(draw);
385:         } else {
386:           w    = 0.5*PetscMin(1.0-x,x);
387:           PetscDrawPushCurrentPoint(draw,x-w,bottom);
388:           if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
389:           PetscDrawPopCurrentPoint(draw);
390:           PetscDrawPushCurrentPoint(draw,x+w,bottom);
391:           if (curfas->smoothu) SNESView(curfas->smoothu,viewer);
392:           PetscDrawPopCurrentPoint(draw);
393:         }
394:         /* this is totally bogus but we have no way of knowing how low the previous one was draw to */
395:         bottom -= 5*th;
396:         if (curfas->next) curfas = (SNES_FAS*)curfas->next->data;
397:         else curfas = NULL;
398:       }
399:     }
400:   }
401:   return(0);
402: }

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

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

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


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

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

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

476:    Collective

478:    Input Arguments:
479: .  snes - SNESFAS

481:    Output Arguments:
482: .  Xcoarse - vector on level one coarser than snes

484:    Level: developer

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

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

502: /*@
503:    SNESFASRestrict - restrict a Vec to the next coarser level

505:    Collective

507:    Input Arguments:
508: +  fine - SNES from which to restrict
509: -  Xfine - vector to restrict

511:    Output Arguments:
512: .  Xcoarse - result of restriction

514:    Level: developer

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

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

536: /*

538: Performs the FAS coarse correction as:

540: fine problem:   F(x) = b
541: coarse problem: F^c(x^c) = b^c

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

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

556:   SNESFASCycleGetCorrection(snes, &next);
557:   if (next) {
558:     fasc = (SNES_FAS*)next->data;

560:     SNESFASCycleGetRestriction(snes, &restrct);
561:     SNESFASCycleGetInterpolation(snes, &interpolate);

563:     X_c  = next->vec_sol;
564:     Xo_c = next->work[0];
565:     F_c  = next->vec_func;
566:     B_c  = next->vec_rhs;

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

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

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

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

597:     if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,snes,0,0,0);}
598:     MatInterpolateAdd(interpolate, X_c, X, X_new);
599:     PetscObjectSetName((PetscObject) X_c, "Coarse correction");
600:     VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view");
601:     PetscObjectSetName((PetscObject) X_new, "Updated Fine solution");
602:     VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view");
603:     if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,snes,0,0,0);}
604:   }
605:   return(0);
606: }

608: /*

610: The additive cycle looks like:

612: xhat = x
613: xhat = dS(x, b)
614: x = coarsecorrection(xhat, b_d)
615: x = x + nu*(xhat - x);
616: (optional) x = uS(x, b)

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

620:  */
621: static PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X)
622: {
623:   Vec                  F, B, Xhat;
624:   Vec                  X_c, Xo_c, F_c, B_c;
625:   PetscErrorCode       ierr;
626:   SNESConvergedReason  reason;
627:   PetscReal            xnorm, fnorm, ynorm;
628:   SNESLineSearchReason lsresult;
629:   SNES                 next;
630:   Mat                  restrct, interpolate;
631:   SNES_FAS             *fas = (SNES_FAS*)snes->data,*fasc;

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

653:     SNESFASRestrict(snes,Xhat,Xo_c);
654:     /* restrict the defect */
655:     MatRestrict(restrct, F, B_c);

657:     /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */
658:     if (fasc->eventresidual) {PetscLogEventBegin(fasc->eventresidual,next,0,0,0);}
659:     SNESComputeFunction(next, Xo_c, F_c);
660:     if (fasc->eventresidual) {PetscLogEventEnd(fasc->eventresidual,next,0,0,0);}
661:     VecCopy(B_c, X_c);
662:     VecCopy(F_c, B_c);
663:     VecCopy(X_c, F_c);
664:     /* set initial guess of the coarse problem to the projected fine solution */
665:     VecCopy(Xo_c, X_c);

667:     /* recurse */
668:     SNESSetInitialFunction(next, F_c);
669:     SNESSolve(next, B_c, X_c);

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

674:     SNESGetConvergedReason(next,&reason);
675:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
676:       snes->reason = SNES_DIVERGED_INNER;
677:       return(0);
678:     }

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

684:     /* additive correction of the coarse direction*/
685:     SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat);
686:     SNESLineSearchGetReason(snes->linesearch, &lsresult);
687:     SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm);
688:     if (lsresult) {
689:       if (++snes->numFailures >= snes->maxFailures) {
690:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
691:         return(0);
692:       }
693:     }
694:   } else {
695:     SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
696:   }
697:   return(0);
698: }

700: /*

702: Defines the FAS cycle as:

704: fine problem: F(x) = b
705: coarse problem: F^c(x) = b^c

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

709: correction:

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

713:  */
714: static PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X)
715: {

718:   Vec            F,B;
719:   SNES           next;

722:   F = snes->vec_func;
723:   B = snes->vec_rhs;
724:   /* pre-smooth -- just update using the pre-smoother */
725:   SNESFASCycleGetCorrection(snes,&next);
726:   SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
727:   if (next) {
728:     SNESFASCoarseCorrection(snes, X, F, X);
729:     SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
730:   }
731:   return(0);
732: }

734: static PetscErrorCode SNESFASCycleSetupPhase_Full(SNES snes)
735: {
736:   SNES           next;
737:   SNES_FAS       *fas = (SNES_FAS*)snes->data;
738:   PetscBool      isFine;

742:   /* pre-smooth -- just update using the pre-smoother */
743:   SNESFASCycleIsFine(snes,&isFine);
744:   SNESFASCycleGetCorrection(snes,&next);
745:   fas->full_stage = 0;
746:   if (next) {SNESFASCycleSetupPhase_Full(next);}
747:   return(0);
748: }

750: static PetscErrorCode SNESFASCycle_Full(SNES snes, Vec X)
751: {
753:   Vec            F,B;
754:   SNES_FAS       *fas = (SNES_FAS*)snes->data;
755:   PetscBool      isFine;
756:   SNES           next;

759:   F = snes->vec_func;
760:   B = snes->vec_rhs;
761:   SNESFASCycleIsFine(snes,&isFine);
762:   SNESFASCycleGetCorrection(snes,&next);

764:   if (isFine) {
765:     SNESFASCycleSetupPhase_Full(snes);
766:   }

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) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
773:       fas->full_downsweep = PETSC_TRUE;
774:       SNESFASCoarseCorrection(snes,X,F,X);
775:       SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
776:       if (fas->level != 1) next->max_its -= 1;
777:     } else {
778:       /* The smoother on the coarse level is the coarse solver */
779:       SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
780:     }
781:     fas->full_stage = 1;
782:   } else if (fas->full_stage == 1) {
783:     if (snes->iter == 0) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
784:     if (next) {
785:       SNESFASCoarseCorrection(snes,X,F,X);
786:       SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
787:     }
788:   }
789:   /* final v-cycle */
790:   if (isFine) {
791:     if (next) {
792:       SNESFASCoarseCorrection(snes,X,F,X);
793:       SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
794:     }
795:   }
796:   return(0);
797: }

799: static PetscErrorCode SNESFASCycle_Kaskade(SNES snes, Vec X)
800: {
802:   Vec            F,B;
803:   SNES           next;

806:   F = snes->vec_func;
807:   B = snes->vec_rhs;
808:   SNESFASCycleGetCorrection(snes,&next);
809:   if (next) {
810:     SNESFASCoarseCorrection(snes,X,F,X);
811:     SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
812:   } else {
813:     SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
814:   }
815:   return(0);
816: }

818: PetscBool SNEScite = PETSC_FALSE;
819: const char SNESCitation[] = "@techreport{pbmkbsxt2012,\n"
820:                             "  title = {Composing Scalable Nonlinear Algebraic Solvers},\n"
821:                             "  author = {Peter Brune and Mathew Knepley and Barry Smith and Xuemin Tu},\n"
822:                             "  year = 2013,\n"
823:                             "  type = Preprint,\n"
824:                             "  number = {ANL/MCS-P2010-0112},\n"
825:                             "  institution = {Argonne National Laboratory}\n}\n";

827: static PetscErrorCode SNESSolve_FAS(SNES snes)
828: {
830:   PetscInt       i, maxits;
831:   Vec            X, F;
832:   PetscReal      fnorm;
833:   SNES_FAS       *fas = (SNES_FAS*)snes->data,*ffas;
834:   DM             dm;
835:   PetscBool      isFine;


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

841:   PetscCitationsRegister(SNESCitation,&SNEScite);
842:   maxits       = snes->max_its;      /* maximum number of iterations */
843:   snes->reason = SNES_CONVERGED_ITERATING;
844:   X            = snes->vec_sol;
845:   F            = snes->vec_func;

847:   SNESFASCycleIsFine(snes, &isFine);
848:   /*norm setup */
849:   PetscObjectSAWsTakeAccess((PetscObject)snes);
850:   snes->iter = 0;
851:   snes->norm = 0.;
852:   PetscObjectSAWsGrantAccess((PetscObject)snes);
853:   if (!snes->vec_func_init_set) {
854:     if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,snes,0,0,0);}
855:     SNESComputeFunction(snes,X,F);
856:     if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,snes,0,0,0);}
857:   } else snes->vec_func_init_set = PETSC_FALSE;

859:   VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
860:   SNESCheckFunctionNorm(snes,fnorm);
861:   PetscObjectSAWsTakeAccess((PetscObject)snes);
862:   snes->norm = fnorm;
863:   PetscObjectSAWsGrantAccess((PetscObject)snes);
864:   SNESLogConvergenceHistory(snes,fnorm,0);
865:   SNESMonitor(snes,0,fnorm);

867:   /* test convergence */
868:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
869:   if (snes->reason) return(0);


872:   if (isFine) {
873:     /* propagate scale-dependent data up the hierarchy */
874:     SNESGetDM(snes,&dm);
875:     for (ffas=fas; ffas->next; ffas=(SNES_FAS*)ffas->next->data) {
876:       DM dmcoarse;
877:       SNESGetDM(ffas->next,&dmcoarse);
878:       DMRestrict(dm,ffas->restrct,ffas->rscale,ffas->inject,dmcoarse);
879:       dm   = dmcoarse;
880:     }
881:   }

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

886:     if (snes->ops->update) {
887:       (*snes->ops->update)(snes, snes->iter);
888:     }
889:     if (fas->fastype == SNES_FAS_MULTIPLICATIVE) {
890:       SNESFASCycle_Multiplicative(snes, X);
891:     } else if (fas->fastype == SNES_FAS_ADDITIVE) {
892:       SNESFASCycle_Additive(snes, X);
893:     } else if (fas->fastype == SNES_FAS_FULL) {
894:       SNESFASCycle_Full(snes, X);
895:     } else if (fas->fastype ==SNES_FAS_KASKADE) {
896:       SNESFASCycle_Kaskade(snes, X);
897:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported FAS type");

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

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

921: /*MC

923: SNESFAS - Full Approximation Scheme nonlinear multigrid solver.

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

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

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

949: Level: beginner

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

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

958: PETSC_EXTERN PetscErrorCode SNESCreate_FAS(SNES snes)
959: {
960:   SNES_FAS       *fas;

964:   snes->ops->destroy        = SNESDestroy_FAS;
965:   snes->ops->setup          = SNESSetUp_FAS;
966:   snes->ops->setfromoptions = SNESSetFromOptions_FAS;
967:   snes->ops->view           = SNESView_FAS;
968:   snes->ops->solve          = SNESSolve_FAS;
969:   snes->ops->reset          = SNESReset_FAS;

971:   snes->usesksp = PETSC_FALSE;
972:   snes->usesnpc = PETSC_FALSE;

974:   if (!snes->tolerancesset) {
975:     snes->max_funcs = 30000;
976:     snes->max_its   = 10000;
977:   }

979:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

981:   PetscNewLog(snes,&fas);

983:   snes->data                  = (void*) fas;
984:   fas->level                  = 0;
985:   fas->levels                 = 1;
986:   fas->n_cycles               = 1;
987:   fas->max_up_it              = 1;
988:   fas->max_down_it            = 1;
989:   fas->smoothu                = NULL;
990:   fas->smoothd                = NULL;
991:   fas->next                   = NULL;
992:   fas->previous               = NULL;
993:   fas->fine                   = snes;
994:   fas->interpolate            = NULL;
995:   fas->restrct                = NULL;
996:   fas->inject                 = NULL;
997:   fas->usedmfornumberoflevels = PETSC_FALSE;
998:   fas->fastype                = SNES_FAS_MULTIPLICATIVE;
999:   fas->full_downsweep         = PETSC_FALSE;

1001:   fas->eventsmoothsetup    = 0;
1002:   fas->eventsmoothsolve    = 0;
1003:   fas->eventresidual       = 0;
1004:   fas->eventinterprestrict = 0;
1005:   return(0);
1006: }