Actual source code: fas.c

  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;

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

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

 26:   /* recursively resets and then destroys */
 27:   SNESReset_FAS(snes);
 28:   SNESDestroy(&fas->next);
 29:   PetscFree(fas);
 30:   return 0;
 31: }

 33: static PetscErrorCode SNESFASSetUpLineSearch_Private(SNES snes, SNES smooth)
 34: {
 35:   SNESLineSearch linesearch;
 36:   SNESLineSearch slinesearch;
 37:   void          *lsprectx, *lspostctx;
 38:   PetscErrorCode (*precheck)(SNESLineSearch, Vec, Vec, PetscBool *, void *);
 39:   PetscErrorCode (*postcheck)(SNESLineSearch, Vec, Vec, Vec, PetscBool *, PetscBool *, void *);

 41:   if (!snes->linesearch) return 0;
 42:   SNESGetLineSearch(snes, &linesearch);
 43:   SNESGetLineSearch(smooth, &slinesearch);
 44:   SNESLineSearchGetPreCheck(linesearch, &precheck, &lsprectx);
 45:   SNESLineSearchGetPostCheck(linesearch, &postcheck, &lspostctx);
 46:   SNESLineSearchSetPreCheck(slinesearch, precheck, lsprectx);
 47:   SNESLineSearchSetPostCheck(slinesearch, postcheck, lspostctx);
 48:   PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);
 49:   return 0;
 50: }

 52: static PetscErrorCode SNESFASCycleSetUpSmoother_Private(SNES snes, SNES smooth)
 53: {
 54:   SNES_FAS *fas = (SNES_FAS *)snes->data;

 56:   PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)smooth);
 57:   SNESSetFromOptions(smooth);
 58:   SNESFASSetUpLineSearch_Private(snes, smooth);

 60:   PetscObjectReference((PetscObject)snes->vec_sol);
 61:   PetscObjectReference((PetscObject)snes->vec_sol_update);
 62:   PetscObjectReference((PetscObject)snes->vec_func);
 63:   smooth->vec_sol        = snes->vec_sol;
 64:   smooth->vec_sol_update = snes->vec_sol_update;
 65:   smooth->vec_func       = snes->vec_func;

 67:   if (fas->eventsmoothsetup) PetscLogEventBegin(fas->eventsmoothsetup, smooth, 0, 0, 0);
 68:   SNESSetUp(smooth);
 69:   if (fas->eventsmoothsetup) PetscLogEventEnd(fas->eventsmoothsetup, smooth, 0, 0, 0);
 70:   return 0;
 71: }

 73: static PetscErrorCode SNESSetUp_FAS(SNES snes)
 74: {
 75:   SNES_FAS *fas = (SNES_FAS *)snes->data;
 76:   PetscInt  dm_levels;
 77:   SNES      next;
 78:   PetscBool isFine, hasCreateRestriction, hasCreateInjection;

 80:   SNESFASCycleIsFine(snes, &isFine);
 81:   if (fas->usedmfornumberoflevels && isFine) {
 82:     DMGetRefineLevel(snes->dm, &dm_levels);
 83:     dm_levels++;
 84:     if (dm_levels > fas->levels) {
 85:       /* reset the number of levels */
 86:       SNESFASSetLevels(snes, dm_levels, NULL);
 87:       SNESSetFromOptions(snes);
 88:     }
 89:   }
 90:   SNESFASCycleGetCorrection(snes, &next);
 91:   if (!isFine) snes->gridsequence = 0; /* no grid sequencing inside the multigrid hierarchy! */

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

 95:   /* set up the smoothers if they haven't already been set up */
 96:   if (!fas->smoothd) SNESFASCycleCreateSmoother_Private(snes, &fas->smoothd);

 98:   if (snes->dm) {
 99:     /* set the smoother DMs properly */
100:     if (fas->smoothu) SNESSetDM(fas->smoothu, snes->dm);
101:     SNESSetDM(fas->smoothd, snes->dm);
102:     /* construct EVERYTHING from the DM -- including the progressive set of smoothers */
103:     if (next) {
104:       /* for now -- assume the DM and the evaluation functions have been set externally */
105:       if (!next->dm) {
106:         DMCoarsen(snes->dm, PetscObjectComm((PetscObject)next), &next->dm);
107:         SNESSetDM(next, next->dm);
108:       }
109:       /* set the interpolation and restriction from the DM */
110:       if (!fas->interpolate) {
111:         DMCreateInterpolation(next->dm, snes->dm, &fas->interpolate, &fas->rscale);
112:         if (!fas->restrct) {
113:           DMHasCreateRestriction(next->dm, &hasCreateRestriction);
114:           /* DM can create restrictions, use that */
115:           if (hasCreateRestriction) {
116:             DMCreateRestriction(next->dm, snes->dm, &fas->restrct);
117:           } else {
118:             PetscObjectReference((PetscObject)fas->interpolate);
119:             fas->restrct = fas->interpolate;
120:           }
121:         }
122:       }
123:       /* set the injection from the DM */
124:       if (!fas->inject) {
125:         DMHasCreateInjection(next->dm, &hasCreateInjection);
126:         if (hasCreateInjection) DMCreateInjection(next->dm, snes->dm, &fas->inject);
127:       }
128:     }
129:   }

131:   /*pass the smoother, function, and jacobian up to the next level if it's not user set already */
132:   if (fas->galerkin) {
133:     if (next) SNESSetFunction(next, NULL, SNESFASGalerkinFunctionDefault, next);
134:     if (fas->smoothd && fas->level != fas->levels - 1) SNESSetFunction(fas->smoothd, NULL, SNESFASGalerkinFunctionDefault, snes);
135:     if (fas->smoothu && fas->level != fas->levels - 1) SNESSetFunction(fas->smoothu, NULL, SNESFASGalerkinFunctionDefault, snes);
136:   }

138:   /* sets the down (pre) smoother's default norm and sets it from options */
139:   if (fas->smoothd) {
140:     if (fas->level == 0 && fas->levels != 1) {
141:       SNESSetNormSchedule(fas->smoothd, SNES_NORM_NONE);
142:     } else {
143:       SNESSetNormSchedule(fas->smoothd, SNES_NORM_FINAL_ONLY);
144:     }
145:     SNESFASCycleSetUpSmoother_Private(snes, fas->smoothd);
146:   }

148:   /* sets the up (post) smoother's default norm and sets it from options */
149:   if (fas->smoothu) {
150:     if (fas->level != fas->levels - 1) {
151:       SNESSetNormSchedule(fas->smoothu, SNES_NORM_NONE);
152:     } else {
153:       SNESSetNormSchedule(fas->smoothu, SNES_NORM_FINAL_ONLY);
154:     }
155:     SNESFASCycleSetUpSmoother_Private(snes, fas->smoothu);
156:   }

158:   if (next) {
159:     /* gotta set up the solution vector for this to work */
160:     if (!next->vec_sol) SNESFASCreateCoarseVec(snes, &next->vec_sol);
161:     if (!next->vec_rhs) SNESFASCreateCoarseVec(snes, &next->vec_rhs);
162:     PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next);
163:     SNESFASSetUpLineSearch_Private(snes, next);
164:     SNESSetUp(next);
165:   }

