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:   PetscFunctionBegin;
 11:   PetscCall(SNESDestroy(&fas->smoothu));
 12:   PetscCall(SNESDestroy(&fas->smoothd));
 13:   PetscCall(MatDestroy(&fas->inject));
 14:   PetscCall(MatDestroy(&fas->interpolate));
 15:   PetscCall(MatDestroy(&fas->restrct));
 16:   PetscCall(VecDestroy(&fas->rscale));
 17:   PetscCall(VecDestroy(&fas->Xg));
 18:   PetscCall(VecDestroy(&fas->Fg));
 19:   if (fas->next) PetscCall(SNESReset(fas->next));
 20:   PetscFunctionReturn(PETSC_SUCCESS);
 21: }

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

 27:   PetscFunctionBegin;
 28:   /* recursively resets and then destroys */
 29:   PetscCall(SNESReset_FAS(snes));
 30:   PetscCall(SNESDestroy(&fas->next));
 31:   PetscCall(PetscFree(fas));
 32:   PetscFunctionReturn(PETSC_SUCCESS);
 33: }

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

 43:   PetscFunctionBegin;
 44:   if (!snes->linesearch) PetscFunctionReturn(PETSC_SUCCESS);
 45:   PetscCall(SNESGetLineSearch(snes, &linesearch));
 46:   PetscCall(SNESGetLineSearch(smooth, &slinesearch));
 47:   PetscCall(SNESLineSearchGetPreCheck(linesearch, &precheck, &lsprectx));
 48:   PetscCall(SNESLineSearchGetPostCheck(linesearch, &postcheck, &lspostctx));
 49:   PetscCall(SNESLineSearchSetPreCheck(slinesearch, precheck, lsprectx));
 50:   PetscCall(SNESLineSearchSetPostCheck(slinesearch, postcheck, lspostctx));
 51:   PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch));
 52:   PetscFunctionReturn(PETSC_SUCCESS);
 53: }

 55: static PetscErrorCode SNESFASCycleSetUpSmoother_Private(SNES snes, SNES smooth)
 56: {
 57:   SNES_FAS *fas = (SNES_FAS *)snes->data;

 59:   PetscFunctionBegin;
 60:   PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)smooth));
 61:   PetscCall(SNESSetFromOptions(smooth));
 62:   PetscCall(SNESFASSetUpLineSearch_Private(snes, smooth));

 64:   PetscCall(PetscObjectReference((PetscObject)snes->vec_sol));
 65:   PetscCall(PetscObjectReference((PetscObject)snes->vec_sol_update));
 66:   PetscCall(PetscObjectReference((PetscObject)snes->vec_func));
 67:   smooth->vec_sol        = snes->vec_sol;
 68:   smooth->vec_sol_update = snes->vec_sol_update;
 69:   smooth->vec_func       = snes->vec_func;

 71:   if (fas->eventsmoothsetup) PetscCall(PetscLogEventBegin(fas->eventsmoothsetup, smooth, 0, 0, 0));
 72:   PetscCall(SNESSetUp(smooth));
 73:   if (fas->eventsmoothsetup) PetscCall(PetscLogEventEnd(fas->eventsmoothsetup, smooth, 0, 0, 0));
 74:   PetscFunctionReturn(PETSC_SUCCESS);
 75: }

 77: static PetscErrorCode SNESSetUp_FAS(SNES snes)
 78: {
 79:   SNES_FAS *fas = (SNES_FAS *)snes->data;
 80:   PetscInt  dm_levels;
 81:   SNES      next;
 82:   PetscBool isFine, hasCreateRestriction, hasCreateInjection;

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

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

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

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

136:   /*pass the smoother, function, and jacobian up to the next level if it's not user set already */
137:   if (fas->galerkin) {
138:     if (next) PetscCall(SNESSetFunction(next, NULL, SNESFASGalerkinFunctionDefault, next));
139:     if (fas->smoothd && fas->level != fas->levels - 1) PetscCall(SNESSetFunction(fas->smoothd, NULL, SNESFASGalerkinFunctionDefault, snes));
140:     if (fas->smoothu && fas->level != fas->levels - 1) PetscCall(SNESSetFunction(fas->smoothu, NULL, SNESFASGalerkinFunctionDefault, snes));
141:   }

143:   /* sets the down (pre) smoother's default norm and sets it from options */
144:   if (fas->smoothd) {
145:     if (fas->level == 0) {
146:       PetscCall(SNESSetNormSchedule(fas->smoothd, SNES_NORM_ALWAYS));
147:     } else {
148:       PetscCall(SNESSetNormSchedule(fas->smoothd, SNES_NORM_FINAL_ONLY));
149:     }
150:     PetscCall(SNESFASCycleSetUpSmoother_Private(snes, fas->smoothd));
151:   }

153:   /* sets the up (post) smoother's default norm and sets it from options */
154:   if (fas->smoothu) {
155:     if (fas->level != fas->levels - 1) {
156:       PetscCall(SNESSetNormSchedule(fas->smoothu, SNES_NORM_NONE));
157:     } else {
158:       PetscCall(SNESSetNormSchedule(fas->smoothu, SNES_NORM_FINAL_ONLY));
159:     }
160:     PetscCall(SNESFASCycleSetUpSmoother_Private(snes, fas->smoothu));
161:   }

163:   if (next) {
164:     /* gotta set up the solution vector for this to work */
165:     if (!next->vec_sol) PetscCall(SNESFASCreateCoarseVec(snes, &next->vec_sol));
166:     if (!next->vec_rhs) PetscCall(SNESFASCreateCoarseVec(snes, &next->vec_rhs));
167:     PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next));
168:     PetscCall(SNESFASSetUpLineSearch_Private(snes, next));
169:     PetscCall(SNESSetUp(next));
170:   }

