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