167:   /* setup FAS work vectors */
168:   if (fas->galerkin) {
169:     VecDuplicate(snes->vec_sol, &fas->Xg);
170:     VecDuplicate(snes->vec_sol, &fas->Fg);
171:   }
172:   return 0;
173: }

175: static PetscErrorCode SNESSetFromOptions_FAS(SNES snes, PetscOptionItems *PetscOptionsObject)
176: {
177:   SNES_FAS      *fas    = (SNES_FAS *)snes->data;
178:   PetscInt       levels = 1;
179:   PetscBool      flg = PETSC_FALSE, upflg = PETSC_FALSE, downflg = PETSC_FALSE, monflg = PETSC_FALSE, galerkinflg = PETSC_FALSE, continuationflg = PETSC_FALSE;
180:   SNESFASType    fastype;
181:   const char    *optionsprefix;
182:   SNESLineSearch linesearch;
183:   PetscInt       m, n_up, n_down;
184:   SNES           next;
185:   PetscBool      isFine;

187:   SNESFASCycleIsFine(snes, &isFine);
188:   PetscOptionsHeadBegin(PetscOptionsObject, "SNESFAS Options-----------------------------------");

190:   /* number of levels -- only process most options on the finest level */
191:   if (isFine) {
192:     PetscOptionsInt("-snes_fas_levels", "Number of Levels", "SNESFASSetLevels", levels, &levels, &flg);
193:     if (!flg && snes->dm) {
194:       DMGetRefineLevel(snes->dm, &levels);
195:       levels++;
196:       fas->usedmfornumberoflevels = PETSC_TRUE;
197:     }
198:     SNESFASSetLevels(snes, levels, NULL);
199:     fastype = fas->fastype;
200:     PetscOptionsEnum("-snes_fas_type", "FAS correction type", "SNESFASSetType", SNESFASTypes, (PetscEnum)fastype, (PetscEnum *)&fastype, &flg);
201:     if (flg) SNESFASSetType(snes, fastype);

203:     SNESGetOptionsPrefix(snes, &optionsprefix);
204:     PetscOptionsInt("-snes_fas_cycles", "Number of cycles", "SNESFASSetCycles", fas->n_cycles, &m, &flg);
205:     if (flg) SNESFASSetCycles(snes, m);
206:     PetscOptionsBool("-snes_fas_continuation", "Corrected grid-sequence continuation", "SNESFASSetContinuation", fas->continuation, &continuationflg, &flg);
207:     if (flg) SNESFASSetContinuation(snes, continuationflg);

209:     PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin", "SNESFASSetGalerkin", fas->galerkin, &galerkinflg, &flg);
210:     if (flg) SNESFASSetGalerkin(snes, galerkinflg);

212:     if (fas->fastype == SNES_FAS_FULL) {
213:       PetscOptionsBool("-snes_fas_full_downsweep", "Smooth on the initial down sweep for full FAS cycles", "SNESFASFullSetDownSweep", fas->full_downsweep, &fas->full_downsweep, &flg);
214:       if (flg) SNESFASFullSetDownSweep(snes, fas->full_downsweep);
215:       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);
216:       if (flg) SNESFASFullSetTotal(snes, fas->full_total);
217:     }

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

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

223:     {
224:       PetscViewer       viewer;
225:       PetscViewerFormat format;
226:       PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_fas_monitor", &viewer, &format, &monflg);
227:       if (monflg) {
228:         PetscViewerAndFormat *vf;
229:         PetscViewerAndFormatCreate(viewer, format, &vf);
230:         PetscObjectDereference((PetscObject)viewer);
231:         SNESFASSetMonitor(snes, vf, PETSC_TRUE);
232:       }
233:     }
234:     flg    = PETSC_FALSE;
235:     monflg = PETSC_TRUE;
236:     PetscOptionsBool("-snes_fas_log", "Log times for each FAS level", "SNESFASSetLog", monflg, &monflg, &flg);
237:     if (flg) SNESFASSetLog(snes, monflg);
238:   }