172:   /* setup FAS work vectors */
173:   if (fas->galerkin) {
174:     PetscCall(VecDuplicate(snes->vec_sol, &fas->Xg));
175:     PetscCall(VecDuplicate(snes->vec_sol, &fas->Fg));
176:   }
177:   PetscFunctionReturn(PETSC_SUCCESS);
178: }

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

192:   PetscFunctionBegin;
193:   PetscCall(SNESFASCycleIsFine(snes, &isFine));
194:   PetscOptionsHeadBegin(PetscOptionsObject, "SNESFAS Options-----------------------------------");

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

209:     PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
210:     PetscCall(PetscOptionsInt("-snes_fas_cycles", "Number of cycles", "SNESFASSetCycles", fas->n_cycles, &m, &flg));
211:     if (flg) PetscCall(SNESFASSetCycles(snes, m));
212:     PetscCall(PetscOptionsBool("-snes_fas_continuation", "Corrected grid-sequence continuation", "SNESFASSetContinuation", fas->continuation, &continuationflg, &flg));
213:     if (flg) PetscCall(SNESFASSetContinuation(snes, continuationflg));

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

218:     if (fas->fastype == SNES_FAS_FULL) {
219:       PetscCall(PetscOptionsBool("-snes_fas_full_downsweep", "Smooth on the initial down sweep for full FAS cycles", "SNESFASFullSetDownSweep", fas->full_downsweep, &fas->full_downsweep, &flg));
220:       if (flg) PetscCall(SNESFASFullSetDownSweep(snes, fas->full_downsweep));
221:       PetscCall(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));
222:       if (flg) PetscCall(SNESFASFullSetTotal(snes, fas->full_total));
223:     }

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

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

229:     {
230:       PetscViewer       viewer;
231:       PetscViewerFormat format;
232:       PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_fas_monitor", &viewer, &format, &monflg));
233:       if (monflg) {
234:         PetscViewerAndFormat *vf;
235:         PetscCall(PetscViewerAndFormatCreate(viewer, format, &vf));
236:         PetscCall(PetscViewerDestroy(&viewer));
237:         PetscCall(SNESFASSetMonitor(snes, vf, PETSC_TRUE));
238:       }
239:     }
240:     flg    = PETSC_FALSE;
241:     monflg = PETSC_TRUE;
242:     PetscCall(PetscOptionsBool("-snes_fas_log", "Log times for each FAS level", "SNESFASSetLog", monflg, &monflg, &flg));
243:     if (flg) PetscCall(SNESFASSetLog(snes, monflg));
244:   }
245:   PetscCall(PetscOptionsBool("-snes_fas_monitor_correction", "View the tau correction at each iteration", "SNESFASCoarseCorrection", fas->monitorCorrection, &fas->monitorCorrection, &flg));

247:   PetscOptionsHeadEnd();

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

253:   /* set up the default line search for coarse grid corrections */
254:   if (fas->fastype == SNES_FAS_ADDITIVE) {
255:     if (!snes->linesearch) {
256:       PetscCall(SNESGetLineSearch(snes, &linesearch));
257:       PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHSECANT));
258:     }
259:   }

261:   /* recursive option setting for the smoothers */
262:   PetscCall(SNESFASCycleGetCorrection(snes, &next));
263:   if (next) PetscCall(SNESSetFromOptions(next));
264:   PetscFunctionReturn(PETSC_SUCCESS);
265: }

