1 /* Defines the basic SNES object */ 2 #include <../src/snes/impls/fas/fasimpls.h> /*I "petscsnesfas.h" I*/ 3 4 const char *const SNESFASTypes[] = {"MULTIPLICATIVE","ADDITIVE","SNESFASType","SNES_FAS",0}; 5 6 extern PetscErrorCode SNESDestroy_FAS(SNES snes); 7 extern PetscErrorCode SNESSetUp_FAS(SNES snes); 8 extern PetscErrorCode SNESSetFromOptions_FAS(SNES snes); 9 extern PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer); 10 extern PetscErrorCode SNESSolve_FAS(SNES snes); 11 extern PetscErrorCode SNESReset_FAS(SNES snes); 12 extern PetscErrorCode SNESFASGalerkinDefaultFunction(SNES, Vec, Vec, void *); 13 extern PetscErrorCode SNESFASCycleCreateSmoother_Private(SNES, SNES *); 14 15 /*MC 16 17 SNESFAS - Full Approximation Scheme nonlinear multigrid solver. 18 19 The nonlinear problem is solved by correction using coarse versions 20 of the nonlinear problem. This problem is perturbed so that a projected 21 solution of the fine problem elicits no correction from the coarse problem. 22 23 Options Database: 24 + -snes_fas_levels - The number of levels 25 . -snes_fas_cycles<1> - The number of cycles -- 1 for V, 2 for W 26 . -snes_fas_type<additive, multiplicative> - Additive or multiplicative cycle 27 . -snes_fas_galerkin<PETSC_FALSE> - Form coarse problems by projection back upon the fine problem 28 . -snes_fas_smoothup<1> - The number of iterations of the post-smoother 29 . -snes_fas_smoothdown<1> - The number of iterations of the pre-smoother 30 . -snes_fas_monitor - Monitor progress of all of the levels 31 . -fas_levels_snes_ - SNES options for all smoothers 32 . -fas_levels_cycle_snes_ - SNES options for all cycles 33 . -fas_levels_i_snes_ - SNES options for the smoothers on level i 34 . -fas_levels_i_cycle_snes_ - SNES options for the cycle on level i 35 - -fas_coarse_snes_ - SNES options for the coarsest smoother 36 37 Notes: 38 The organization of the FAS solver is slightly different from the organization of PCMG 39 As each level has smoother SNES instances(down and potentially up) and a cycle SNES instance. 40 The cycle SNES instance may be used for monitoring convergence on a particular level. 41 42 Level: beginner 43 44 .seealso: PCMG, SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 45 M*/ 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "SNESCreate_FAS" 49 PETSC_EXTERN_C PetscErrorCode SNESCreate_FAS(SNES snes) 50 { 51 SNES_FAS * fas; 52 PetscErrorCode ierr; 53 54 PetscFunctionBegin; 55 snes->ops->destroy = SNESDestroy_FAS; 56 snes->ops->setup = SNESSetUp_FAS; 57 snes->ops->setfromoptions = SNESSetFromOptions_FAS; 58 snes->ops->view = SNESView_FAS; 59 snes->ops->solve = SNESSolve_FAS; 60 snes->ops->reset = SNESReset_FAS; 61 62 snes->usesksp = PETSC_FALSE; 63 snes->usespc = PETSC_FALSE; 64 65 if (!snes->tolerancesset) { 66 snes->max_funcs = 30000; 67 snes->max_its = 10000; 68 } 69 70 ierr = PetscNewLog(snes, SNES_FAS, &fas);CHKERRQ(ierr); 71 snes->data = (void*) fas; 72 fas->level = 0; 73 fas->levels = 1; 74 fas->n_cycles = 1; 75 fas->max_up_it = 1; 76 fas->max_down_it = 1; 77 fas->smoothu = PETSC_NULL; 78 fas->smoothd = PETSC_NULL; 79 fas->next = PETSC_NULL; 80 fas->previous = PETSC_NULL; 81 fas->fine = snes; 82 fas->interpolate = PETSC_NULL; 83 fas->restrct = PETSC_NULL; 84 fas->inject = PETSC_NULL; 85 fas->monitor = PETSC_NULL; 86 fas->usedmfornumberoflevels = PETSC_FALSE; 87 fas->fastype = SNES_FAS_MULTIPLICATIVE; 88 PetscFunctionReturn(0); 89 } 90 91 #undef __FUNCT__ 92 #define __FUNCT__ "SNESReset_FAS" 93 PetscErrorCode SNESReset_FAS(SNES snes) 94 { 95 PetscErrorCode ierr = 0; 96 SNES_FAS * fas = (SNES_FAS *)snes->data; 97 98 PetscFunctionBegin; 99 ierr = SNESDestroy(&fas->smoothu);CHKERRQ(ierr); 100 ierr = SNESDestroy(&fas->smoothd);CHKERRQ(ierr); 101 ierr = MatDestroy(&fas->inject);CHKERRQ(ierr); 102 ierr = MatDestroy(&fas->interpolate);CHKERRQ(ierr); 103 ierr = MatDestroy(&fas->restrct);CHKERRQ(ierr); 104 ierr = VecDestroy(&fas->rscale);CHKERRQ(ierr); 105 if (fas->next) ierr = SNESReset(fas->next);CHKERRQ(ierr); 106 107 PetscFunctionReturn(0); 108 } 109 110 #undef __FUNCT__ 111 #define __FUNCT__ "SNESDestroy_FAS" 112 PetscErrorCode SNESDestroy_FAS(SNES snes) 113 { 114 SNES_FAS * fas = (SNES_FAS *)snes->data; 115 PetscErrorCode ierr = 0; 116 117 PetscFunctionBegin; 118 /* recursively resets and then destroys */ 119 ierr = SNESReset(snes);CHKERRQ(ierr); 120 if (fas->next) ierr = SNESDestroy(&fas->next);CHKERRQ(ierr); 121 ierr = PetscFree(fas);CHKERRQ(ierr); 122 123 PetscFunctionReturn(0); 124 } 125 126 #undef __FUNCT__ 127 #define __FUNCT__ "SNESSetUp_FAS" 128 PetscErrorCode SNESSetUp_FAS(SNES snes) 129 { 130 SNES_FAS *fas = (SNES_FAS *) snes->data; 131 PetscErrorCode ierr; 132 VecScatter injscatter; 133 PetscInt dm_levels; 134 Vec vec_sol, vec_func, vec_sol_update, vec_rhs; /* preserve these if they're set through the reset */ 135 SNES next; 136 PetscBool isFine; 137 SNESLineSearch linesearch; 138 SNESLineSearch slinesearch; 139 void *lsprectx,*lspostctx; 140 PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*); 141 PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool *,PetscBool *,void*); 142 143 PetscFunctionBegin; 144 ierr = SNESFASCycleIsFine(snes, &isFine);CHKERRQ(ierr); 145 if (fas->usedmfornumberoflevels && isFine) { 146 ierr = DMGetRefineLevel(snes->dm,&dm_levels);CHKERRQ(ierr); 147 dm_levels++; 148 if (dm_levels > fas->levels) { 149 /* we don't want the solution and func vectors to be destroyed in the SNESReset when it's called in SNESFASSetLevels_FAS*/ 150 vec_sol = snes->vec_sol; 151 vec_func = snes->vec_func; 152 vec_sol_update = snes->vec_sol_update; 153 vec_rhs = snes->vec_rhs; 154 snes->vec_sol = PETSC_NULL; 155 snes->vec_func = PETSC_NULL; 156 snes->vec_sol_update = PETSC_NULL; 157 snes->vec_rhs = PETSC_NULL; 158 159 /* reset the number of levels */ 160 ierr = SNESFASSetLevels(snes,dm_levels,PETSC_NULL);CHKERRQ(ierr); 161 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 162 163 snes->vec_sol = vec_sol; 164 snes->vec_func = vec_func; 165 snes->vec_rhs = vec_rhs; 166 snes->vec_sol_update = vec_sol_update; 167 } 168 } 169 ierr = SNESFASCycleGetCorrection(snes, &next);CHKERRQ(ierr); 170 if (!isFine) snes->gridsequence = 0; /* no grid sequencing inside the multigrid hierarchy! */ 171 172 if (fas->fastype == SNES_FAS_MULTIPLICATIVE) { 173 ierr = SNESDefaultGetWork(snes, 2);CHKERRQ(ierr); /* work vectors used for intergrid transfers */ 174 } else { 175 ierr = SNESDefaultGetWork(snes, 2);CHKERRQ(ierr); /* work vectors used for intergrid transfers */ 176 } 177 178 /* set up the smoothers if they haven't already been set up */ 179 if (!fas->smoothd) { 180 ierr = SNESFASCycleCreateSmoother_Private(snes, &fas->smoothd);CHKERRQ(ierr); 181 } 182 183 if (snes->dm) { 184 /* set the smoother DMs properly */ 185 if (fas->smoothu) ierr = SNESSetDM(fas->smoothu, snes->dm);CHKERRQ(ierr); 186 ierr = SNESSetDM(fas->smoothd, snes->dm);CHKERRQ(ierr); 187 /* construct EVERYTHING from the DM -- including the progressive set of smoothers */ 188 if (next) { 189 /* for now -- assume the DM and the evaluation functions have been set externally */ 190 if (!next->dm) { 191 ierr = DMCoarsen(snes->dm, ((PetscObject)next)->comm, &next->dm);CHKERRQ(ierr); 192 ierr = SNESSetDM(next, next->dm);CHKERRQ(ierr); 193 } 194 /* set the interpolation and restriction from the DM */ 195 if (!fas->interpolate) { 196 ierr = DMCreateInterpolation(next->dm, snes->dm, &fas->interpolate, &fas->rscale);CHKERRQ(ierr); 197 if (!fas->restrct) { 198 ierr = PetscObjectReference((PetscObject)fas->interpolate);CHKERRQ(ierr); 199 fas->restrct = fas->interpolate; 200 } 201 } 202 /* set the injection from the DM */ 203 if (!fas->inject) { 204 ierr = DMCreateInjection(next->dm, snes->dm, &injscatter);CHKERRQ(ierr); 205 ierr = MatCreateScatter(((PetscObject)snes)->comm, injscatter, &fas->inject);CHKERRQ(ierr); 206 ierr = VecScatterDestroy(&injscatter);CHKERRQ(ierr); 207 } 208 } 209 } 210 /*pass the smoother, function, and jacobian up to the next level if it's not user set already */ 211 if (fas->galerkin) { 212 if (next) 213 ierr = SNESSetFunction(next, PETSC_NULL, SNESFASGalerkinDefaultFunction, next);CHKERRQ(ierr); 214 if (fas->smoothd && fas->level != fas->levels - 1) ierr = SNESSetFunction(fas->smoothd, PETSC_NULL, SNESFASGalerkinDefaultFunction, snes);CHKERRQ(ierr); 215 if (fas->smoothu && fas->level != fas->levels - 1) ierr = SNESSetFunction(fas->smoothu, PETSC_NULL, SNESFASGalerkinDefaultFunction, snes);CHKERRQ(ierr); 216 } 217 218 /* sets the down (pre) smoother's default norm and sets it from options */ 219 if (fas->smoothd) { 220 if (fas->level == 0 && fas->levels != 1) { 221 ierr = SNESSetNormType(fas->smoothd, SNES_NORM_NONE);CHKERRQ(ierr); 222 } else { 223 ierr = SNESSetNormType(fas->smoothd, SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 224 } 225 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothd);CHKERRQ(ierr); 226 ierr = SNESSetFromOptions(fas->smoothd);CHKERRQ(ierr); 227 ierr = SNESGetSNESLineSearch(snes,&linesearch);CHKERRQ(ierr); 228 ierr = SNESGetSNESLineSearch(fas->smoothd,&slinesearch);CHKERRQ(ierr); 229 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 230 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 231 ierr = SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);CHKERRQ(ierr); 232 ierr = SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);CHKERRQ(ierr); 233 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);CHKERRQ(ierr); 234 } 235 236 /* sets the up (post) smoother's default norm and sets it from options */ 237 if (fas->smoothu) { 238 if (fas->level != fas->levels - 1) { 239 ierr = SNESSetNormType(fas->smoothu, SNES_NORM_NONE);CHKERRQ(ierr); 240 } else { 241 ierr = SNESSetNormType(fas->smoothu, SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 242 } 243 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothu);CHKERRQ(ierr); 244 ierr = SNESSetFromOptions(fas->smoothu);CHKERRQ(ierr); 245 ierr = SNESGetSNESLineSearch(snes,&linesearch);CHKERRQ(ierr); 246 ierr = SNESGetSNESLineSearch(fas->smoothu,&slinesearch);CHKERRQ(ierr); 247 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 248 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 249 ierr = SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);CHKERRQ(ierr); 250 ierr = SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);CHKERRQ(ierr); 251 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);CHKERRQ(ierr); 252 } 253 254 if (next) { 255 /* gotta set up the solution vector for this to work */ 256 if (!next->vec_sol) {ierr = SNESFASCreateCoarseVec(snes,&next->vec_sol);CHKERRQ(ierr);} 257 if (!next->vec_rhs) {ierr = SNESFASCreateCoarseVec(snes,&next->vec_rhs);CHKERRQ(ierr);} 258 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next);CHKERRQ(ierr); 259 ierr = SNESGetSNESLineSearch(snes,&linesearch);CHKERRQ(ierr); 260 ierr = SNESGetSNESLineSearch(fas->next,&slinesearch);CHKERRQ(ierr); 261 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 262 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 263 ierr = SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);CHKERRQ(ierr); 264 ierr = SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);CHKERRQ(ierr); 265 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);CHKERRQ(ierr); 266 ierr = SNESSetUp(next);CHKERRQ(ierr); 267 } 268 /* setup FAS work vectors */ 269 if (fas->galerkin) { 270 ierr = VecDuplicate(snes->vec_sol, &fas->Xg);CHKERRQ(ierr); 271 ierr = VecDuplicate(snes->vec_sol, &fas->Fg);CHKERRQ(ierr); 272 } 273 PetscFunctionReturn(0); 274 } 275 276 #undef __FUNCT__ 277 #define __FUNCT__ "SNESSetFromOptions_FAS" 278 PetscErrorCode SNESSetFromOptions_FAS(SNES snes) 279 { 280 SNES_FAS *fas = (SNES_FAS *) snes->data; 281 PetscInt levels = 1; 282 PetscBool flg = PETSC_FALSE, upflg = PETSC_FALSE, downflg = PETSC_FALSE, monflg = PETSC_FALSE, galerkinflg = PETSC_FALSE; 283 PetscErrorCode ierr; 284 char monfilename[PETSC_MAX_PATH_LEN]; 285 SNESFASType fastype; 286 const char *optionsprefix; 287 SNESLineSearch linesearch; 288 PetscInt m, n_up, n_down; 289 SNES next; 290 PetscBool isFine; 291 292 PetscFunctionBegin; 293 ierr = SNESFASCycleIsFine(snes, &isFine);CHKERRQ(ierr); 294 ierr = PetscOptionsHead("SNESFAS Options-----------------------------------");CHKERRQ(ierr); 295 296 /* number of levels -- only process most options on the finest level */ 297 if (isFine) { 298 ierr = PetscOptionsInt("-snes_fas_levels", "Number of Levels", "SNESFASSetLevels", levels, &levels, &flg);CHKERRQ(ierr); 299 if (!flg && snes->dm) { 300 ierr = DMGetRefineLevel(snes->dm,&levels);CHKERRQ(ierr); 301 levels++; 302 fas->usedmfornumberoflevels = PETSC_TRUE; 303 } 304 ierr = SNESFASSetLevels(snes, levels, PETSC_NULL);CHKERRQ(ierr); 305 fastype = fas->fastype; 306 ierr = PetscOptionsEnum("-snes_fas_type","FAS correction type","SNESFASSetType",SNESFASTypes,(PetscEnum)fastype,(PetscEnum*)&fastype,&flg);CHKERRQ(ierr); 307 if (flg) { 308 ierr = SNESFASSetType(snes, fastype);CHKERRQ(ierr); 309 } 310 311 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 312 ierr = PetscOptionsInt("-snes_fas_cycles","Number of cycles","SNESFASSetCycles",fas->n_cycles,&m,&flg);CHKERRQ(ierr); 313 if (flg) { 314 ierr = SNESFASSetCycles(snes, m);CHKERRQ(ierr); 315 } 316 317 ierr = PetscOptionsBool("-snes_fas_galerkin", "Form coarse problems with Galerkin","SNESFASSetGalerkin",fas->galerkin,&galerkinflg,&flg);CHKERRQ(ierr); 318 if (flg) { 319 ierr = SNESFASSetGalerkin(snes, galerkinflg);CHKERRQ(ierr); 320 } 321 322 ierr = PetscOptionsInt("-snes_fas_smoothup","Number of post-smoothing steps","SNESFASSetNumberSmoothUp",fas->max_up_it,&n_up,&upflg);CHKERRQ(ierr); 323 324 ierr = PetscOptionsInt("-snes_fas_smoothdown","Number of pre-smoothing steps","SNESFASSetNumberSmoothDown",fas->max_down_it,&n_down,&downflg);CHKERRQ(ierr); 325 326 ierr = PetscOptionsString("-snes_fas_monitor","Monitor FAS progress","SNESFASSetMonitor","stdout",monfilename,PETSC_MAX_PATH_LEN,&monflg);CHKERRQ(ierr); 327 if (monflg) ierr = SNESFASSetMonitor(snes, PETSC_TRUE);CHKERRQ(ierr); 328 } 329 330 ierr = PetscOptionsTail();CHKERRQ(ierr); 331 /* setup from the determined types if there is no pointwise procedure or smoother defined */ 332 if (upflg) { 333 ierr = SNESFASSetNumberSmoothUp(snes,n_up);CHKERRQ(ierr); 334 } 335 if (downflg) { 336 ierr = SNESFASSetNumberSmoothDown(snes,n_down);CHKERRQ(ierr); 337 } 338 339 /* set up the default line search for coarse grid corrections */ 340 if (fas->fastype == SNES_FAS_ADDITIVE) { 341 if (!snes->linesearch) { 342 ierr = SNESGetSNESLineSearch(snes, &linesearch);CHKERRQ(ierr); 343 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);CHKERRQ(ierr); 344 } 345 } 346 347 ierr = SNESFASCycleGetCorrection(snes, &next);CHKERRQ(ierr); 348 /* recursive option setting for the smoothers */ 349 if (next) {ierr = SNESSetFromOptions(next);CHKERRQ(ierr);} 350 PetscFunctionReturn(0); 351 } 352 353 #undef __FUNCT__ 354 #define __FUNCT__ "SNESView_FAS" 355 PetscErrorCode SNESView_FAS(SNES snes, PetscViewer viewer) 356 { 357 SNES_FAS *fas = (SNES_FAS *) snes->data; 358 PetscBool isFine,iascii,isdraw; 359 PetscInt i; 360 PetscErrorCode ierr; 361 SNES smoothu, smoothd, levelsnes; 362 363 PetscFunctionBegin; 364 ierr = SNESFASCycleIsFine(snes, &isFine);CHKERRQ(ierr); 365 if (isFine) { 366 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 367 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 368 if (iascii) { 369 ierr = PetscViewerASCIIPrintf(viewer, "FAS: type is %s, levels=%D, cycles=%D\n", SNESFASTypes[fas->fastype], fas->levels, fas->n_cycles);CHKERRQ(ierr); 370 if (fas->galerkin) { 371 ierr = PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid function evaluation\n");CHKERRQ(ierr); 372 } else { 373 ierr = PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid function evaluation\n");CHKERRQ(ierr); 374 } 375 for (i=0; i<fas->levels; i++) { 376 ierr = SNESFASGetCycleSNES(snes, i, &levelsnes);CHKERRQ(ierr); 377 ierr = SNESFASCycleGetSmootherUp(levelsnes, &smoothu);CHKERRQ(ierr); 378 ierr = SNESFASCycleGetSmootherDown(levelsnes, &smoothd);CHKERRQ(ierr); 379 if (!i) { 380 ierr = PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level %D -------------------------------\n",i);CHKERRQ(ierr); 381 } else { 382 ierr = PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);CHKERRQ(ierr); 383 } 384 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 385 ierr = SNESView(smoothd,viewer);CHKERRQ(ierr); 386 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 387 if (i && (smoothd == smoothu)) { 388 ierr = PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");CHKERRQ(ierr); 389 } else if (i) { 390 ierr = PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);CHKERRQ(ierr); 391 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 392 ierr = SNESView(smoothu,viewer);CHKERRQ(ierr); 393 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 394 } 395 } 396 } else if (isdraw) { 397 PetscDraw draw; 398 PetscReal x,w,y,bottom,th,wth; 399 SNES_FAS *curfas = fas; 400 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 401 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 402 ierr = PetscDrawStringGetSize(draw,&wth,&th);CHKERRQ(ierr); 403 bottom = y - th; 404 while (curfas) { 405 if (!curfas->smoothu) { 406 ierr = PetscDrawPushCurrentPoint(draw,x,bottom);CHKERRQ(ierr); 407 if (curfas->smoothd) ierr = SNESView(curfas->smoothd,viewer);CHKERRQ(ierr); 408 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 409 } else { 410 w = 0.5*PetscMin(1.0-x,x); 411 ierr = PetscDrawPushCurrentPoint(draw,x-w,bottom);CHKERRQ(ierr); 412 if (curfas->smoothd) ierr = SNESView(curfas->smoothd,viewer);CHKERRQ(ierr); 413 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 414 ierr = PetscDrawPushCurrentPoint(draw,x+w,bottom);CHKERRQ(ierr); 415 if (curfas->smoothu) ierr = SNESView(curfas->smoothu,viewer);CHKERRQ(ierr); 416 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 417 } 418 /* this is totally bogus but we have no way of knowing how low the previous one was draw to */ 419 bottom -= 5*th; 420 if (curfas->next) { 421 curfas = (SNES_FAS*)curfas->next->data; 422 } else { 423 curfas = PETSC_NULL; 424 } 425 } 426 } else { 427 SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for SNESFAS",((PetscObject)viewer)->type_name); 428 } 429 } 430 PetscFunctionReturn(0); 431 } 432 433 #undef __FUNCT__ 434 #define __FUNCT__ "SNESFASDownSmooth_Private" 435 /* 436 Defines the action of the downsmoother 437 */ 438 PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm) 439 { 440 PetscErrorCode ierr = 0; 441 SNESConvergedReason reason; 442 Vec FPC; 443 SNES smoothd; 444 445 PetscFunctionBegin; 446 ierr = SNESFASCycleGetSmootherDown(snes, &smoothd);CHKERRQ(ierr); 447 ierr = SNESSetInitialFunction(smoothd, F);CHKERRQ(ierr); 448 ierr = SNESSetInitialFunctionNorm(smoothd, *fnorm);CHKERRQ(ierr); 449 ierr = SNESSolve(smoothd, B, X);CHKERRQ(ierr); 450 /* check convergence reason for the smoother */ 451 ierr = SNESGetConvergedReason(smoothd,&reason);CHKERRQ(ierr); 452 if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN)) { 453 snes->reason = SNES_DIVERGED_INNER; 454 PetscFunctionReturn(0); 455 } 456 ierr = SNESGetFunction(smoothd, &FPC, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 457 ierr = VecCopy(FPC, F);CHKERRQ(ierr); 458 ierr = SNESGetFunctionNorm(smoothd, fnorm);CHKERRQ(ierr); 459 PetscFunctionReturn(0); 460 } 461 462 463 #undef __FUNCT__ 464 #define __FUNCT__ "SNESFASUpSmooth_Private" 465 /* 466 Defines the action of the upsmoother 467 */ 468 PetscErrorCode SNESFASUpSmooth_Private (SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm) 469 { 470 PetscErrorCode ierr = 0; 471 SNESConvergedReason reason; 472 Vec FPC; 473 SNES smoothu; 474 475 PetscFunctionBegin; 476 ierr = SNESFASCycleGetSmootherUp(snes, &smoothu);CHKERRQ(ierr); 477 ierr = SNESSolve(smoothu, B, X);CHKERRQ(ierr); 478 /* check convergence reason for the smoother */ 479 ierr = SNESGetConvergedReason(smoothu,&reason);CHKERRQ(ierr); 480 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 481 snes->reason = SNES_DIVERGED_INNER; 482 PetscFunctionReturn(0); 483 } 484 ierr = SNESGetFunction(smoothu, &FPC, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 485 ierr = VecCopy(FPC, F);CHKERRQ(ierr); 486 ierr = SNESGetFunctionNorm(smoothu, fnorm);CHKERRQ(ierr); 487 PetscFunctionReturn(0); 488 } 489 490 #undef __FUNCT__ 491 #define __FUNCT__ "SNESFASCreateCoarseVec" 492 /*@ 493 SNESFASCreateCoarseVec - create Vec corresponding to a state vector on one level coarser than current level 494 495 Collective 496 497 Input Arguments: 498 . snes - SNESFAS 499 500 Output Arguments: 501 . Xcoarse - vector on level one coarser than snes 502 503 Level: developer 504 505 .seealso: SNESFASSetRestriction(), SNESFASRestrict() 506 @*/ 507 PetscErrorCode SNESFASCreateCoarseVec(SNES snes,Vec *Xcoarse) 508 { 509 PetscErrorCode ierr; 510 SNES_FAS *fas = (SNES_FAS*)snes->data; 511 512 PetscFunctionBegin; 513 if (fas->rscale) {ierr = VecDuplicate(fas->rscale,Xcoarse);CHKERRQ(ierr);} 514 else if (fas->inject) {ierr = MatGetVecs(fas->inject,Xcoarse,PETSC_NULL);CHKERRQ(ierr);} 515 else SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or injection");CHKERRQ(ierr); 516 PetscFunctionReturn(0); 517 } 518 519 #undef __FUNCT__ 520 #define __FUNCT__ "SNESFASRestrict" 521 /*@ 522 SNESFASRestrict - restrict a Vec to the next coarser level 523 524 Collective 525 526 Input Arguments: 527 + fine - SNES from which to restrict 528 - Xfine - vector to restrict 529 530 Output Arguments: 531 . Xcoarse - result of restriction 532 533 Level: developer 534 535 .seealso: SNESFASSetRestriction(), SNESFASSetInjection() 536 @*/ 537 PetscErrorCode SNESFASRestrict(SNES fine,Vec Xfine,Vec Xcoarse) 538 { 539 PetscErrorCode ierr; 540 SNES_FAS *fas = (SNES_FAS*)fine->data; 541 542 PetscFunctionBegin; 543 PetscValidHeaderSpecific(fine,SNES_CLASSID,1); 544 PetscValidHeaderSpecific(Xfine,VEC_CLASSID,2); 545 PetscValidHeaderSpecific(Xcoarse,VEC_CLASSID,3); 546 if (fas->inject) { 547 ierr = MatRestrict(fas->inject,Xfine,Xcoarse);CHKERRQ(ierr); 548 } else { 549 ierr = MatRestrict(fas->restrct,Xfine,Xcoarse);CHKERRQ(ierr); 550 ierr = VecPointwiseMult(Xcoarse,fas->rscale,Xcoarse);CHKERRQ(ierr); 551 } 552 PetscFunctionReturn(0); 553 } 554 555 #undef __FUNCT__ 556 #define __FUNCT__ "SNESFASCoarseCorrection" 557 /* 558 559 Performs the FAS coarse correction as: 560 561 fine problem: F(x) = 0 562 coarse problem: F^c(x) = b^c 563 564 b^c = F^c(I^c_fx^f - I^c_fF(x)) 565 566 */ 567 PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new) 568 { 569 PetscErrorCode ierr; 570 Vec X_c, Xo_c, F_c, B_c; 571 SNESConvergedReason reason; 572 SNES next; 573 Mat restrct, interpolate; 574 PetscFunctionBegin; 575 ierr = SNESFASCycleGetCorrection(snes, &next);CHKERRQ(ierr); 576 if (next) { 577 ierr = SNESFASCycleGetRestriction(snes, &restrct);CHKERRQ(ierr); 578 ierr = SNESFASCycleGetInterpolation(snes, &interpolate);CHKERRQ(ierr); 579 580 X_c = next->vec_sol; 581 Xo_c = next->work[0]; 582 F_c = next->vec_func; 583 B_c = next->vec_rhs; 584 585 ierr = SNESFASRestrict(snes,X,Xo_c);CHKERRQ(ierr); 586 /* restrict the defect */ 587 ierr = MatRestrict(restrct, F, B_c);CHKERRQ(ierr); 588 /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */ 589 ierr = SNESComputeFunction(next, Xo_c, F_c);CHKERRQ(ierr); 590 ierr = VecCopy(B_c, X_c);CHKERRQ(ierr); 591 ierr = VecCopy(F_c, B_c);CHKERRQ(ierr); 592 ierr = VecCopy(X_c, F_c);CHKERRQ(ierr); 593 /* set initial guess of the coarse problem to the projected fine solution */ 594 ierr = VecCopy(Xo_c, X_c);CHKERRQ(ierr); 595 596 /* recurse to the next level */ 597 ierr = SNESSetInitialFunction(next, F_c);CHKERRQ(ierr); 598 ierr = SNESSolve(next, B_c, X_c);CHKERRQ(ierr); 599 ierr = SNESGetConvergedReason(next,&reason);CHKERRQ(ierr); 600 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 601 snes->reason = SNES_DIVERGED_INNER; 602 PetscFunctionReturn(0); 603 } 604 /* correct as x <- x + I(x^c - Rx)*/ 605 ierr = VecAXPY(X_c, -1.0, Xo_c);CHKERRQ(ierr); 606 ierr = MatInterpolateAdd(interpolate, X_c, X, X_new);CHKERRQ(ierr); 607 } 608 PetscFunctionReturn(0); 609 } 610 611 #undef __FUNCT__ 612 #define __FUNCT__ "SNESFASCycle_Additive" 613 /* 614 615 The additive cycle looks like: 616 617 xhat = x 618 xhat = dS(x, b) 619 x = coarsecorrection(xhat, b_d) 620 x = x + nu*(xhat - x); 621 (optional) x = uS(x, b) 622 623 With the coarse RHS (defect correction) as below. 624 625 */ 626 PetscErrorCode SNESFASCycle_Additive(SNES snes, Vec X) 627 { 628 Vec F, B, Xhat; 629 Vec X_c, Xo_c, F_c, B_c; 630 PetscErrorCode ierr; 631 SNESConvergedReason reason; 632 PetscReal xnorm, fnorm, ynorm; 633 PetscBool lssuccess; 634 SNES next; 635 Mat restrct, interpolate; 636 PetscFunctionBegin; 637 ierr = SNESFASCycleGetCorrection(snes, &next);CHKERRQ(ierr); 638 F = snes->vec_func; 639 B = snes->vec_rhs; 640 Xhat = snes->work[1]; 641 ierr = VecCopy(X, Xhat);CHKERRQ(ierr); 642 /* recurse first */ 643 if (next) { 644 ierr = SNESFASCycleGetRestriction(snes, &restrct);CHKERRQ(ierr); 645 ierr = SNESFASCycleGetInterpolation(snes, &interpolate);CHKERRQ(ierr); 646 ierr = SNESComputeFunction(snes, Xhat, F);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 ierr = SNESComputeFunction(next, Xo_c, F_c);CHKERRQ(ierr); 659 ierr = VecCopy(B_c, X_c);CHKERRQ(ierr); 660 ierr = VecCopy(F_c, B_c);CHKERRQ(ierr); 661 ierr = VecCopy(X_c, F_c);CHKERRQ(ierr); 662 /* set initial guess of the coarse problem to the projected fine solution */ 663 ierr = VecCopy(Xo_c, X_c);CHKERRQ(ierr); 664 665 /* recurse */ 666 ierr = SNESSetInitialFunction(next, F_c);CHKERRQ(ierr); 667 ierr = SNESSolve(next, B_c, X_c);CHKERRQ(ierr); 668 669 /* smooth on this level */ 670 ierr = SNESFASDownSmooth_Private(snes, B, X, F, &fnorm);CHKERRQ(ierr); 671 672 ierr = SNESGetConvergedReason(next,&reason);CHKERRQ(ierr); 673 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 674 snes->reason = SNES_DIVERGED_INNER; 675 PetscFunctionReturn(0); 676 } 677 678 /* correct as x <- x + I(x^c - Rx)*/ 679 ierr = VecAYPX(X_c, -1.0, Xo_c);CHKERRQ(ierr); 680 ierr = MatInterpolate(interpolate, X_c, Xhat);CHKERRQ(ierr); 681 682 /* additive correction of the coarse direction*/ 683 ierr = SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Xhat);CHKERRQ(ierr); 684 ierr = SNESLineSearchGetSuccess(snes->linesearch, &lssuccess);CHKERRQ(ierr); 685 if (!lssuccess) { 686 if (++snes->numFailures >= snes->maxFailures) { 687 snes->reason = SNES_DIVERGED_LINE_SEARCH; 688 PetscFunctionReturn(0); 689 } 690 } 691 ierr = SNESLineSearchGetNorms(snes->linesearch, &xnorm, &snes->norm, &ynorm);CHKERRQ(ierr); 692 } else { 693 ierr = SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);CHKERRQ(ierr); 694 } 695 PetscFunctionReturn(0); 696 } 697 698 #undef __FUNCT__ 699 #define __FUNCT__ "SNESFASCycle_Multiplicative" 700 /* 701 702 Defines the FAS cycle as: 703 704 fine problem: F(x) = 0 705 coarse problem: F^c(x) = b^c 706 707 b^c = F^c(I^c_fx^f - I^c_fF(x)) 708 709 correction: 710 711 x = x + I(x^c - Rx) 712 713 */ 714 PetscErrorCode SNESFASCycle_Multiplicative(SNES snes, Vec X) 715 { 716 717 PetscErrorCode ierr; 718 Vec F,B; 719 SNES_FAS *fas = (SNES_FAS *)snes->data; 720 721 PetscFunctionBegin; 722 F = snes->vec_func; 723 B = snes->vec_rhs; 724 /* pre-smooth -- just update using the pre-smoother */ 725 ierr = SNESFASDownSmooth_Private(snes, B, X, F, &snes->norm);CHKERRQ(ierr); 726 if (fas->level != 0) { 727 ierr = SNESFASCoarseCorrection(snes, X, F, X);CHKERRQ(ierr); 728 ierr = SNESFASUpSmooth_Private(snes, B, X, F, &snes->norm);CHKERRQ(ierr); 729 } 730 PetscFunctionReturn(0); 731 } 732 733 #undef __FUNCT__ 734 #define __FUNCT__ "SNESSolve_FAS" 735 736 PetscErrorCode SNESSolve_FAS(SNES snes) 737 { 738 PetscErrorCode ierr; 739 PetscInt i, maxits; 740 Vec X, F; 741 PetscReal fnorm; 742 SNES_FAS *fas = (SNES_FAS *)snes->data,*ffas; 743 DM dm; 744 PetscBool isFine; 745 746 PetscFunctionBegin; 747 maxits = snes->max_its; /* maximum number of iterations */ 748 snes->reason = SNES_CONVERGED_ITERATING; 749 X = snes->vec_sol; 750 F = snes->vec_func; 751 752 ierr = SNESFASCycleIsFine(snes, &isFine);CHKERRQ(ierr); 753 /*norm setup */ 754 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 755 snes->iter = 0; 756 snes->norm = 0.; 757 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 758 if (!snes->vec_func_init_set) { 759 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 760 if (snes->domainerror) { 761 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 762 PetscFunctionReturn(0); 763 } 764 } else { 765 snes->vec_func_init_set = PETSC_FALSE; 766 } 767 768 if (!snes->norm_init_set) { 769 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 770 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 771 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 772 } else { 773 fnorm = snes->norm_init; 774 snes->norm_init_set = PETSC_FALSE; 775 } 776 777 snes->norm = fnorm; 778 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 779 SNESLogConvHistory(snes,fnorm,0); 780 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 781 782 /* set parameter for default relative tolerance convergence test */ 783 snes->ttol = fnorm*snes->rtol; 784 /* test convergence */ 785 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 786 if (snes->reason) PetscFunctionReturn(0); 787 788 789 if (isFine) { 790 /* propagate scale-dependent data up the hierarchy */ 791 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 792 for (ffas=fas; ffas->next; ffas=(SNES_FAS*)ffas->next->data) { 793 DM dmcoarse; 794 ierr = SNESGetDM(ffas->next,&dmcoarse);CHKERRQ(ierr); 795 ierr = DMRestrict(dm,ffas->restrct,ffas->rscale,ffas->inject,dmcoarse);CHKERRQ(ierr); 796 dm = dmcoarse; 797 } 798 } 799 800 for (i = 0; i < maxits; i++) { 801 /* Call general purpose update function */ 802 803 if (snes->ops->update) { 804 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 805 } 806 if (fas->fastype == SNES_FAS_MULTIPLICATIVE) { 807 ierr = SNESFASCycle_Multiplicative(snes, X);CHKERRQ(ierr); 808 } else { 809 ierr = SNESFASCycle_Additive(snes, X);CHKERRQ(ierr); 810 } 811 812 /* check for FAS cycle divergence */ 813 if (snes->reason != SNES_CONVERGED_ITERATING) { 814 PetscFunctionReturn(0); 815 } 816 817 /* Monitor convergence */ 818 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 819 snes->iter = i+1; 820 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 821 SNESLogConvHistory(snes,snes->norm,0); 822 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 823 /* Test for convergence */ 824 if (isFine) { 825 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,snes->norm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 826 if (snes->reason) break; 827 } 828 } 829 if (i == maxits) { 830 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);CHKERRQ(ierr); 831 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 832 } 833 PetscFunctionReturn(0); 834 } 835