240:   PetscOptionsHeadEnd();

242:   /* setup from the determined types if there is no pointwise procedure or smoother defined */
243:   if (upflg) SNESFASSetNumberSmoothUp(snes, n_up);
244:   if (downflg) SNESFASSetNumberSmoothDown(snes, n_down);

246:   /* set up the default line search for coarse grid corrections */
247:   if (fas->fastype == SNES_FAS_ADDITIVE) {
248:     if (!snes->linesearch) {
249:       SNESGetLineSearch(snes, &linesearch);
250:       SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
251:     }
252:   }

254:   /* recursive option setting for the smoothers */
255:   SNESFASCycleGetCorrection(snes, &next);
256:   if (next) SNESSetFromOptions(next);
257:   return 0;
258: }

260: #include <petscdraw.h>
261: static PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer)
262: {
263:   SNES_FAS *fas = (SNES_FAS *)snes->data;
264:   PetscBool isFine, iascii, isdraw;
265:   PetscInt  i;
266:   SNES      smoothu, smoothd, levelsnes;

268:   SNESFASCycleIsFine(snes, &isFine);
269:   if (isFine) {
270:     PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
271:     PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
272:     if (iascii) {
273:       PetscViewerASCIIPrintf(viewer, "  type is %s, levels=%" PetscInt_FMT ", cycles=%" PetscInt_FMT "\n", SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles);
274:       if (fas->galerkin) {
275:         PetscViewerASCIIPrintf(viewer, "  Using Galerkin computed coarse grid function evaluation\n");
276:       } else {
277:         PetscViewerASCIIPrintf(viewer, "  Not using Galerkin computed coarse grid function evaluation\n");
278:       }
279:       for (i = 0; i < fas->levels; i++) {
280:         SNESFASGetCycleSNES(snes, i, &levelsnes);
281:         SNESFASCycleGetSmootherUp(levelsnes, &smoothu);
282:         SNESFASCycleGetSmootherDown(levelsnes, &smoothd);
283:         if (!i) {
284:           PetscViewerASCIIPrintf(viewer, "  Coarse grid solver -- level %" PetscInt_FMT " -------------------------------\n", i);
285:         } else {
286:           PetscViewerASCIIPrintf(viewer, "  Down solver (pre-smoother) on level %" PetscInt_FMT " -------------------------------\n", i);
287:         }
288:         PetscViewerASCIIPushTab(viewer);
289:         if (smoothd) {
290:           SNESView(smoothd, viewer);
291:         } else {
292:           PetscViewerASCIIPrintf(viewer, "Not yet available\n");
293:         }
294:         PetscViewerASCIIPopTab(viewer);
295:         if (i && (smoothd == smoothu)) {
296:           PetscViewerASCIIPrintf(viewer, "  Up solver (post-smoother) same as down solver (pre-smoother)\n");
297:         } else if (i) {
298:           PetscViewerASCIIPrintf(viewer, "  Up solver (post-smoother) on level %" PetscInt_FMT " -------------------------------\n", i);
299:           PetscViewerASCIIPushTab(viewer);
300:           if (smoothu) {
301:             SNESView(smoothu, viewer);
302:           } else {
303:             PetscViewerASCIIPrintf(viewer, "Not yet available\n");
304:           }
305:           PetscViewerASCIIPopTab(viewer);
306:         }
307:       }
308:     } else if (isdraw) {
309:       PetscDraw draw;
310:       PetscReal x, w, y, bottom, th, wth;
311:       SNES_FAS *curfas = fas;
312:       PetscViewerDrawGetDraw(viewer, 0, &draw);
313:       PetscDrawGetCurrentPoint(draw, &x, &y);
314:       PetscDrawStringGetSize(draw, &wth, &th);
315:       bottom = y - th;
316:       while (curfas) {
317:         if (!curfas->smoothu) {
318:           PetscDrawPushCurrentPoint(draw, x, bottom);
319:           if (curfas->smoothd) SNESView(curfas->smoothd, viewer);
320:           PetscDrawPopCurrentPoint(draw);
321:         } else {
322:           w = 0.5 * PetscMin(1.0 - x, x);
323:           PetscDrawPushCurrentPoint(draw, x - w, bottom);
324:           if (curfas->smoothd) SNESView(curfas->smoothd, viewer);
325:           PetscDrawPopCurrentPoint(draw);
326:           PetscDrawPushCurrentPoint(draw, x + w, bottom);
327:           if (curfas->smoothu) SNESView(curfas->smoothu, viewer);
328:           PetscDrawPopCurrentPoint(draw);
329:         }
330:         /* this is totally bogus but we have no way of knowing how low the previous one was draw to */
331:         bottom -= 5 * th;
332:         if (curfas->next) curfas = (SNES_FAS *)curfas->next->data;
333:         else curfas = NULL;
334:       }
335:     }
336:   }
337:   return 0;
338: }

