Actual source code: fas.c

petsc-3.15.0 2021-04-05
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",NULL};

  6: static PetscErrorCode SNESReset_FAS(SNES snes)
  7: {
  8:   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:   VecDestroy(&fas->Xg);
 19:   VecDestroy(&fas->Fg);
 20:   if (fas->next) {SNESReset(fas->next);}
 21:   return(0);
 22: }

 24: static PetscErrorCode SNESDestroy_FAS(SNES snes)
 25: {
 26:   SNES_FAS       *fas = (SNES_FAS*)snes->data;

 30:   /* recursively resets and then destroys */
 31:   SNESReset_FAS(snes);
 32:   SNESDestroy(&fas->next);
 33:   PetscFree(fas);
 34:   return(0);
 35: }

 37: static PetscErrorCode SNESFASSetUpLineSearch_Private(SNES snes, SNES smooth)
 38: {
 39:   SNESLineSearch linesearch;
 40:   SNESLineSearch slinesearch;
 41:   void           *lsprectx,*lspostctx;
 42:   PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*);
 43:   PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);

 47:   if (!snes->linesearch) return(0);
 48:   SNESGetLineSearch(snes,&linesearch);
 49:   SNESGetLineSearch(smooth,&slinesearch);
 50:   SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
 51:   SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
 52:   SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);
 53:   SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);
 54:   PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
 55:   return(0);
 56: }

 58: static PetscErrorCode SNESFASCycleSetUpSmoother_Private(SNES snes, SNES smooth)
 59: {
 60:   SNES_FAS       *fas = (SNES_FAS*) snes->data;

 64:   PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)smooth);
 65:   SNESSetFromOptions(smooth);
 66:   SNESFASSetUpLineSearch_Private(snes, smooth);

 68:   PetscObjectReference((PetscObject)snes->vec_sol);
 69:   PetscObjectReference((PetscObject)snes->vec_sol_update);
 70:   PetscObjectReference((PetscObject)snes->vec_func);
 71:   smooth->vec_sol        = snes->vec_sol;
 72:   smooth->vec_sol_update = snes->vec_sol_update;
 73:   smooth->vec_func       = snes->vec_func;

 75:   if (fas->eventsmoothsetup) {PetscLogEventBegin(fas->eventsmoothsetup,smooth,0,0,0);}
 76:   SNESSetUp(smooth);
 77:   if (fas->eventsmoothsetup) {PetscLogEventEnd(fas->eventsmoothsetup,smooth,0,0,0);}
 78:   return(0);
 79: }

 81: static PetscErrorCode SNESSetUp_FAS(SNES snes)
 82: {
 83:   SNES_FAS       *fas = (SNES_FAS*) snes->data;
 85:   PetscInt       dm_levels;
 86:   SNES           next;
 87:   PetscBool      isFine, hasCreateRestriction, hasCreateInjection;

 90:   SNESFASCycleIsFine(snes, &isFine);
 91:   if (fas->usedmfornumberoflevels && isFine) {
 92:     DMGetRefineLevel(snes->dm,&dm_levels);
 93:     dm_levels++;
 94:     if (dm_levels > fas->levels) {
 95:       /* reset the number of levels */
 96:       SNESFASSetLevels(snes,dm_levels,NULL);
 97:       SNESSetFromOptions(snes);
 98:     }
 99:   }
100:   SNESFASCycleGetCorrection(snes, &next);
101:   if (!isFine) snes->gridsequence = 0; /* no grid sequencing inside the multigrid hierarchy! */

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

105:   /* set up the smoothers if they haven't already been set up */
106:   if (!fas->smoothd) {
107:     SNESFASCycleCreateSmoother_Private(snes, &fas->smoothd);
108:   }

110:   if (snes->dm) {
111:     /* set the smoother DMs properly */
112:     if (fas->smoothu) {SNESSetDM(fas->smoothu, snes->dm);}
113:     SNESSetDM(fas->smoothd, snes->dm);
114:     /* construct EVERYTHING from the DM -- including the progressive set of smoothers */
115:     if (next) {
116:       /* for now -- assume the DM and the evaluation functions have been set externally */
117:       if (!next->dm) {
118:         DMCoarsen(snes->dm, PetscObjectComm((PetscObject)next), &next->dm);
119:         SNESSetDM(next, next->dm);
120:       }
121:       /* set the interpolation and restriction from the DM */
122:       if (!fas->interpolate) {
123:         DMCreateInterpolation(next->dm, snes->dm, &fas->interpolate, &fas->rscale);
124:         if (!fas->restrct) {
125:           DMHasCreateRestriction(next->dm, &hasCreateRestriction);
126:           /* DM can create restrictions, use that */
127:           if (hasCreateRestriction) {
128:             DMCreateRestriction(next->dm, snes->dm, &fas->restrct);
129:           } else {
130:             PetscObjectReference((PetscObject)fas->interpolate);
131:             fas->restrct = fas->interpolate;
132:           }
133:         }
134:       }
135:       /* set the injection from the DM */
136:       if (!fas->inject) {
137:         DMHasCreateInjection(next->dm, &hasCreateInjection);
138:         if (hasCreateInjection) {
139:           DMCreateInjection(next->dm, snes->dm, &fas->inject);
140:         }
141:       }
142:     }
143:   }

145:   /*pass the smoother, function, and jacobian up to the next level if it's not user set already */
146:   if (fas->galerkin) {
147:     if (next) {
148:       SNESSetFunction(next, NULL, SNESFASGalerkinFunctionDefault, next);
149:     }
150:     if (fas->smoothd && fas->level != fas->levels - 1) {
151:       SNESSetFunction(fas->smoothd, NULL, SNESFASGalerkinFunctionDefault, snes);
152:     }
153:     if (fas->smoothu && fas->level != fas->levels - 1) {
154:       SNESSetFunction(fas->smoothu, NULL, SNESFASGalerkinFunctionDefault, snes);
155:     }
156:   }

158:   /* sets the down (pre) smoother's default norm and sets it from options */
159:   if (fas->smoothd) {
160:     if (fas->level == 0 && fas->levels != 1) {
161:       SNESSetNormSchedule(fas->smoothd, SNES_NORM_NONE);
162:     } else {
163:       SNESSetNormSchedule(fas->smoothd, SNES_NORM_FINAL_ONLY);
164:     }
165:     SNESFASCycleSetUpSmoother_Private(snes, fas->smoothd);
166:   }

168:   /* sets the up (post) smoother's default norm and sets it from options */
169:   if (fas->smoothu) {
170:     if (fas->level != fas->levels - 1) {
171:       SNESSetNormSchedule(fas->smoothu, SNES_NORM_NONE);
172:     } else {
173:       SNESSetNormSchedule(fas->smoothu, SNES_NORM_FINAL_ONLY);
174:     }
175:     SNESFASCycleSetUpSmoother_Private(snes, fas->smoothu);
176:   }

178:   if (next) {
179:     /* gotta set up the solution vector for this to work */
180:     if (!next->vec_sol) {SNESFASCreateCoarseVec(snes,&next->vec_sol);}
181:     if (!next->vec_rhs) {SNESFASCreateCoarseVec(snes,&next->vec_rhs);}
182:     PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next);
183:     SNESFASSetUpLineSearch_Private(snes, next);
184:     SNESSetUp(next);
185:   }

