/* Defines the basic SNES object */ #include <../src/snes/impls/fas/fasimpls.h> /*I "petscsnes.h" I*/ const char *const SNESFASTypes[] = {"MULTIPLICATIVE", "ADDITIVE", "FULL", "KASKADE", "SNESFASType", "SNES_FAS", NULL}; static PetscErrorCode SNESReset_FAS(SNES snes) { SNES_FAS *fas = (SNES_FAS *)snes->data; PetscFunctionBegin; PetscCall(SNESDestroy(&fas->smoothu)); PetscCall(SNESDestroy(&fas->smoothd)); PetscCall(MatDestroy(&fas->inject)); PetscCall(MatDestroy(&fas->interpolate)); PetscCall(MatDestroy(&fas->restrct)); PetscCall(VecDestroy(&fas->rscale)); PetscCall(VecDestroy(&fas->Xg)); PetscCall(VecDestroy(&fas->Fg)); if (fas->next) PetscCall(SNESReset(fas->next)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode SNESDestroy_FAS(SNES snes) { SNES_FAS *fas = (SNES_FAS *)snes->data; PetscFunctionBegin; /* recursively resets and then destroys */ PetscCall(SNESReset_FAS(snes)); PetscCall(SNESDestroy(&fas->next)); PetscCall(PetscFree(fas)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode SNESFASSetUpLineSearch_Private(SNES snes, SNES smooth) { SNESLineSearch linesearch; SNESLineSearch slinesearch; void *lsprectx, *lspostctx; PetscErrorCode (*precheck)(SNESLineSearch, Vec, Vec, PetscBool *, void *); PetscErrorCode (*postcheck)(SNESLineSearch, Vec, Vec, Vec, PetscBool *, PetscBool *, void *); PetscFunctionBegin; if (!snes->linesearch) PetscFunctionReturn(PETSC_SUCCESS); PetscCall(SNESGetLineSearch(snes, &linesearch)); PetscCall(SNESGetLineSearch(smooth, &slinesearch)); PetscCall(SNESLineSearchGetPreCheck(linesearch, &precheck, &lsprectx)); PetscCall(SNESLineSearchGetPostCheck(linesearch, &postcheck, &lspostctx)); PetscCall(SNESLineSearchSetPreCheck(slinesearch, precheck, lsprectx)); PetscCall(SNESLineSearchSetPostCheck(slinesearch, postcheck, lspostctx)); PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode SNESFASCycleSetUpSmoother_Private(SNES snes, SNES smooth) { SNES_FAS *fas = (SNES_FAS *)snes->data; PetscFunctionBegin; PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)smooth)); PetscCall(SNESSetFromOptions(smooth)); PetscCall(SNESFASSetUpLineSearch_Private(snes, smooth)); PetscCall(PetscObjectReference((PetscObject)snes->vec_sol)); PetscCall(PetscObjectReference((PetscObject)snes->vec_sol_update)); PetscCall(PetscObjectReference((PetscObject)snes->vec_func)); smooth->vec_sol = snes->vec_sol; smooth->vec_sol_update = snes->vec_sol_update; smooth->vec_func = snes->vec_func; if (fas->eventsmoothsetup) PetscCall(PetscLogEventBegin(fas->eventsmoothsetup, smooth, 0, 0, 0)); PetscCall(SNESSetUp(smooth)); if (fas->eventsmoothsetup) PetscCall(PetscLogEventEnd(fas->eventsmoothsetup, smooth, 0, 0, 0)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode SNESSetUp_FAS(SNES snes) { SNES_FAS *fas = (SNES_FAS *)snes->data; PetscInt dm_levels; SNES next; PetscBool isFine, hasCreateRestriction, hasCreateInjection; PetscFunctionBegin; PetscCall(SNESFASCycleIsFine(snes, &isFine)); if (fas->usedmfornumberoflevels && isFine) { PetscCall(DMGetRefineLevel(snes->dm, &dm_levels)); dm_levels++; if (dm_levels > fas->levels) { /* reset the number of levels */ PetscCall(SNESFASSetLevels(snes, dm_levels, NULL)); PetscCall(SNESSetFromOptions(snes)); } } PetscCall(SNESFASCycleGetCorrection(snes, &next)); if (!isFine) snes->gridsequence = 0; /* no grid sequencing inside the multigrid hierarchy! */ PetscCall(SNESSetWorkVecs(snes, 2)); /* work vectors used for intergrid transfers */ /* set up the smoothers if they haven't already been set up */ if (!fas->smoothd) PetscCall(SNESFASCycleCreateSmoother_Private(snes, &fas->smoothd)); if (snes->dm) { /* set the smoother DMs properly */ if (fas->smoothu) PetscCall(SNESSetDM(fas->smoothu, snes->dm)); PetscCall(SNESSetDM(fas->smoothd, snes->dm)); /* construct EVERYTHING from the DM -- including the progressive set of smoothers */ if (next) { /* for now -- assume the DM and the evaluation functions have been set externally */ if (!next->dm) { PetscCall(DMCoarsen(snes->dm, PetscObjectComm((PetscObject)next), &next->dm)); PetscCall(SNESSetDM(next, next->dm)); } /* set the interpolation and restriction from the DM */ if (!fas->interpolate) { PetscCall(DMCreateInterpolation(next->dm, snes->dm, &fas->interpolate, &fas->rscale)); if (!fas->restrct) { PetscCall(DMHasCreateRestriction(next->dm, &hasCreateRestriction)); /* DM can create restrictions, use that */ if (hasCreateRestriction) { PetscCall(DMCreateRestriction(next->dm, snes->dm, &fas->restrct)); } else { PetscCall(PetscObjectReference((PetscObject)fas->interpolate)); fas->restrct = fas->interpolate; } } } /* set the injection from the DM */ if (!fas->inject) { PetscCall(DMHasCreateInjection(next->dm, &hasCreateInjection)); if (hasCreateInjection) PetscCall(DMCreateInjection(next->dm, snes->dm, &fas->inject)); } } } /*pass the smoother, function, and jacobian up to the next level if it's not user set already */ if (fas->galerkin) { if (next) PetscCall(SNESSetFunction(next, NULL, SNESFASGalerkinFunctionDefault, next)); if (fas->smoothd && fas->level != fas->levels - 1) PetscCall(SNESSetFunction(fas->smoothd, NULL, SNESFASGalerkinFunctionDefault, snes)); if (fas->smoothu && fas->level != fas->levels - 1) PetscCall(SNESSetFunction(fas->smoothu, NULL, SNESFASGalerkinFunctionDefault, snes)); } /* sets the down (pre) smoother's default norm and sets it from options */ if (fas->smoothd) { if (fas->level == 0) { PetscCall(SNESSetNormSchedule(fas->smoothd, SNES_NORM_ALWAYS)); } else { PetscCall(SNESSetNormSchedule(fas->smoothd, SNES_NORM_FINAL_ONLY)); } PetscCall(SNESFASCycleSetUpSmoother_Private(snes, fas->smoothd)); } /* sets the up (post) smoother's default norm and sets it from options */ if (fas->smoothu) { if (fas->level != fas->levels - 1) { PetscCall(SNESSetNormSchedule(fas->smoothu, SNES_NORM_NONE)); } else { PetscCall(SNESSetNormSchedule(fas->smoothu, SNES_NORM_FINAL_ONLY)); } PetscCall(SNESFASCycleSetUpSmoother_Private(snes, fas->smoothu)); } if (next) { /* gotta set up the solution vector for this to work */ if (!next->vec_sol) PetscCall(SNESFASCreateCoarseVec(snes, &next->vec_sol)); if (!next->vec_rhs) PetscCall(SNESFASCreateCoarseVec(snes, &next->vec_rhs)); PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next)); PetscCall(SNESFASSetUpLineSearch_Private(snes, next)); PetscCall(SNESSetUp(next)); } /* setup FAS work vectors */ if (fas->galerkin) { PetscCall(VecDuplicate(snes->vec_sol, &fas->Xg)); PetscCall(VecDuplicate(snes->vec_sol, &fas->Fg)); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode SNESSetFromOptions_FAS(SNES snes, PetscOptionItems *PetscOptionsObject) { SNES_FAS *fas = (SNES_FAS *)snes->data; PetscInt levels = 1; PetscBool flg = PETSC_FALSE, upflg = PETSC_FALSE, downflg = PETSC_FALSE, monflg = PETSC_FALSE, galerkinflg = PETSC_FALSE, continuationflg = PETSC_FALSE; SNESFASType fastype; const char *optionsprefix; SNESLineSearch linesearch; PetscInt m, n_up, n_down; SNES next; PetscBool isFine; PetscFunctionBegin; PetscCall(SNESFASCycleIsFine(snes, &isFine)); PetscOptionsHeadBegin(PetscOptionsObject, "SNESFAS Options-----------------------------------"); /* number of levels -- only process most options on the finest level */ if (isFine) { PetscCall(PetscOptionsInt("-snes_fas_levels", "Number of Levels", "SNESFASSetLevels", levels, &levels, &flg)); if (!flg && snes->dm) { PetscCall(DMGetRefineLevel(snes->dm, &levels)); levels++; fas->usedmfornumberoflevels = PETSC_TRUE; } PetscCall(SNESFASSetLevels(snes, levels, NULL)); fastype = fas->fastype; PetscCall(PetscOptionsEnum("-snes_fas_type", "FAS correction type", "SNESFASSetType", SNESFASTypes, (PetscEnum)fastype, (PetscEnum *)&fastype, &flg)); if (flg) PetscCall(SNESFASSetType(snes, fastype)); PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix)); PetscCall(PetscOptionsInt("-snes_fas_cycles", "Number of cycles", "SNESFASSetCycles", fas->n_cycles, &m, &flg)); if (flg) PetscCall(SNESFASSetCycles(snes, m)); PetscCall(PetscOptionsBool("-snes_fas_continuation", "Corrected grid-sequence continuation", "SNESFASSetContinuation", fas->continuation, &continuationflg, &flg)); if (flg) PetscCall(SNESFASSetContinuation(snes, continuationflg)); PetscCall(PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin", "SNESFASSetGalerkin", fas->galerkin, &galerkinflg, &flg)); if (flg) PetscCall(SNESFASSetGalerkin(snes, galerkinflg)); if (fas->fastype == SNES_FAS_FULL) { PetscCall(PetscOptionsBool("-snes_fas_full_downsweep", "Smooth on the initial down sweep for full FAS cycles", "SNESFASFullSetDownSweep", fas->full_downsweep, &fas->full_downsweep, &flg)); if (flg) PetscCall(SNESFASFullSetDownSweep(snes, fas->full_downsweep)); 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)); if (flg) PetscCall(SNESFASFullSetTotal(snes, fas->full_total)); } PetscCall(PetscOptionsInt("-snes_fas_smoothup", "Number of post-smoothing steps", "SNESFASSetNumberSmoothUp", fas->max_up_it, &n_up, &upflg)); PetscCall(PetscOptionsInt("-snes_fas_smoothdown", "Number of pre-smoothing steps", "SNESFASSetNumberSmoothDown", fas->max_down_it, &n_down, &downflg)); { PetscViewer viewer; PetscViewerFormat format; PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_fas_monitor", &viewer, &format, &monflg)); if (monflg) { PetscViewerAndFormat *vf; PetscCall(PetscViewerAndFormatCreate(viewer, format, &vf)); PetscCall(PetscViewerDestroy(&viewer)); PetscCall(SNESFASSetMonitor(snes, vf, PETSC_TRUE)); } } flg = PETSC_FALSE; monflg = PETSC_TRUE; PetscCall(PetscOptionsBool("-snes_fas_log", "Log times for each FAS level", "SNESFASSetLog", monflg, &monflg, &flg)); if (flg) PetscCall(SNESFASSetLog(snes, monflg)); } PetscOptionsHeadEnd(); /* setup from the determined types if there is no pointwise procedure or smoother defined */ if (upflg) PetscCall(SNESFASSetNumberSmoothUp(snes, n_up)); if (downflg) PetscCall(SNESFASSetNumberSmoothDown(snes, n_down)); /* set up the default line search for coarse grid corrections */ if (fas->fastype == SNES_FAS_ADDITIVE) { if (!snes->linesearch) { PetscCall(SNESGetLineSearch(snes, &linesearch)); PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHL2)); } } /* recursive option setting for the smoothers */ PetscCall(SNESFASCycleGetCorrection(snes, &next)); if (next) PetscCall(SNESSetFromOptions(next)); PetscFunctionReturn(PETSC_SUCCESS); } #include static PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer) { SNES_FAS *fas = (SNES_FAS *)snes->data; PetscBool isFine, iascii, isdraw; PetscInt i; SNES smoothu, smoothd, levelsnes; PetscFunctionBegin; PetscCall(SNESFASCycleIsFine(snes, &isFine)); if (isFine) { PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); if (iascii) { PetscCall(PetscViewerASCIIPrintf(viewer, " type is %s, levels=%" PetscInt_FMT ", cycles=%" PetscInt_FMT "\n", SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles)); if (fas->galerkin) { PetscCall(PetscViewerASCIIPrintf(viewer, " Using Galerkin computed coarse grid function evaluation\n")); } else { PetscCall(PetscViewerASCIIPrintf(viewer, " Not using Galerkin computed coarse grid function evaluation\n")); } for (i = 0; i < fas->levels; i++) { PetscCall(SNESFASGetCycleSNES(snes, i, &levelsnes)); PetscCall(SNESFASCycleGetSmootherUp(levelsnes, &smoothu)); PetscCall(SNESFASCycleGetSmootherDown(levelsnes, &smoothd)); if (!i) { PetscCall(PetscViewerASCIIPrintf(viewer, " Coarse grid solver -- level %" PetscInt_FMT " -------------------------------\n", i)); } else { PetscCall(PetscViewerASCIIPrintf(viewer, " Down solver (pre-smoother) on level %" PetscInt_FMT " -------------------------------\n", i)); } PetscCall(PetscViewerASCIIPushTab(viewer)); if (smoothd) { PetscCall(SNESView(smoothd, viewer)); } else { PetscCall(PetscViewerASCIIPrintf(viewer, "Not yet available\n")); } PetscCall(PetscViewerASCIIPopTab(viewer)); if (i && (smoothd == smoothu)) { PetscCall(PetscViewerASCIIPrintf(viewer, " Up solver (post-smoother) same as down solver (pre-smoother)\n")); } else if (i) { PetscCall(PetscViewerASCIIPrintf(viewer, " Up solver (post-smoother) on level %" PetscInt_FMT " -------------------------------\n", i)); PetscCall(PetscViewerASCIIPushTab(viewer)); if (smoothu) { PetscCall(SNESView(smoothu, viewer)); } else { PetscCall(PetscViewerASCIIPrintf(viewer, "Not yet available\n")); } PetscCall(PetscViewerASCIIPopTab(viewer)); } } } else if (isdraw) { PetscDraw draw; PetscReal x, w, y, bottom, th, wth; SNES_FAS *curfas = fas; PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y)); PetscCall(PetscDrawStringGetSize(draw, &wth, &th)); bottom = y - th; while (curfas) { if (!curfas->smoothu) { PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom)); if (curfas->smoothd) PetscCall(SNESView(curfas->smoothd, viewer)); PetscCall(PetscDrawPopCurrentPoint(draw)); } else { w = 0.5 * PetscMin(1.0 - x, x); PetscCall(PetscDrawPushCurrentPoint(draw, x - w, bottom)); if (curfas->smoothd) PetscCall(SNESView(curfas->smoothd, viewer)); PetscCall(PetscDrawPopCurrentPoint(draw)); PetscCall(PetscDrawPushCurrentPoint(draw, x + w, bottom)); if (curfas->smoothu) PetscCall(SNESView(curfas->smoothu, viewer)); PetscCall(PetscDrawPopCurrentPoint(draw)); } /* this is totally bogus but we have no way of knowing how low the previous one was draw to */ bottom -= 5 * th; if (curfas->next) curfas = (SNES_FAS *)curfas->next->data; else curfas = NULL; } } } PetscFunctionReturn(PETSC_SUCCESS); } /* Defines the action of the downsmoother */ static PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm) { SNESConvergedReason reason; Vec FPC; SNES smoothd; PetscBool flg; SNES_FAS *fas = (SNES_FAS *)snes->data; PetscFunctionBegin; PetscCall(SNESFASCycleGetSmootherDown(snes, &smoothd)); PetscCall(SNESSetInitialFunction(smoothd, F)); if (fas->eventsmoothsolve) PetscCall(PetscLogEventBegin(fas->eventsmoothsolve, smoothd, B, X, 0)); PetscCall(SNESSolve(smoothd, B, X)); if (fas->eventsmoothsolve) PetscCall(PetscLogEventEnd(fas->eventsmoothsolve, smoothd, B, X, 0)); /* check convergence reason for the smoother */ PetscCall(SNESGetConvergedReason(smoothd, &reason)); if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) { snes->reason = SNES_DIVERGED_INNER; PetscFunctionReturn(PETSC_SUCCESS); } PetscCall(SNESGetFunction(smoothd, &FPC, NULL, NULL)); PetscCall(SNESGetAlwaysComputesFinalResidual(smoothd, &flg)); if (!flg) PetscCall(SNESComputeFunction(smoothd, X, FPC)); PetscCall(VecCopy(FPC, F)); if (fnorm) PetscCall(VecNorm(F, NORM_2, fnorm)); PetscFunctionReturn(PETSC_SUCCESS); } /* Defines the action of the upsmoother */ static PetscErrorCode SNESFASUpSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm) { SNESConvergedReason reason; Vec FPC; SNES smoothu; PetscBool flg; SNES_FAS *fas = (SNES_FAS *)snes->data; PetscFunctionBegin; PetscCall(SNESFASCycleGetSmootherUp(snes, &smoothu)); if (fas->eventsmoothsolve) PetscCall(PetscLogEventBegin(fas->eventsmoothsolve, smoothu, 0, 0, 0)); PetscCall(SNESSolve(smoothu, B, X)); if (fas->eventsmoothsolve) PetscCall(PetscLogEventEnd(fas->eventsmoothsolve, smoothu, 0, 0, 0)); /* check convergence reason for the smoother */ PetscCall(SNESGetConvergedReason(smoothu, &reason)); if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) { snes->reason = SNES_DIVERGED_INNER; PetscFunctionReturn(PETSC_SUCCESS); } PetscCall(SNESGetFunction(smoothu, &FPC, NULL, NULL)); PetscCall(SNESGetAlwaysComputesFinalResidual(smoothu, &flg)); if (!flg) PetscCall(SNESComputeFunction(smoothu, X, FPC)); PetscCall(VecCopy(FPC, F)); if (fnorm) PetscCall(VecNorm(F, NORM_2, fnorm)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESFASCreateCoarseVec - create a `Vec` corresponding to a state vector on one level coarser than the current level Collective Input Parameter: . snes - `SNESFAS` object Output Parameter: . Xcoarse - vector on level one coarser than the current level Level: developer .seealso: [](ch_snes), `SNESFASSetRestriction()`, `SNESFASRestrict()`, `SNESFAS` @*/ PetscErrorCode SNESFASCreateCoarseVec(SNES snes, Vec *Xcoarse) { SNES_FAS *fas; PetscFunctionBegin; PetscValidHeaderSpecificType(snes, SNES_CLASSID, 1, SNESFAS); PetscAssertPointer(Xcoarse, 2); fas = (SNES_FAS *)snes->data; if (fas->rscale) { PetscCall(VecDuplicate(fas->rscale, Xcoarse)); } else if (fas->interpolate) { PetscCall(MatCreateVecs(fas->interpolate, Xcoarse, NULL)); } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must set rscale or interpolation"); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESFASRestrict - restrict a `Vec` to the next coarser level Collective Input Parameters: + fine - `SNES` from which to restrict - Xfine - vector to restrict Output Parameter: . Xcoarse - result of restriction Level: developer .seealso: [](ch_snes), `SNES`, `SNESFAS`, `SNESFASSetRestriction()`, `SNESFASSetInjection()`, `SNESFASCreateCoarseVec()` @*/ PetscErrorCode SNESFASRestrict(SNES fine, Vec Xfine, Vec Xcoarse) { SNES_FAS *fas; PetscFunctionBegin; PetscValidHeaderSpecificType(fine, SNES_CLASSID, 1, SNESFAS); PetscValidHeaderSpecific(Xfine, VEC_CLASSID, 2); PetscValidHeaderSpecific(Xcoarse, VEC_CLASSID, 3); fas = (SNES_FAS *)fine->data; if (fas->inject) { PetscCall(MatRestrict(fas->inject, Xfine, Xcoarse)); } else { PetscCall(MatRestrict(fas->restrct, Xfine, Xcoarse)); PetscCall(VecPointwiseMult(Xcoarse, fas->rscale, Xcoarse)); } PetscFunctionReturn(PETSC_SUCCESS); } /* Performs a variant of FAS using the interpolated total coarse solution fine problem: F(x) = b coarse problem: F^c(x^c) = Rb, Initial guess Rx interpolated solution: x^f = I x^c (total solution interpolation */ static PetscErrorCode SNESFASInterpolatedCoarseSolution(SNES snes, Vec X, Vec X_new) { Vec X_c, B_c; SNESConvergedReason reason; SNES next; Mat restrct, interpolate; SNES_FAS *fasc; PetscFunctionBegin; PetscCall(SNESFASCycleGetCorrection(snes, &next)); if (next) { fasc = (SNES_FAS *)next->data; PetscCall(SNESFASCycleGetRestriction(snes, &restrct)); PetscCall(SNESFASCycleGetInterpolation(snes, &interpolate)); X_c = next->vec_sol; if (fasc->eventinterprestrict) PetscCall(PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0)); /* restrict the total solution: Rb */ PetscCall(SNESFASRestrict(snes, X, X_c)); B_c = next->vec_rhs; if (snes->vec_rhs) { /* restrict the total rhs defect: Rb */ PetscCall(MatRestrict(restrct, snes->vec_rhs, B_c)); } else { PetscCall(VecSet(B_c, 0.)); } if (fasc->eventinterprestrict) PetscCall(PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0)); PetscCall(SNESSolve(next, B_c, X_c)); PetscCall(SNESGetConvergedReason(next, &reason)); if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { snes->reason = SNES_DIVERGED_INNER; PetscFunctionReturn(PETSC_SUCCESS); } /* x^f <- Ix^c*/ DM dmc, dmf; PetscCall(SNESGetDM(next, &dmc)); PetscCall(SNESGetDM(snes, &dmf)); if (fasc->eventinterprestrict) PetscCall(PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0)); PetscCall(DMInterpolateSolution(dmc, dmf, interpolate, X_c, X_new)); if (fasc->eventinterprestrict) PetscCall(PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0)); PetscCall(PetscObjectSetName((PetscObject)X_c, "Coarse solution")); PetscCall(VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view")); PetscCall(PetscObjectSetName((PetscObject)X_new, "Updated Fine solution")); PetscCall(VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view")); } PetscFunctionReturn(PETSC_SUCCESS); } /* Performs the FAS coarse correction as: fine problem: F(x) = b coarse problem: F^c(x^c) = b^c b^c = F^c(Rx) - R(F(x) - b) */ static PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new) { Vec X_c, Xo_c, F_c, B_c; SNESConvergedReason reason; SNES next; Mat restrct, interpolate; SNES_FAS *fasc; PetscFunctionBegin; PetscCall(SNESFASCycleGetCorrection(snes, &next)); if (next) { fasc = (SNES_FAS *)next->data; PetscCall(SNESFASCycleGetRestriction(snes, &restrct)); PetscCall(SNESFASCycleGetInterpolation(snes, &interpolate)); X_c = next->vec_sol; Xo_c = next->work[0]; F_c = next->vec_func; B_c = next->vec_rhs; if (fasc->eventinterprestrict) PetscCall(PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0)); PetscCall(SNESFASRestrict(snes, X, Xo_c)); /* restrict the defect: R(F(x) - b) */ PetscCall(MatRestrict(restrct, F, B_c)); if (fasc->eventinterprestrict) PetscCall(PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0)); if (fasc->eventresidual) PetscCall(PetscLogEventBegin(fasc->eventresidual, next, 0, 0, 0)); /* F_c = F^c(Rx) - R(F(x) - b) since the second term was sitting in next->vec_rhs */ PetscCall(SNESComputeFunction(next, Xo_c, F_c)); if (fasc->eventresidual) PetscCall(PetscLogEventEnd(fasc->eventresidual, next, 0, 0, 0)); /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */ PetscCall(VecCopy(B_c, X_c)); PetscCall(VecCopy(F_c, B_c)); PetscCall(VecCopy(X_c, F_c)); /* set initial guess of the coarse problem to the projected fine solution */ PetscCall(VecCopy(Xo_c, X_c)); /* recurse to the next level */ PetscCall(SNESSetInitialFunction(next, F_c)); PetscCall(SNESSolve(next, B_c, X_c)); PetscCall(SNESGetConvergedReason(next, &reason)); if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { snes->reason = SNES_DIVERGED_INNER; PetscFunctionReturn(PETSC_SUCCESS); } /* correct as x <- x + I(x^c - Rx)*/ PetscCall(VecAXPY(X_c, -1.0, Xo_c)); if (fasc->eventinterprestrict) PetscCall(PetscLogEventBegin(fasc->eventinterprestrict, snes, 0, 0, 0)); PetscCall(MatInterpolateAdd(interpolate, X_c, X, X_new)); if (fasc->eventinterprestrict) PetscCall(PetscLogEventEnd(fasc->eventinterprestrict, snes, 0, 0, 0)); PetscCall(PetscObjectSetName((PetscObject)X_c, "Coarse correction")); PetscCall(VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view")); PetscCall(PetscObjectSetName((PetscObject)X_new, "Updated Fine solution")); PetscCall(VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view")); } PetscFunctionReturn(PETSC_SUCCESS); } /* The additive cycle is: xhat = x xhat = dS(x, b) x = coarsecorrection(xhat, b_d) x = x + nu*(xhat - x); (optional) x = uS(x, b) With the coarse RHS (defect correction) as below. */ static PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X) { Vec F, B, Xhat; Vec X_c, Xo_c, F_c, B_c; SNESConvergedReason reason; PetscReal xnorm, fnorm, ynorm; SNESLineSearchReason lsresult; SNES next; Mat restrct, interpolate; SNES_FAS *fas = (SNES_FAS *)snes->data, *fasc; PetscFunctionBegin; PetscCall(SNESFASCycleGetCorrection(snes, &next)); F = snes->vec_func; B = snes->vec_rhs; Xhat = snes->work[1]; PetscCall(VecCopy(X, Xhat)); /* recurse first */ if (next) { fasc = (SNES_FAS *)next->data; PetscCall(SNESFASCycleGetRestriction(snes, &restrct)); PetscCall(SNESFASCycleGetInterpolation(snes, &interpolate)); if (fas->eventresidual) PetscCall(PetscLogEventBegin(fas->eventresidual, snes, 0, 0, 0)); PetscCall(SNESComputeFunction(snes, Xhat, F)); if (fas->eventresidual) PetscCall(PetscLogEventEnd(fas->eventresidual, snes, 0, 0, 0)); PetscCall(VecNorm(F, NORM_2, &fnorm)); X_c = next->vec_sol; Xo_c = next->work[0]; F_c = next->vec_func; B_c = next->vec_rhs; PetscCall(SNESFASRestrict(snes, Xhat, Xo_c)); /* restrict the defect */ PetscCall(MatRestrict(restrct, F, B_c)); /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */ if (fasc->eventresidual) PetscCall(PetscLogEventBegin(fasc->eventresidual, next, 0, 0, 0)); PetscCall(SNESComputeFunction(next, Xo_c, F_c)); if (fasc->eventresidual) PetscCall(PetscLogEventEnd(fasc->eventresidual, next, 0, 0, 0)); PetscCall(VecCopy(B_c, X_c)); PetscCall(VecCopy(F_c, B_c)); PetscCall(VecCopy(X_c, F_c)); /* set initial guess of the coarse problem to the projected fine solution */ PetscCall(VecCopy(Xo_c, X_c)); /* recurse */ PetscCall(SNESSetInitialFunction(next, F_c)); PetscCall(SNESSolve(next, B_c, X_c)); /* smooth on this level */ PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &fnorm)); PetscCall(SNESGetConvergedReason(next, &reason)); if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { snes->reason = SNES_DIVERGED_INNER; PetscFunctionReturn(PETSC_SUCCESS); } /* correct as x <- x + I(x^c - Rx)*/ PetscCall(VecAYPX(X_c, -1.0, Xo_c)); PetscCall(MatInterpolate(interpolate, X_c, Xhat)); /* additive correction of the coarse direction*/ PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat)); PetscCall(SNESLineSearchGetReason(snes->linesearch, &lsresult)); PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm)); if (lsresult) { if (++snes->numFailures >= snes->maxFailures) { snes->reason = SNES_DIVERGED_LINE_SEARCH; PetscFunctionReturn(PETSC_SUCCESS); } } } else { PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm)); } PetscFunctionReturn(PETSC_SUCCESS); } /* Defines the FAS cycle as: fine problem: F(x) = b coarse problem: F^c(x) = b^c b^c = F^c(Rx) - R(F(x) - b) correction: x = x + I(x^c - Rx) */ static PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X) { Vec F, B; SNES next; PetscFunctionBegin; F = snes->vec_func; B = snes->vec_rhs; /* pre-smooth -- just update using the pre-smoother */ PetscCall(SNESFASCycleGetCorrection(snes, &next)); PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm)); if (next) { PetscCall(SNESFASCoarseCorrection(snes, X, F, X)); PetscCall(SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm)); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode SNESFASCycleSetupPhase_Full(SNES snes) { SNES next; SNES_FAS *fas = (SNES_FAS *)snes->data; PetscBool isFine; PetscFunctionBegin; /* pre-smooth -- just update using the pre-smoother */ PetscCall(SNESFASCycleIsFine(snes, &isFine)); PetscCall(SNESFASCycleGetCorrection(snes, &next)); fas->full_stage = 0; if (next) PetscCall(SNESFASCycleSetupPhase_Full(next)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode SNESFASCycle_Full(SNES snes, Vec X) { Vec F, B; SNES_FAS *fas = (SNES_FAS *)snes->data; PetscBool isFine; SNES next; PetscFunctionBegin; F = snes->vec_func; B = snes->vec_rhs; PetscCall(SNESFASCycleIsFine(snes, &isFine)); PetscCall(SNESFASCycleGetCorrection(snes, &next)); if (isFine) PetscCall(SNESFASCycleSetupPhase_Full(snes)); if (fas->full_stage == 0) { /* downsweep */ if (next) { if (fas->level != 1) next->max_its += 1; if (fas->full_downsweep) PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm)); fas->full_downsweep = PETSC_TRUE; if (fas->full_total) PetscCall(SNESFASInterpolatedCoarseSolution(snes, X, X)); else PetscCall(SNESFASCoarseCorrection(snes, X, F, X)); fas->full_total = PETSC_FALSE; PetscCall(SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm)); if (fas->level != 1) next->max_its -= 1; } else { /* The smoother on the coarse level is the coarse solver */ PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm)); } fas->full_stage = 1; } else if (fas->full_stage == 1) { if (snes->iter == 0) PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm)); if (next) { PetscCall(SNESFASCoarseCorrection(snes, X, F, X)); PetscCall(SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm)); } } /* final v-cycle */ if (isFine) { if (next) { PetscCall(SNESFASCoarseCorrection(snes, X, F, X)); PetscCall(SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm)); } } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode SNESFASCycle_Kaskade(SNES snes, Vec X) { Vec F, B; SNES next; PetscFunctionBegin; F = snes->vec_func; B = snes->vec_rhs; PetscCall(SNESFASCycleGetCorrection(snes, &next)); if (next) { PetscCall(SNESFASCoarseCorrection(snes, X, F, X)); PetscCall(SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm)); } else { PetscCall(SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm)); } PetscFunctionReturn(PETSC_SUCCESS); } PetscBool SNEScite = PETSC_FALSE; const char SNESCitation[] = "@Article{bruneknepleysmithtu15," " title = {Composing Scalable Nonlinear Algebraic Solvers}," " author = {Peter R. Brune and Matthew G. Knepley and Barry F. Smith and Xuemin Tu}," " journal = {SIAM Review}," " volume = {57}," " number = {4}," " pages = {535--565}," " doi = {10.1137/130936725}," " year = {2015}\n"; static PetscErrorCode SNESSolve_FAS(SNES snes) { PetscInt i; Vec X, F; PetscReal fnorm; SNES_FAS *fas = (SNES_FAS *)snes->data, *ffas; DM dm; PetscBool isFine; PetscFunctionBegin; 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); PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite)); snes->reason = SNES_CONVERGED_ITERATING; X = snes->vec_sol; F = snes->vec_func; PetscCall(SNESFASCycleIsFine(snes, &isFine)); /* norm setup */ PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); snes->iter = 0; snes->norm = 0; PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); if (!snes->vec_func_init_set) { if (fas->eventresidual) PetscCall(PetscLogEventBegin(fas->eventresidual, snes, 0, 0, 0)); PetscCall(SNESComputeFunction(snes, X, F)); if (fas->eventresidual) PetscCall(PetscLogEventEnd(fas->eventresidual, snes, 0, 0, 0)); } else snes->vec_func_init_set = PETSC_FALSE; PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F|| */ SNESCheckFunctionNorm(snes, fnorm); PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); snes->norm = fnorm; PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0)); /* test convergence */ PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm)); PetscCall(SNESMonitor(snes, snes->iter, fnorm)); if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS); if (isFine) { /* propagate scale-dependent data up the hierarchy */ PetscCall(SNESGetDM(snes, &dm)); for (ffas = fas; ffas->next; ffas = (SNES_FAS *)ffas->next->data) { DM dmcoarse; PetscCall(SNESGetDM(ffas->next, &dmcoarse)); PetscCall(DMRestrict(dm, ffas->restrct, ffas->rscale, ffas->inject, dmcoarse)); dm = dmcoarse; } } for (i = 0; i < snes->max_its; i++) { /* Call general purpose update function */ PetscTryTypeMethod(snes, update, snes->iter); if (fas->fastype == SNES_FAS_MULTIPLICATIVE) { PetscCall(SNESFASCycle_Multiplicative(snes, X)); } else if (fas->fastype == SNES_FAS_ADDITIVE) { PetscCall(SNESFASCycle_Additive(snes, X)); } else if (fas->fastype == SNES_FAS_FULL) { PetscCall(SNESFASCycle_Full(snes, X)); } else if (fas->fastype == SNES_FAS_KASKADE) { PetscCall(SNESFASCycle_Kaskade(snes, X)); } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Unsupported FAS type"); /* check for FAS cycle divergence */ if (snes->reason != SNES_CONVERGED_ITERATING) PetscFunctionReturn(PETSC_SUCCESS); /* Monitor convergence */ PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); snes->iter = i + 1; PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0)); PetscCall(SNESConverged(snes, snes->iter, 0.0, 0.0, snes->norm)); PetscCall(SNESMonitor(snes, snes->iter, snes->norm)); if (snes->reason) break; } PetscFunctionReturn(PETSC_SUCCESS); } /*MC SNESFAS - An implementation of the Full Approximation Scheme nonlinear multigrid solver, FAS, or nonlinear multigrid {cite}`bruneknepleysmithtu15` for solving nonlinear systems of equations with `SNES`. The nonlinear problem is solved by correction using coarse versions of the nonlinear problem. This problem is perturbed so that a projected solution of the fine problem elicits no correction from the coarse problem. Options Database Keys and Prefixes: + -snes_fas_levels - The number of levels . -snes_fas_cycles - The number of cycles -- 1 for V, 2 for W . -snes_fas_type - Additive or multiplicative cycle . -snes_fas_galerkin - Form coarse problems by projection back upon the fine problem . -snes_fas_smoothup - The number of iterations of the post-smoother . -snes_fas_smoothdown - The number of iterations of the pre-smoother . -snes_fas_monitor - Monitor progress of all of the levels . -snes_fas_full_downsweep - call the downsmooth on the initial downsweep of full FAS . -fas_levels_snes_ - prefix for `SNES` options for all smoothers . -fas_levels_cycle_snes_ - prefix for `SNES` options for all cycles . -fas_levels_i_snes_ - prefix `SNES` options for the smoothers on level i . -fas_levels_i_cycle_snes_ - prefix for `SNES` options for the cycle on level i - -fas_coarse_snes_ - prefix for `SNES` options for the coarsest smoother Level: beginner Note: The organization of the `SNESFAS` solver is slightly different from the organization of `PCMG` As each level has smoother `SNES` instances(down and potentially up) and a cycle `SNES` instance. The cycle `SNES` instance may be used for monitoring convergence on a particular level. .seealso: [](ch_snes), `PCMG`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESType`, `SNESFASSetRestriction()`, `SNESFASSetInjection()`, `SNESFASFullGetTotal()`, `SNESFASSetType()`, `SNESFASGetType()`, `SNESFASSetLevels()`, `SNESFASGetLevels()`, `SNESFASGetCycleSNES()`, `SNESFASSetNumberSmoothUp()`, `SNESFASSetNumberSmoothDown()`, `SNESFASSetContinuation()`, `SNESFASSetCycles()`, `SNESFASSetMonitor()`, `SNESFASSetLog()`, `SNESFASCycleSetCycles()`, `SNESFASCycleGetSmoother()`, `SNESFASCycleGetSmootherUp()`, `SNESFASCycleGetSmootherDown()`, `SNESFASCycleGetCorrection()`, `SNESFASCycleGetInterpolation()`, `SNESFASCycleGetRestriction()`, `SNESFASCycleGetInjection()`, `SNESFASCycleGetRScale()`, `SNESFASCycleIsFine()`, `SNESFASSetInterpolation()`, `SNESFASGetInterpolation()`, `SNESFASSetRestriction()`, `SNESFASGetRestriction()`, `SNESFASSetInjection()`, `SNESFASGetInjection()`, `SNESFASSetRScale()`,`SNESFASGetSmoother()`, `SNESFASGetSmootherDown()`, `SNESFASGetSmootherUp()`, `SNESFASGetCoarseSolve()`, `SNESFASFullSetDownSweep()`, `SNESFASFullSetTotal()`, `SNESFASFullGetTotal()` M*/ PETSC_EXTERN PetscErrorCode SNESCreate_FAS(SNES snes) { SNES_FAS *fas; PetscFunctionBegin; snes->ops->destroy = SNESDestroy_FAS; snes->ops->setup = SNESSetUp_FAS; snes->ops->setfromoptions = SNESSetFromOptions_FAS; snes->ops->view = SNESView_FAS; snes->ops->solve = SNESSolve_FAS; snes->ops->reset = SNESReset_FAS; snes->usesksp = PETSC_FALSE; snes->usesnpc = PETSC_FALSE; PetscCall(SNESParametersInitialize(snes)); PetscObjectParameterSetDefault(snes, max_funcs, 30000); PetscObjectParameterSetDefault(snes, max_its, 10000); snes->alwayscomputesfinalresidual = PETSC_TRUE; PetscCall(PetscNew(&fas)); snes->data = (void *)fas; fas->level = 0; fas->levels = 1; fas->n_cycles = 1; fas->max_up_it = 1; fas->max_down_it = 1; fas->smoothu = NULL; fas->smoothd = NULL; fas->next = NULL; fas->previous = NULL; fas->fine = snes; fas->interpolate = NULL; fas->restrct = NULL; fas->inject = NULL; fas->usedmfornumberoflevels = PETSC_FALSE; fas->fastype = SNES_FAS_MULTIPLICATIVE; fas->full_downsweep = PETSC_FALSE; fas->full_total = PETSC_FALSE; fas->eventsmoothsetup = 0; fas->eventsmoothsolve = 0; fas->eventresidual = 0; fas->eventinterprestrict = 0; PetscFunctionReturn(PETSC_SUCCESS); }