267: #include <petscdraw.h>
268: static PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer)
269: {
270:   SNES_FAS *fas = (SNES_FAS *)snes->data;
271:   PetscBool isFine, isascii, isdraw;
272:   PetscInt  i;
273:   SNES      smoothu, smoothd, levelsnes;

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

348: /*
349: Defines the action of the downsmoother
350:  */
351: static PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
352: {
353:   SNESConvergedReason reason;
354:   Vec                 FPC;
355:   SNES                smoothd;
356:   PetscBool           flg;
357:   SNES_FAS           *fas = (SNES_FAS *)snes->data;

359:   PetscFunctionBegin;
360:   PetscCall(SNESFASCycleGetSmootherDown(snes, &smoothd));
361:   PetscCall(SNESSetInitialFunction(smoothd, F));
362:   if (fas->eventsmoothsolve) PetscCall(PetscLogEventBegin(fas->eventsmoothsolve, smoothd, B, X, 0));
363:   PetscCall(SNESSolve(smoothd, B, X));
364:   if (fas->eventsmoothsolve) PetscCall(PetscLogEventEnd(fas->eventsmoothsolve, smoothd, B, X, 0));
365:   /* check convergence reason for the smoother */
366:   PetscCall(SNESGetConvergedReason(smoothd, &reason));
367:   if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
368:     snes->reason = SNES_DIVERGED_INNER;
369:     PetscFunctionReturn(PETSC_SUCCESS);
370:   }

372:   PetscCall(SNESGetFunction(smoothd, &FPC, NULL, NULL));
373:   PetscCall(SNESGetAlwaysComputesFinalResidual(smoothd, &flg));
374:   if (!flg) PetscCall(SNESComputeFunction(smoothd, X, FPC));
375:   PetscCall(VecCopy(FPC, F));
376:   if (fnorm) {
377:     PetscCall(VecNorm(F, NORM_2, fnorm));
378:     SNESCheckFunctionDomainError(snes, *fnorm);
379:   }
380:   PetscFunctionReturn(PETSC_SUCCESS);
381: }

383: /*
384: Defines the action of the upsmoother
385:  */
386: static PetscErrorCode SNESFASUpSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
387: {
388:   SNESConvergedReason reason;
389:   Vec                 FPC;
390:   SNES                smoothu;
391:   PetscBool           flg;
392:   SNES_FAS           *fas = (SNES_FAS *)snes->data;

394:   PetscFunctionBegin;
395:   PetscCall(SNESFASCycleGetSmootherUp(snes, &smoothu));
396:   if (fas->eventsmoothsolve) PetscCall(PetscLogEventBegin(fas->eventsmoothsolve, smoothu, 0, 0, 0));
397:   PetscCall(SNESSolve(smoothu, B, X));
398:   if (fas->eventsmoothsolve) PetscCall(PetscLogEventEnd(fas->eventsmoothsolve, smoothu, 0, 0, 0));
399:   /* check convergence reason for the smoother */
400:   PetscCall(SNESGetConvergedReason(smoothu, &reason));
401:   if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
402:     snes->reason = SNES_DIVERGED_INNER;
403:     PetscFunctionReturn(PETSC_SUCCESS);
404:   }
405:   PetscCall(SNESGetFunction(smoothu, &FPC, NULL, NULL));
406:   PetscCall(SNESGetAlwaysComputesFinalResidual(smoothu, &flg));
407:   if (!flg) PetscCall(SNESComputeFunction(smoothu, X, FPC));
408:   PetscCall(VecCopy(FPC, F));
409:   if (fnorm) {
410:     PetscCall(VecNorm(F, NORM_2, fnorm));
411:     SNESCheckFunctionDomainError(snes, *fnorm);
412:   }
413:   PetscFunctionReturn(PETSC_SUCCESS);
414: }

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

419:   Collective

421:   Input Parameter:
422: . snes - `SNESFAS` object

424:   Output Parameter:
425: . Xcoarse - vector on level one coarser than the current level

427:   Level: developer

429: .seealso: [](ch_snes), `SNESFASSetRestriction()`, `SNESFASRestrict()`, `SNESFAS`
430: @*/
431: PetscErrorCode SNESFASCreateCoarseVec(SNES snes, Vec *Xcoarse)
432: {
433:   SNES_FAS *fas;

435:   PetscFunctionBegin;
437:   PetscAssertPointer(Xcoarse, 2);
438:   fas = (SNES_FAS *)snes->data;
439:   if (fas->rscale) {
440:     PetscCall(VecDuplicate(fas->rscale, Xcoarse));
441:   } else if (fas->interpolate) {
442:     PetscCall(MatCreateVecs(fas->interpolate, Xcoarse, NULL));
443:   } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must set rscale or interpolation");
444:   PetscFunctionReturn(PETSC_SUCCESS);
445: }

447: /*@
448:   SNESFASRestrict - restrict a `Vec` to the next coarser level

450:   Collective

452:   Input Parameters:
453: + fine  - `SNES` from which to restrict
454: - Xfine - vector to restrict

456:   Output Parameter:
457: . Xcoarse - result of restriction

459:   Level: developer

461: .seealso: [](ch_snes), `SNES`, `SNESFAS`, `SNESFASSetRestriction()`, `SNESFASSetInjection()`, `SNESFASCreateCoarseVec()`
462: @*/
463: PetscErrorCode SNESFASRestrict(SNES fine, Vec Xfine, Vec Xcoarse)
464: {
465:   SNES_FAS *fas;

467:   PetscFunctionBegin;
471:   fas = (SNES_FAS *)fine->data;
472:   if (fas->inject) {
473:     PetscCall(MatRestrict(fas->inject, Xfine, Xcoarse));
474:   } else {
475:     PetscCall(MatRestrict(fas->restrct, Xfine, Xcoarse));
476:     PetscCall(VecPointwiseMult(Xcoarse, fas->rscale, Xcoarse));
477:   }
478:   PetscFunctionReturn(PETSC_SUCCESS);
479: }