340: /*
341: Defines the action of the downsmoother
342:  */
343: static PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
344: {
345:   SNESConvergedReason reason;
346:   Vec                 FPC;
347:   SNES                smoothd;
348:   PetscBool           flg;
349:   SNES_FAS           *fas = (SNES_FAS *)snes->data;

351:   SNESFASCycleGetSmootherDown(snes, &smoothd);
352:   SNESSetInitialFunction(smoothd, F);
353:   if (fas->eventsmoothsolve) PetscLogEventBegin(fas->eventsmoothsolve, smoothd, B, X, 0);
354:   SNESSolve(smoothd, B, X);
355:   if (fas->eventsmoothsolve) PetscLogEventEnd(fas->eventsmoothsolve, smoothd, B, X, 0);
356:   /* check convergence reason for the smoother */
357:   SNESGetConvergedReason(smoothd, &reason);
358:   if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
359:     snes->reason = SNES_DIVERGED_INNER;
360:     return 0;
361:   }

363:   SNESGetFunction(smoothd, &FPC, NULL, NULL);
364:   SNESGetAlwaysComputesFinalResidual(smoothd, &flg);
365:   if (!flg) SNESComputeFunction(smoothd, X, FPC);
366:   VecCopy(FPC, F);
367:   if (fnorm) VecNorm(F, NORM_2, fnorm);
368:   return 0;
369: }

371: /*
372: Defines the action of the upsmoother
373:  */
374: static PetscErrorCode SNESFASUpSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
375: {
376:   SNESConvergedReason reason;
377:   Vec                 FPC;
378:   SNES                smoothu;
379:   PetscBool           flg;
380:   SNES_FAS           *fas = (SNES_FAS *)snes->data;

382:   SNESFASCycleGetSmootherUp(snes, &smoothu);
383:   if (fas->eventsmoothsolve) PetscLogEventBegin(fas->eventsmoothsolve, smoothu, 0, 0, 0);
384:   SNESSolve(smoothu, B, X);
385:   if (fas->eventsmoothsolve) PetscLogEventEnd(fas->eventsmoothsolve, smoothu, 0, 0, 0);
386:   /* check convergence reason for the smoother */
387:   SNESGetConvergedReason(smoothu, &reason);
388:   if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
389:     snes->reason = SNES_DIVERGED_INNER;
390:     return 0;
391:   }
392:   SNESGetFunction(smoothu, &FPC, NULL, NULL);
393:   SNESGetAlwaysComputesFinalResidual(smoothu, &flg);
394:   if (!flg) SNESComputeFunction(smoothu, X, FPC);
395:   VecCopy(FPC, F);
396:   if (fnorm) VecNorm(F, NORM_2, fnorm);
397:   return 0;
398: }

400: /*@
401:    SNESFASCreateCoarseVec - create `Vec` corresponding to a state vector on one level coarser than current level

403:    Collective

405:    Input Parameter:
406: .  snes - `SNESFAS` object

408:    Output Parameter:
409: .  Xcoarse - vector on level one coarser than snes

411:    Level: developer

413: .seealso: `SNESFASSetRestriction()`, `SNESFASRestrict()`
414: @*/
415: PetscErrorCode SNESFASCreateCoarseVec(SNES snes, Vec *Xcoarse)
416: {
417:   SNES_FAS *fas;

421:   fas = (SNES_FAS *)snes->data;
422:   if (fas->rscale) {
423:     VecDuplicate(fas->rscale, Xcoarse);
424:   } else if (fas->inject) {
425:     MatCreateVecs(fas->inject, Xcoarse, NULL);
426:   } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must set restriction or injection");
427:   return 0;
428: }