187:   /* setup FAS work vectors */
188:   if (fas->galerkin) {
189:     VecDuplicate(snes->vec_sol, &fas->Xg);
190:     VecDuplicate(snes->vec_sol, &fas->Fg);
191:   }
192:   return(0);
193: }

195: static PetscErrorCode SNESSetFromOptions_FAS(PetscOptionItems *PetscOptionsObject,SNES snes)
196: {
197:   SNES_FAS       *fas   = (SNES_FAS*) snes->data;
198:   PetscInt       levels = 1;
199:   PetscBool      flg    = PETSC_FALSE, upflg = PETSC_FALSE, downflg = PETSC_FALSE, monflg = PETSC_FALSE, galerkinflg = PETSC_FALSE,continuationflg = PETSC_FALSE;
201:   SNESFASType    fastype;
202:   const char     *optionsprefix;
203:   SNESLineSearch linesearch;
204:   PetscInt       m, n_up, n_down;
205:   SNES           next;
206:   PetscBool      isFine;

209:   SNESFASCycleIsFine(snes, &isFine);
210:   PetscOptionsHead(PetscOptionsObject,"SNESFAS Options-----------------------------------");

212:   /* number of levels -- only process most options on the finest level */
213:   if (isFine) {
214:     PetscOptionsInt("-snes_fas_levels", "Number of Levels", "SNESFASSetLevels", levels, &levels, &flg);
215:     if (!flg && snes->dm) {
216:       DMGetRefineLevel(snes->dm,&levels);
217:       levels++;
218:       fas->usedmfornumberoflevels = PETSC_TRUE;
219:     }
220:     SNESFASSetLevels(snes, levels, NULL);
221:     fastype = fas->fastype;
222:     PetscOptionsEnum("-snes_fas_type","FAS correction type","SNESFASSetType",SNESFASTypes,(PetscEnum)fastype,(PetscEnum*)&fastype,&flg);
223:     if (flg) {
224:       SNESFASSetType(snes, fastype);
225:     }

227:     SNESGetOptionsPrefix(snes, &optionsprefix);
228:     PetscOptionsInt("-snes_fas_cycles","Number of cycles","SNESFASSetCycles",fas->n_cycles,&m,&flg);
229:     if (flg) {
230:       SNESFASSetCycles(snes, m);
231:     }
232:     PetscOptionsBool("-snes_fas_continuation","Corrected grid-sequence continuation","SNESFASSetContinuation",fas->continuation,&continuationflg,&flg);
233:     if (flg) {
234:       SNESFASSetContinuation(snes,continuationflg);
235:     }

237:     PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin","SNESFASSetGalerkin",fas->galerkin,&galerkinflg,&flg);
238:     if (flg) {
239:       SNESFASSetGalerkin(snes, galerkinflg);
240:     }

242:     if (fas->fastype == SNES_FAS_FULL) {
243:       PetscOptionsBool("-snes_fas_full_downsweep","Smooth on the initial down sweep for full FAS cycles","SNESFASFullSetDownSweep",fas->full_downsweep,&fas->full_downsweep,&flg);
244:       if (flg) {SNESFASFullSetDownSweep(snes,fas->full_downsweep);}
245:       PetscOptionsBool("-snes_fas_full_total","Use total restriction and interpolaton on the indial down and up sweeps for the full FAS cycle","SNESFASFullSetUseTotal",fas->full_total,&fas->full_total,&flg);
246:       if (flg) {SNESFASFullSetTotal(snes,fas->full_total);}
247:     }

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

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

253:     {
254:       PetscViewer viewer;
255:       PetscViewerFormat format;
256:       PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,"-snes_fas_monitor",&viewer,&format,&monflg);
257:       if (monflg) {
258:         PetscViewerAndFormat *vf;
259:         PetscViewerAndFormatCreate(viewer,format,&vf);
260:         PetscObjectDereference((PetscObject)viewer);
261:         SNESFASSetMonitor(snes,vf,PETSC_TRUE);
262:       }
263:     }
264:     flg    = PETSC_FALSE;
265:     monflg = PETSC_TRUE;
266:     PetscOptionsBool("-snes_fas_log","Log times for each FAS level","SNESFASSetLog",monflg,&monflg,&flg);
267:     if (flg) {SNESFASSetLog(snes,monflg);}
268:   }