481: /*
482: Performs a variant of FAS using the interpolated total coarse solution

484: fine problem:   F(x) = b
485: coarse problem: F^c(x^c) = Rb, Initial guess Rx
486: interpolated solution: x^f = I x^c (total solution interpolation
487:  */
488: static PetscErrorCode SNESFASInterpolatedCoarseSolution(SNES snes, Vec X, Vec X_new)
489: {
490:   Vec                 X_c, B_c;
491:   SNESConvergedReason reason;
492:   SNES                next;
493:   Mat                 restrct, interpolate;
494:   SNES_FAS           *fasc;

496:   PetscFunctionBegin;
497:   PetscCall(SNESFASCycleGetCorrection(snes, &next));
498:   if (next) {
499:     fasc = (SNES_FAS *)next->data;

501:     PetscCall(SNESFASCycleGetRestriction(snes, &restrct));
502:     PetscCall(SNESFASCycleGetInterpolation(snes, &interpolate));

504:     X_c = next->vec_sol;

506:     if (fasc->eventinterprestrict) PetscCall(PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0));
507:     /* restrict the total solution: Rb */
508:     PetscCall(SNESFASRestrict(snes, X, X_c));
509:     B_c = next->vec_rhs;
510:     if (snes->vec_rhs) {
511:       /* restrict the total rhs defect: Rb */
512:       PetscCall(MatRestrict(restrct, snes->vec_rhs, B_c));
513:     } else {
514:       PetscCall(VecSet(B_c, 0.));
515:     }
516:     if (fasc->eventinterprestrict) PetscCall(PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0));

518:     PetscCall(SNESSolve(next, B_c, X_c));
519:     PetscCall(SNESGetConvergedReason(next, &reason));
520:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
521:       snes->reason = SNES_DIVERGED_INNER;
522:       PetscFunctionReturn(PETSC_SUCCESS);
523:     }
524:     /* x^f <- Ix^c*/
525:     DM dmc, dmf;

527:     PetscCall(SNESGetDM(next, &dmc));
528:     PetscCall(SNESGetDM(snes, &dmf));
529:     if (fasc->eventinterprestrict) PetscCall(PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0));
530:     PetscCall(DMInterpolateSolution(dmc, dmf, interpolate, X_c, X_new));
531:     if (fasc->eventinterprestrict) PetscCall(PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0));
532:     PetscCall(PetscObjectSetName((PetscObject)X_c, "Coarse solution"));
533:     PetscCall(VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view"));
534:     PetscCall(PetscObjectSetName((PetscObject)X_new, "Updated Fine solution"));
535:     PetscCall(VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view"));
536:   }
537:   PetscFunctionReturn(PETSC_SUCCESS);
538: }

540: /*
541: Performs the FAS coarse correction as:

543: fine problem:   F(x) = b
544: coarse problem: F^c(x^c) = b^c

546: b^c = F^c(Rx) - R(F(x) - b)
547:     = tau + R b
548:  */
549: static PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new)
550: {
551:   PetscBool           monitorCorrection = ((SNES_FAS *)snes->data)->monitorCorrection;
552:   Vec                 X_c, Xo_c, F_c, B_c;
553:   SNESConvergedReason reason;
554:   SNES                next;
555:   Mat                 restrct, interpolate;
556:   SNES_FAS           *fasc;

558:   PetscFunctionBegin;
559:   PetscCall(SNESFASCycleGetCorrection(snes, &next));
560:   if (next) {
561:     fasc = (SNES_FAS *)next->data;

563:     PetscCall(SNESFASCycleGetRestriction(snes, &restrct));
564:     PetscCall(SNESFASCycleGetInterpolation(snes, &interpolate));

566:     X_c  = next->vec_sol;
567:     Xo_c = next->work[0];
568:     F_c  = next->vec_func;
569:     B_c  = next->vec_rhs;

571:     if (fasc->eventinterprestrict) PetscCall(PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0));
572:     PetscCall(SNESFASRestrict(snes, X, Xo_c));
573:     /* restrict the defect: R(F(x) - b) */
574:     PetscCall(MatRestrict(restrct, F, B_c));
575:     if (monitorCorrection) {
576:       PetscReal fnorm, rnorm;

578:       PetscCall(VecNorm(F, NORM_2, &fnorm));
579:       PetscCall(VecNorm(B_c, NORM_2, &rnorm));
580:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), "||F(X)|| %g\n||R (F(X) - b)|| %g\n", (double)fnorm, (double)rnorm));
581:     }
582:     if (fasc->eventinterprestrict) PetscCall(PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0));

584:     if (fasc->eventresidual) PetscCall(PetscLogEventBegin(fasc->eventresidual, next, 0, 0, 0));
585:     /* F_c = F^c(Rx) - R(F(x) - b) since the second term was sitting in next->vec_rhs */
586:     PetscCall(SNESComputeFunction(next, Xo_c, F_c));
587:     if (monitorCorrection) {
588:       PetscReal xnorm, fnorm;

590:       PetscCall(VecNorm(Xo_c, NORM_2, &xnorm));
591:       PetscCall(VecNorm(F_c, NORM_2, &fnorm));
592:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), "||R(X)|| %g\n||tau + Rb = F_c (R X) - R(F(X) - b)|| %g\n", (double)xnorm, (double)fnorm));
593:       PetscCall(PetscObjectSetName((PetscObject)F_c, "tau"));
594:       PetscCall(VecViewFromOptions(F_c, (PetscObject)next, "-tau_view"));
595:     }
596:     if (fasc->eventresidual) PetscCall(PetscLogEventEnd(fasc->eventresidual, next, 0, 0, 0));