430: /*@
431:    SNESFASRestrict - restrict a `Vec` to the next coarser level

433:    Collective

435:    Input Parameters:
436: +  fine - `SNES` from which to restrict
437: -  Xfine - vector to restrict

439:    Output Parameter:
440: .  Xcoarse - result of restriction

442:    Level: developer

444: .seealso: `SNES`, `SNESFAS`, `SNESFASSetRestriction()`, `SNESFASSetInjection()`
445: @*/
446: PetscErrorCode SNESFASRestrict(SNES fine, Vec Xfine, Vec Xcoarse)
447: {
448:   SNES_FAS *fas;

453:   fas = (SNES_FAS *)fine->data;
454:   if (fas->inject) {
455:     MatRestrict(fas->inject, Xfine, Xcoarse);
456:   } else {
457:     MatRestrict(fas->restrct, Xfine, Xcoarse);
458:     VecPointwiseMult(Xcoarse, fas->rscale, Xcoarse);
459:   }
460:   return 0;
461: }

463: /*

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

467: fine problem:   F(x) = b
468: coarse problem: F^c(x^c) = Rb, Initial guess Rx
469: interpolated solution: x^f = I x^c (total solution interpolation

471:  */
472: static PetscErrorCode SNESFASInterpolatedCoarseSolution(SNES snes, Vec X, Vec X_new)
473: {
474:   Vec                 X_c, B_c;
475:   SNESConvergedReason reason;
476:   SNES                next;
477:   Mat                 restrct, interpolate;
478:   SNES_FAS           *fasc;

480:   SNESFASCycleGetCorrection(snes, &next);
481:   if (next) {
482:     fasc = (SNES_FAS *)next->data;

484:     SNESFASCycleGetRestriction(snes, &restrct);
485:     SNESFASCycleGetInterpolation(snes, &interpolate);

487:     X_c = next->vec_sol;

489:     if (fasc->eventinterprestrict) PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0);
490:     /* restrict the total solution: Rb */
491:     SNESFASRestrict(snes, X, X_c);
492:     B_c = next->vec_rhs;
493:     if (snes->vec_rhs) {
494:       /* restrict the total rhs defect: Rb */
495:       MatRestrict(restrct, snes->vec_rhs, B_c);
496:     } else {
497:       VecSet(B_c, 0.);
498:     }
499:     if (fasc->eventinterprestrict) PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0);

501:     SNESSolve(next, B_c, X_c);
502:     SNESGetConvergedReason(next, &reason);
503:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
504:       snes->reason = SNES_DIVERGED_INNER;
505:       return 0;
506:     }
507:     /* x^f <- Ix^c*/
508:     DM dmc, dmf;

510:     SNESGetDM(next, &dmc);
511:     SNESGetDM(snes, &dmf);
512:     if (fasc->eventinterprestrict) PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0);
513:     DMInterpolateSolution(dmc, dmf, interpolate, X_c, X_new);
514:     if (fasc->eventinterprestrict) PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0);
515:     PetscObjectSetName((PetscObject)X_c, "Coarse solution");
516:     VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view");
517:     PetscObjectSetName((PetscObject)X_new, "Updated Fine solution");
518:     VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view");
519:   }
520:   return 0;
521: }

523: /*

525: Performs the FAS coarse correction as:

527: fine problem:   F(x) = b
528: coarse problem: F^c(x^c) = b^c

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

532:  */
533: PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new)
534: {
535:   Vec                 X_c, Xo_c, F_c, B_c;
536:   SNESConvergedReason reason;
537:   SNES                next;
538:   Mat                 restrct, interpolate;
539:   SNES_FAS           *fasc;

541:   SNESFASCycleGetCorrection(snes, &next);
542:   if (next) {
543:     fasc = (SNES_FAS *)next->data;

545:     SNESFASCycleGetRestriction(snes, &restrct);
546:     SNESFASCycleGetInterpolation(snes, &interpolate);

548:     X_c  = next->vec_sol;
549:     Xo_c = next->work[0];
550:     F_c  = next->vec_func;
551:     B_c  = next->vec_rhs;

553:     if (fasc->eventinterprestrict) PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0);
554:     SNESFASRestrict(snes, X, Xo_c);
555:     /* restrict the defect: R(F(x) - b) */
556:     MatRestrict(restrct, F, B_c);
557:     if (fasc->eventinterprestrict) PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0);

559:     if (fasc->eventresidual) PetscLogEventBegin(fasc->eventresidual, next, 0, 0, 0);
560:     /* F_c = F^c(Rx) - R(F(x) - b) since the second term was sitting in next->vec_rhs */
561:     SNESComputeFunction(next, Xo_c, F_c);
562:     if (fasc->eventresidual) PetscLogEventEnd(fasc->eventresidual, next, 0, 0, 0);

564:     /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */
565:     VecCopy(B_c, X_c);
566:     VecCopy(F_c, B_c);
567:     VecCopy(X_c, F_c);
568:     /* set initial guess of the coarse problem to the projected fine solution */
569:     VecCopy(Xo_c, X_c);

571:     /* recurse to the next level */
572:     SNESSetInitialFunction(next, F_c);
573:     SNESSolve(next, B_c, X_c);
574:     SNESGetConvergedReason(next, &reason);
575:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
576:       snes->reason = SNES_DIVERGED_INNER;
577:       return 0;
578:     }
579:     /* correct as x <- x + I(x^c - Rx)*/
580:     VecAXPY(X_c, -1.0, Xo_c);