270:   PetscOptionsTail();

272:   /* setup from the determined types if there is no pointwise procedure or smoother defined */
273:   if (upflg) {
274:     SNESFASSetNumberSmoothUp(snes,n_up);
275:   }
276:   if (downflg) {
277:     SNESFASSetNumberSmoothDown(snes,n_down);
278:   }

280:   /* set up the default line search for coarse grid corrections */
281:   if (fas->fastype == SNES_FAS_ADDITIVE) {
282:     if (!snes->linesearch) {
283:       SNESGetLineSearch(snes, &linesearch);
284:       SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
285:     }
286:   }

288:   /* recursive option setting for the smoothers */
289:   SNESFASCycleGetCorrection(snes, &next);
290:   if (next) {SNESSetFromOptions(next);}
291:   return(0);
292: }

294: #include <petscdraw.h>
295: static PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer)
296: {
297:   SNES_FAS       *fas = (SNES_FAS*) snes->data;
298:   PetscBool      isFine,iascii,isdraw;
299:   PetscInt       i;
301:   SNES           smoothu, smoothd, levelsnes;

304:   SNESFASCycleIsFine(snes, &isFine);
305:   if (isFine) {
306:     PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
307:     PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
308:     if (iascii) {
309:       PetscViewerASCIIPrintf(viewer, "  type is %s, levels=%D, cycles=%D\n",  SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles);
310:       if (fas->galerkin) {
311:         PetscViewerASCIIPrintf(viewer,"  Using Galerkin computed coarse grid function evaluation\n");
312:       } else {
313:         PetscViewerASCIIPrintf(viewer,"  Not using Galerkin computed coarse grid function evaluation\n");
314:       }
315:       for (i=0; i<fas->levels; i++) {
316:         SNESFASGetCycleSNES(snes, i, &levelsnes);
317:         SNESFASCycleGetSmootherUp(levelsnes, &smoothu);
318:         SNESFASCycleGetSmootherDown(levelsnes, &smoothd);
319:         if (!i) {
320:           PetscViewerASCIIPrintf(viewer,"  Coarse grid solver -- level %D -------------------------------\n",i);
321:         } else {
322:           PetscViewerASCIIPrintf(viewer,"  Down solver (pre-smoother) on level %D -------------------------------\n",i);
323:         }
324:         PetscViewerASCIIPushTab(viewer);
325:         if (smoothd) {
326:           SNESView(smoothd,viewer);
327:         } else {
328:           PetscViewerASCIIPrintf(viewer,"Not yet available\n");
329:         }
330:         PetscViewerASCIIPopTab(viewer);
331:         if (i && (smoothd == smoothu)) {
332:           PetscViewerASCIIPrintf(viewer,"  Up solver (post-smoother) same as down solver (pre-smoother)\n");
333:         } else if (i) {
334:           PetscViewerASCIIPrintf(viewer,"  Up solver (post-smoother) on level %D -------------------------------\n",i);
335:           PetscViewerASCIIPushTab(viewer);
336:           if (smoothu) {
337:             SNESView(smoothu,viewer);
338:           } else {
339:             PetscViewerASCIIPrintf(viewer,"Not yet available\n");
340:           }
341:           PetscViewerASCIIPopTab(viewer);
342:         }
343:       }
344:     } else if (isdraw) {
345:       PetscDraw draw;
346:       PetscReal x,w,y,bottom,th,wth;
347:       SNES_FAS  *curfas = fas;
348:       PetscViewerDrawGetDraw(viewer,0,&draw);
349:       PetscDrawGetCurrentPoint(draw,&x,&y);
350:       PetscDrawStringGetSize(draw,&wth,&th);
351:       bottom = y - th;
352:       while (curfas) {
353:         if (!curfas->smoothu) {
354:           PetscDrawPushCurrentPoint(draw,x,bottom);
355:           if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
356:           PetscDrawPopCurrentPoint(draw);
357:         } else {
358:           w    = 0.5*PetscMin(1.0-x,x);
359:           PetscDrawPushCurrentPoint(draw,x-w,bottom);
360:           if (curfas->smoothd) SNESView(curfas->smoothd,viewer);
361:           PetscDrawPopCurrentPoint(draw);
362:           PetscDrawPushCurrentPoint(draw,x+w,bottom);
363:           if (curfas->smoothu) SNESView(curfas->smoothu,viewer);
364:           PetscDrawPopCurrentPoint(draw);
365:         }
366:         /* this is totally bogus but we have no way of knowing how low the previous one was draw to */
367:         bottom -= 5*th;
368:         if (curfas->next) curfas = (SNES_FAS*)curfas->next->data;
369:         else curfas = NULL;
370:       }
371:     }
372:   }
373:   return(0);
374: }