598:     /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */
599:     PetscCall(VecCopy(B_c, X_c)); // Now x^c holds R(F(x) - b)
600:     PetscCall(VecCopy(F_c, B_c)); // Now b^c holds F^c(Rx) - R(F(x) - b)
601:     PetscCall(VecCopy(X_c, F_c)); // Now F^c holds R(F(x) - b)
602:     /* set initial guess of the coarse problem to the projected fine solution */
603:     PetscCall(VecCopy(Xo_c, X_c));

605:     /* recurse to the next level */
606:     // So x^c_0 = R x and F^c_0 = F^c(x^c_0) - b^c = F^c(x^c_0) - (F^c(x^c_0) - R(F(x) - b)) = R(F(x) - b)
607:     PetscCall(SNESSetInitialFunction(next, F_c));
608:     PetscCall(SNESSolve(next, B_c, X_c));
609:     PetscCall(SNESGetConvergedReason(next, &reason));
610:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
611:       snes->reason = SNES_DIVERGED_INNER;
612:       PetscFunctionReturn(PETSC_SUCCESS);
613:     }
614:     /* correct as x <- x + I(x^c - Rx)*/
615:     PetscCall(VecAXPY(X_c, -1.0, Xo_c));

617:     if (fasc->eventinterprestrict) PetscCall(PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0));
618:     PetscCall(MatInterpolateAdd(interpolate, X_c, X, X_new));
619:     if (fasc->eventinterprestrict) PetscCall(PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0));
620:     if (monitorCorrection) {
621:       PetscReal xnorm, xonorm, inorm;

623:       PetscCall(VecNorm(X_c, NORM_2, &xnorm));
624:       PetscCall(VecNorm(X, NORM_2, &xonorm));
625:       PetscCall(VecNorm(X_new, NORM_2, &inorm));
626:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), "||X_c - Xo_c|| %g\n||X|| %g\n||X + I (X_c - X_co)|| %g\n", (double)xnorm, (double)xonorm, (double)inorm));
627:     }
628:     // TODO Check for snes->b, Technically we should strip out R b here if it is nonzero
629:     PetscCall(PetscObjectSetName((PetscObject)B_c, "Tau correction"));
630:     PetscCall(VecViewFromOptions(B_c, NULL, "-fas_tau_correction_view"));
631:     PetscCall(PetscObjectSetName((PetscObject)X_c, "Coarse correction"));
632:     PetscCall(VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view"));
633:     PetscCall(PetscObjectSetName((PetscObject)X_new, "Updated Fine solution"));
634:     PetscCall(VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view"));
635:   }
636:   PetscFunctionReturn(PETSC_SUCCESS);
637: }

639: /*
640: The additive cycle is:

642: xhat = x
643: xhat = dS(x, b)
644: x = coarsecorrection(xhat, b_d)
645: x = x + nu*(xhat - x);
646: (optional) x = uS(x, b)

648: With the coarse RHS (defect correction) as below.
649:  */
650: static PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X)
651: {
652:   Vec                 F, B, Xhat;
653:   Vec                 X_c, Xo_c, F_c, B_c;
654:   SNESConvergedReason reason;
655:   PetscReal           xnorm, fnorm, ynorm;
656:   SNES                next;
657:   Mat                 restrct, interpolate;
658:   SNES_FAS           *fas = (SNES_FAS *)snes->data, *fasc;

660:   PetscFunctionBegin;
661:   PetscCall(SNESFASCycleGetCorrection(snes, &next));
662:   F    = snes->vec_func;
663:   B    = snes->vec_rhs;
664:   Xhat = snes->work[1];
665:   PetscCall(VecCopy(X, Xhat));
666:   /* recurse first */
667:   if (next) {
668:     fasc = (SNES_FAS *)next->data;
669:     PetscCall(SNESFASCycleGetRestriction(snes, &restrct));
670:     PetscCall(SNESFASCycleGetInterpolation(snes, &interpolate));
671:     if (fas->eventresidual) PetscCall(PetscLogEventBegin(fas->eventresidual, snes, 0, 0, 0));
672:     PetscCall(SNESComputeFunction(snes, Xhat, F));
673:     if (fas->eventresidual) PetscCall(PetscLogEventEnd(fas->eventresidual, snes, 0, 0, 0));
674:     PetscCall(VecNorm(F, NORM_2, &fnorm));
675:     SNESCheckFunctionDomainError(snes, fnorm);
676:     X_c  = next->vec_sol;
677:     Xo_c = next->work[0];
678:     F_c  = next->vec_func;
679:     B_c  = next->vec_rhs;

681:     PetscCall(SNESFASRestrict(snes, Xhat, Xo_c));
682:     /* restrict the defect */
683:     PetscCall(MatRestrict(restrct, F, B_c));

685:     /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */
686:     if (fasc->eventresidual) PetscCall(PetscLogEventBegin(fasc->eventresidual, next, 0, 0, 0));
687:     PetscCall(SNESComputeFunction(next, Xo_c, F_c));
688:     if (fasc->eventresidual) PetscCall(PetscLogEventEnd(fasc->eventresidual, next, 0, 0, 0));
689:     PetscCall(VecCopy(B_c, X_c));
690:     PetscCall(VecCopy(F_c, B_c));
691:     PetscCall(VecCopy(X_c, F_c));
692:     /* set initial guess of the coarse problem to the projected fine solution */
693:     PetscCall(VecCopy(Xo_c, X_c));

695:     /* recurse */
696:     PetscCall(SNESSetInitialFunction(next, F_c));
697:     PetscCall(SNESSolve(next, B_c, X_c));

699:     /* smooth on this level */
700:     PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &fnorm));