582:     if (fasc->eventinterprestrict) PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0);
583:     MatInterpolateAdd(interpolate, X_c, X, X_new);
584:     if (fasc->eventinterprestrict) PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0);
585:     PetscObjectSetName((PetscObject)X_c, "Coarse correction");
586:     VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view");
587:     PetscObjectSetName((PetscObject)X_new, "Updated Fine solution");
588:     VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view");
589:   }
590:   return 0;
591: }

593: /*

595: The additive cycle looks like:

597: xhat = x
598: xhat = dS(x, b)
599: x = coarsecorrection(xhat, b_d)
600: x = x + nu*(xhat - x);
601: (optional) x = uS(x, b)

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

605:  */
606: static PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X)
607: {
608:   Vec                  F, B, Xhat;
609:   Vec                  X_c, Xo_c, F_c, B_c;
610:   SNESConvergedReason  reason;
611:   PetscReal            xnorm, fnorm, ynorm;
612:   SNESLineSearchReason lsresult;
613:   SNES                 next;
614:   Mat                  restrct, interpolate;
615:   SNES_FAS            *fas = (SNES_FAS *)snes->data, *fasc;

617:   SNESFASCycleGetCorrection(snes, &next);
618:   F    = snes->vec_func;
619:   B    = snes->vec_rhs;
620:   Xhat = snes->work[1];
621:   VecCopy(X, Xhat);
622:   /* recurse first */
623:   if (next) {
624:     fasc = (SNES_FAS *)next->data;
625:     SNESFASCycleGetRestriction(snes, &restrct);
626:     SNESFASCycleGetInterpolation(snes, &interpolate);
627:     if (fas->eventresidual) PetscLogEventBegin(fas->eventresidual, snes, 0, 0, 0);
628:     SNESComputeFunction(snes, Xhat, F);
629:     if (fas->eventresidual) PetscLogEventEnd(fas->eventresidual, snes, 0, 0, 0);
630:     VecNorm(F, NORM_2, &fnorm);
631:     X_c  = next->vec_sol;
632:     Xo_c = next->work[0];
633:     F_c  = next->vec_func;
634:     B_c  = next->vec_rhs;

636:     SNESFASRestrict(snes, Xhat, Xo_c);
637:     /* restrict the defect */
638:     MatRestrict(restrct, F, B_c);

640:     /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */
641:     if (fasc->eventresidual) PetscLogEventBegin(fasc->eventresidual, next, 0, 0, 0);
642:     SNESComputeFunction(next, Xo_c, F_c);
643:     if (fasc->eventresidual) PetscLogEventEnd(fasc->eventresidual, next, 0, 0, 0);
644:     VecCopy(B_c, X_c);
645:     VecCopy(F_c, B_c);
646:     VecCopy(X_c, F_c);
647:     /* set initial guess of the coarse problem to the projected fine solution */
648:     VecCopy(Xo_c, X_c);

650:     /* recurse */
651:     SNESSetInitialFunction(next, F_c);
652:     SNESSolve(next, B_c, X_c);

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

657:     SNESGetConvergedReason(next, &reason);
658:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
659:       snes->reason = SNES_DIVERGED_INNER;
660:       return 0;
661:     }

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

667:     /* additive correction of the coarse direction*/
668:     SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat);
669:     SNESLineSearchGetReason(snes->linesearch, &lsresult);
670:     SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm);
671:     if (lsresult) {
672:       if (++snes->numFailures >= snes->maxFailures) {
673:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
674:         return 0;
675:       }
676:     }
677:   } else {
678:     SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
679:   }
680:   return 0;
681: }

683: /*

685: Defines the FAS cycle as:

687: fine problem: F(x) = b
688: coarse problem: F^c(x) = b^c

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

692: correction:

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

696:  */
697: static PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X)
698: {
699:   Vec  F, B;
700:   SNES next;

702:   F = snes->vec_func;
703:   B = snes->vec_rhs;
704:   /* pre-smooth -- just update using the pre-smoother */
705:   SNESFASCycleGetCorrection(snes, &next);
706:   SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
707:   if (next) {
708:     SNESFASCoarseCorrection(snes, X, F, X);
709:     SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
710:   }
711:   return 0;
712: }

714: static PetscErrorCode SNESFASCycleSetupPhase_Full(SNES snes)
715: {
716:   SNES      next;
717:   SNES_FAS *fas = (SNES_FAS *)snes->data;
718:   PetscBool isFine;

720:   /* pre-smooth -- just update using the pre-smoother */
721:   SNESFASCycleIsFine(snes, &isFine);
722:   SNESFASCycleGetCorrection(snes, &next);
723:   fas->full_stage = 0;
724:   if (next) SNESFASCycleSetupPhase_Full(next);
725:   return 0;
726: }