376: /*
377: Defines the action of the downsmoother
378:  */
379: static PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
380: {
381:   PetscErrorCode      ierr;
382:   SNESConvergedReason reason;
383:   Vec                 FPC;
384:   SNES                smoothd;
385:   PetscBool           flg;
386:   SNES_FAS            *fas = (SNES_FAS*) snes->data;

389:   SNESFASCycleGetSmootherDown(snes, &smoothd);
390:   SNESSetInitialFunction(smoothd, F);
391:   if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,smoothd,B,X,0);}
392:   SNESSolve(smoothd, B, X);
393:   if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,smoothd,B,X,0);}
394:   /* check convergence reason for the smoother */
395:   SNESGetConvergedReason(smoothd,&reason);
396:   if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
397:     snes->reason = SNES_DIVERGED_INNER;
398:     return(0);
399:   }

401:   SNESGetFunction(smoothd, &FPC, NULL, NULL);
402:   SNESGetAlwaysComputesFinalResidual(smoothd, &flg);
403:   if (!flg) {
404:     SNESComputeFunction(smoothd, X, FPC);
405:   }
406:   VecCopy(FPC, F);
407:   if (fnorm) {VecNorm(F,NORM_2,fnorm);}
408:   return(0);
409: }


412: /*
413: Defines the action of the upsmoother
414:  */
415: static PetscErrorCode SNESFASUpSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
416: {
417:   PetscErrorCode      ierr;
418:   SNESConvergedReason reason;
419:   Vec                 FPC;
420:   SNES                smoothu;
421:   PetscBool           flg;
422:   SNES_FAS            *fas = (SNES_FAS*) snes->data;

425:   SNESFASCycleGetSmootherUp(snes, &smoothu);
426:   if (fas->eventsmoothsolve) {PetscLogEventBegin(fas->eventsmoothsolve,smoothu,0,0,0);}
427:   SNESSolve(smoothu, B, X);
428:   if (fas->eventsmoothsolve) {PetscLogEventEnd(fas->eventsmoothsolve,smoothu,0,0,0);}
429:   /* check convergence reason for the smoother */
430:   SNESGetConvergedReason(smoothu,&reason);
431:   if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
432:     snes->reason = SNES_DIVERGED_INNER;
433:     return(0);
434:   }
435:   SNESGetFunction(smoothu, &FPC, NULL, NULL);
436:   SNESGetAlwaysComputesFinalResidual(smoothu, &flg);
437:   if (!flg) {
438:     SNESComputeFunction(smoothu, X, FPC);
439:   }
440:   VecCopy(FPC, F);
441:   if (fnorm) {VecNorm(F,NORM_2,fnorm);}
442:   return(0);
443: }

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

448:    Collective

450:    Input Arguments:
451: .  snes - SNESFAS

453:    Output Arguments:
454: .  Xcoarse - vector on level one coarser than snes

456:    Level: developer

458: .seealso: SNESFASSetRestriction(), SNESFASRestrict()
459: @*/
460: PetscErrorCode SNESFASCreateCoarseVec(SNES snes,Vec *Xcoarse)
461: {
463:   SNES_FAS       *fas;

468:   fas = (SNES_FAS*)snes->data;
469:   if (fas->rscale) {
470:     VecDuplicate(fas->rscale,Xcoarse);
471:   } else if (fas->inject) {
472:     MatCreateVecs(fas->inject,Xcoarse,NULL);
473:   } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or injection");
474:   return(0);
475: }

