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