728: static PetscErrorCode SNESFASCycle_Full(SNES snes, Vec X)
729: {
730:   Vec       F, B;
731:   SNES_FAS *fas = (SNES_FAS *)snes->data;
732:   PetscBool isFine;
733:   SNES      next;

735:   F = snes->vec_func;
736:   B = snes->vec_rhs;
737:   SNESFASCycleIsFine(snes, &isFine);
738:   SNESFASCycleGetCorrection(snes, &next);

740:   if (isFine) SNESFASCycleSetupPhase_Full(snes);

742:   if (fas->full_stage == 0) {
743:     /* downsweep */
744:     if (next) {
745:       if (fas->level != 1) next->max_its += 1;
746:       if (fas->full_downsweep) SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
747:       fas->full_downsweep = PETSC_TRUE;
748:       if (fas->full_total) SNESFASInterpolatedCoarseSolution(snes, X, X);
749:       else SNESFASCoarseCorrection(snes, X, F, X);
750:       fas->full_total = PETSC_FALSE;
751:       SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
752:       if (fas->level != 1) next->max_its -= 1;
753:     } else {
754:       /* The smoother on the coarse level is the coarse solver */
755:       SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
756:     }
757:     fas->full_stage = 1;
758:   } else if (fas->full_stage == 1) {
759:     if (snes->iter == 0) SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
760:     if (next) {
761:       SNESFASCoarseCorrection(snes, X, F, X);
762:       SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
763:     }
764:   }
765:   /* final v-cycle */
766:   if (isFine) {
767:     if (next) {
768:       SNESFASCoarseCorrection(snes, X, F, X);
769:       SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
770:     }
771:   }
772:   return 0;
773: }

775: static PetscErrorCode SNESFASCycle_Kaskade(SNES snes, Vec X)
776: {
777:   Vec  F, B;
778:   SNES next;

780:   F = snes->vec_func;
781:   B = snes->vec_rhs;
782:   SNESFASCycleGetCorrection(snes, &next);
783:   if (next) {
784:     SNESFASCoarseCorrection(snes, X, F, X);
785:     SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);
786:   } else {
787:     SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);
788:   }
789:   return 0;
790: }

792: PetscBool  SNEScite       = PETSC_FALSE;
793: const char SNESCitation[] = "@techreport{pbmkbsxt2012,\n"
794:                             "  title = {Composing Scalable Nonlinear Algebraic Solvers},\n"
795:                             "  author = {Peter Brune and Mathew Knepley and Barry Smith and Xuemin Tu},\n"
796:                             "  year = 2013,\n"
797:                             "  type = Preprint,\n"
798:                             "  number = {ANL/MCS-P2010-0112},\n"
799:                             "  institution = {Argonne National Laboratory}\n}\n";

801: static PetscErrorCode SNESSolve_FAS(SNES snes)
802: {
803:   PetscInt  i;
804:   Vec       X, F;
805:   PetscReal fnorm;
806:   SNES_FAS *fas = (SNES_FAS *)snes->data, *ffas;
807:   DM        dm;
808:   PetscBool isFine;


812:   PetscCitationsRegister(SNESCitation, &SNEScite);
813:   snes->reason = SNES_CONVERGED_ITERATING;
814:   X            = snes->vec_sol;
815:   F            = snes->vec_func;

817:   SNESFASCycleIsFine(snes, &isFine);
818:   /* norm setup */
819:   PetscObjectSAWsTakeAccess((PetscObject)snes);
820:   snes->iter = 0;
821:   snes->norm = 0;
822:   PetscObjectSAWsGrantAccess((PetscObject)snes);
823:   if (!snes->vec_func_init_set) {
824:     if (fas->eventresidual) PetscLogEventBegin(fas->eventresidual, snes, 0, 0, 0);
825:     SNESComputeFunction(snes, X, F);
826:     if (fas->eventresidual) PetscLogEventEnd(fas->eventresidual, snes, 0, 0, 0);
827:   } else snes->vec_func_init_set = PETSC_FALSE;

829:   VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
830:   SNESCheckFunctionNorm(snes, fnorm);
831:   PetscObjectSAWsTakeAccess((PetscObject)snes);
832:   snes->norm = fnorm;
833:   PetscObjectSAWsGrantAccess((PetscObject)snes);
834:   SNESLogConvergenceHistory(snes, fnorm, 0);
835:   SNESMonitor(snes, snes->iter, fnorm);

837:   /* test convergence */
838:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
839:   if (snes->reason) return 0;

841:   if (isFine) {
842:     /* propagate scale-dependent data up the hierarchy */
843:     SNESGetDM(snes, &dm);
844:     for (ffas = fas; ffas->next; ffas = (SNES_FAS *)ffas->next->data) {
845:       DM dmcoarse;
846:       SNESGetDM(ffas->next, &dmcoarse);
847:       DMRestrict(dm, ffas->restrct, ffas->rscale, ffas->inject, dmcoarse);
848:       dm = dmcoarse;
849:     }
850:   }

852:   for (i = 0; i < snes->max_its; i++) {
853:     /* Call general purpose update function */
854:     PetscTryTypeMethod(snes, update, snes->iter);

856:     if (fas->fastype == SNES_FAS_MULTIPLICATIVE) {
857:       SNESFASCycle_Multiplicative(snes, X);
858:     } else if (fas->fastype == SNES_FAS_ADDITIVE) {
859:       SNESFASCycle_Additive(snes, X);
860:     } else if (fas->fastype == SNES_FAS_FULL) {
861:       SNESFASCycle_Full(snes, X);
862:     } else if (fas->fastype == SNES_FAS_KASKADE) {
863:       SNESFASCycle_Kaskade(snes, X);
864:     } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Unsupported FAS type");

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

869:     /* Monitor convergence */
870:     PetscObjectSAWsTakeAccess((PetscObject)snes);
871:     snes->iter = i + 1;
872:     PetscObjectSAWsGrantAccess((PetscObject)snes);
873:     SNESLogConvergenceHistory(snes, snes->norm, 0);
874:     SNESMonitor(snes, snes->iter, snes->norm);
875:     /* Test for convergence */
876:     if (isFine) {
877:       PetscUseTypeMethod(snes, converged, snes->iter, 0.0, 0.0, snes->norm, &snes->reason, snes->cnvP);
878:       if (snes->reason) break;
879:     }
880:   }
881:   if (i == snes->max_its) {
882:     PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", i);
883:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
884:   }
885:   return 0;
886: }