477: /*@
478:    SNESFASRestrict - restrict a Vec to the next coarser level

480:    Collective

482:    Input Arguments:
483: +  fine - SNES from which to restrict
484: -  Xfine - vector to restrict

486:    Output Arguments:
487: .  Xcoarse - result of restriction

489:    Level: developer

491: .seealso: SNESFASSetRestriction(), SNESFASSetInjection()
492: @*/
493: PetscErrorCode SNESFASRestrict(SNES fine,Vec Xfine,Vec Xcoarse)
494: {
496:   SNES_FAS       *fas;

502:   fas = (SNES_FAS*)fine->data;
503:   if (fas->inject) {
504:     MatRestrict(fas->inject,Xfine,Xcoarse);
505:   } else {
506:     MatRestrict(fas->restrct,Xfine,Xcoarse);
507:     VecPointwiseMult(Xcoarse,fas->rscale,Xcoarse);
508:   }
509:   return(0);
510: }

512: /*

514: Performs a variant of FAS using the interpolated total coarse solution

516: fine problem:   F(x) = b
517: coarse problem: F^c(x^c) = Rb, Initial guess Rx
518: interpolated solution: x^f = I x^c (total solution interpolation

520:  */
521: static PetscErrorCode SNESFASInterpolatedCoarseSolution(SNES snes, Vec X, Vec X_new)
522: {
523:   PetscErrorCode      ierr;
524:   Vec                 X_c, B_c;
525:   SNESConvergedReason reason;
526:   SNES                next;
527:   Mat                 restrct, interpolate;
528:   SNES_FAS            *fasc;

531:   SNESFASCycleGetCorrection(snes, &next);
532:   if (next) {
533:     fasc = (SNES_FAS*)next->data;

535:     SNESFASCycleGetRestriction(snes, &restrct);
536:     SNESFASCycleGetInterpolation(snes, &interpolate);

538:     X_c  = next->vec_sol;

540:     if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,snes,0,0,0);}
541:     /* restrict the total solution: Rb */
542:     SNESFASRestrict(snes, X, X_c);
543:     B_c = next->vec_rhs;
544:     if (snes->vec_rhs) {
545:       /* restrict the total rhs defect: Rb */
546:       MatRestrict(restrct, snes->vec_rhs, B_c);
547:     } else {
548:       VecSet(B_c, 0.);
549:     }
550:     if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,snes,0,0,0);}

552:     SNESSolve(next, B_c, X_c);
553:     SNESGetConvergedReason(next,&reason);
554:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
555:       snes->reason = SNES_DIVERGED_INNER;
556:       return(0);
557:     }
558:     /* x^f <- Ix^c*/
559:     DM dmc,dmf;

561:     SNESGetDM(next, &dmc);
562:     SNESGetDM(snes, &dmf);
563:     if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,snes,0,0,0);}
564:     DMInterpolateSolution(dmc, dmf, interpolate, X_c, X_new);
565:     if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,snes,0,0,0);}
566:     PetscObjectSetName((PetscObject) X_c, "Coarse solution");
567:     VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view");
568:     PetscObjectSetName((PetscObject) X_new, "Updated Fine solution");
569:     VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view");
570:   }
571:   return(0);
572: }

574: /*

576: Performs the FAS coarse correction as:

578: fine problem:   F(x) = b
579: coarse problem: F^c(x^c) = b^c

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

583:  */
584: PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new)
585: {
586:   PetscErrorCode      ierr;
587:   Vec                 X_c, Xo_c, F_c, B_c;
588:   SNESConvergedReason reason;
589:   SNES                next;
590:   Mat                 restrct, interpolate;
591:   SNES_FAS            *fasc;

594:   SNESFASCycleGetCorrection(snes, &next);
595:   if (next) {
596:     fasc = (SNES_FAS*)next->data;

598:     SNESFASCycleGetRestriction(snes, &restrct);
599:     SNESFASCycleGetInterpolation(snes, &interpolate);

601:     X_c  = next->vec_sol;
602:     Xo_c = next->work[0];
603:     F_c  = next->vec_func;
604:     B_c  = next->vec_rhs;

606:     if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,snes,0,0,0);}
607:     SNESFASRestrict(snes, X, Xo_c);
608:     /* restrict the defect: R(F(x) - b) */
609:     MatRestrict(restrct, F, B_c);
610:     if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,snes,0,0,0);}

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

617:     /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */
618:     VecCopy(B_c, X_c);
619:     VecCopy(F_c, B_c);
620:     VecCopy(X_c, F_c);
621:     /* set initial guess of the coarse problem to the projected fine solution */
622:     VecCopy(Xo_c, X_c);

624:     /* recurse to the next level */
625:     SNESSetInitialFunction(next, F_c);
626:     SNESSolve(next, B_c, X_c);
627:     SNESGetConvergedReason(next,&reason);
628:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
629:       snes->reason = SNES_DIVERGED_INNER;
630:       return(0);
631:     }
632:     /* correct as x <- x + I(x^c - Rx)*/
633:     VecAXPY(X_c, -1.0, Xo_c);

