1 /* Defines the basic SNES object */ 2 #include <../src/snes/impls/fas/fasimpls.h> /*I "petscsnes.h" I*/ 3 4 const char *const SNESFASTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","SNESFASType","SNES_FAS",0}; 5 6 static PetscErrorCode SNESReset_FAS(SNES snes) 7 { 8 PetscErrorCode ierr = 0; 9 SNES_FAS * fas = (SNES_FAS*)snes->data; 10 11 PetscFunctionBegin; 12 ierr = SNESDestroy(&fas->smoothu);CHKERRQ(ierr); 13 ierr = SNESDestroy(&fas->smoothd);CHKERRQ(ierr); 14 ierr = MatDestroy(&fas->inject);CHKERRQ(ierr); 15 ierr = MatDestroy(&fas->interpolate);CHKERRQ(ierr); 16 ierr = MatDestroy(&fas->restrct);CHKERRQ(ierr); 17 ierr = VecDestroy(&fas->rscale);CHKERRQ(ierr); 18 if (fas->galerkin) { 19 ierr = VecDestroy(&fas->Xg);CHKERRQ(ierr); 20 ierr = VecDestroy(&fas->Fg);CHKERRQ(ierr); 21 } 22 if (fas->next) {ierr = SNESReset(fas->next);CHKERRQ(ierr);} 23 PetscFunctionReturn(0); 24 } 25 26 static PetscErrorCode SNESDestroy_FAS(SNES snes) 27 { 28 SNES_FAS * fas = (SNES_FAS*)snes->data; 29 PetscErrorCode ierr = 0; 30 31 PetscFunctionBegin; 32 /* recursively resets and then destroys */ 33 ierr = SNESReset(snes);CHKERRQ(ierr); 34 if (fas->next) { 35 ierr = SNESDestroy(&fas->next);CHKERRQ(ierr); 36 } 37 ierr = PetscFree(fas);CHKERRQ(ierr); 38 PetscFunctionReturn(0); 39 } 40 41 static PetscErrorCode SNESSetUp_FAS(SNES snes) 42 { 43 SNES_FAS *fas = (SNES_FAS*) snes->data; 44 PetscErrorCode ierr; 45 PetscInt dm_levels; 46 Vec vec_sol, vec_func, vec_sol_update, vec_rhs; /* preserve these if they're set through the reset */ 47 SNES next; 48 PetscBool isFine, hasCreateRestriction; 49 SNESLineSearch linesearch; 50 SNESLineSearch slinesearch; 51 void *lsprectx,*lspostctx; 52 PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*); 53 PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*); 54 55 PetscFunctionBegin; 56 ierr = SNESFASCycleIsFine(snes, &isFine);CHKERRQ(ierr); 57 if (fas->usedmfornumberoflevels && isFine) { 58 ierr = DMGetRefineLevel(snes->dm,&dm_levels);CHKERRQ(ierr); 59 dm_levels++; 60 if (dm_levels > fas->levels) { 61 /* we don't want the solution and func vectors to be destroyed in the SNESReset when it's called in SNESFASSetLevels_FAS*/ 62 vec_sol = snes->vec_sol; 63 vec_func = snes->vec_func; 64 vec_sol_update = snes->vec_sol_update; 65 vec_rhs = snes->vec_rhs; 66 snes->vec_sol = NULL; 67 snes->vec_func = NULL; 68 snes->vec_sol_update = NULL; 69 snes->vec_rhs = NULL; 70 71 /* reset the number of levels */ 72 ierr = SNESFASSetLevels(snes,dm_levels,NULL);CHKERRQ(ierr); 73 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 74 75 snes->vec_sol = vec_sol; 76 snes->vec_func = vec_func; 77 snes->vec_rhs = vec_rhs; 78 snes->vec_sol_update = vec_sol_update; 79 } 80 } 81 ierr = SNESFASCycleGetCorrection(snes, &next);CHKERRQ(ierr); 82 if (!isFine) snes->gridsequence = 0; /* no grid sequencing inside the multigrid hierarchy! */ 83 84 ierr = SNESSetWorkVecs(snes, 2);CHKERRQ(ierr); /* work vectors used for intergrid transfers */ 85 86 /* set up the smoothers if they haven't already been set up */ 87 if (!fas->smoothd) { 88 ierr = SNESFASCycleCreateSmoother_Private(snes, &fas->smoothd);CHKERRQ(ierr); 89 } 90 91 if (snes->dm) { 92 /* set the smoother DMs properly */ 93 if (fas->smoothu) ierr = SNESSetDM(fas->smoothu, snes->dm);CHKERRQ(ierr); 94 ierr = SNESSetDM(fas->smoothd, snes->dm);CHKERRQ(ierr); 95 /* construct EVERYTHING from the DM -- including the progressive set of smoothers */ 96 if (next) { 97 /* for now -- assume the DM and the evaluation functions have been set externally */ 98 if (!next->dm) { 99 ierr = DMCoarsen(snes->dm, PetscObjectComm((PetscObject)next), &next->dm);CHKERRQ(ierr); 100 ierr = SNESSetDM(next, next->dm);CHKERRQ(ierr); 101 } 102 /* set the interpolation and restriction from the DM */ 103 if (!fas->interpolate) { 104 ierr = DMCreateInterpolation(next->dm, snes->dm, &fas->interpolate, &fas->rscale);CHKERRQ(ierr); 105 if (!fas->restrct) { 106 ierr = DMHasCreateRestriction(next->dm, &hasCreateRestriction);CHKERRQ(ierr); 107 /* DM can create restrictions, use that */ 108 if (hasCreateRestriction) { 109 ierr = DMCreateRestriction(next->dm, snes->dm, &fas->restrct);CHKERRQ(ierr); 110 } else { 111 ierr = PetscObjectReference((PetscObject)fas->interpolate);CHKERRQ(ierr); 112 fas->restrct = fas->interpolate; 113 } 114 } 115 } 116 /* set the injection from the DM */ 117 if (!fas->inject) { 118 PetscBool flg; 119 ierr = DMHasCreateInjection(next->dm, &flg);CHKERRQ(ierr); 120 if (flg) { 121 ierr = DMCreateInjection(next->dm, snes->dm, &fas->inject);CHKERRQ(ierr); 122 } 123 } 124 } 125 } 126 /*pass the smoother, function, and jacobian up to the next level if it's not user set already */ 127 if (fas->galerkin) { 128 if (next) { 129 ierr = SNESSetFunction(next, NULL, SNESFASGalerkinFunctionDefault, next);CHKERRQ(ierr); 130 } 131 if (fas->smoothd && fas->level != fas->levels - 1) { 132 ierr = SNESSetFunction(fas->smoothd, NULL, SNESFASGalerkinFunctionDefault, snes);CHKERRQ(ierr); 133 } 134 if (fas->smoothu && fas->level != fas->levels - 1) { 135 ierr = SNESSetFunction(fas->smoothu, NULL, SNESFASGalerkinFunctionDefault, snes);CHKERRQ(ierr); 136 } 137 } 138 139 /* sets the down (pre) smoother's default norm and sets it from options */ 140 if (fas->smoothd) { 141 if (fas->level == 0 && fas->levels != 1) { 142 ierr = SNESSetNormSchedule(fas->smoothd, SNES_NORM_NONE);CHKERRQ(ierr); 143 } else { 144 ierr = SNESSetNormSchedule(fas->smoothd, SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 145 } 146 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothd);CHKERRQ(ierr); 147 ierr = SNESSetFromOptions(fas->smoothd);CHKERRQ(ierr); 148 if (snes->linesearch) { 149 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 150 ierr = SNESGetLineSearch(fas->smoothd,&slinesearch);CHKERRQ(ierr); 151 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 152 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 153 ierr = SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);CHKERRQ(ierr); 154 ierr = SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);CHKERRQ(ierr); 155 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);CHKERRQ(ierr); 156 } 157 158 fas->smoothd->vec_sol = snes->vec_sol; 159 ierr = PetscObjectReference((PetscObject)snes->vec_sol);CHKERRQ(ierr); 160 fas->smoothd->vec_sol_update = snes->vec_sol_update; 161 ierr = PetscObjectReference((PetscObject)snes->vec_sol_update);CHKERRQ(ierr); 162 fas->smoothd->vec_func = snes->vec_func; 163 ierr = PetscObjectReference((PetscObject)snes->vec_func);CHKERRQ(ierr); 164 165 if (fas->eventsmoothsetup) {ierr = PetscLogEventBegin(fas->eventsmoothsetup,fas->smoothd,0,0,0);CHKERRQ(ierr);} 166 ierr = SNESSetUp(fas->smoothd);CHKERRQ(ierr); 167 if (fas->eventsmoothsetup) {ierr = PetscLogEventEnd(fas->eventsmoothsetup,fas->smoothd,0,0,0);CHKERRQ(ierr);} 168 } 169 170 /* sets the up (post) smoother's default norm and sets it from options */ 171 if (fas->smoothu) { 172 if (fas->level != fas->levels - 1) { 173 ierr = SNESSetNormSchedule(fas->smoothu, SNES_NORM_NONE);CHKERRQ(ierr); 174 } else { 175 ierr = SNESSetNormSchedule(fas->smoothu, SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 176 } 177 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothu);CHKERRQ(ierr); 178 ierr = SNESSetFromOptions(fas->smoothu);CHKERRQ(ierr); 179 if (snes->linesearch) { 180 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 181 ierr = SNESGetLineSearch(fas->smoothu,&slinesearch);CHKERRQ(ierr); 182 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 183 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 184 ierr = SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);CHKERRQ(ierr); 185 ierr = SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);CHKERRQ(ierr); 186 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);CHKERRQ(ierr); 187 } 188 189 fas->smoothu->vec_sol = snes->vec_sol; 190 ierr = PetscObjectReference((PetscObject)snes->vec_sol);CHKERRQ(ierr); 191 fas->smoothu->vec_sol_update = snes->vec_sol_update; 192 ierr = PetscObjectReference((PetscObject)snes->vec_sol_update);CHKERRQ(ierr); 193 fas->smoothu->vec_func = snes->vec_func; 194 ierr = PetscObjectReference((PetscObject)snes->vec_func);CHKERRQ(ierr); 195 196 if (fas->eventsmoothsetup) {ierr = PetscLogEventBegin(fas->eventsmoothsetup,fas->smoothu,0,0,0);CHKERRQ(ierr);} 197 ierr = SNESSetUp(fas->smoothu);CHKERRQ(ierr); 198 if (fas->eventsmoothsetup) {ierr = PetscLogEventEnd(fas->eventsmoothsetup,fas->smoothu,0,0,0);CHKERRQ(ierr);} 199 200 } 201 202 if (next) { 203 /* gotta set up the solution vector for this to work */ 204 if (!next->vec_sol) {ierr = SNESFASCreateCoarseVec(snes,&next->vec_sol);CHKERRQ(ierr);} 205 if (!next->vec_rhs) {ierr = SNESFASCreateCoarseVec(snes,&next->vec_rhs);CHKERRQ(ierr);} 206 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next);CHKERRQ(ierr); 207 if (snes->linesearch) { 208 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 209 ierr = SNESGetLineSearch(fas->next,&slinesearch);CHKERRQ(ierr); 210 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 211 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 212 ierr = SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);CHKERRQ(ierr); 213 ierr = SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);CHKERRQ(ierr); 214 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);CHKERRQ(ierr); 215 } 216 ierr = SNESSetUp(next);CHKERRQ(ierr); 217 } 218 /* setup FAS work vectors */ 219 if (fas->galerkin) { 220 ierr = VecDuplicate(snes->vec_sol, &fas->Xg);CHKERRQ(ierr); 221 ierr = VecDuplicate(snes->vec_sol, &fas->Fg);CHKERRQ(ierr); 222 } 223 PetscFunctionReturn(0); 224 } 225 226 static PetscErrorCode SNESSetFromOptions_FAS(PetscOptionItems *PetscOptionsObject,SNES snes) 227 { 228 SNES_FAS *fas = (SNES_FAS*) snes->data; 229 PetscInt levels = 1; 230 PetscBool flg = PETSC_FALSE, upflg = PETSC_FALSE, downflg = PETSC_FALSE, monflg = PETSC_FALSE, galerkinflg = PETSC_FALSE,continuationflg = PETSC_FALSE; 231 PetscErrorCode ierr; 232 SNESFASType fastype; 233 const char *optionsprefix; 234 SNESLineSearch linesearch; 235 PetscInt m, n_up, n_down; 236 SNES next; 237 PetscBool isFine; 238 239 PetscFunctionBegin; 240 ierr = SNESFASCycleIsFine(snes, &isFine);CHKERRQ(ierr); 241 ierr = PetscOptionsHead(PetscOptionsObject,"SNESFAS Options-----------------------------------");CHKERRQ(ierr); 242 243 /* number of levels -- only process most options on the finest level */ 244 if (isFine) { 245 ierr = PetscOptionsInt("-snes_fas_levels", "Number of Levels", "SNESFASSetLevels", levels, &levels, &flg);CHKERRQ(ierr); 246 if (!flg && snes->dm) { 247 ierr = DMGetRefineLevel(snes->dm,&levels);CHKERRQ(ierr); 248 levels++; 249 fas->usedmfornumberoflevels = PETSC_TRUE; 250 } 251 ierr = SNESFASSetLevels(snes, levels, NULL);CHKERRQ(ierr); 252 fastype = fas->fastype; 253 ierr = PetscOptionsEnum("-snes_fas_type","FAS correction type","SNESFASSetType",SNESFASTypes,(PetscEnum)fastype,(PetscEnum*)&fastype,&flg);CHKERRQ(ierr); 254 if (flg) { 255 ierr = SNESFASSetType(snes, fastype);CHKERRQ(ierr); 256 } 257 258 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 259 ierr = PetscOptionsInt("-snes_fas_cycles","Number of cycles","SNESFASSetCycles",fas->n_cycles,&m,&flg);CHKERRQ(ierr); 260 if (flg) { 261 ierr = SNESFASSetCycles(snes, m);CHKERRQ(ierr); 262 } 263 ierr = PetscOptionsBool("-snes_fas_continuation","Corrected grid-sequence continuation","SNESFASSetContinuation",fas->continuation,&continuationflg,&flg);CHKERRQ(ierr); 264 if (flg) { 265 ierr = SNESFASSetContinuation(snes,continuationflg);CHKERRQ(ierr); 266 } 267 268 ierr = PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin","SNESFASSetGalerkin",fas->galerkin,&galerkinflg,&flg);CHKERRQ(ierr); 269 if (flg) { 270 ierr = SNESFASSetGalerkin(snes, galerkinflg);CHKERRQ(ierr); 271 } 272 273 if (fas->fastype == SNES_FAS_FULL) { 274 ierr = PetscOptionsBool("-snes_fas_full_downsweep","Smooth on the initial down sweep for full FAS cycles","SNESFASFullSetDownSweep",fas->full_downsweep,&fas->full_downsweep,&flg);CHKERRQ(ierr); 275 if (flg) {ierr = SNESFASFullSetDownSweep(snes,fas->full_downsweep);CHKERRQ(ierr);} 276 } 277 278 ierr = PetscOptionsInt("-snes_fas_smoothup","Number of post-smoothing steps","SNESFASSetNumberSmoothUp",fas->max_up_it,&n_up,&upflg);CHKERRQ(ierr); 279 280 ierr = PetscOptionsInt("-snes_fas_smoothdown","Number of pre-smoothing steps","SNESFASSetNumberSmoothDown",fas->max_down_it,&n_down,&downflg);CHKERRQ(ierr); 281 282 { 283 PetscViewer viewer; 284 PetscViewerFormat format; 285 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,"-snes_fas_monitor",&viewer,&format,&monflg);CHKERRQ(ierr); 286 if (monflg) { 287 PetscViewerAndFormat *vf; 288 ierr = PetscViewerAndFormatCreate(viewer,format,&vf);CHKERRQ(ierr); 289 ierr = PetscObjectDereference((PetscObject)viewer);CHKERRQ(ierr); 290 ierr = SNESFASSetMonitor(snes,vf,PETSC_TRUE);CHKERRQ(ierr); 291 } 292 } 293 flg = PETSC_FALSE; 294 monflg = PETSC_TRUE; 295 ierr = PetscOptionsBool("-snes_fas_log","Log times for each FAS level","SNESFASSetLog",monflg,&monflg,&flg);CHKERRQ(ierr); 296 if (flg) {ierr = SNESFASSetLog(snes,monflg);CHKERRQ(ierr);} 297 } 298 299 ierr = PetscOptionsTail();CHKERRQ(ierr); 300 /* setup from the determined types if there is no pointwise procedure or smoother defined */ 301 if (upflg) { 302 ierr = SNESFASSetNumberSmoothUp(snes,n_up);CHKERRQ(ierr); 303 } 304 if (downflg) { 305 ierr = SNESFASSetNumberSmoothDown(snes,n_down);CHKERRQ(ierr); 306 } 307 308 /* set up the default line search for coarse grid corrections */ 309 if (fas->fastype == SNES_FAS_ADDITIVE) { 310 if (!snes->linesearch) { 311 ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr); 312 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);CHKERRQ(ierr); 313 } 314 } 315 316 ierr = SNESFASCycleGetCorrection(snes, &next);CHKERRQ(ierr); 317 /* recursive option setting for the smoothers */ 318 if (next) {ierr = SNESSetFromOptions(next);CHKERRQ(ierr);} 319 PetscFunctionReturn(0); 320 } 321 322 #include <petscdraw.h> 323 static PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer) 324 { 325 SNES_FAS *fas = (SNES_FAS*) snes->data; 326 PetscBool isFine,iascii,isdraw; 327 PetscInt i; 328 PetscErrorCode ierr; 329 SNES smoothu, smoothd, levelsnes; 330 331 PetscFunctionBegin; 332 ierr = SNESFASCycleIsFine(snes, &isFine);CHKERRQ(ierr); 333 if (isFine) { 334 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 335 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 336 if (iascii) { 337 ierr = PetscViewerASCIIPrintf(viewer, " type is %s, levels=%D, cycles=%D\n", SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles);CHKERRQ(ierr); 338 if (fas->galerkin) { 339 ierr = PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid function evaluation\n");CHKERRQ(ierr); 340 } else { 341 ierr = PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid function evaluation\n");CHKERRQ(ierr); 342 } 343 for (i=0; i<fas->levels; i++) { 344 ierr = SNESFASGetCycleSNES(snes, i, &levelsnes);CHKERRQ(ierr); 345 ierr = SNESFASCycleGetSmootherUp(levelsnes, &smoothu);CHKERRQ(ierr); 346 ierr = SNESFASCycleGetSmootherDown(levelsnes, &smoothd);CHKERRQ(ierr); 347 if (!i) { 348 ierr = PetscViewerASCIIPrintf(viewer," Coarse grid solver -- level %D -------------------------------\n",i);CHKERRQ(ierr); 349 } else { 350 ierr = PetscViewerASCIIPrintf(viewer," Down solver (pre-smoother) on level %D -------------------------------\n",i);CHKERRQ(ierr); 351 } 352 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 353 if (smoothd) { 354 ierr = SNESView(smoothd,viewer);CHKERRQ(ierr); 355 } else { 356 ierr = PetscViewerASCIIPrintf(viewer,"Not yet available\n");CHKERRQ(ierr); 357 } 358 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 359 if (i && (smoothd == smoothu)) { 360 ierr = PetscViewerASCIIPrintf(viewer," Up solver (post-smoother) same as down solver (pre-smoother)\n");CHKERRQ(ierr); 361 } else if (i) { 362 ierr = PetscViewerASCIIPrintf(viewer," Up solver (post-smoother) on level %D -------------------------------\n",i);CHKERRQ(ierr); 363 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 364 if (smoothu) { 365 ierr = SNESView(smoothu,viewer);CHKERRQ(ierr); 366 } else { 367 ierr = PetscViewerASCIIPrintf(viewer,"Not yet available\n");CHKERRQ(ierr); 368 } 369 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 370 } 371 } 372 } else if (isdraw) { 373 PetscDraw draw; 374 PetscReal x,w,y,bottom,th,wth; 375 SNES_FAS *curfas = fas; 376 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 377 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 378 ierr = PetscDrawStringGetSize(draw,&wth,&th);CHKERRQ(ierr); 379 bottom = y - th; 380 while (curfas) { 381 if (!curfas->smoothu) { 382 ierr = PetscDrawPushCurrentPoint(draw,x,bottom);CHKERRQ(ierr); 383 if (curfas->smoothd) ierr = SNESView(curfas->smoothd,viewer);CHKERRQ(ierr); 384 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 385 } else { 386 w = 0.5*PetscMin(1.0-x,x); 387 ierr = PetscDrawPushCurrentPoint(draw,x-w,bottom);CHKERRQ(ierr); 388 if (curfas->smoothd) ierr = SNESView(curfas->smoothd,viewer);CHKERRQ(ierr); 389 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 390 ierr = PetscDrawPushCurrentPoint(draw,x+w,bottom);CHKERRQ(ierr); 391 if (curfas->smoothu) ierr = SNESView(curfas->smoothu,viewer);CHKERRQ(ierr); 392 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 393 } 394 /* this is totally bogus but we have no way of knowing how low the previous one was draw to */ 395 bottom -= 5*th; 396 if (curfas->next) curfas = (SNES_FAS*)curfas->next->data; 397 else curfas = NULL; 398 } 399 } 400 } 401 PetscFunctionReturn(0); 402 } 403 404 /* 405 Defines the action of the downsmoother 406 */ 407 static PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm) 408 { 409 PetscErrorCode ierr = 0; 410 SNESConvergedReason reason; 411 Vec FPC; 412 SNES smoothd; 413 PetscBool flg; 414 SNES_FAS *fas = (SNES_FAS*) snes->data; 415 416 PetscFunctionBegin; 417 ierr = SNESFASCycleGetSmootherDown(snes, &smoothd);CHKERRQ(ierr); 418 ierr = SNESSetInitialFunction(smoothd, F);CHKERRQ(ierr); 419 if (fas->eventsmoothsolve) {ierr = PetscLogEventBegin(fas->eventsmoothsolve,smoothd,B,X,0);CHKERRQ(ierr);} 420 ierr = SNESSolve(smoothd, B, X);CHKERRQ(ierr); 421 if (fas->eventsmoothsolve) {ierr = PetscLogEventEnd(fas->eventsmoothsolve,smoothd,B,X,0);CHKERRQ(ierr);} 422 /* check convergence reason for the smoother */ 423 ierr = SNESGetConvergedReason(smoothd,&reason);CHKERRQ(ierr); 424 if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) { 425 snes->reason = SNES_DIVERGED_INNER; 426 PetscFunctionReturn(0); 427 } 428 429 ierr = SNESGetFunction(smoothd, &FPC, NULL, NULL);CHKERRQ(ierr); 430 ierr = SNESGetAlwaysComputesFinalResidual(smoothd, &flg);CHKERRQ(ierr); 431 if (!flg) { 432 ierr = SNESComputeFunction(smoothd, X, FPC);CHKERRQ(ierr); 433 } 434 ierr = VecCopy(FPC, F);CHKERRQ(ierr); 435 if (fnorm) {ierr = VecNorm(F,NORM_2,fnorm);CHKERRQ(ierr);} 436 PetscFunctionReturn(0); 437 } 438 439 440 /* 441 Defines the action of the upsmoother 442 */ 443 static PetscErrorCode SNESFASUpSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm) 444 { 445 PetscErrorCode ierr = 0; 446 SNESConvergedReason reason; 447 Vec FPC; 448 SNES smoothu; 449 PetscBool flg; 450 SNES_FAS *fas = (SNES_FAS*) snes->data; 451 452 PetscFunctionBegin; 453 ierr = SNESFASCycleGetSmootherUp(snes, &smoothu);CHKERRQ(ierr); 454 if (fas->eventsmoothsolve) {ierr = PetscLogEventBegin(fas->eventsmoothsolve,smoothu,0,0,0);CHKERRQ(ierr);} 455 ierr = SNESSolve(smoothu, B, X);CHKERRQ(ierr); 456 if (fas->eventsmoothsolve) {ierr = PetscLogEventEnd(fas->eventsmoothsolve,smoothu,0,0,0);CHKERRQ(ierr);} 457 /* check convergence reason for the smoother */ 458 ierr = SNESGetConvergedReason(smoothu,&reason);CHKERRQ(ierr); 459 if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) { 460 snes->reason = SNES_DIVERGED_INNER; 461 PetscFunctionReturn(0); 462 } 463 ierr = SNESGetFunction(smoothu, &FPC, NULL, NULL);CHKERRQ(ierr); 464 ierr = SNESGetAlwaysComputesFinalResidual(smoothu, &flg);CHKERRQ(ierr); 465 if (!flg) { 466 ierr = SNESComputeFunction(smoothu, X, FPC);CHKERRQ(ierr); 467 } 468 ierr = VecCopy(FPC, F);CHKERRQ(ierr); 469 if (fnorm) {ierr = VecNorm(F,NORM_2,fnorm);CHKERRQ(ierr);} 470 PetscFunctionReturn(0); 471 } 472 473 /*@ 474 SNESFASCreateCoarseVec - create Vec corresponding to a state vector on one level coarser than current level 475 476 Collective 477 478 Input Arguments: 479 . snes - SNESFAS 480 481 Output Arguments: 482 . Xcoarse - vector on level one coarser than snes 483 484 Level: developer 485 486 .seealso: SNESFASSetRestriction(), SNESFASRestrict() 487 @*/ 488 PetscErrorCode SNESFASCreateCoarseVec(SNES snes,Vec *Xcoarse) 489 { 490 PetscErrorCode ierr; 491 SNES_FAS *fas = (SNES_FAS*)snes->data; 492 493 PetscFunctionBegin; 494 if (fas->rscale) { 495 ierr = VecDuplicate(fas->rscale,Xcoarse);CHKERRQ(ierr); 496 } else if (fas->inject) { 497 ierr = MatCreateVecs(fas->inject,Xcoarse,NULL);CHKERRQ(ierr); 498 } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or injection"); 499 PetscFunctionReturn(0); 500 } 501 502 /*@ 503 SNESFASRestrict - restrict a Vec to the next coarser level 504 505 Collective 506 507 Input Arguments: 508 + fine - SNES from which to restrict 509 - Xfine - vector to restrict 510 511 Output Arguments: 512 . Xcoarse - result of restriction 513 514 Level: developer 515 516 .seealso: SNESFASSetRestriction(), SNESFASSetInjection() 517 @*/ 518 PetscErrorCode SNESFASRestrict(SNES fine,Vec Xfine,Vec Xcoarse) 519 { 520 PetscErrorCode ierr; 521 SNES_FAS *fas = (SNES_FAS*)fine->data; 522 523 PetscFunctionBegin; 524 PetscValidHeaderSpecific(fine,SNES_CLASSID,1); 525 PetscValidHeaderSpecific(Xfine,VEC_CLASSID,2); 526 PetscValidHeaderSpecific(Xcoarse,VEC_CLASSID,3); 527 if (fas->inject) { 528 ierr = MatRestrict(fas->inject,Xfine,Xcoarse);CHKERRQ(ierr); 529 } else { 530 ierr = MatRestrict(fas->restrct,Xfine,Xcoarse);CHKERRQ(ierr); 531 ierr = VecPointwiseMult(Xcoarse,fas->rscale,Xcoarse);CHKERRQ(ierr); 532 } 533 PetscFunctionReturn(0); 534 } 535 536 /* 537 538 Performs the FAS coarse correction as: 539 540 fine problem: F(x) = b 541 coarse problem: F^c(x^c) = b^c 542 543 b^c = F^c(Rx) - R(F(x) - b) 544 545 */ 546 PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new) 547 { 548 PetscErrorCode ierr; 549 Vec X_c, Xo_c, F_c, B_c; 550 SNESConvergedReason reason; 551 SNES next; 552 Mat restrct, interpolate; 553 SNES_FAS *fasc; 554 555 PetscFunctionBegin; 556 ierr = SNESFASCycleGetCorrection(snes, &next);CHKERRQ(ierr); 557 if (next) { 558 fasc = (SNES_FAS*)next->data; 559 560 ierr = SNESFASCycleGetRestriction(snes, &restrct);CHKERRQ(ierr); 561 ierr = SNESFASCycleGetInterpolation(snes, &interpolate);CHKERRQ(ierr); 562 563 X_c = next->vec_sol; 564 Xo_c = next->work[0]; 565 F_c = next->vec_func; 566 B_c = next->vec_rhs; 567 568 if (fasc->eventinterprestrict) {ierr = PetscLogEventBegin(fasc->eventinterprestrict,snes,0,0,0);CHKERRQ(ierr);} 569 ierr = SNESFASRestrict(snes, X, Xo_c);CHKERRQ(ierr); 570 /* restrict the defect: R(F(x) - b) */ 571 ierr = MatRestrict(restrct, F, B_c);CHKERRQ(ierr); 572 if (fasc->eventinterprestrict) {ierr = PetscLogEventEnd(fasc->eventinterprestrict,snes,0,0,0);CHKERRQ(ierr);} 573 574 if (fasc->eventresidual) {ierr = PetscLogEventBegin(fasc->eventresidual,next,0,0,0);CHKERRQ(ierr);} 575 /* F_c = F^c(Rx) - R(F(x) - b) since the second term was sitting in next->vec_rhs */ 576 ierr = SNESComputeFunction(next, Xo_c, F_c);CHKERRQ(ierr); 577 if (fasc->eventresidual) {ierr = PetscLogEventEnd(fasc->eventresidual,next,0,0,0);CHKERRQ(ierr);} 578 579 /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */ 580 ierr = VecCopy(B_c, X_c);CHKERRQ(ierr); 581 ierr = VecCopy(F_c, B_c);CHKERRQ(ierr); 582 ierr = VecCopy(X_c, F_c);CHKERRQ(ierr); 583 /* set initial guess of the coarse problem to the projected fine solution */ 584 ierr = VecCopy(Xo_c, X_c);CHKERRQ(ierr); 585 586 /* recurse to the next level */ 587 ierr = SNESSetInitialFunction(next, F_c);CHKERRQ(ierr); 588 ierr = SNESSolve(next, B_c, X_c);CHKERRQ(ierr); 589 ierr = SNESGetConvergedReason(next,&reason);CHKERRQ(ierr); 590 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 591 snes->reason = SNES_DIVERGED_INNER; 592 PetscFunctionReturn(0); 593 } 594 /* correct as x <- x + I(x^c - Rx)*/ 595 ierr = VecAXPY(X_c, -1.0, Xo_c);CHKERRQ(ierr); 596 597 if (fasc->eventinterprestrict) {ierr = PetscLogEventBegin(fasc->eventinterprestrict,snes,0,0,0);CHKERRQ(ierr);} 598 ierr = MatInterpolateAdd(interpolate, X_c, X, X_new);CHKERRQ(ierr); 599 ierr = PetscObjectSetName((PetscObject) X_c, "Coarse correction");CHKERRQ(ierr); 600 ierr = VecViewFromOptions(X_c, NULL, "-fas_coarse_solution_view");CHKERRQ(ierr); 601 ierr = PetscObjectSetName((PetscObject) X_new, "Updated Fine solution");CHKERRQ(ierr); 602 ierr = VecViewFromOptions(X_new, NULL, "-fas_levels_1_solution_view");CHKERRQ(ierr); 603 if (fasc->eventinterprestrict) {ierr = PetscLogEventEnd(fasc->eventinterprestrict,snes,0,0,0);CHKERRQ(ierr);} 604 } 605 PetscFunctionReturn(0); 606 } 607 608 /* 609 610 The additive cycle looks like: 611 612 xhat = x 613 xhat = dS(x, b) 614 x = coarsecorrection(xhat, b_d) 615 x = x + nu*(xhat - x); 616 (optional) x = uS(x, b) 617 618 With the coarse RHS (defect correction) as below. 619 620 */ 621 static PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X) 622 { 623 Vec F, B, Xhat; 624 Vec X_c, Xo_c, F_c, B_c; 625 PetscErrorCode ierr; 626 SNESConvergedReason reason; 627 PetscReal xnorm, fnorm, ynorm; 628 SNESLineSearchReason lsresult; 629 SNES next; 630 Mat restrct, interpolate; 631 SNES_FAS *fas = (SNES_FAS*)snes->data,*fasc; 632 633 PetscFunctionBegin; 634 ierr = SNESFASCycleGetCorrection(snes, &next);CHKERRQ(ierr); 635 F = snes->vec_func; 636 B = snes->vec_rhs; 637 Xhat = snes->work[1]; 638 ierr = VecCopy(X, Xhat);CHKERRQ(ierr); 639 /* recurse first */ 640 if (next) { 641 fasc = (SNES_FAS*)next->data; 642 ierr = SNESFASCycleGetRestriction(snes, &restrct);CHKERRQ(ierr); 643 ierr = SNESFASCycleGetInterpolation(snes, &interpolate);CHKERRQ(ierr); 644 if (fas->eventresidual) {ierr = PetscLogEventBegin(fas->eventresidual,snes,0,0,0);CHKERRQ(ierr);} 645 ierr = SNESComputeFunction(snes, Xhat, F);CHKERRQ(ierr); 646 if (fas->eventresidual) {ierr = PetscLogEventEnd(fas->eventresidual,snes,0,0,0);CHKERRQ(ierr);} 647 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); 648 X_c = next->vec_sol; 649 Xo_c = next->work[0]; 650 F_c = next->vec_func; 651 B_c = next->vec_rhs; 652 653 ierr = SNESFASRestrict(snes,Xhat,Xo_c);CHKERRQ(ierr); 654 /* restrict the defect */ 655 ierr = MatRestrict(restrct, F, B_c);CHKERRQ(ierr); 656 657 /* solve the coarse problem corresponding to F^c(x^c) = b^c = Rb + F^c(Rx) - RF(x) */ 658 if (fasc->eventresidual) {ierr = PetscLogEventBegin(fasc->eventresidual,next,0,0,0);CHKERRQ(ierr);} 659 ierr = SNESComputeFunction(next, Xo_c, F_c);CHKERRQ(ierr); 660 if (fasc->eventresidual) {ierr = PetscLogEventEnd(fasc->eventresidual,next,0,0,0);CHKERRQ(ierr);} 661 ierr = VecCopy(B_c, X_c);CHKERRQ(ierr); 662 ierr = VecCopy(F_c, B_c);CHKERRQ(ierr); 663 ierr = VecCopy(X_c, F_c);CHKERRQ(ierr); 664 /* set initial guess of the coarse problem to the projected fine solution */ 665 ierr = VecCopy(Xo_c, X_c);CHKERRQ(ierr); 666 667 /* recurse */ 668 ierr = SNESSetInitialFunction(next, F_c);CHKERRQ(ierr); 669 ierr = SNESSolve(next, B_c, X_c);CHKERRQ(ierr); 670 671 /* smooth on this level */ 672 ierr = SNESFASDownSmooth_Private(snes, B, X, F, &fnorm);CHKERRQ(ierr); 673 674 ierr = SNESGetConvergedReason(next,&reason);CHKERRQ(ierr); 675 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 676 snes->reason = SNES_DIVERGED_INNER; 677 PetscFunctionReturn(0); 678 } 679 680 /* correct as x <- x + I(x^c - Rx)*/ 681 ierr = VecAYPX(X_c, -1.0, Xo_c);CHKERRQ(ierr); 682 ierr = MatInterpolate(interpolate, X_c, Xhat);CHKERRQ(ierr); 683 684 /* additive correction of the coarse direction*/ 685 ierr = SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat);CHKERRQ(ierr); 686 ierr = SNESLineSearchGetReason(snes->linesearch, &lsresult);CHKERRQ(ierr); 687 ierr = SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm);CHKERRQ(ierr); 688 if (lsresult) { 689 if (++snes->numFailures >= snes->maxFailures) { 690 snes->reason = SNES_DIVERGED_LINE_SEARCH; 691 PetscFunctionReturn(0); 692 } 693 } 694 } else { 695 ierr = SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);CHKERRQ(ierr); 696 } 697 PetscFunctionReturn(0); 698 } 699 700 /* 701 702 Defines the FAS cycle as: 703 704 fine problem: F(x) = b 705 coarse problem: F^c(x) = b^c 706 707 b^c = F^c(Rx) - R(F(x) - b) 708 709 correction: 710 711 x = x + I(x^c - Rx) 712 713 */ 714 static PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X) 715 { 716 717 PetscErrorCode ierr; 718 Vec F,B; 719 SNES next; 720 721 PetscFunctionBegin; 722 F = snes->vec_func; 723 B = snes->vec_rhs; 724 /* pre-smooth -- just update using the pre-smoother */ 725 ierr = SNESFASCycleGetCorrection(snes,&next);CHKERRQ(ierr); 726 ierr = SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);CHKERRQ(ierr); 727 if (next) { 728 ierr = SNESFASCoarseCorrection(snes, X, F, X);CHKERRQ(ierr); 729 ierr = SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);CHKERRQ(ierr); 730 } 731 PetscFunctionReturn(0); 732 } 733 734 static PetscErrorCode SNESFASCycleSetupPhase_Full(SNES snes) 735 { 736 SNES next; 737 SNES_FAS *fas = (SNES_FAS*)snes->data; 738 PetscBool isFine; 739 PetscErrorCode ierr; 740 741 PetscFunctionBegin; 742 /* pre-smooth -- just update using the pre-smoother */ 743 ierr = SNESFASCycleIsFine(snes,&isFine);CHKERRQ(ierr); 744 ierr = SNESFASCycleGetCorrection(snes,&next);CHKERRQ(ierr); 745 fas->full_stage = 0; 746 if (next) {ierr = SNESFASCycleSetupPhase_Full(next);CHKERRQ(ierr);} 747 PetscFunctionReturn(0); 748 } 749 750 static PetscErrorCode SNESFASCycle_Full(SNES snes, Vec X) 751 { 752 PetscErrorCode ierr; 753 Vec F,B; 754 SNES_FAS *fas = (SNES_FAS*)snes->data; 755 PetscBool isFine; 756 SNES next; 757 758 PetscFunctionBegin; 759 F = snes->vec_func; 760 B = snes->vec_rhs; 761 ierr = SNESFASCycleIsFine(snes,&isFine);CHKERRQ(ierr); 762 ierr = SNESFASCycleGetCorrection(snes,&next);CHKERRQ(ierr); 763 764 if (isFine) { 765 ierr = SNESFASCycleSetupPhase_Full(snes);CHKERRQ(ierr); 766 } 767 768 if (fas->full_stage == 0) { 769 /* downsweep */ 770 if (next) { 771 if (fas->level != 1) next->max_its += 1; 772 if (fas->full_downsweep) {ierr = SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);CHKERRQ(ierr);} 773 fas->full_downsweep = PETSC_TRUE; 774 ierr = SNESFASCoarseCorrection(snes,X,F,X);CHKERRQ(ierr); 775 ierr = SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);CHKERRQ(ierr); 776 if (fas->level != 1) next->max_its -= 1; 777 } else { 778 /* The smoother on the coarse level is the coarse solver */ 779 ierr = SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);CHKERRQ(ierr); 780 } 781 fas->full_stage = 1; 782 } else if (fas->full_stage == 1) { 783 if (snes->iter == 0) {ierr = SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);CHKERRQ(ierr);} 784 if (next) { 785 ierr = SNESFASCoarseCorrection(snes,X,F,X);CHKERRQ(ierr); 786 ierr = SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);CHKERRQ(ierr); 787 } 788 } 789 /* final v-cycle */ 790 if (isFine) { 791 if (next) { 792 ierr = SNESFASCoarseCorrection(snes,X,F,X);CHKERRQ(ierr); 793 ierr = SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);CHKERRQ(ierr); 794 } 795 } 796 PetscFunctionReturn(0); 797 } 798 799 static PetscErrorCode SNESFASCycle_Kaskade(SNES snes, Vec X) 800 { 801 PetscErrorCode ierr; 802 Vec F,B; 803 SNES next; 804 805 PetscFunctionBegin; 806 F = snes->vec_func; 807 B = snes->vec_rhs; 808 ierr = SNESFASCycleGetCorrection(snes,&next);CHKERRQ(ierr); 809 if (next) { 810 ierr = SNESFASCoarseCorrection(snes,X,F,X);CHKERRQ(ierr); 811 ierr = SNESFASUpSmooth_Private(snes,B,X,F,&snes->norm);CHKERRQ(ierr); 812 } else { 813 ierr = SNESFASDownSmooth_Private(snes,B,X,F,&snes->norm);CHKERRQ(ierr); 814 } 815 PetscFunctionReturn(0); 816 } 817 818 PetscBool SNEScite = PETSC_FALSE; 819 const char SNESCitation[] = "@techreport{pbmkbsxt2012,\n" 820 " title = {Composing Scalable Nonlinear Algebraic Solvers},\n" 821 " author = {Peter Brune and Mathew Knepley and Barry Smith and Xuemin Tu},\n" 822 " year = 2013,\n" 823 " type = Preprint,\n" 824 " number = {ANL/MCS-P2010-0112},\n" 825 " institution = {Argonne National Laboratory}\n}\n"; 826 827 static PetscErrorCode SNESSolve_FAS(SNES snes) 828 { 829 PetscErrorCode ierr; 830 PetscInt i, maxits; 831 Vec X, F; 832 PetscReal fnorm; 833 SNES_FAS *fas = (SNES_FAS*)snes->data,*ffas; 834 DM dm; 835 PetscBool isFine; 836 837 PetscFunctionBegin; 838 839 if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name); 840 841 ierr = PetscCitationsRegister(SNESCitation,&SNEScite);CHKERRQ(ierr); 842 maxits = snes->max_its; /* maximum number of iterations */ 843 snes->reason = SNES_CONVERGED_ITERATING; 844 X = snes->vec_sol; 845 F = snes->vec_func; 846 847 ierr = SNESFASCycleIsFine(snes, &isFine);CHKERRQ(ierr); 848 /*norm setup */ 849 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 850 snes->iter = 0; 851 snes->norm = 0.; 852 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 853 if (!snes->vec_func_init_set) { 854 if (fas->eventresidual) {ierr = PetscLogEventBegin(fas->eventresidual,snes,0,0,0);CHKERRQ(ierr);} 855 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 856 if (fas->eventresidual) {ierr = PetscLogEventEnd(fas->eventresidual,snes,0,0,0);CHKERRQ(ierr);} 857 } else snes->vec_func_init_set = PETSC_FALSE; 858 859 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 860 SNESCheckFunctionNorm(snes,fnorm); 861 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 862 snes->norm = fnorm; 863 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 864 ierr = SNESLogConvergenceHistory(snes,fnorm,0);CHKERRQ(ierr); 865 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 866 867 /* test convergence */ 868 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 869 if (snes->reason) PetscFunctionReturn(0); 870 871 872 if (isFine) { 873 /* propagate scale-dependent data up the hierarchy */ 874 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 875 for (ffas=fas; ffas->next; ffas=(SNES_FAS*)ffas->next->data) { 876 DM dmcoarse; 877 ierr = SNESGetDM(ffas->next,&dmcoarse);CHKERRQ(ierr); 878 ierr = DMRestrict(dm,ffas->restrct,ffas->rscale,ffas->inject,dmcoarse);CHKERRQ(ierr); 879 dm = dmcoarse; 880 } 881 } 882 883 for (i = 0; i < maxits; i++) { 884 /* Call general purpose update function */ 885 886 if (snes->ops->update) { 887 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 888 } 889 if (fas->fastype == SNES_FAS_MULTIPLICATIVE) { 890 ierr = SNESFASCycle_Multiplicative(snes, X);CHKERRQ(ierr); 891 } else if (fas->fastype == SNES_FAS_ADDITIVE) { 892 ierr = SNESFASCycle_Additive(snes, X);CHKERRQ(ierr); 893 } else if (fas->fastype == SNES_FAS_FULL) { 894 ierr = SNESFASCycle_Full(snes, X);CHKERRQ(ierr); 895 } else if (fas->fastype ==SNES_FAS_KASKADE) { 896 ierr = SNESFASCycle_Kaskade(snes, X);CHKERRQ(ierr); 897 } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported FAS type"); 898 899 /* check for FAS cycle divergence */ 900 if (snes->reason != SNES_CONVERGED_ITERATING) PetscFunctionReturn(0); 901 902 /* Monitor convergence */ 903 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 904 snes->iter = i+1; 905 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 906 ierr = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr); 907 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 908 /* Test for convergence */ 909 if (isFine) { 910 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,snes->norm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 911 if (snes->reason) break; 912 } 913 } 914 if (i == maxits) { 915 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);CHKERRQ(ierr); 916 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 917 } 918 PetscFunctionReturn(0); 919 } 920 921 /*MC 922 923 SNESFAS - Full Approximation Scheme nonlinear multigrid solver. 924 925 The nonlinear problem is solved by correction using coarse versions 926 of the nonlinear problem. This problem is perturbed so that a projected 927 solution of the fine problem elicits no correction from the coarse problem. 928 929 Options Database: 930 + -snes_fas_levels - The number of levels 931 . -snes_fas_cycles<1> - The number of cycles -- 1 for V, 2 for W 932 . -snes_fas_type<additive,multiplicative,full,kaskade> - Additive or multiplicative cycle 933 . -snes_fas_galerkin<PETSC_FALSE> - Form coarse problems by projection back upon the fine problem 934 . -snes_fas_smoothup<1> - The number of iterations of the post-smoother 935 . -snes_fas_smoothdown<1> - The number of iterations of the pre-smoother 936 . -snes_fas_monitor - Monitor progress of all of the levels 937 . -snes_fas_full_downsweep<PETSC_FALSE> - call the downsmooth on the initial downsweep of full FAS 938 . -fas_levels_snes_ - SNES options for all smoothers 939 . -fas_levels_cycle_snes_ - SNES options for all cycles 940 . -fas_levels_i_snes_ - SNES options for the smoothers on level i 941 . -fas_levels_i_cycle_snes_ - SNES options for the cycle on level i 942 - -fas_coarse_snes_ - SNES options for the coarsest smoother 943 944 Notes: 945 The organization of the FAS solver is slightly different from the organization of PCMG 946 As each level has smoother SNES instances(down and potentially up) and a cycle SNES instance. 947 The cycle SNES instance may be used for monitoring convergence on a particular level. 948 949 Level: beginner 950 951 References: 952 . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers", 953 SIAM Review, 57(4), 2015 954 955 .seealso: PCMG, SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 956 M*/ 957 958 PETSC_EXTERN PetscErrorCode SNESCreate_FAS(SNES snes) 959 { 960 SNES_FAS *fas; 961 PetscErrorCode ierr; 962 963 PetscFunctionBegin; 964 snes->ops->destroy = SNESDestroy_FAS; 965 snes->ops->setup = SNESSetUp_FAS; 966 snes->ops->setfromoptions = SNESSetFromOptions_FAS; 967 snes->ops->view = SNESView_FAS; 968 snes->ops->solve = SNESSolve_FAS; 969 snes->ops->reset = SNESReset_FAS; 970 971 snes->usesksp = PETSC_FALSE; 972 snes->usesnpc = PETSC_FALSE; 973 974 if (!snes->tolerancesset) { 975 snes->max_funcs = 30000; 976 snes->max_its = 10000; 977 } 978 979 snes->alwayscomputesfinalresidual = PETSC_TRUE; 980 981 ierr = PetscNewLog(snes,&fas);CHKERRQ(ierr); 982 983 snes->data = (void*) fas; 984 fas->level = 0; 985 fas->levels = 1; 986 fas->n_cycles = 1; 987 fas->max_up_it = 1; 988 fas->max_down_it = 1; 989 fas->smoothu = NULL; 990 fas->smoothd = NULL; 991 fas->next = NULL; 992 fas->previous = NULL; 993 fas->fine = snes; 994 fas->interpolate = NULL; 995 fas->restrct = NULL; 996 fas->inject = NULL; 997 fas->usedmfornumberoflevels = PETSC_FALSE; 998 fas->fastype = SNES_FAS_MULTIPLICATIVE; 999 fas->full_downsweep = PETSC_FALSE; 1000 1001 fas->eventsmoothsetup = 0; 1002 fas->eventsmoothsolve = 0; 1003 fas->eventresidual = 0; 1004 fas->eventinterprestrict = 0; 1005 PetscFunctionReturn(0); 1006 } 1007