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