635:     if (fasc->eventinterprestrict) {PetscLogEventBegin(fasc->eventinterprestrict,snes,0,0,0);}
636:     MatInterpolateAdd(interpolate, X_c, X, X_new);
637:     if (fasc->eventinterprestrict) {PetscLogEventEnd(fasc->eventinterprestrict,snes,0,0,0);}
638:     PetscObjectSetName((PetscObject) X_c, "Coarse correction");
639:     VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view");
640:     PetscObjectSetName((PetscObject) X_new, "Updated Fine solution");
641:     VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view");
642:   }
643:   return(0);
644: }

646: /*

648: The additive cycle looks like:

650: xhat = x
651: xhat = dS(x, b)
652: x = coarsecorrection(xhat, b_d)
653: x = x + nu*(xhat - x);
654: (optional) x = uS(x, b)

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

658:  */
659: static PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X)
660: {
661:   Vec                  F, B, Xhat;
662:   Vec                  X_c, Xo_c, F_c, B_c;
663:   PetscErrorCode       ierr;
664:   SNESConvergedReason  reason;
665:   PetscReal            xnorm, fnorm, ynorm;
666:   SNESLineSearchReason lsresult;
667:   SNES                 next;
668:   Mat                  restrct, interpolate;
669:   SNES_FAS             *fas = (SNES_FAS*)snes->data,*fasc;

672:   SNESFASCycleGetCorrection(snes, &next);
673:   F    = snes->vec_func;
674:   B    = snes->vec_rhs;
675:   Xhat = snes->work[1];
676:   VecCopy(X, Xhat);
677:   /* recurse first */
678:   if (next) {
679:     fasc = (SNES_FAS*)next->data;
680:     SNESFASCycleGetRestriction(snes, &restrct);
681:     SNESFASCycleGetInterpolation(snes, &interpolate);
682:     if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,snes,0,0,0);}
683:     SNESComputeFunction(snes, Xhat, F);
684:     if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,snes,0,0,0);}
685:     VecNorm(F, NORM_2, &fnorm);
686:     X_c  = next->vec_sol;
687:     Xo_c = next->work[0];
688:     F_c  = next->vec_func;
689:     B_c  = next->vec_rhs;

691:     SNESFASRestrict(snes,Xhat,Xo_c);
692:     /* restrict the defect */
693:     MatRestrict(restrct, F, B_c);

695:     /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */
696:     if (fasc->eventresidual) {PetscLogEventBegin(fasc->eventresidual,next,0,0,0);}
697:     SNESComputeFunction(next, Xo_c, F_c);
698:     if (fasc->eventresidual) {PetscLogEventEnd(fasc->eventresidual,next,0,0,0);}
699:     VecCopy(B_c, X_c);
700:     VecCopy(F_c, B_c);
701:     VecCopy(X_c, F_c);
702:     /* set initial guess of the coarse problem to the projected fine solution */
703:     VecCopy(Xo_c, X_c);

705:     /* recurse */
706:     SNESSetInitialFunction(next, F_c);
707:     SNESSolve(next, B_c, X_c);

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

712:     SNESGetConvergedReason(next,&reason);
713:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
714:       snes->reason = SNES_DIVERGED_INNER;
715:       return(0);
716:     }

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

722:     /* additive correction of the coarse direction*/
723:     SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat);
724:     SNESLineSearchGetReason(snes->linesearch, &lsresult);
725:     SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm);
726:     if (lsresult) {
727:       if (++snes->numFailures >= snes->maxFailures) {
728:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
729:         return(0);
730:       }
731:     }
732:   } else {
733:     SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
734:   }
735:   return(0);
736: }

738: /*

740: Defines the FAS cycle as:

742: fine problem: F(x) = b
743: coarse problem: F^c(x) = b^c

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

747: correction:

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

751:  */
752: static PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X)
753: {

756:   Vec            F,B;
757:   SNES           next;

760:   F = snes->vec_func;
761:   B = snes->vec_rhs;
762:   /* pre-smooth -- just update using the pre-smoother */
763:   SNESFASCycleGetCorrection(snes, &next);
764:   SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
765:   if (next) {
766:     SNESFASCoarseCorrection(snes, X, F, X);
767:     SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
768:   }
769:   return(0);
770: }

772: static PetscErrorCode SNESFASCycleSetupPhase_Full(SNES snes)
773: {
774:   SNES           next;
775:   SNES_FAS       *fas = (SNES_FAS*)snes->data;
776:   PetscBool      isFine;

780:   /* pre-smooth -- just update using the pre-smoother */
781:   SNESFASCycleIsFine(snes,&isFine);
782:   SNESFASCycleGetCorrection(snes,&next);
783:   fas->full_stage = 0;
784:   if (next) {SNESFASCycleSetupPhase_Full(next);}
785:   return(0);
786: }

