xref: /petsc/src/snes/impls/fas/fas.c (revision bd89dbf26d8a5efecb980364933175da61864cd7)
1421d9b32SPeter Brune /* Defines the basic SNES object */
2a055c8caSBarry Smith #include <../src/snes/impls/fas/fasimpls.h> /*I  "petscsnes.h"  I*/
3421d9b32SPeter Brune 
49e5d0892SLisandro Dalcin const char *const SNESFASTypes[] = {"MULTIPLICATIVE", "ADDITIVE", "FULL", "KASKADE", "SNESFASType", "SNES_FAS", NULL};
507144faaSPeter Brune 
SNESReset_FAS(SNES snes)6d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESReset_FAS(SNES snes)
7d71ae5a4SJacob Faibussowitsch {
8421d9b32SPeter Brune   SNES_FAS *fas = (SNES_FAS *)snes->data;
9421d9b32SPeter Brune 
10421d9b32SPeter Brune   PetscFunctionBegin;
119566063dSJacob Faibussowitsch   PetscCall(SNESDestroy(&fas->smoothu));
129566063dSJacob Faibussowitsch   PetscCall(SNESDestroy(&fas->smoothd));
139566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&fas->inject));
149566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&fas->interpolate));
159566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&fas->restrct));
169566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&fas->rscale));
179566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&fas->Xg));
189566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&fas->Fg));
199566063dSJacob Faibussowitsch   if (fas->next) PetscCall(SNESReset(fas->next));
203ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
21421d9b32SPeter Brune }
22421d9b32SPeter Brune 
SNESDestroy_FAS(SNES snes)23d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESDestroy_FAS(SNES snes)
24d71ae5a4SJacob Faibussowitsch {
25421d9b32SPeter Brune   SNES_FAS *fas = (SNES_FAS *)snes->data;
26421d9b32SPeter Brune 
27421d9b32SPeter Brune   PetscFunctionBegin;
28421d9b32SPeter Brune   /* recursively resets and then destroys */
299566063dSJacob Faibussowitsch   PetscCall(SNESReset_FAS(snes));
309566063dSJacob Faibussowitsch   PetscCall(SNESDestroy(&fas->next));
319566063dSJacob Faibussowitsch   PetscCall(PetscFree(fas));
323ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
33421d9b32SPeter Brune }
34421d9b32SPeter Brune 
SNESFASSetUpLineSearch_Private(SNES snes,SNES smooth)35d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESFASSetUpLineSearch_Private(SNES snes, SNES smooth)
36d71ae5a4SJacob Faibussowitsch {
37f833ba53SLisandro Dalcin   SNESLineSearch linesearch;
38f833ba53SLisandro Dalcin   SNESLineSearch slinesearch;
39f833ba53SLisandro Dalcin   void          *lsprectx, *lspostctx;
40f833ba53SLisandro Dalcin   PetscErrorCode (*precheck)(SNESLineSearch, Vec, Vec, PetscBool *, void *);
41f833ba53SLisandro Dalcin   PetscErrorCode (*postcheck)(SNESLineSearch, Vec, Vec, Vec, PetscBool *, PetscBool *, void *);
42f833ba53SLisandro Dalcin 
43f833ba53SLisandro Dalcin   PetscFunctionBegin;
443ba16761SJacob Faibussowitsch   if (!snes->linesearch) PetscFunctionReturn(PETSC_SUCCESS);
459566063dSJacob Faibussowitsch   PetscCall(SNESGetLineSearch(snes, &linesearch));
469566063dSJacob Faibussowitsch   PetscCall(SNESGetLineSearch(smooth, &slinesearch));
479566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetPreCheck(linesearch, &precheck, &lsprectx));
489566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetPostCheck(linesearch, &postcheck, &lspostctx));
499566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetPreCheck(slinesearch, precheck, lsprectx));
509566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetPostCheck(slinesearch, postcheck, lspostctx));
519566063dSJacob Faibussowitsch   PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch));
523ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
53f833ba53SLisandro Dalcin }
54f833ba53SLisandro Dalcin 
SNESFASCycleSetUpSmoother_Private(SNES snes,SNES smooth)55d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESFASCycleSetUpSmoother_Private(SNES snes, SNES smooth)
56d71ae5a4SJacob Faibussowitsch {
57f833ba53SLisandro Dalcin   SNES_FAS *fas = (SNES_FAS *)snes->data;
58f833ba53SLisandro Dalcin 
59f833ba53SLisandro Dalcin   PetscFunctionBegin;
609566063dSJacob Faibussowitsch   PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)smooth));
619566063dSJacob Faibussowitsch   PetscCall(SNESSetFromOptions(smooth));
629566063dSJacob Faibussowitsch   PetscCall(SNESFASSetUpLineSearch_Private(snes, smooth));
63f833ba53SLisandro Dalcin 
649566063dSJacob Faibussowitsch   PetscCall(PetscObjectReference((PetscObject)snes->vec_sol));
659566063dSJacob Faibussowitsch   PetscCall(PetscObjectReference((PetscObject)snes->vec_sol_update));
669566063dSJacob Faibussowitsch   PetscCall(PetscObjectReference((PetscObject)snes->vec_func));
67f833ba53SLisandro Dalcin   smooth->vec_sol        = snes->vec_sol;
68f833ba53SLisandro Dalcin   smooth->vec_sol_update = snes->vec_sol_update;
69f833ba53SLisandro Dalcin   smooth->vec_func       = snes->vec_func;
70f833ba53SLisandro Dalcin 
719566063dSJacob Faibussowitsch   if (fas->eventsmoothsetup) PetscCall(PetscLogEventBegin(fas->eventsmoothsetup, smooth, 0, 0, 0));
729566063dSJacob Faibussowitsch   PetscCall(SNESSetUp(smooth));
739566063dSJacob Faibussowitsch   if (fas->eventsmoothsetup) PetscCall(PetscLogEventEnd(fas->eventsmoothsetup, smooth, 0, 0, 0));
743ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
75f833ba53SLisandro Dalcin }
76f833ba53SLisandro Dalcin 
SNESSetUp_FAS(SNES snes)77d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESSetUp_FAS(SNES snes)
78d71ae5a4SJacob Faibussowitsch {
7948bfdf8aSPeter Brune   SNES_FAS *fas = (SNES_FAS *)snes->data;
80d1adcc6fSPeter Brune   PetscInt  dm_levels;
81ab8d36c9SPeter Brune   SNES      next;
82f833ba53SLisandro Dalcin   PetscBool isFine, hasCreateRestriction, hasCreateInjection;
83eff52c0eSPeter Brune 
846b2b7091SBarry Smith   PetscFunctionBegin;
859566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleIsFine(snes, &isFine));
86ab8d36c9SPeter Brune   if (fas->usedmfornumberoflevels && isFine) {
879566063dSJacob Faibussowitsch     PetscCall(DMGetRefineLevel(snes->dm, &dm_levels));
88d1adcc6fSPeter Brune     dm_levels++;
89cc05f883SPeter Brune     if (dm_levels > fas->levels) {
903dccd265SPeter Brune       /* reset the number of levels */
919566063dSJacob Faibussowitsch       PetscCall(SNESFASSetLevels(snes, dm_levels, NULL));
929566063dSJacob Faibussowitsch       PetscCall(SNESSetFromOptions(snes));
93d1adcc6fSPeter Brune     }
94d1adcc6fSPeter Brune   }
959566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleGetCorrection(snes, &next));
96ab8d36c9SPeter Brune   if (!isFine) snes->gridsequence = 0; /* no grid sequencing inside the multigrid hierarchy! */
973dccd265SPeter Brune 
989566063dSJacob Faibussowitsch   PetscCall(SNESSetWorkVecs(snes, 2)); /* work vectors used for intergrid transfers */
99cc05f883SPeter Brune 
100ab8d36c9SPeter Brune   /* set up the smoothers if they haven't already been set up */
10148a46eb9SPierre Jolivet   if (!fas->smoothd) PetscCall(SNESFASCycleCreateSmoother_Private(snes, &fas->smoothd));
102ab8d36c9SPeter Brune 
10379d9a41aSPeter Brune   if (snes->dm) {
104ab8d36c9SPeter Brune     /* set the smoother DMs properly */
1059566063dSJacob Faibussowitsch     if (fas->smoothu) PetscCall(SNESSetDM(fas->smoothu, snes->dm));
1069566063dSJacob Faibussowitsch     PetscCall(SNESSetDM(fas->smoothd, snes->dm));
10779d9a41aSPeter Brune     /* construct EVERYTHING from the DM -- including the progressive set of smoothers */
108ab8d36c9SPeter Brune     if (next) {
10979d9a41aSPeter Brune       /* for now -- assume the DM and the evaluation functions have been set externally */
110ab8d36c9SPeter Brune       if (!next->dm) {
1119566063dSJacob Faibussowitsch         PetscCall(DMCoarsen(snes->dm, PetscObjectComm((PetscObject)next), &next->dm));
1129566063dSJacob Faibussowitsch         PetscCall(SNESSetDM(next, next->dm));
11379d9a41aSPeter Brune       }
11479d9a41aSPeter Brune       /* set the interpolation and restriction from the DM */
11579d9a41aSPeter Brune       if (!fas->interpolate) {
1169566063dSJacob Faibussowitsch         PetscCall(DMCreateInterpolation(next->dm, snes->dm, &fas->interpolate, &fas->rscale));
117bccf9bb3SJed Brown         if (!fas->restrct) {
1189566063dSJacob Faibussowitsch           PetscCall(DMHasCreateRestriction(next->dm, &hasCreateRestriction));
1190a7266b2SPatrick Farrell           /* DM can create restrictions, use that */
1200a7266b2SPatrick Farrell           if (hasCreateRestriction) {
1219566063dSJacob Faibussowitsch             PetscCall(DMCreateRestriction(next->dm, snes->dm, &fas->restrct));
1220a7266b2SPatrick Farrell           } else {
1239566063dSJacob Faibussowitsch             PetscCall(PetscObjectReference((PetscObject)fas->interpolate));
12479d9a41aSPeter Brune             fas->restrct = fas->interpolate;
12579d9a41aSPeter Brune           }
126bccf9bb3SJed Brown         }
1270a7266b2SPatrick Farrell       }
12879d9a41aSPeter Brune       /* set the injection from the DM */
12979d9a41aSPeter Brune       if (!fas->inject) {
1309566063dSJacob Faibussowitsch         PetscCall(DMHasCreateInjection(next->dm, &hasCreateInjection));
13148a46eb9SPierre Jolivet         if (hasCreateInjection) PetscCall(DMCreateInjection(next->dm, snes->dm, &fas->inject));
13279d9a41aSPeter Brune       }
13379d9a41aSPeter Brune     }
13423e68893SLawrence Mitchell   }
135f833ba53SLisandro Dalcin 
13679d9a41aSPeter Brune   /*pass the smoother, function, and jacobian up to the next level if it's not user set already */
13779d9a41aSPeter Brune   if (fas->galerkin) {
1381baa6e33SBarry Smith     if (next) PetscCall(SNESSetFunction(next, NULL, SNESFASGalerkinFunctionDefault, next));
13948a46eb9SPierre Jolivet     if (fas->smoothd && fas->level != fas->levels - 1) PetscCall(SNESSetFunction(fas->smoothd, NULL, SNESFASGalerkinFunctionDefault, snes));
14048a46eb9SPierre Jolivet     if (fas->smoothu && fas->level != fas->levels - 1) PetscCall(SNESSetFunction(fas->smoothu, NULL, SNESFASGalerkinFunctionDefault, snes));
14179d9a41aSPeter Brune   }
14279d9a41aSPeter Brune 
143534ebe21SPeter Brune   /* sets the down (pre) smoother's default norm and sets it from options */
144534ebe21SPeter Brune   if (fas->smoothd) {
1452d157150SStefano Zampini     if (fas->level == 0) {
1462d157150SStefano Zampini       PetscCall(SNESSetNormSchedule(fas->smoothd, SNES_NORM_ALWAYS));
147534ebe21SPeter Brune     } else {
1489566063dSJacob Faibussowitsch       PetscCall(SNESSetNormSchedule(fas->smoothd, SNES_NORM_FINAL_ONLY));
149534ebe21SPeter Brune     }
1509566063dSJacob Faibussowitsch     PetscCall(SNESFASCycleSetUpSmoother_Private(snes, fas->smoothd));
151534ebe21SPeter Brune   }
152534ebe21SPeter Brune 
153534ebe21SPeter Brune   /* sets the up (post) smoother's default norm and sets it from options */
154534ebe21SPeter Brune   if (fas->smoothu) {
155534ebe21SPeter Brune     if (fas->level != fas->levels - 1) {
1569566063dSJacob Faibussowitsch       PetscCall(SNESSetNormSchedule(fas->smoothu, SNES_NORM_NONE));
157534ebe21SPeter Brune     } else {
1589566063dSJacob Faibussowitsch       PetscCall(SNESSetNormSchedule(fas->smoothu, SNES_NORM_FINAL_ONLY));
159534ebe21SPeter Brune     }
1609566063dSJacob Faibussowitsch     PetscCall(SNESFASCycleSetUpSmoother_Private(snes, fas->smoothu));
161534ebe21SPeter Brune   }
162d06165b7SPeter Brune 
163ab8d36c9SPeter Brune   if (next) {
16479d9a41aSPeter Brune     /* gotta set up the solution vector for this to work */
1659566063dSJacob Faibussowitsch     if (!next->vec_sol) PetscCall(SNESFASCreateCoarseVec(snes, &next->vec_sol));
1669566063dSJacob Faibussowitsch     if (!next->vec_rhs) PetscCall(SNESFASCreateCoarseVec(snes, &next->vec_rhs));
1679566063dSJacob Faibussowitsch     PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next));
1689566063dSJacob Faibussowitsch     PetscCall(SNESFASSetUpLineSearch_Private(snes, next));
1699566063dSJacob Faibussowitsch     PetscCall(SNESSetUp(next));
17079d9a41aSPeter Brune   }
171f833ba53SLisandro Dalcin 
1726273346dSPeter Brune   /* setup FAS work vectors */
1736273346dSPeter Brune   if (fas->galerkin) {
1749566063dSJacob Faibussowitsch     PetscCall(VecDuplicate(snes->vec_sol, &fas->Xg));
1759566063dSJacob Faibussowitsch     PetscCall(VecDuplicate(snes->vec_sol, &fas->Fg));
1766273346dSPeter Brune   }
1773ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
178421d9b32SPeter Brune }
179421d9b32SPeter Brune 
SNESSetFromOptions_FAS(SNES snes,PetscOptionItems PetscOptionsObject)180ce78bad3SBarry Smith static PetscErrorCode SNESSetFromOptions_FAS(SNES snes, PetscOptionItems PetscOptionsObject)
181d71ae5a4SJacob Faibussowitsch {
182ee78dd50SPeter Brune   SNES_FAS      *fas    = (SNES_FAS *)snes->data;
183ee78dd50SPeter Brune   PetscInt       levels = 1;
18487f44e3fSPeter Brune   PetscBool      flg = PETSC_FALSE, upflg = PETSC_FALSE, downflg = PETSC_FALSE, monflg = PETSC_FALSE, galerkinflg = PETSC_FALSE, continuationflg = PETSC_FALSE;
18507144faaSPeter Brune   SNESFASType    fastype;
186fde0ff24SPeter Brune   const char    *optionsprefix;
187f1c6b773SPeter Brune   SNESLineSearch linesearch;
18866585501SPeter Brune   PetscInt       m, n_up, n_down;
189ab8d36c9SPeter Brune   SNES           next;
190ab8d36c9SPeter Brune   PetscBool      isFine;
191421d9b32SPeter Brune 
192421d9b32SPeter Brune   PetscFunctionBegin;
1939566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleIsFine(snes, &isFine));
194d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "SNESFAS Options-----------------------------------");
195ee78dd50SPeter Brune 
196ab8d36c9SPeter Brune   /* number of levels -- only process most options on the finest level */
197ab8d36c9SPeter Brune   if (isFine) {
1989566063dSJacob Faibussowitsch     PetscCall(PetscOptionsInt("-snes_fas_levels", "Number of Levels", "SNESFASSetLevels", levels, &levels, &flg));
199c732cbdbSBarry Smith     if (!flg && snes->dm) {
2009566063dSJacob Faibussowitsch       PetscCall(DMGetRefineLevel(snes->dm, &levels));
201c732cbdbSBarry Smith       levels++;
202d1adcc6fSPeter Brune       fas->usedmfornumberoflevels = PETSC_TRUE;
203c732cbdbSBarry Smith     }
2049566063dSJacob Faibussowitsch     PetscCall(SNESFASSetLevels(snes, levels, NULL));
20507144faaSPeter Brune     fastype = fas->fastype;
2069566063dSJacob Faibussowitsch     PetscCall(PetscOptionsEnum("-snes_fas_type", "FAS correction type", "SNESFASSetType", SNESFASTypes, (PetscEnum)fastype, (PetscEnum *)&fastype, &flg));
2071baa6e33SBarry Smith     if (flg) PetscCall(SNESFASSetType(snes, fastype));
208ee78dd50SPeter Brune 
2099566063dSJacob Faibussowitsch     PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
2109566063dSJacob Faibussowitsch     PetscCall(PetscOptionsInt("-snes_fas_cycles", "Number of cycles", "SNESFASSetCycles", fas->n_cycles, &m, &flg));
2111baa6e33SBarry Smith     if (flg) PetscCall(SNESFASSetCycles(snes, m));
2129566063dSJacob Faibussowitsch     PetscCall(PetscOptionsBool("-snes_fas_continuation", "Corrected grid-sequence continuation", "SNESFASSetContinuation", fas->continuation, &continuationflg, &flg));
2131baa6e33SBarry Smith     if (flg) PetscCall(SNESFASSetContinuation(snes, continuationflg));
214fde0ff24SPeter Brune 
2159566063dSJacob Faibussowitsch     PetscCall(PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin", "SNESFASSetGalerkin", fas->galerkin, &galerkinflg, &flg));
2161baa6e33SBarry Smith     if (flg) PetscCall(SNESFASSetGalerkin(snes, galerkinflg));
217ee78dd50SPeter Brune 
218928e959bSPeter Brune     if (fas->fastype == SNES_FAS_FULL) {
2199566063dSJacob Faibussowitsch       PetscCall(PetscOptionsBool("-snes_fas_full_downsweep", "Smooth on the initial down sweep for full FAS cycles", "SNESFASFullSetDownSweep", fas->full_downsweep, &fas->full_downsweep, &flg));
2209566063dSJacob Faibussowitsch       if (flg) PetscCall(SNESFASFullSetDownSweep(snes, fas->full_downsweep));
2219566063dSJacob Faibussowitsch       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));
2229566063dSJacob Faibussowitsch       if (flg) PetscCall(SNESFASFullSetTotal(snes, fas->full_total));
223928e959bSPeter Brune     }
224928e959bSPeter Brune 
2259566063dSJacob Faibussowitsch     PetscCall(PetscOptionsInt("-snes_fas_smoothup", "Number of post-smoothing steps", "SNESFASSetNumberSmoothUp", fas->max_up_it, &n_up, &upflg));
226162d76ddSPeter Brune 
2279566063dSJacob Faibussowitsch     PetscCall(PetscOptionsInt("-snes_fas_smoothdown", "Number of pre-smoothing steps", "SNESFASSetNumberSmoothDown", fas->max_down_it, &n_down, &downflg));
228162d76ddSPeter Brune 
229d142ab34SLawrence Mitchell     {
230d142ab34SLawrence Mitchell       PetscViewer       viewer;
231d142ab34SLawrence Mitchell       PetscViewerFormat format;
232648c30bcSBarry Smith       PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_fas_monitor", &viewer, &format, &monflg));
233d142ab34SLawrence Mitchell       if (monflg) {
234d142ab34SLawrence Mitchell         PetscViewerAndFormat *vf;
2359566063dSJacob Faibussowitsch         PetscCall(PetscViewerAndFormatCreate(viewer, format, &vf));
236648c30bcSBarry Smith         PetscCall(PetscViewerDestroy(&viewer));
2379566063dSJacob Faibussowitsch         PetscCall(SNESFASSetMonitor(snes, vf, PETSC_TRUE));
238d142ab34SLawrence Mitchell       }
239d142ab34SLawrence Mitchell     }
2400dd27c6cSPeter Brune     flg    = PETSC_FALSE;
2410dd27c6cSPeter Brune     monflg = PETSC_TRUE;
2429566063dSJacob Faibussowitsch     PetscCall(PetscOptionsBool("-snes_fas_log", "Log times for each FAS level", "SNESFASSetLog", monflg, &monflg, &flg));
2439566063dSJacob Faibussowitsch     if (flg) PetscCall(SNESFASSetLog(snes, monflg));
244ab8d36c9SPeter Brune   }
245ee78dd50SPeter Brune 
246d0609cedSBarry Smith   PetscOptionsHeadEnd();
247f833ba53SLisandro Dalcin 
2488cc86e31SPeter Brune   /* setup from the determined types if there is no pointwise procedure or smoother defined */
2491baa6e33SBarry Smith   if (upflg) PetscCall(SNESFASSetNumberSmoothUp(snes, n_up));
2501baa6e33SBarry Smith   if (downflg) PetscCall(SNESFASSetNumberSmoothDown(snes, n_down));
251eff52c0eSPeter Brune 
2529e764e56SPeter Brune   /* set up the default line search for coarse grid corrections */
2539e764e56SPeter Brune   if (fas->fastype == SNES_FAS_ADDITIVE) {
2549e764e56SPeter Brune     if (!snes->linesearch) {
2559566063dSJacob Faibussowitsch       PetscCall(SNESGetLineSearch(snes, &linesearch));
256a99ef635SJonas Heinzmann       PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHSECANT));
2579e764e56SPeter Brune     }
2589e764e56SPeter Brune   }
2599e764e56SPeter Brune 
260ee78dd50SPeter Brune   /* recursive option setting for the smoothers */
2619566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleGetCorrection(snes, &next));
2629566063dSJacob Faibussowitsch   if (next) PetscCall(SNESSetFromOptions(next));
2633ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
264421d9b32SPeter Brune }
265421d9b32SPeter Brune 
2669804daf3SBarry Smith #include <petscdraw.h>
SNESView_FAS(SNES snes,PetscViewer viewer)267d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer)
268d71ae5a4SJacob Faibussowitsch {
269421d9b32SPeter Brune   SNES_FAS *fas = (SNES_FAS *)snes->data;
2709f196a02SMartin Diehl   PetscBool isFine, isascii, isdraw;
271ab8d36c9SPeter Brune   PetscInt  i;
272ab8d36c9SPeter Brune   SNES      smoothu, smoothd, levelsnes;
273421d9b32SPeter Brune 
274421d9b32SPeter Brune   PetscFunctionBegin;
2759566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleIsFine(snes, &isFine));
276ab8d36c9SPeter Brune   if (isFine) {
2779f196a02SMartin Diehl     PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2789566063dSJacob Faibussowitsch     PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
2799f196a02SMartin Diehl     if (isascii) {
28063a3b9bcSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  type is %s, levels=%" PetscInt_FMT ", cycles=%" PetscInt_FMT "\n", SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles));
281ab8d36c9SPeter Brune       if (fas->galerkin) {
2829566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(viewer, "  Using Galerkin computed coarse grid function evaluation\n"));
283421d9b32SPeter Brune       } else {
2849566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(viewer, "  Not using Galerkin computed coarse grid function evaluation\n"));
285421d9b32SPeter Brune       }
286ab8d36c9SPeter Brune       for (i = 0; i < fas->levels; i++) {
2879566063dSJacob Faibussowitsch         PetscCall(SNESFASGetCycleSNES(snes, i, &levelsnes));
2889566063dSJacob Faibussowitsch         PetscCall(SNESFASCycleGetSmootherUp(levelsnes, &smoothu));
2899566063dSJacob Faibussowitsch         PetscCall(SNESFASCycleGetSmootherDown(levelsnes, &smoothd));
290ab8d36c9SPeter Brune         if (!i) {
29163a3b9bcSJacob Faibussowitsch           PetscCall(PetscViewerASCIIPrintf(viewer, "  Coarse grid solver -- level %" PetscInt_FMT " -------------------------------\n", i));
292421d9b32SPeter Brune         } else {
29363a3b9bcSJacob Faibussowitsch           PetscCall(PetscViewerASCIIPrintf(viewer, "  Down solver (pre-smoother) on level %" PetscInt_FMT " -------------------------------\n", i));
294421d9b32SPeter Brune         }
2959566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPushTab(viewer));
296166b3ea4SJed Brown         if (smoothd) {
2979566063dSJacob Faibussowitsch           PetscCall(SNESView(smoothd, viewer));
298166b3ea4SJed Brown         } else {
2999566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIPrintf(viewer, "Not yet available\n"));
300166b3ea4SJed Brown         }
3019566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPopTab(viewer));
302ab8d36c9SPeter Brune         if (i && (smoothd == smoothu)) {
3039566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIPrintf(viewer, "  Up solver (post-smoother) same as down solver (pre-smoother)\n"));
304ab8d36c9SPeter Brune         } else if (i) {
30563a3b9bcSJacob Faibussowitsch           PetscCall(PetscViewerASCIIPrintf(viewer, "  Up solver (post-smoother) on level %" PetscInt_FMT " -------------------------------\n", i));
3069566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIPushTab(viewer));
307166b3ea4SJed Brown           if (smoothu) {
3089566063dSJacob Faibussowitsch             PetscCall(SNESView(smoothu, viewer));
309166b3ea4SJed Brown           } else {
3109566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIPrintf(viewer, "Not yet available\n"));
311166b3ea4SJed Brown           }
3129566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIPopTab(viewer));
313ab8d36c9SPeter Brune         }
314ab8d36c9SPeter Brune       }
315656ede7eSPeter Brune     } else if (isdraw) {
316656ede7eSPeter Brune       PetscDraw draw;
317b4375e8dSPeter Brune       PetscReal x, w, y, bottom, th, wth;
318656ede7eSPeter Brune       SNES_FAS *curfas = fas;
3199566063dSJacob Faibussowitsch       PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
3209566063dSJacob Faibussowitsch       PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
3219566063dSJacob Faibussowitsch       PetscCall(PetscDrawStringGetSize(draw, &wth, &th));
322656ede7eSPeter Brune       bottom = y - th;
323656ede7eSPeter Brune       while (curfas) {
324b4375e8dSPeter Brune         if (!curfas->smoothu) {
3259566063dSJacob Faibussowitsch           PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
3269566063dSJacob Faibussowitsch           if (curfas->smoothd) PetscCall(SNESView(curfas->smoothd, viewer));
3279566063dSJacob Faibussowitsch           PetscCall(PetscDrawPopCurrentPoint(draw));
328b4375e8dSPeter Brune         } else {
329b4375e8dSPeter Brune           w = 0.5 * PetscMin(1.0 - x, x);
3309566063dSJacob Faibussowitsch           PetscCall(PetscDrawPushCurrentPoint(draw, x - w, bottom));
3319566063dSJacob Faibussowitsch           if (curfas->smoothd) PetscCall(SNESView(curfas->smoothd, viewer));
3329566063dSJacob Faibussowitsch           PetscCall(PetscDrawPopCurrentPoint(draw));
3339566063dSJacob Faibussowitsch           PetscCall(PetscDrawPushCurrentPoint(draw, x + w, bottom));
3349566063dSJacob Faibussowitsch           if (curfas->smoothu) PetscCall(SNESView(curfas->smoothu, viewer));
3359566063dSJacob Faibussowitsch           PetscCall(PetscDrawPopCurrentPoint(draw));
336b4375e8dSPeter Brune         }
337656ede7eSPeter Brune         /* this is totally bogus but we have no way of knowing how low the previous one was draw to */
338656ede7eSPeter Brune         bottom -= 5 * th;
3391aa26658SKarl Rupp         if (curfas->next) curfas = (SNES_FAS *)curfas->next->data;
3400298fd71SBarry Smith         else curfas = NULL;
341656ede7eSPeter Brune       }
342421d9b32SPeter Brune     }
343ab8d36c9SPeter Brune   }
3443ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
345421d9b32SPeter Brune }
346421d9b32SPeter Brune 
34739bd7f45SPeter Brune /*
34839bd7f45SPeter Brune Defines the action of the downsmoother
34939bd7f45SPeter Brune  */
SNESFASDownSmooth_Private(SNES snes,Vec B,Vec X,Vec F,PetscReal * fnorm)350d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
351d71ae5a4SJacob Faibussowitsch {
352742fe5e2SPeter Brune   SNESConvergedReason reason;
353ab8d36c9SPeter Brune   Vec                 FPC;
354ab8d36c9SPeter Brune   SNES                smoothd;
3556cbb2f26SLawrence Mitchell   PetscBool           flg;
3560dd27c6cSPeter Brune   SNES_FAS           *fas = (SNES_FAS *)snes->data;
3576e111a19SKarl Rupp 
358421d9b32SPeter Brune   PetscFunctionBegin;
3599566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleGetSmootherDown(snes, &smoothd));
3609566063dSJacob Faibussowitsch   PetscCall(SNESSetInitialFunction(smoothd, F));
3619566063dSJacob Faibussowitsch   if (fas->eventsmoothsolve) PetscCall(PetscLogEventBegin(fas->eventsmoothsolve, smoothd, B, X, 0));
3629566063dSJacob Faibussowitsch   PetscCall(SNESSolve(smoothd, B, X));
3639566063dSJacob Faibussowitsch   if (fas->eventsmoothsolve) PetscCall(PetscLogEventEnd(fas->eventsmoothsolve, smoothd, B, X, 0));
364742fe5e2SPeter Brune   /* check convergence reason for the smoother */
3659566063dSJacob Faibussowitsch   PetscCall(SNESGetConvergedReason(smoothd, &reason));
3661ff805c3SPeter Brune   if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
367742fe5e2SPeter Brune     snes->reason = SNES_DIVERGED_INNER;
3683ba16761SJacob Faibussowitsch     PetscFunctionReturn(PETSC_SUCCESS);
369742fe5e2SPeter Brune   }
3706cbb2f26SLawrence Mitchell 
3719566063dSJacob Faibussowitsch   PetscCall(SNESGetFunction(smoothd, &FPC, NULL, NULL));
3729566063dSJacob Faibussowitsch   PetscCall(SNESGetAlwaysComputesFinalResidual(smoothd, &flg));
37348a46eb9SPierre Jolivet   if (!flg) PetscCall(SNESComputeFunction(smoothd, X, FPC));
3749566063dSJacob Faibussowitsch   PetscCall(VecCopy(FPC, F));
375*76c63389SBarry Smith   if (fnorm) {
376*76c63389SBarry Smith     PetscCall(VecNorm(F, NORM_2, fnorm));
377*76c63389SBarry Smith     SNESCheckFunctionDomainError(snes, *fnorm);
378*76c63389SBarry Smith   }
3793ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
38039bd7f45SPeter Brune }
38139bd7f45SPeter Brune 
38239bd7f45SPeter Brune /*
38307144faaSPeter Brune Defines the action of the upsmoother
38439bd7f45SPeter Brune  */
SNESFASUpSmooth_Private(SNES snes,Vec B,Vec X,Vec F,PetscReal * fnorm)385d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESFASUpSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
386d71ae5a4SJacob Faibussowitsch {
38739bd7f45SPeter Brune   SNESConvergedReason reason;
388ab8d36c9SPeter Brune   Vec                 FPC;
389ab8d36c9SPeter Brune   SNES                smoothu;
3906cbb2f26SLawrence Mitchell   PetscBool           flg;
3910dd27c6cSPeter Brune   SNES_FAS           *fas = (SNES_FAS *)snes->data;
392ab8d36c9SPeter Brune 
3936e111a19SKarl Rupp   PetscFunctionBegin;
3949566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleGetSmootherUp(snes, &smoothu));
3959566063dSJacob Faibussowitsch   if (fas->eventsmoothsolve) PetscCall(PetscLogEventBegin(fas->eventsmoothsolve, smoothu, 0, 0, 0));
3969566063dSJacob Faibussowitsch   PetscCall(SNESSolve(smoothu, B, X));
3979566063dSJacob Faibussowitsch   if (fas->eventsmoothsolve) PetscCall(PetscLogEventEnd(fas->eventsmoothsolve, smoothu, 0, 0, 0));
39839bd7f45SPeter Brune   /* check convergence reason for the smoother */
3999566063dSJacob Faibussowitsch   PetscCall(SNESGetConvergedReason(smoothu, &reason));
4001ff805c3SPeter Brune   if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
40139bd7f45SPeter Brune     snes->reason = SNES_DIVERGED_INNER;
4023ba16761SJacob Faibussowitsch     PetscFunctionReturn(PETSC_SUCCESS);
40339bd7f45SPeter Brune   }
4049566063dSJacob Faibussowitsch   PetscCall(SNESGetFunction(smoothu, &FPC, NULL, NULL));
4059566063dSJacob Faibussowitsch   PetscCall(SNESGetAlwaysComputesFinalResidual(smoothu, &flg));
40648a46eb9SPierre Jolivet   if (!flg) PetscCall(SNESComputeFunction(smoothu, X, FPC));
4079566063dSJacob Faibussowitsch   PetscCall(VecCopy(FPC, F));
408*76c63389SBarry Smith   if (fnorm) {
409*76c63389SBarry Smith     PetscCall(VecNorm(F, NORM_2, fnorm));
410*76c63389SBarry Smith     SNESCheckFunctionDomainError(snes, *fnorm);
411*76c63389SBarry Smith   }
4123ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
41339bd7f45SPeter Brune }
41439bd7f45SPeter Brune 
415938e4a01SJed Brown /*@
416420bcc1bSBarry Smith   SNESFASCreateCoarseVec - create a `Vec` corresponding to a state vector on one level coarser than the current level
417938e4a01SJed Brown 
418938e4a01SJed Brown   Collective
419938e4a01SJed Brown 
4204165533cSJose E. Roman   Input Parameter:
421f6dfbefdSBarry Smith . snes - `SNESFAS` object
422938e4a01SJed Brown 
4234165533cSJose E. Roman   Output Parameter:
424420bcc1bSBarry Smith . Xcoarse - vector on level one coarser than the current level
425938e4a01SJed Brown 
426938e4a01SJed Brown   Level: developer
427938e4a01SJed Brown 
428420bcc1bSBarry Smith .seealso: [](ch_snes), `SNESFASSetRestriction()`, `SNESFASRestrict()`, `SNESFAS`
429938e4a01SJed Brown @*/
SNESFASCreateCoarseVec(SNES snes,Vec * Xcoarse)430d71ae5a4SJacob Faibussowitsch PetscErrorCode SNESFASCreateCoarseVec(SNES snes, Vec *Xcoarse)
431d71ae5a4SJacob Faibussowitsch {
432f833ba53SLisandro Dalcin   SNES_FAS *fas;
433938e4a01SJed Brown 
434938e4a01SJed Brown   PetscFunctionBegin;
435f833ba53SLisandro Dalcin   PetscValidHeaderSpecificType(snes, SNES_CLASSID, 1, SNESFAS);
4364f572ea9SToby Isaac   PetscAssertPointer(Xcoarse, 2);
437f833ba53SLisandro Dalcin   fas = (SNES_FAS *)snes->data;
4381aa26658SKarl Rupp   if (fas->rscale) {
4399566063dSJacob Faibussowitsch     PetscCall(VecDuplicate(fas->rscale, Xcoarse));
44090b9e4c9SPablo Brubeck   } else if (fas->interpolate) {
44190b9e4c9SPablo Brubeck     PetscCall(MatCreateVecs(fas->interpolate, Xcoarse, NULL));
44290b9e4c9SPablo Brubeck   } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must set rscale or interpolation");
4433ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
444938e4a01SJed Brown }
445938e4a01SJed Brown 
446e9923e8dSJed Brown /*@
447f6dfbefdSBarry Smith   SNESFASRestrict - restrict a `Vec` to the next coarser level
448e9923e8dSJed Brown 
449e9923e8dSJed Brown   Collective
450e9923e8dSJed Brown 
4514165533cSJose E. Roman   Input Parameters:
452f6dfbefdSBarry Smith + fine  - `SNES` from which to restrict
453e9923e8dSJed Brown - Xfine - vector to restrict
454e9923e8dSJed Brown 
4554165533cSJose E. Roman   Output Parameter:
456e9923e8dSJed Brown . Xcoarse - result of restriction
457e9923e8dSJed Brown 
458e9923e8dSJed Brown   Level: developer
459e9923e8dSJed Brown 
460420bcc1bSBarry Smith .seealso: [](ch_snes), `SNES`, `SNESFAS`, `SNESFASSetRestriction()`, `SNESFASSetInjection()`, `SNESFASCreateCoarseVec()`
461e9923e8dSJed Brown @*/
SNESFASRestrict(SNES fine,Vec Xfine,Vec Xcoarse)462d71ae5a4SJacob Faibussowitsch PetscErrorCode SNESFASRestrict(SNES fine, Vec Xfine, Vec Xcoarse)
463d71ae5a4SJacob Faibussowitsch {
464f833ba53SLisandro Dalcin   SNES_FAS *fas;
465e9923e8dSJed Brown 
466e9923e8dSJed Brown   PetscFunctionBegin;
467f833ba53SLisandro Dalcin   PetscValidHeaderSpecificType(fine, SNES_CLASSID, 1, SNESFAS);
468e9923e8dSJed Brown   PetscValidHeaderSpecific(Xfine, VEC_CLASSID, 2);
469e9923e8dSJed Brown   PetscValidHeaderSpecific(Xcoarse, VEC_CLASSID, 3);
470f833ba53SLisandro Dalcin   fas = (SNES_FAS *)fine->data;
471e9923e8dSJed Brown   if (fas->inject) {
4729566063dSJacob Faibussowitsch     PetscCall(MatRestrict(fas->inject, Xfine, Xcoarse));
473e9923e8dSJed Brown   } else {
4749566063dSJacob Faibussowitsch     PetscCall(MatRestrict(fas->restrct, Xfine, Xcoarse));
4759566063dSJacob Faibussowitsch     PetscCall(VecPointwiseMult(Xcoarse, fas->rscale, Xcoarse));
476e9923e8dSJed Brown   }
4773ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
478e9923e8dSJed Brown }
479e9923e8dSJed Brown 
48039bd7f45SPeter Brune /*
4816dfbafefSToby Isaac Performs a variant of FAS using the interpolated total coarse solution
4826dfbafefSToby Isaac 
4836dfbafefSToby Isaac fine problem:   F(x) = b
4846dfbafefSToby Isaac coarse problem: F^c(x^c) = Rb, Initial guess Rx
4856dfbafefSToby Isaac interpolated solution: x^f = I x^c (total solution interpolation
4866dfbafefSToby Isaac  */
SNESFASInterpolatedCoarseSolution(SNES snes,Vec X,Vec X_new)487d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESFASInterpolatedCoarseSolution(SNES snes, Vec X, Vec X_new)
488d71ae5a4SJacob Faibussowitsch {
4896dfbafefSToby Isaac   Vec                 X_c, B_c;
4906dfbafefSToby Isaac   SNESConvergedReason reason;
4916dfbafefSToby Isaac   SNES                next;
4926dfbafefSToby Isaac   Mat                 restrct, interpolate;
4936dfbafefSToby Isaac   SNES_FAS           *fasc;
4946dfbafefSToby Isaac 
4956dfbafefSToby Isaac   PetscFunctionBegin;
4969566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleGetCorrection(snes, &next));
4976dfbafefSToby Isaac   if (next) {
4986dfbafefSToby Isaac     fasc = (SNES_FAS *)next->data;
4996dfbafefSToby Isaac 
5009566063dSJacob Faibussowitsch     PetscCall(SNESFASCycleGetRestriction(snes, &restrct));
5019566063dSJacob Faibussowitsch     PetscCall(SNESFASCycleGetInterpolation(snes, &interpolate));
5026dfbafefSToby Isaac 
5036dfbafefSToby Isaac     X_c = next->vec_sol;
5046dfbafefSToby Isaac 
5059566063dSJacob Faibussowitsch     if (fasc->eventinterprestrict) PetscCall(PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0));
5066dfbafefSToby Isaac     /* restrict the total solution: Rb */
5079566063dSJacob Faibussowitsch     PetscCall(SNESFASRestrict(snes, X, X_c));
5086dfbafefSToby Isaac     B_c = next->vec_rhs;
5096dfbafefSToby Isaac     if (snes->vec_rhs) {
5106dfbafefSToby Isaac       /* restrict the total rhs defect: Rb */
5119566063dSJacob Faibussowitsch       PetscCall(MatRestrict(restrct, snes->vec_rhs, B_c));
5126dfbafefSToby Isaac     } else {
5139566063dSJacob Faibussowitsch       PetscCall(VecSet(B_c, 0.));
5146dfbafefSToby Isaac     }
5159566063dSJacob Faibussowitsch     if (fasc->eventinterprestrict) PetscCall(PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0));
5166dfbafefSToby Isaac 
5179566063dSJacob Faibussowitsch     PetscCall(SNESSolve(next, B_c, X_c));
5189566063dSJacob Faibussowitsch     PetscCall(SNESGetConvergedReason(next, &reason));
5196dfbafefSToby Isaac     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
5206dfbafefSToby Isaac       snes->reason = SNES_DIVERGED_INNER;
5213ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
5226dfbafefSToby Isaac     }
5236dfbafefSToby Isaac     /* x^f <- Ix^c*/
5246dfbafefSToby Isaac     DM dmc, dmf;
5256dfbafefSToby Isaac 
5269566063dSJacob Faibussowitsch     PetscCall(SNESGetDM(next, &dmc));
5279566063dSJacob Faibussowitsch     PetscCall(SNESGetDM(snes, &dmf));
5289566063dSJacob Faibussowitsch     if (fasc->eventinterprestrict) PetscCall(PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0));
5299566063dSJacob Faibussowitsch     PetscCall(DMInterpolateSolution(dmc, dmf, interpolate, X_c, X_new));
5309566063dSJacob Faibussowitsch     if (fasc->eventinterprestrict) PetscCall(PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0));
5319566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetName((PetscObject)X_c, "Coarse solution"));
5329566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view"));
5339566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetName((PetscObject)X_new, "Updated Fine solution"));
5349566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view"));
5356dfbafefSToby Isaac   }
5363ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
5376dfbafefSToby Isaac }
5386dfbafefSToby Isaac 
5396dfbafefSToby Isaac /*
54039bd7f45SPeter Brune Performs the FAS coarse correction as:
54139bd7f45SPeter Brune 
542b5527d98SMatthew G. Knepley fine problem:   F(x) = b
543b5527d98SMatthew G. Knepley coarse problem: F^c(x^c) = b^c
54439bd7f45SPeter Brune 
545b5527d98SMatthew G. Knepley b^c = F^c(Rx) - R(F(x) - b)
54639bd7f45SPeter Brune  */
SNESFASCoarseCorrection(SNES snes,Vec X,Vec F,Vec X_new)54766976f2fSJacob Faibussowitsch static PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new)
548d71ae5a4SJacob Faibussowitsch {
54939bd7f45SPeter Brune   Vec                 X_c, Xo_c, F_c, B_c;
55039bd7f45SPeter Brune   SNESConvergedReason reason;
551ab8d36c9SPeter Brune   SNES                next;
552ab8d36c9SPeter Brune   Mat                 restrct, interpolate;
5530dd27c6cSPeter Brune   SNES_FAS           *fasc;
5545fd66863SKarl Rupp 
55539bd7f45SPeter Brune   PetscFunctionBegin;
5569566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleGetCorrection(snes, &next));
557ab8d36c9SPeter Brune   if (next) {
5580dd27c6cSPeter Brune     fasc = (SNES_FAS *)next->data;
5590dd27c6cSPeter Brune 
5609566063dSJacob Faibussowitsch     PetscCall(SNESFASCycleGetRestriction(snes, &restrct));
5619566063dSJacob Faibussowitsch     PetscCall(SNESFASCycleGetInterpolation(snes, &interpolate));
562ab8d36c9SPeter Brune 
563ab8d36c9SPeter Brune     X_c  = next->vec_sol;
564ab8d36c9SPeter Brune     Xo_c = next->work[0];
565ab8d36c9SPeter Brune     F_c  = next->vec_func;
566ab8d36c9SPeter Brune     B_c  = next->vec_rhs;
567efe1f98aSPeter Brune 
5689566063dSJacob Faibussowitsch     if (fasc->eventinterprestrict) PetscCall(PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0));
5699566063dSJacob Faibussowitsch     PetscCall(SNESFASRestrict(snes, X, Xo_c));
5705609cbf2SMatthew G. Knepley     /* restrict the defect: R(F(x) - b) */
5719566063dSJacob Faibussowitsch     PetscCall(MatRestrict(restrct, F, B_c));
5729566063dSJacob Faibussowitsch     if (fasc->eventinterprestrict) PetscCall(PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0));
5730dd27c6cSPeter Brune 
5749566063dSJacob Faibussowitsch     if (fasc->eventresidual) PetscCall(PetscLogEventBegin(fasc->eventresidual, next, 0, 0, 0));
5755609cbf2SMatthew G. Knepley     /* F_c = F^c(Rx) - R(F(x) - b) since the second term was sitting in next->vec_rhs */
5769566063dSJacob Faibussowitsch     PetscCall(SNESComputeFunction(next, Xo_c, F_c));
5779566063dSJacob Faibussowitsch     if (fasc->eventresidual) PetscCall(PetscLogEventEnd(fasc->eventresidual, next, 0, 0, 0));
5780dd27c6cSPeter Brune 
5790dd27c6cSPeter Brune     /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */
5809566063dSJacob Faibussowitsch     PetscCall(VecCopy(B_c, X_c));
5819566063dSJacob Faibussowitsch     PetscCall(VecCopy(F_c, B_c));
5829566063dSJacob Faibussowitsch     PetscCall(VecCopy(X_c, F_c));
583ee78dd50SPeter Brune     /* set initial guess of the coarse problem to the projected fine solution */
5849566063dSJacob Faibussowitsch     PetscCall(VecCopy(Xo_c, X_c));
585c90fad12SPeter Brune 
586c90fad12SPeter Brune     /* recurse to the next level */
5879566063dSJacob Faibussowitsch     PetscCall(SNESSetInitialFunction(next, F_c));
5889566063dSJacob Faibussowitsch     PetscCall(SNESSolve(next, B_c, X_c));
5899566063dSJacob Faibussowitsch     PetscCall(SNESGetConvergedReason(next, &reason));
590742fe5e2SPeter Brune     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
591742fe5e2SPeter Brune       snes->reason = SNES_DIVERGED_INNER;
5923ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
593742fe5e2SPeter Brune     }
594fa9694d7SPeter Brune     /* correct as x <- x + I(x^c - Rx)*/
5959566063dSJacob Faibussowitsch     PetscCall(VecAXPY(X_c, -1.0, Xo_c));
5960dd27c6cSPeter Brune 
5979566063dSJacob Faibussowitsch     if (fasc->eventinterprestrict) PetscCall(PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0));
5989566063dSJacob Faibussowitsch     PetscCall(MatInterpolateAdd(interpolate, X_c, X, X_new));
5999566063dSJacob Faibussowitsch     if (fasc->eventinterprestrict) PetscCall(PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0));
6009566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetName((PetscObject)X_c, "Coarse correction"));
6019566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view"));
6029566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetName((PetscObject)X_new, "Updated Fine solution"));
6039566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view"));
604293a7e31SPeter Brune   }
6053ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
60639bd7f45SPeter Brune }
60739bd7f45SPeter Brune 
60839bd7f45SPeter Brune /*
609420bcc1bSBarry Smith The additive cycle is:
61039bd7f45SPeter Brune 
61107144faaSPeter Brune xhat = x
61207144faaSPeter Brune xhat = dS(x, b)
61307144faaSPeter Brune x = coarsecorrection(xhat, b_d)
61407144faaSPeter Brune x = x + nu*(xhat - x);
61539bd7f45SPeter Brune (optional) x = uS(x, b)
61639bd7f45SPeter Brune 
61739bd7f45SPeter Brune With the coarse RHS (defect correction) as below.
61839bd7f45SPeter Brune  */
SNESFASCycle_Additive(SNES snes,Vec X)619d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X)
620d71ae5a4SJacob Faibussowitsch {
62107144faaSPeter Brune   Vec                 F, B, Xhat;
62222c1e704SPeter Brune   Vec                 X_c, Xo_c, F_c, B_c;
62307144faaSPeter Brune   SNESConvergedReason reason;
62422c1e704SPeter Brune   PetscReal           xnorm, fnorm, ynorm;
625ab8d36c9SPeter Brune   SNES                next;
626ab8d36c9SPeter Brune   Mat                 restrct, interpolate;
6270dd27c6cSPeter Brune   SNES_FAS           *fas = (SNES_FAS *)snes->data, *fasc;
6280dd27c6cSPeter Brune 
62939bd7f45SPeter Brune   PetscFunctionBegin;
6309566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleGetCorrection(snes, &next));
63139bd7f45SPeter Brune   F    = snes->vec_func;
63239bd7f45SPeter Brune   B    = snes->vec_rhs;
633e7f468e7SPeter Brune   Xhat = snes->work[1];
6349566063dSJacob Faibussowitsch   PetscCall(VecCopy(X, Xhat));
63507144faaSPeter Brune   /* recurse first */
636ab8d36c9SPeter Brune   if (next) {
6370dd27c6cSPeter Brune     fasc = (SNES_FAS *)next->data;
6389566063dSJacob Faibussowitsch     PetscCall(SNESFASCycleGetRestriction(snes, &restrct));
6399566063dSJacob Faibussowitsch     PetscCall(SNESFASCycleGetInterpolation(snes, &interpolate));
6409566063dSJacob Faibussowitsch     if (fas->eventresidual) PetscCall(PetscLogEventBegin(fas->eventresidual, snes, 0, 0, 0));
6419566063dSJacob Faibussowitsch     PetscCall(SNESComputeFunction(snes, Xhat, F));
6429566063dSJacob Faibussowitsch     if (fas->eventresidual) PetscCall(PetscLogEventEnd(fas->eventresidual, snes, 0, 0, 0));
6439566063dSJacob Faibussowitsch     PetscCall(VecNorm(F, NORM_2, &fnorm));
644*76c63389SBarry Smith     SNESCheckFunctionDomainError(snes, fnorm);
645ab8d36c9SPeter Brune     X_c  = next->vec_sol;
646ab8d36c9SPeter Brune     Xo_c = next->work[0];
647ab8d36c9SPeter Brune     F_c  = next->vec_func;
648ab8d36c9SPeter Brune     B_c  = next->vec_rhs;
64939bd7f45SPeter Brune 
6509566063dSJacob Faibussowitsch     PetscCall(SNESFASRestrict(snes, Xhat, Xo_c));
65107144faaSPeter Brune     /* restrict the defect */
6529566063dSJacob Faibussowitsch     PetscCall(MatRestrict(restrct, F, B_c));
65307144faaSPeter Brune 
65407144faaSPeter Brune     /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */
6559566063dSJacob Faibussowitsch     if (fasc->eventresidual) PetscCall(PetscLogEventBegin(fasc->eventresidual, next, 0, 0, 0));
6569566063dSJacob Faibussowitsch     PetscCall(SNESComputeFunction(next, Xo_c, F_c));
6579566063dSJacob Faibussowitsch     if (fasc->eventresidual) PetscCall(PetscLogEventEnd(fasc->eventresidual, next, 0, 0, 0));
6589566063dSJacob Faibussowitsch     PetscCall(VecCopy(B_c, X_c));
6599566063dSJacob Faibussowitsch     PetscCall(VecCopy(F_c, B_c));
6609566063dSJacob Faibussowitsch     PetscCall(VecCopy(X_c, F_c));
66107144faaSPeter Brune     /* set initial guess of the coarse problem to the projected fine solution */
6629566063dSJacob Faibussowitsch     PetscCall(VecCopy(Xo_c, X_c));
66307144faaSPeter Brune 
66407144faaSPeter Brune     /* recurse */
6659566063dSJacob Faibussowitsch     PetscCall(SNESSetInitialFunction(next, F_c));
6669566063dSJacob Faibussowitsch     PetscCall(SNESSolve(next, B_c, X_c));
66707144faaSPeter Brune 
66807144faaSPeter Brune     /* smooth on this level */
6699566063dSJacob Faibussowitsch     PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &fnorm));
67007144faaSPeter Brune 
6719566063dSJacob Faibussowitsch     PetscCall(SNESGetConvergedReason(next, &reason));
67207144faaSPeter Brune     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
67307144faaSPeter Brune       snes->reason = SNES_DIVERGED_INNER;
6743ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
67507144faaSPeter Brune     }
67607144faaSPeter Brune 
67707144faaSPeter Brune     /* correct as x <- x + I(x^c - Rx)*/
6789566063dSJacob Faibussowitsch     PetscCall(VecAYPX(X_c, -1.0, Xo_c));
6799566063dSJacob Faibussowitsch     PetscCall(MatInterpolate(interpolate, X_c, Xhat));
68007144faaSPeter Brune 
681ddebd997SPeter Brune     /* additive correction of the coarse direction*/
6829566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat));
6839566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm));
684*76c63389SBarry Smith     SNESCheckLineSearchFailure(snes);
68507144faaSPeter Brune   } else {
6869566063dSJacob Faibussowitsch     PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm));
68707144faaSPeter Brune   }
6883ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
68939bd7f45SPeter Brune }
69039bd7f45SPeter Brune 
69139bd7f45SPeter Brune /*
69239bd7f45SPeter Brune Defines the FAS cycle as:
69339bd7f45SPeter Brune 
694b5527d98SMatthew G. Knepley fine problem: F(x) = b
69539bd7f45SPeter Brune coarse problem: F^c(x) = b^c
69639bd7f45SPeter Brune 
697b5527d98SMatthew G. Knepley b^c = F^c(Rx) - R(F(x) - b)
69839bd7f45SPeter Brune 
69939bd7f45SPeter Brune correction:
70039bd7f45SPeter Brune 
70139bd7f45SPeter Brune x = x + I(x^c - Rx)
70239bd7f45SPeter Brune  */
SNESFASCycle_Multiplicative(SNES snes,Vec X)703d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X)
704d71ae5a4SJacob Faibussowitsch {
70539bd7f45SPeter Brune   Vec  F, B;
70634d65b3cSPeter Brune   SNES next;
70739bd7f45SPeter Brune 
70839bd7f45SPeter Brune   PetscFunctionBegin;
70939bd7f45SPeter Brune   F = snes->vec_func;
71039bd7f45SPeter Brune   B = snes->vec_rhs;
71139bd7f45SPeter Brune   /* pre-smooth -- just update using the pre-smoother */
7129566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleGetCorrection(snes, &next));
7139566063dSJacob Faibussowitsch   PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm));
71434d65b3cSPeter Brune   if (next) {
7159566063dSJacob Faibussowitsch     PetscCall(SNESFASCoarseCorrection(snes, X, F, X));
7169566063dSJacob Faibussowitsch     PetscCall(SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm));
717fe6f9142SPeter Brune   }
7183ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
719421d9b32SPeter Brune }
720421d9b32SPeter Brune 
SNESFASCycleSetupPhase_Full(SNES snes)721d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESFASCycleSetupPhase_Full(SNES snes)
722d71ae5a4SJacob Faibussowitsch {
7238c94862eSPeter Brune   SNES      next;
7248c94862eSPeter Brune   SNES_FAS *fas = (SNES_FAS *)snes->data;
7258c94862eSPeter Brune   PetscBool isFine;
7268c94862eSPeter Brune 
7278c94862eSPeter Brune   PetscFunctionBegin;
7288c94862eSPeter Brune   /* pre-smooth -- just update using the pre-smoother */
7299566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleIsFine(snes, &isFine));
7309566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleGetCorrection(snes, &next));
7318c94862eSPeter Brune   fas->full_stage = 0;
7329566063dSJacob Faibussowitsch   if (next) PetscCall(SNESFASCycleSetupPhase_Full(next));
7333ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
7348c94862eSPeter Brune }
7358c94862eSPeter Brune 
SNESFASCycle_Full(SNES snes,Vec X)736d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESFASCycle_Full(SNES snes, Vec X)
737d71ae5a4SJacob Faibussowitsch {
7388c94862eSPeter Brune   Vec       F, B;
7398c94862eSPeter Brune   SNES_FAS *fas = (SNES_FAS *)snes->data;
7408c94862eSPeter Brune   PetscBool isFine;
7418c94862eSPeter Brune   SNES      next;
7428c94862eSPeter Brune 
7438c94862eSPeter Brune   PetscFunctionBegin;
7448c94862eSPeter Brune   F = snes->vec_func;
7458c94862eSPeter Brune   B = snes->vec_rhs;
7469566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleIsFine(snes, &isFine));
7479566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleGetCorrection(snes, &next));
7488c94862eSPeter Brune 
7491baa6e33SBarry Smith   if (isFine) PetscCall(SNESFASCycleSetupPhase_Full(snes));
7508c94862eSPeter Brune 
7518c94862eSPeter Brune   if (fas->full_stage == 0) {
752928e959bSPeter Brune     /* downsweep */
7538c94862eSPeter Brune     if (next) {
7548c94862eSPeter Brune       if (fas->level != 1) next->max_its += 1;
7559566063dSJacob Faibussowitsch       if (fas->full_downsweep) PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm));
756e140b65aSPatrick Farrell       fas->full_downsweep = PETSC_TRUE;
7579566063dSJacob Faibussowitsch       if (fas->full_total) PetscCall(SNESFASInterpolatedCoarseSolution(snes, X, X));
7589566063dSJacob Faibussowitsch       else PetscCall(SNESFASCoarseCorrection(snes, X, F, X));
7596dfbafefSToby Isaac       fas->full_total = PETSC_FALSE;
7609566063dSJacob Faibussowitsch       PetscCall(SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm));
7618c94862eSPeter Brune       if (fas->level != 1) next->max_its -= 1;
7628c94862eSPeter Brune     } else {
763a3a80b83SMatthew G. Knepley       /* The smoother on the coarse level is the coarse solver */
7649566063dSJacob Faibussowitsch       PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm));
7658c94862eSPeter Brune     }
7668c94862eSPeter Brune     fas->full_stage = 1;
7678c94862eSPeter Brune   } else if (fas->full_stage == 1) {
7689566063dSJacob Faibussowitsch     if (snes->iter == 0) PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm));
7698c94862eSPeter Brune     if (next) {
7709566063dSJacob Faibussowitsch       PetscCall(SNESFASCoarseCorrection(snes, X, F, X));
7719566063dSJacob Faibussowitsch       PetscCall(SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm));
7728c94862eSPeter Brune     }
7738c94862eSPeter Brune   }
7748c94862eSPeter Brune   /* final v-cycle */
7758c94862eSPeter Brune   if (isFine) {
7768c94862eSPeter Brune     if (next) {
7779566063dSJacob Faibussowitsch       PetscCall(SNESFASCoarseCorrection(snes, X, F, X));
7789566063dSJacob Faibussowitsch       PetscCall(SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm));
7798c94862eSPeter Brune     }
7808c94862eSPeter Brune   }
7813ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
7828c94862eSPeter Brune }
7838c94862eSPeter Brune 
SNESFASCycle_Kaskade(SNES snes,Vec X)784d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESFASCycle_Kaskade(SNES snes, Vec X)
785d71ae5a4SJacob Faibussowitsch {
78634d65b3cSPeter Brune   Vec  F, B;
78734d65b3cSPeter Brune   SNES next;
78834d65b3cSPeter Brune 
78934d65b3cSPeter Brune   PetscFunctionBegin;
79034d65b3cSPeter Brune   F = snes->vec_func;
79134d65b3cSPeter Brune   B = snes->vec_rhs;
7929566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleGetCorrection(snes, &next));
79334d65b3cSPeter Brune   if (next) {
7949566063dSJacob Faibussowitsch     PetscCall(SNESFASCoarseCorrection(snes, X, F, X));
7959566063dSJacob Faibussowitsch     PetscCall(SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm));
79634d65b3cSPeter Brune   } else {
7979566063dSJacob Faibussowitsch     PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm));
79834d65b3cSPeter Brune   }
7993ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
80034d65b3cSPeter Brune }
80134d65b3cSPeter Brune 
802fffbeea8SBarry Smith PetscBool  SNEScite       = PETSC_FALSE;
8031d27aa22SBarry Smith const char SNESCitation[] = "@Article{bruneknepleysmithtu15,"
8041d27aa22SBarry Smith                             "  title         = {Composing Scalable Nonlinear Algebraic Solvers},"
8051d27aa22SBarry Smith                             "  author        = {Peter R. Brune and Matthew G. Knepley and Barry F. Smith and Xuemin Tu},"
8061d27aa22SBarry Smith                             "  journal       = {SIAM Review},"
8071d27aa22SBarry Smith                             "  volume        = {57},"
8081d27aa22SBarry Smith                             "  number        = {4},"
8091d27aa22SBarry Smith                             "  pages         = {535--565},"
8101d27aa22SBarry Smith                             "  doi           = {10.1137/130936725},"
8111d27aa22SBarry Smith                             "  year          = {2015}\n";
812fffbeea8SBarry Smith 
SNESSolve_FAS(SNES snes)813d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESSolve_FAS(SNES snes)
814d71ae5a4SJacob Faibussowitsch {
815f833ba53SLisandro Dalcin   PetscInt  i;
816ddb5aff1SPeter Brune   Vec       X, F;
817fe6f9142SPeter Brune   PetscReal fnorm;
818b17ce1afSJed Brown   SNES_FAS *fas = (SNES_FAS *)snes->data, *ffas;
819b17ce1afSJed Brown   DM        dm;
820e70c42e5SPeter Brune   PetscBool isFine;
821b17ce1afSJed Brown 
822421d9b32SPeter Brune   PetscFunctionBegin;
8230b121fc5SBarry Smith   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);
824c579b300SPatrick Farrell 
8259566063dSJacob Faibussowitsch   PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
826fe6f9142SPeter Brune   snes->reason = SNES_CONVERGED_ITERATING;
827fa9694d7SPeter Brune   X            = snes->vec_sol;
828f5a6d4f9SBarry Smith   F            = snes->vec_func;
829293a7e31SPeter Brune 
8309566063dSJacob Faibussowitsch   PetscCall(SNESFASCycleIsFine(snes, &isFine));
831293a7e31SPeter Brune   /* norm setup */
8329566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
833fe6f9142SPeter Brune   snes->iter = 0;
834f833ba53SLisandro Dalcin   snes->norm = 0;
8359566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
836e4ed7901SPeter Brune   if (!snes->vec_func_init_set) {
8379566063dSJacob Faibussowitsch     if (fas->eventresidual) PetscCall(PetscLogEventBegin(fas->eventresidual, snes, 0, 0, 0));
8389566063dSJacob Faibussowitsch     PetscCall(SNESComputeFunction(snes, X, F));
8399566063dSJacob Faibussowitsch     if (fas->eventresidual) PetscCall(PetscLogEventEnd(fas->eventresidual, snes, 0, 0, 0));
8401aa26658SKarl Rupp   } else snes->vec_func_init_set = PETSC_FALSE;
841e4ed7901SPeter Brune 
8429566063dSJacob Faibussowitsch   PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F||  */
843*76c63389SBarry Smith   SNESCheckFunctionDomainError(snes, fnorm);
8449566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
845fe6f9142SPeter Brune   snes->norm = fnorm;
8469566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
8479566063dSJacob Faibussowitsch   PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
848fe6f9142SPeter Brune 
849fe6f9142SPeter Brune   /* test convergence */
8502d157150SStefano Zampini   PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
8512d157150SStefano Zampini   PetscCall(SNESMonitor(snes, snes->iter, fnorm));
8523ba16761SJacob Faibussowitsch   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
853e4ed7901SPeter Brune 
854b9c2fdf1SPeter Brune   if (isFine) {
855b9c2fdf1SPeter Brune     /* propagate scale-dependent data up the hierarchy */
8569566063dSJacob Faibussowitsch     PetscCall(SNESGetDM(snes, &dm));
857b17ce1afSJed Brown     for (ffas = fas; ffas->next; ffas = (SNES_FAS *)ffas->next->data) {
858b17ce1afSJed Brown       DM dmcoarse;
8599566063dSJacob Faibussowitsch       PetscCall(SNESGetDM(ffas->next, &dmcoarse));
8609566063dSJacob Faibussowitsch       PetscCall(DMRestrict(dm, ffas->restrct, ffas->rscale, ffas->inject, dmcoarse));
861b17ce1afSJed Brown       dm = dmcoarse;
862b17ce1afSJed Brown     }
863b9c2fdf1SPeter Brune   }
864b17ce1afSJed Brown 
865f833ba53SLisandro Dalcin   for (i = 0; i < snes->max_its; i++) {
866fe6f9142SPeter Brune     /* Call general purpose update function */
867dbbe0bcdSBarry Smith     PetscTryTypeMethod(snes, update, snes->iter);
868f833ba53SLisandro Dalcin 
86907144faaSPeter Brune     if (fas->fastype == SNES_FAS_MULTIPLICATIVE) {
8709566063dSJacob Faibussowitsch       PetscCall(SNESFASCycle_Multiplicative(snes, X));
8718c94862eSPeter Brune     } else if (fas->fastype == SNES_FAS_ADDITIVE) {
8729566063dSJacob Faibussowitsch       PetscCall(SNESFASCycle_Additive(snes, X));
8738c94862eSPeter Brune     } else if (fas->fastype == SNES_FAS_FULL) {
8749566063dSJacob Faibussowitsch       PetscCall(SNESFASCycle_Full(snes, X));
87534d65b3cSPeter Brune     } else if (fas->fastype == SNES_FAS_KASKADE) {
8769566063dSJacob Faibussowitsch       PetscCall(SNESFASCycle_Kaskade(snes, X));
8776c4ed002SBarry Smith     } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Unsupported FAS type");
878742fe5e2SPeter Brune 
879742fe5e2SPeter Brune     /* check for FAS cycle divergence */
8803ba16761SJacob Faibussowitsch     if (snes->reason != SNES_CONVERGED_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
881b9c2fdf1SPeter Brune 
882c90fad12SPeter Brune     /* Monitor convergence */
8839566063dSJacob Faibussowitsch     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
884c90fad12SPeter Brune     snes->iter = i + 1;
8859566063dSJacob Faibussowitsch     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
8869566063dSJacob Faibussowitsch     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
8872d157150SStefano Zampini     PetscCall(SNESConverged(snes, snes->iter, 0.0, 0.0, snes->norm));
8889566063dSJacob Faibussowitsch     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
889c90fad12SPeter Brune     if (snes->reason) break;
890fe6f9142SPeter Brune   }
8913ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
892421d9b32SPeter Brune }
89340244768SBarry Smith 
89440244768SBarry Smith /*MC
8950b4b7b1cSBarry Smith    SNESFAS - An implementation of the Full Approximation Scheme nonlinear multigrid solver, FAS, or nonlinear multigrid {cite}`bruneknepleysmithtu15` for
8960b4b7b1cSBarry Smith              solving nonlinear systems of equations with `SNES`.
89740244768SBarry Smith 
89840244768SBarry Smith    The nonlinear problem is solved by correction using coarse versions
89940244768SBarry Smith    of the nonlinear problem.  This problem is perturbed so that a projected
90040244768SBarry Smith    solution of the fine problem elicits no correction from the coarse problem.
90140244768SBarry Smith 
902f6dfbefdSBarry Smith    Options Database Keys and Prefixes:
9031d27aa22SBarry Smith +   -snes_fas_levels <l>                                  - The number of levels
9041d27aa22SBarry Smith .   -snes_fas_cycles <c>                                  - The number of cycles -- 1 for V, 2 for W
90540244768SBarry Smith .   -snes_fas_type<additive,multiplicative,full,kaskade>  - Additive or multiplicative cycle
9061d27aa22SBarry Smith .   -snes_fas_galerkin <false,true>                       - Form coarse problems by projection back upon the fine problem
9071d27aa22SBarry Smith .   -snes_fas_smoothup <u>                                - The number of iterations of the post-smoother
9081d27aa22SBarry Smith .   -snes_fas_smoothdown <d>                              - The number of iterations of the pre-smoother
90940244768SBarry Smith .   -snes_fas_monitor                                     - Monitor progress of all of the levels
9101d27aa22SBarry Smith .   -snes_fas_full_downsweep <false,true>                 - call the downsmooth on the initial downsweep of full FAS
9111d27aa22SBarry Smith .   -fas_levels_snes_                                     - prefix for `SNES` options for all smoothers
9121d27aa22SBarry Smith .   -fas_levels_cycle_snes_                               - prefix for `SNES` options for all cycles
9131d27aa22SBarry Smith .   -fas_levels_i_snes_                                   - prefix `SNES` options for the smoothers on level i
9141d27aa22SBarry Smith .   -fas_levels_i_cycle_snes_                             - prefix for `SNES` options for the cycle on level i
9151d27aa22SBarry Smith -   -fas_coarse_snes_                                     - prefix for `SNES` options for the coarsest smoother
91640244768SBarry Smith 
917420bcc1bSBarry Smith    Level: beginner
918420bcc1bSBarry Smith 
919f6dfbefdSBarry Smith    Note:
9200b4b7b1cSBarry Smith    The organization of the `SNESFAS` solver is slightly different from the organization of `PCMG`
921f6dfbefdSBarry Smith    As each level has smoother `SNES` instances(down and potentially up) and a cycle `SNES` instance.
922f6dfbefdSBarry Smith    The cycle `SNES` instance may be used for monitoring convergence on a particular level.
92340244768SBarry Smith 
924420bcc1bSBarry Smith .seealso: [](ch_snes), `PCMG`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESType`, `SNESFASSetRestriction()`, `SNESFASSetInjection()`,
925420bcc1bSBarry Smith           `SNESFASFullGetTotal()`, `SNESFASSetType()`, `SNESFASGetType()`, `SNESFASSetLevels()`, `SNESFASGetLevels()`, `SNESFASGetCycleSNES()`,
926420bcc1bSBarry Smith           `SNESFASSetNumberSmoothUp()`, `SNESFASSetNumberSmoothDown()`, `SNESFASSetContinuation()`, `SNESFASSetCycles()`, `SNESFASSetMonitor()`,
927420bcc1bSBarry Smith           `SNESFASSetLog()`, `SNESFASCycleSetCycles()`, `SNESFASCycleGetSmoother()`, `SNESFASCycleGetSmootherUp()`, `SNESFASCycleGetSmootherDown()`,
928420bcc1bSBarry Smith           `SNESFASCycleGetCorrection()`, `SNESFASCycleGetInterpolation()`, `SNESFASCycleGetRestriction()`, `SNESFASCycleGetInjection()`,
929420bcc1bSBarry Smith           `SNESFASCycleGetRScale()`, `SNESFASCycleIsFine()`, `SNESFASSetInterpolation()`, `SNESFASGetInterpolation()`, `SNESFASSetRestriction()`,
930420bcc1bSBarry Smith           `SNESFASGetRestriction()`, `SNESFASSetInjection()`, `SNESFASGetInjection()`, `SNESFASSetRScale()`,`SNESFASGetSmoother()`,
931420bcc1bSBarry Smith           `SNESFASGetSmootherDown()`, `SNESFASGetSmootherUp()`, `SNESFASGetCoarseSolve()`, `SNESFASFullSetDownSweep()`, `SNESFASFullSetTotal()`, `SNESFASFullGetTotal()`
93240244768SBarry Smith M*/
93340244768SBarry Smith 
SNESCreate_FAS(SNES snes)934d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode SNESCreate_FAS(SNES snes)
935d71ae5a4SJacob Faibussowitsch {
93640244768SBarry Smith   SNES_FAS *fas;
93740244768SBarry Smith 
93840244768SBarry Smith   PetscFunctionBegin;
93940244768SBarry Smith   snes->ops->destroy        = SNESDestroy_FAS;
94040244768SBarry Smith   snes->ops->setup          = SNESSetUp_FAS;
94140244768SBarry Smith   snes->ops->setfromoptions = SNESSetFromOptions_FAS;
94240244768SBarry Smith   snes->ops->view           = SNESView_FAS;
94340244768SBarry Smith   snes->ops->solve          = SNESSolve_FAS;
94440244768SBarry Smith   snes->ops->reset          = SNESReset_FAS;
94540244768SBarry Smith 
94640244768SBarry Smith   snes->usesksp = PETSC_FALSE;
947efd4aadfSBarry Smith   snes->usesnpc = PETSC_FALSE;
94840244768SBarry Smith 
94977e5a1f9SBarry Smith   PetscCall(SNESParametersInitialize(snes));
95077e5a1f9SBarry Smith   PetscObjectParameterSetDefault(snes, max_funcs, 30000);
95177e5a1f9SBarry Smith   PetscObjectParameterSetDefault(snes, max_its, 10000);
95240244768SBarry Smith 
9534fc747eaSLawrence Mitchell   snes->alwayscomputesfinalresidual = PETSC_TRUE;
9544fc747eaSLawrence Mitchell 
9554dfa11a4SJacob Faibussowitsch   PetscCall(PetscNew(&fas));
95640244768SBarry Smith 
95740244768SBarry Smith   snes->data                  = (void *)fas;
95840244768SBarry Smith   fas->level                  = 0;
95940244768SBarry Smith   fas->levels                 = 1;
96040244768SBarry Smith   fas->n_cycles               = 1;
96140244768SBarry Smith   fas->max_up_it              = 1;
96240244768SBarry Smith   fas->max_down_it            = 1;
96340244768SBarry Smith   fas->smoothu                = NULL;
96440244768SBarry Smith   fas->smoothd                = NULL;
96540244768SBarry Smith   fas->next                   = NULL;
96640244768SBarry Smith   fas->previous               = NULL;
96740244768SBarry Smith   fas->fine                   = snes;
96840244768SBarry Smith   fas->interpolate            = NULL;
96940244768SBarry Smith   fas->restrct                = NULL;
97040244768SBarry Smith   fas->inject                 = NULL;
97140244768SBarry Smith   fas->usedmfornumberoflevels = PETSC_FALSE;
97240244768SBarry Smith   fas->fastype                = SNES_FAS_MULTIPLICATIVE;
97340244768SBarry Smith   fas->full_downsweep         = PETSC_FALSE;
9746dfbafefSToby Isaac   fas->full_total             = PETSC_FALSE;
97540244768SBarry Smith 
97640244768SBarry Smith   fas->eventsmoothsetup    = 0;
97740244768SBarry Smith   fas->eventsmoothsolve    = 0;
97840244768SBarry Smith   fas->eventresidual       = 0;
97940244768SBarry Smith   fas->eventinterprestrict = 0;
9803ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
98140244768SBarry Smith }
982