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