788: static PetscErrorCode SNESFASCycle_Full(SNES snes, Vec X)
789: {
791:   Vec            F,B;
792:   SNES_FAS       *fas = (SNES_FAS*)snes->data;
793:   PetscBool      isFine;
794:   SNES           next;

797:   F = snes->vec_func;
798:   B = snes->vec_rhs;
799:   SNESFASCycleIsFine(snes,&isFine);
800:   SNESFASCycleGetCorrection(snes,&next);

802:   if (isFine) {
803:     SNESFASCycleSetupPhase_Full(snes);
804:   }

806:   if (fas->full_stage == 0) {
807:     /* downsweep */
808:     if (next) {
809:       if (fas->level != 1) next->max_its += 1;
810:       if (fas->full_downsweep) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
811:       fas->full_downsweep = PETSC_TRUE;
812:       if (fas->full_total) {SNESFASInterpolatedCoarseSolution(snes,X,X);}
813:       else                 {SNESFASCoarseCorrection(snes,X,F,X);}
814:       fas->full_total = PETSC_FALSE;
815:       SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
816:       if (fas->level != 1) next->max_its -= 1;
817:     } else {
818:       /* The smoother on the coarse level is the coarse solver */
819:       SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
820:     }
821:     fas->full_stage = 1;
822:   } else if (fas->full_stage == 1) {
823:     if (snes->iter == 0) {SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);}
824:     if (next) {
825:       SNESFASCoarseCorrection(snes,X,F,X);
826:       SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
827:     }
828:   }
829:   /* final v-cycle */
830:   if (isFine) {
831:     if (next) {
832:       SNESFASCoarseCorrection(snes,X,F,X);
833:       SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
834:     }
835:   }
836:   return(0);
837: }

839: static PetscErrorCode SNESFASCycle_Kaskade(SNES snes, Vec X)
840: {
842:   Vec            F,B;
843:   SNES           next;

846:   F = snes->vec_func;
847:   B = snes->vec_rhs;
848:   SNESFASCycleGetCorrection(snes,&next);
849:   if (next) {
850:     SNESFASCoarseCorrection(snes,X,F,X);
851:     SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);
852:   } else {
853:     SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);
854:   }
855:   return(0);
856: }

858: PetscBool SNEScite = PETSC_FALSE;
859: const char SNESCitation[] = "@techreport{pbmkbsxt2012,\n"
860:                             "  title = {Composing Scalable Nonlinear Algebraic Solvers},\n"
861:                             "  author = {Peter Brune and Mathew Knepley and Barry Smith and Xuemin Tu},\n"
862:                             "  year = 2013,\n"
863:                             "  type = Preprint,\n"
864:                             "  number = {ANL/MCS-P2010-0112},\n"
865:                             "  institution = {Argonne National Laboratory}\n}\n";