702:     PetscCall(SNESGetConvergedReason(next, &reason));
703:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
704:       snes->reason = SNES_DIVERGED_INNER;
705:       PetscFunctionReturn(PETSC_SUCCESS);
706:     }

708:     /* correct as x <- x + I(x^c - Rx)*/
709:     PetscCall(VecAYPX(X_c, -1.0, Xo_c));
710:     PetscCall(MatInterpolate(interpolate, X_c, Xhat));

712:     /* additive correction of the coarse direction*/
713:     PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat));
714:     PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm));
715:     SNESCheckLineSearchFailure(snes);
716:   } else {
717:     PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm));
718:   }
719:   PetscFunctionReturn(PETSC_SUCCESS);
720: }

722: /*
723: Defines the FAS cycle as:

725: fine problem: F(x) = b
726: coarse problem: F^c(x) = b^c

728: b^c = F^c(Rx) - R(F(x) - b)
729:     = F^c(Rx) - R(F(x)) + R b
730:     = tau + R b

732: correction:

734: x = x + I(x^c - Rx)
735:  */
736: static PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X)
737: {
738:   Vec  F, B;
739:   SNES next;

741:   PetscFunctionBegin;
742:   F = snes->vec_func;
743:   B = snes->vec_rhs;
744:   /* pre-smooth -- just update using the pre-smoother */
745:   PetscCall(SNESFASCycleGetCorrection(snes, &next));
746:   PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm));
747:   if (next) {
748:     PetscCall(SNESFASCoarseCorrection(snes, X, F, X));
749:     PetscCall(SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm));
750:   }
751:   PetscFunctionReturn(PETSC_SUCCESS);
752: }

754: static PetscErrorCode SNESFASCycleSetupPhase_Full(SNES snes)
755: {
756:   SNES      next;
757:   SNES_FAS *fas = (SNES_FAS *)snes->data;
758:   PetscBool isFine;

760:   PetscFunctionBegin;
761:   /* pre-smooth -- just update using the pre-smoother */
762:   PetscCall(SNESFASCycleIsFine(snes, &isFine));
763:   PetscCall(SNESFASCycleGetCorrection(snes, &next));
764:   fas->full_stage = 0;
765:   if (next) PetscCall(SNESFASCycleSetupPhase_Full(next));
766:   PetscFunctionReturn(PETSC_SUCCESS);
767: }

769: static PetscErrorCode SNESFASCycle_Full(SNES snes, Vec X)
770: {
771:   Vec       F, B;
772:   SNES_FAS *fas = (SNES_FAS *)snes->data;
773:   PetscBool isFine;
774:   SNES      next;

776:   PetscFunctionBegin;
777:   F = snes->vec_func;
778:   B = snes->vec_rhs;
779:   PetscCall(SNESFASCycleIsFine(snes, &isFine));
780:   PetscCall(SNESFASCycleGetCorrection(snes, &next));

782:   if (isFine) PetscCall(SNESFASCycleSetupPhase_Full(snes));

784:   if (fas->full_stage == 0) {
785:     /* downsweep */
786:     if (next) {
787:       if (fas->level != 1) next->max_its += 1;
788:       if (fas->full_downsweep) PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm));
789:       fas->full_downsweep = PETSC_TRUE;
790:       if (fas->full_total) PetscCall(SNESFASInterpolatedCoarseSolution(snes, X, X));
791:       else PetscCall(SNESFASCoarseCorrection(snes, X, F, X));
792:       fas->full_total = PETSC_FALSE;
793:       PetscCall(SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm));
794:       if (fas->level != 1) next->max_its -= 1;
795:     } else {
796:       /* The smoother on the coarse level is the coarse solver */
797:       PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm));
798:     }
799:     fas->full_stage = 1;
800:   } else if (fas->full_stage == 1) {
801:     if (snes->iter == 0) PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm));
802:     if (next) {
803:       PetscCall(SNESFASCoarseCorrection(snes, X, F, X));
804:       PetscCall(SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm));
805:     }
806:   }
807:   /* final v-cycle */
808:   if (isFine) {
809:     if (next) {
810:       PetscCall(SNESFASCoarseCorrection(snes, X, F, X));
811:       PetscCall(SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm));
812:     }
813:   }
814:   PetscFunctionReturn(PETSC_SUCCESS);
815: }

