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