867: static PetscErrorCode SNESSolve_FAS(SNES snes)
868: {
870:   PetscInt       i;
871:   Vec            X, F;
872:   PetscReal      fnorm;
873:   SNES_FAS       *fas = (SNES_FAS*)snes->data,*ffas;
874:   DM             dm;
875:   PetscBool      isFine;

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

880:   PetscCitationsRegister(SNESCitation,&SNEScite);
881:   snes->reason = SNES_CONVERGED_ITERATING;
882:   X            = snes->vec_sol;
883:   F            = snes->vec_func;

885:   SNESFASCycleIsFine(snes, &isFine);
886:   /* norm setup */
887:   PetscObjectSAWsTakeAccess((PetscObject)snes);
888:   snes->iter = 0;
889:   snes->norm = 0;
890:   PetscObjectSAWsGrantAccess((PetscObject)snes);
891:   if (!snes->vec_func_init_set) {
892:     if (fas->eventresidual) {PetscLogEventBegin(fas->eventresidual,snes,0,0,0);}
893:     SNESComputeFunction(snes,X,F);
894:     if (fas->eventresidual) {PetscLogEventEnd(fas->eventresidual,snes,0,0,0);}
895:   } else snes->vec_func_init_set = PETSC_FALSE;

897:   VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
898:   SNESCheckFunctionNorm(snes,fnorm);
899:   PetscObjectSAWsTakeAccess((PetscObject)snes);
900:   snes->norm = fnorm;
901:   PetscObjectSAWsGrantAccess((PetscObject)snes);
902:   SNESLogConvergenceHistory(snes,fnorm,0);
903:   SNESMonitor(snes,snes->iter,fnorm);

905:   /* test convergence */
906:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
907:   if (snes->reason) return(0);


910:   if (isFine) {
911:     /* propagate scale-dependent data up the hierarchy */
912:     SNESGetDM(snes,&dm);
913:     for (ffas=fas; ffas->next; ffas=(SNES_FAS*)ffas->next->data) {
914:       DM dmcoarse;
915:       SNESGetDM(ffas->next,&dmcoarse);
916:       DMRestrict(dm,ffas->restrct,ffas->rscale,ffas->inject,dmcoarse);
917:       dm   = dmcoarse;
918:     }
919:   }

921:   for (i = 0; i < snes->max_its; i++) {
922:     /* Call general purpose update function */
923:     if (snes->ops->update) {
924:       (*snes->ops->update)(snes, snes->iter);
925:     }

927:     if (fas->fastype == SNES_FAS_MULTIPLICATIVE) {
928:       SNESFASCycle_Multiplicative(snes, X);
929:     } else if (fas->fastype == SNES_FAS_ADDITIVE) {
930:       SNESFASCycle_Additive(snes, X);
931:     } else if (fas->fastype == SNES_FAS_FULL) {
932:       SNESFASCycle_Full(snes, X);
933:     } else if (fas->fastype == SNES_FAS_KASKADE) {
934:       SNESFASCycle_Kaskade(snes, X);
935:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported FAS type");

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

940:     /* Monitor convergence */
941:     PetscObjectSAWsTakeAccess((PetscObject)snes);
942:     snes->iter = i+1;
943:     PetscObjectSAWsGrantAccess((PetscObject)snes);
944:     SNESLogConvergenceHistory(snes,snes->norm,0);
945:     SNESMonitor(snes,snes->iter,snes->norm);
946:     /* Test for convergence */
947:     if (isFine) {
948:       (*snes->ops->converged)(snes,snes->iter,0.0,0.0,snes->norm,&snes->reason,snes->cnvP);
949:       if (snes->reason) break;
950:     }
951:   }
952:   if (i == snes->max_its) {
953:     PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", i);
954:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
955:   }
956:   return(0);
957: }

959: /*MC

961: SNESFAS - Full Approximation Scheme nonlinear multigrid solver.

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

967: Options Database:
968: +   -snes_fas_levels -  The number of levels
969: .   -snes_fas_cycles<1> -  The number of cycles -- 1 for V, 2 for W
970: .   -snes_fas_type<additive,multiplicative,full,kaskade>  -  Additive or multiplicative cycle
971: .   -snes_fas_galerkin<PETSC_FALSE> -  Form coarse problems by projection back upon the fine problem
972: .   -snes_fas_smoothup<1> -  The number of iterations of the post-smoother
973: .   -snes_fas_smoothdown<1> -  The number of iterations of the pre-smoother
974: .   -snes_fas_monitor -  Monitor progress of all of the levels
975: .   -snes_fas_full_downsweep<PETSC_FALSE> - call the downsmooth on the initial downsweep of full FAS
976: .   -fas_levels_snes_ -  SNES options for all smoothers
977: .   -fas_levels_cycle_snes_ -  SNES options for all cycles
978: .   -fas_levels_i_snes_ -  SNES options for the smoothers on level i
979: .   -fas_levels_i_cycle_snes_ - SNES options for the cycle on level i
980: -   -fas_coarse_snes_ -  SNES options for the coarsest smoother

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

987: Level: beginner

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

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

996: PETSC_EXTERN PetscErrorCode SNESCreate_FAS(SNES snes)
997: {
998:   SNES_FAS       *fas;

1002:   snes->ops->destroy        = SNESDestroy_FAS;
1003:   snes->ops->setup          = SNESSetUp_FAS;
1004:   snes->ops->setfromoptions = SNESSetFromOptions_FAS;
1005:   snes->ops->view           = SNESView_FAS;
1006:   snes->ops->solve          = SNESSolve_FAS;
1007:   snes->ops->reset          = SNESReset_FAS;

1009:   snes->usesksp = PETSC_FALSE;
1010:   snes->usesnpc = PETSC_FALSE;

1012:   if (!snes->tolerancesset) {
1013:     snes->max_funcs = 30000;
1014:     snes->max_its   = 10000;
1015:   }

1017:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

1019:   PetscNewLog(snes,&fas);

1021:   snes->data                  = (void*) fas;
1022:   fas->level                  = 0;
1023:   fas->levels                 = 1;
1024:   fas->n_cycles               = 1;
1025:   fas->max_up_it              = 1;
1026:   fas->max_down_it            = 1;
1027:   fas->smoothu                = NULL;
1028:   fas->smoothd                = NULL;
1029:   fas->next                   = NULL;
1030:   fas->previous               = NULL;
1031:   fas->fine                   = snes;
1032:   fas->interpolate            = NULL;
1033:   fas->restrct                = NULL;
1034:   fas->inject                 = NULL;
1035:   fas->usedmfornumberoflevels = PETSC_FALSE;
1036:   fas->fastype                = SNES_FAS_MULTIPLICATIVE;
1037:   fas->full_downsweep         = PETSC_FALSE;
1038:   fas->full_total             = PETSC_FALSE;

1040:   fas->eventsmoothsetup    = 0;
1041:   fas->eventsmoothsolve    = 0;
1042:   fas->eventresidual       = 0;
1043:   fas->eventinterprestrict = 0;
1044:   return(0);
1045: }