817: static PetscErrorCode SNESFASCycle_Kaskade(SNES snes, Vec X)
818: {
819:   Vec  F, B;
820:   SNES next;

822:   PetscFunctionBegin;
823:   F = snes->vec_func;
824:   B = snes->vec_rhs;
825:   PetscCall(SNESFASCycleGetCorrection(snes, &next));
826:   if (next) {
827:     PetscCall(SNESFASCoarseCorrection(snes, X, F, X));
828:     PetscCall(SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm));
829:   } else {
830:     PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm));
831:   }
832:   PetscFunctionReturn(PETSC_SUCCESS);
833: }

835: PetscBool  SNEScite       = PETSC_FALSE;
836: const char SNESCitation[] = "@Article{bruneknepleysmithtu15,"
837:                             "  title         = {Composing Scalable Nonlinear Algebraic Solvers},"
838:                             "  author        = {Peter R. Brune and Matthew G. Knepley and Barry F. Smith and Xuemin Tu},"
839:                             "  journal       = {SIAM Review},"
840:                             "  volume        = {57},"
841:                             "  number        = {4},"
842:                             "  pages         = {535--565},"
843:                             "  doi           = {10.1137/130936725},"
844:                             "  year          = {2015}\n";

846: static PetscErrorCode SNESSolve_FAS(SNES snes)
847: {
848:   PetscInt  i;
849:   Vec       X, F;
850:   PetscReal fnorm;
851:   SNES_FAS *fas = (SNES_FAS *)snes->data, *ffas;
852:   DM        dm;
853:   PetscBool isFine;

855:   PetscFunctionBegin;
856:   PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);

858:   PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
859:   snes->reason = SNES_CONVERGED_ITERATING;
860:   X            = snes->vec_sol;
861:   F            = snes->vec_func;

863:   PetscCall(SNESFASCycleIsFine(snes, &isFine));
864:   /* norm setup */
865:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
866:   snes->iter = 0;
867:   snes->norm = 0;
868:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
869:   if (!snes->vec_func_init_set) {
870:     if (fas->eventresidual) PetscCall(PetscLogEventBegin(fas->eventresidual, snes, 0, 0, 0));
871:     PetscCall(SNESComputeFunction(snes, X, F));
872:     if (fas->eventresidual) PetscCall(PetscLogEventEnd(fas->eventresidual, snes, 0, 0, 0));
873:   } else snes->vec_func_init_set = PETSC_FALSE;

875:   PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F||  */
876:   SNESCheckFunctionDomainError(snes, fnorm);
877:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
878:   snes->norm = fnorm;
879:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
880:   PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));

882:   /* test convergence */
883:   PetscCall(SNESConverged(snes, snes->iter, 0.0, 0.0, fnorm));
884:   PetscCall(SNESMonitor(snes, snes->iter, fnorm));
885:   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

887:   if (isFine) {
888:     /* propagate scale-dependent data up the hierarchy */
889:     PetscCall(SNESGetDM(snes, &dm));
890:     for (ffas = fas; ffas->next; ffas = (SNES_FAS *)ffas->next->data) {
891:       DM dmcoarse;
892:       PetscCall(SNESGetDM(ffas->next, &dmcoarse));
893:       PetscCall(DMRestrict(dm, ffas->restrct, ffas->rscale, ffas->inject, dmcoarse));
894:       dm = dmcoarse;
895:     }
896:   }

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

902:     if (fas->fastype == SNES_FAS_MULTIPLICATIVE) PetscCall(SNESFASCycle_Multiplicative(snes, X));
903:     else if (fas->fastype == SNES_FAS_ADDITIVE) PetscCall(SNESFASCycle_Additive(snes, X));
904:     else if (fas->fastype == SNES_FAS_FULL) PetscCall(SNESFASCycle_Full(snes, X));
905:     else {
906:       PetscCheck(fas->fastype == SNES_FAS_KASKADE, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Unsupported FAS type");
907:       PetscCall(SNESFASCycle_Kaskade(snes, X));
908:     }

910:     /* check for FAS cycle divergence */
911:     if (snes->reason != SNES_CONVERGED_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);

913:     /* Monitor convergence */
914:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
915:     snes->iter = i + 1;
916:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
917:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
918:     PetscCall(SNESConverged(snes, snes->iter, 0.0, 0.0, snes->norm));
919:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
920:     if (snes->reason) break;
921:   }
922:   PetscFunctionReturn(PETSC_SUCCESS);
923: }