888: /*MC

890: SNESFAS - Full Approximation Scheme nonlinear multigrid solver.

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

896:    Options Database Keys and Prefixes:
897: +   -snes_fas_levels -  The number of levels
898: .   -snes_fas_cycles<1> -  The number of cycles -- 1 for V, 2 for W
899: .   -snes_fas_type<additive,multiplicative,full,kaskade>  -  Additive or multiplicative cycle
900: .   -snes_fas_galerkin<`PETSC_FALSE`> -  Form coarse problems by projection back upon the fine problem
901: .   -snes_fas_smoothup<1> -  The number of iterations of the post-smoother
902: .   -snes_fas_smoothdown<1> -  The number of iterations of the pre-smoother
903: .   -snes_fas_monitor -  Monitor progress of all of the levels
904: .   -snes_fas_full_downsweep<`PETSC_FALSE`> - call the downsmooth on the initial downsweep of full FAS
905: .   -fas_levels_snes_ -  `SNES` options for all smoothers
906: .   -fas_levels_cycle_snes_ -  `SNES` options for all cycles
907: .   -fas_levels_i_snes_ -  `SNES` options for the smoothers on level i
908: .   -fas_levels_i_cycle_snes_ - `SNES` options for the cycle on level i
909: -   -fas_coarse_snes_ -  `SNES` options for the coarsest smoother

911:    Note:
912:    The organization of the FAS solver is slightly different from the organization of `PCMG`
913:    As each level has smoother `SNES` instances(down and potentially up) and a cycle `SNES` instance.
914:    The cycle `SNES` instance may be used for monitoring convergence on a particular level.

916:    Level: beginner

918:    References:
919: .  * - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
920:    SIAM Review, 57(4), 2015

922: .seealso: `PCMG`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESType`, `SNESFASSetRestriction()`, `SNESFASSetInjection()`,
923:           `SNESFASFullGetTotal()`
924: M*/

926: PETSC_EXTERN PetscErrorCode SNESCreate_FAS(SNES snes)
927: {
928:   SNES_FAS *fas;

930:   snes->ops->destroy        = SNESDestroy_FAS;
931:   snes->ops->setup          = SNESSetUp_FAS;
932:   snes->ops->setfromoptions = SNESSetFromOptions_FAS;
933:   snes->ops->view           = SNESView_FAS;
934:   snes->ops->solve          = SNESSolve_FAS;
935:   snes->ops->reset          = SNESReset_FAS;

937:   snes->usesksp = PETSC_FALSE;
938:   snes->usesnpc = PETSC_FALSE;

940:   if (!snes->tolerancesset) {
941:     snes->max_funcs = 30000;
942:     snes->max_its   = 10000;
943:   }

945:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

947:   PetscNew(&fas);

949:   snes->data                  = (void *)fas;
950:   fas->level                  = 0;
951:   fas->levels                 = 1;
952:   fas->n_cycles               = 1;
953:   fas->max_up_it              = 1;
954:   fas->max_down_it            = 1;
955:   fas->smoothu                = NULL;
956:   fas->smoothd                = NULL;
957:   fas->next                   = NULL;
958:   fas->previous               = NULL;
959:   fas->fine                   = snes;
960:   fas->interpolate            = NULL;
961:   fas->restrct                = NULL;
962:   fas->inject                 = NULL;
963:   fas->usedmfornumberoflevels = PETSC_FALSE;
964:   fas->fastype                = SNES_FAS_MULTIPLICATIVE;
965:   fas->full_downsweep         = PETSC_FALSE;
966:   fas->full_total             = PETSC_FALSE;

968:   fas->eventsmoothsetup    = 0;
969:   fas->eventsmoothsolve    = 0;
970:   fas->eventresidual       = 0;
971:   fas->eventinterprestrict = 0;
972:   return 0;
973: }