925: /*MC
926:    SNESFAS - An implementation of the Full Approximation Scheme nonlinear multigrid solver, FAS, or nonlinear multigrid {cite}`bruneknepleysmithtu15` for
927:              solving nonlinear systems of equations with `SNES`.

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

933:    Options Database Keys and Prefixes:
934: +   -snes_fas_levels l                                    - The number of levels
935: .   -snes_fas_cycles (1|2)                                - The number of cycles -- 1 for V, 2 for W
936: .   -snes_fas_type (additive|multiplicative|full|kaskade) - Additive or multiplicative cycle
937: .   -snes_fas_galerkin (true|false)                       - Form coarse problems by projection back upon the fine problem
938: .   -snes_fas_smoothup u                                  - The number of iterations of the post-smoother
939: .   -snes_fas_smoothdown d                                - The number of iterations of the pre-smoother
940: .   -snes_fas_monitor                                     - Monitor progress of all of the levels
941: .   -snes_fas_full_downsweep (true|false)                 - call the downsmooth on the initial downsweep of full FAS
942: .   -fas_levels_snes_                                     - prefix for `SNES` options for all smoothers
943: .   -fas_levels_cycle_snes_                               - prefix for `SNES` options for all cycles
944: .   -fas_levels_i_snes_                                   - prefix `SNES` options for the smoothers on level i
945: .   -fas_levels_i_cycle_snes_                             - prefix for `SNES` options for the cycle on level i
946: -   -fas_coarse_snes_                                     - prefix for `SNES` options for the coarsest smoother

948:    Level: beginner

950:    Note:
951:    The organization of the `SNESFAS` solver is slightly different from the organization of `PCMG`
952:    As each level has smoother `SNES` instances(down and potentially up) and a cycle `SNES` instance.
953:    The cycle `SNES` instance may be used for monitoring convergence on a particular level.

955: .seealso: [](ch_snes), `PCMG`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESType`, `SNESFASSetRestriction()`, `SNESFASSetInjection()`,
956:           `SNESFASFullGetTotal()`, `SNESFASSetType()`, `SNESFASGetType()`, `SNESFASSetLevels()`, `SNESFASGetLevels()`, `SNESFASGetCycleSNES()`,
957:           `SNESFASSetNumberSmoothUp()`, `SNESFASSetNumberSmoothDown()`, `SNESFASSetContinuation()`, `SNESFASSetCycles()`, `SNESFASSetMonitor()`,
958:           `SNESFASSetLog()`, `SNESFASCycleSetCycles()`, `SNESFASCycleGetSmoother()`, `SNESFASCycleGetSmootherUp()`, `SNESFASCycleGetSmootherDown()`,
959:           `SNESFASCycleGetCorrection()`, `SNESFASCycleGetInterpolation()`, `SNESFASCycleGetRestriction()`, `SNESFASCycleGetInjection()`,
960:           `SNESFASCycleGetRScale()`, `SNESFASCycleIsFine()`, `SNESFASSetInterpolation()`, `SNESFASGetInterpolation()`,
961:           `SNESFASGetRestriction()`, `SNESFASGetInjection()`, `SNESFASSetRScale()`, `SNESFASGetSmoother()`,
962:           `SNESFASGetSmootherDown()`, `SNESFASGetSmootherUp()`, `SNESFASGetCoarseSolve()`, `SNESFASFullSetDownSweep()`, `SNESFASFullSetTotal()`
963: M*/

965: PETSC_EXTERN PetscErrorCode SNESCreate_FAS(SNES snes)
966: {
967:   SNES_FAS *fas;

969:   PetscFunctionBegin;
970:   snes->ops->destroy        = SNESDestroy_FAS;
971:   snes->ops->setup          = SNESSetUp_FAS;
972:   snes->ops->setfromoptions = SNESSetFromOptions_FAS;
973:   snes->ops->view           = SNESView_FAS;
974:   snes->ops->solve          = SNESSolve_FAS;
975:   snes->ops->reset          = SNESReset_FAS;

977:   snes->usesksp = PETSC_FALSE;
978:   snes->usesnpc = PETSC_FALSE;

980:   PetscCall(SNESParametersInitialize(snes));
981:   PetscObjectParameterSetDefault(snes, max_funcs, 30000);
982:   PetscObjectParameterSetDefault(snes, max_its, 10000);

984:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

986:   PetscCall(PetscNew(&fas));

988:   snes->data                  = (void *)fas;
989:   fas->level                  = 0;
990:   fas->levels                 = 1;
991:   fas->n_cycles               = 1;
992:   fas->max_up_it              = 1;
993:   fas->max_down_it            = 1;
994:   fas->smoothu                = NULL;
995:   fas->smoothd                = NULL;
996:   fas->next                   = NULL;
997:   fas->previous               = NULL;
998:   fas->fine                   = snes;
999:   fas->interpolate            = NULL;
1000:   fas->restrct                = NULL;
1001:   fas->inject                 = NULL;
1002:   fas->usedmfornumberoflevels = PETSC_FALSE;
1003:   fas->fastype                = SNES_FAS_MULTIPLICATIVE;
1004:   fas->full_downsweep         = PETSC_FALSE;
1005:   fas->full_total             = PETSC_FALSE;

1007:   fas->eventsmoothsetup    = 0;
1008:   fas->eventsmoothsolve    = 0;
1009:   fas->eventresidual       = 0;
1010:   fas->eventinterprestrict = 0;
1011:   PetscFunctionReturn(PETSC_SUCCESS);
1012: }