1 #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/ 2 #include <petscdm.h> 3 4 typedef struct { 5 PetscInt n; /* local subdomains */ 6 SNES *subsnes; /* nonlinear solvers for each subdomain */ 7 Vec *x; /* solution vectors */ 8 Vec *xl; /* solution local vectors */ 9 Vec *y; /* step vectors */ 10 Vec *b; /* rhs vectors */ 11 VecScatter *oscatter; /* scatter from global space to the subdomain global space */ 12 VecScatter *iscatter; /* scatter from global space to the nonoverlapping subdomain space */ 13 VecScatter *gscatter; /* scatter from global space to the subdomain local space */ 14 PCASMType type; /* ASM type */ 15 PetscBool usesdm; /* use the DM for setting up the subproblems */ 16 PetscBool finaljacobian; /* compute the jacobian of the converged solution */ 17 18 /* logging events */ 19 PetscLogEvent eventrestrictinterp; 20 PetscLogEvent eventsubsolve; 21 } SNES_NASM; 22 23 const char *const SNESNASMTypes[] = {"NONE","RESTRICT","INTERPOLATE","BASIC","PCASMType","PC_ASM_",0}; 24 25 #undef __FUNCT__ 26 #define __FUNCT__ "SNESReset_NASM" 27 PetscErrorCode SNESReset_NASM(SNES snes) 28 { 29 SNES_NASM *nasm = (SNES_NASM*)snes->data; 30 PetscErrorCode ierr; 31 PetscInt i; 32 33 PetscFunctionBegin; 34 for (i=0; i<nasm->n; i++) { 35 if (nasm->xl) { ierr = VecDestroy(&nasm->xl[i]);CHKERRQ(ierr); } 36 if (nasm->x) { ierr = VecDestroy(&nasm->x[i]);CHKERRQ(ierr); } 37 if (nasm->y) { ierr = VecDestroy(&nasm->y[i]);CHKERRQ(ierr); } 38 if (nasm->b) { ierr = VecDestroy(&nasm->b[i]);CHKERRQ(ierr); } 39 40 if (nasm->subsnes) { ierr = SNESDestroy(&nasm->subsnes[i]);CHKERRQ(ierr); } 41 if (nasm->oscatter) { ierr = VecScatterDestroy(&nasm->oscatter[i]);CHKERRQ(ierr); } 42 if (nasm->iscatter) { ierr = VecScatterDestroy(&nasm->iscatter[i]);CHKERRQ(ierr); } 43 if (nasm->gscatter) { ierr = VecScatterDestroy(&nasm->gscatter[i]);CHKERRQ(ierr); } 44 } 45 46 if (nasm->x) {ierr = PetscFree(nasm->x);CHKERRQ(ierr);} 47 if (nasm->xl) {ierr = PetscFree(nasm->xl);CHKERRQ(ierr);} 48 if (nasm->y) {ierr = PetscFree(nasm->y);CHKERRQ(ierr);} 49 if (nasm->b) {ierr = PetscFree(nasm->b);CHKERRQ(ierr);} 50 51 if (nasm->subsnes) {ierr = PetscFree(nasm->subsnes);CHKERRQ(ierr);} 52 if (nasm->oscatter) {ierr = PetscFree(nasm->oscatter);CHKERRQ(ierr);} 53 if (nasm->iscatter) {ierr = PetscFree(nasm->iscatter);CHKERRQ(ierr);} 54 if (nasm->gscatter) {ierr = PetscFree(nasm->gscatter);CHKERRQ(ierr);} 55 56 nasm->eventrestrictinterp = 0; 57 nasm->eventsubsolve = 0; 58 PetscFunctionReturn(0); 59 } 60 61 #undef __FUNCT__ 62 #define __FUNCT__ "SNESDestroy_NASM" 63 PetscErrorCode SNESDestroy_NASM(SNES snes) 64 { 65 PetscErrorCode ierr; 66 67 PetscFunctionBegin; 68 ierr = SNESReset_NASM(snes);CHKERRQ(ierr); 69 ierr = PetscFree(snes->data);CHKERRQ(ierr); 70 PetscFunctionReturn(0); 71 } 72 73 #undef __FUNCT__ 74 #define __FUNCT__ "DMGlobalToLocalSubDomainDirichletHook_Private" 75 PetscErrorCode DMGlobalToLocalSubDomainDirichletHook_Private(DM dm,Vec g,InsertMode mode,Vec l,void *ctx) 76 { 77 PetscErrorCode ierr; 78 Vec bcs = (Vec)ctx; 79 80 PetscFunctionBegin; 81 ierr = VecCopy(bcs,l);CHKERRQ(ierr); 82 PetscFunctionReturn(0); 83 } 84 85 #undef __FUNCT__ 86 #define __FUNCT__ "SNESSetUp_NASM" 87 PetscErrorCode SNESSetUp_NASM(SNES snes) 88 { 89 SNES_NASM *nasm = (SNES_NASM*)snes->data; 90 PetscErrorCode ierr; 91 DM dm,subdm; 92 DM *subdms; 93 PetscInt i; 94 const char *optionsprefix; 95 Vec F; 96 97 PetscFunctionBegin; 98 if (!nasm->subsnes) { 99 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 100 if (dm) { 101 nasm->usesdm = PETSC_TRUE; 102 ierr = DMCreateDomainDecomposition(dm,&nasm->n,NULL,NULL,NULL,&subdms);CHKERRQ(ierr); 103 if (!subdms) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM has no default decomposition defined. Set subsolves manually with SNESNASMSetSubdomains()."); 104 ierr = DMCreateDomainDecompositionScatters(dm,nasm->n,subdms,&nasm->iscatter,&nasm->oscatter,&nasm->gscatter);CHKERRQ(ierr); 105 106 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 107 ierr = PetscMalloc(nasm->n*sizeof(SNES),&nasm->subsnes);CHKERRQ(ierr); 108 109 for (i=0; i<nasm->n; i++) { 110 ierr = SNESCreate(PETSC_COMM_SELF,&nasm->subsnes[i]);CHKERRQ(ierr); 111 ierr = SNESAppendOptionsPrefix(nasm->subsnes[i],optionsprefix);CHKERRQ(ierr); 112 ierr = SNESAppendOptionsPrefix(nasm->subsnes[i],"sub_");CHKERRQ(ierr); 113 ierr = SNESSetDM(nasm->subsnes[i],subdms[i]);CHKERRQ(ierr); 114 ierr = SNESSetFromOptions(nasm->subsnes[i]);CHKERRQ(ierr); 115 ierr = DMDestroy(&subdms[i]);CHKERRQ(ierr); 116 } 117 ierr = PetscFree(subdms);CHKERRQ(ierr); 118 } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Cannot construct local problems automatically without a DM!"); 119 } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set subproblems manually if there is no DM!"); 120 /* allocate the global vectors */ 121 if (!nasm->x) { 122 ierr = PetscMalloc(nasm->n*sizeof(Vec),&nasm->x);CHKERRQ(ierr); 123 ierr = PetscMemzero(nasm->x,nasm->n*sizeof(Vec));CHKERRQ(ierr); 124 } 125 if (!nasm->xl) { 126 ierr = PetscMalloc(nasm->n*sizeof(Vec),&nasm->xl);CHKERRQ(ierr); 127 ierr = PetscMemzero(nasm->xl,nasm->n*sizeof(Vec));CHKERRQ(ierr); 128 } 129 if (!nasm->y) { 130 ierr = PetscMalloc(nasm->n*sizeof(Vec),&nasm->y);CHKERRQ(ierr); 131 ierr = PetscMemzero(nasm->y,nasm->n*sizeof(Vec));CHKERRQ(ierr); 132 } 133 if (!nasm->b) { 134 ierr = PetscMalloc(nasm->n*sizeof(Vec),&nasm->b);CHKERRQ(ierr); 135 ierr = PetscMemzero(nasm->b,nasm->n*sizeof(Vec));CHKERRQ(ierr); 136 } 137 138 for (i=0; i<nasm->n; i++) { 139 ierr = SNESGetFunction(nasm->subsnes[i],&F,NULL,NULL);CHKERRQ(ierr); 140 if (!nasm->x[i]) {ierr = VecDuplicate(F,&nasm->x[i]);CHKERRQ(ierr);} 141 if (!nasm->y[i]) {ierr = VecDuplicate(F,&nasm->y[i]);CHKERRQ(ierr);} 142 if (!nasm->b[i]) {ierr = VecDuplicate(F,&nasm->b[i]);CHKERRQ(ierr);} 143 if (!nasm->xl[i]) { 144 ierr = SNESGetDM(nasm->subsnes[i],&subdm);CHKERRQ(ierr); 145 ierr = DMCreateLocalVector(subdm,&nasm->xl[i]);CHKERRQ(ierr); 146 } 147 ierr = DMGlobalToLocalHookAdd(subdm,DMGlobalToLocalSubDomainDirichletHook_Private,NULL,nasm->xl[i]);CHKERRQ(ierr); 148 } 149 if (nasm->finaljacobian) {ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr);} 150 PetscFunctionReturn(0); 151 } 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "SNESSetFromOptions_NASM" 155 PetscErrorCode SNESSetFromOptions_NASM(SNES snes) 156 { 157 PetscErrorCode ierr; 158 PCASMType asmtype; 159 PetscBool flg,monflg; 160 SNES_NASM *nasm = (SNES_NASM*)snes->data; 161 162 PetscFunctionBegin; 163 ierr = PetscOptionsHead("Nonlinear Additive Schwartz options");CHKERRQ(ierr); 164 ierr = PetscOptionsEnum("-snes_nasm_type","Type of restriction/extension","",SNESNASMTypes,(PetscEnum)nasm->type,(PetscEnum*)&asmtype,&flg);CHKERRQ(ierr); 165 if (flg) nasm->type = asmtype; 166 flg = PETSC_FALSE; 167 monflg = PETSC_TRUE; 168 ierr = PetscOptionsBool("-snes_nasm_log","Log times for subSNES solves and restriction","",monflg,&monflg,&flg);CHKERRQ(ierr); 169 ierr = PetscOptionsBool("-snes_nasm_finaljacobian","Compute the global jacobian of the final iterate (for ASPIN)","",nasm->finaljacobian,&nasm->finaljacobian,NULL);CHKERRQ(ierr); 170 if (flg) { 171 ierr = PetscLogEventRegister("SNESNASMSubSolve",((PetscObject)snes)->classid,&nasm->eventsubsolve);CHKERRQ(ierr); 172 ierr = PetscLogEventRegister("SNESNASMRestrict",((PetscObject)snes)->classid,&nasm->eventrestrictinterp);CHKERRQ(ierr); 173 } 174 ierr = PetscOptionsTail();CHKERRQ(ierr); 175 PetscFunctionReturn(0); 176 } 177 178 #undef __FUNCT__ 179 #define __FUNCT__ "SNESView_NASM" 180 PetscErrorCode SNESView_NASM(SNES snes, PetscViewer viewer) 181 { 182 PetscFunctionBegin; 183 SNES_NASM *nasm = (SNES_NASM*)snes->data; 184 PetscErrorCode ierr; 185 PetscMPIInt rank; 186 PetscInt i,N; 187 PetscBool iascii,isstring; 188 PetscViewer sviewer; 189 MPI_Comm comm; 190 191 PetscFunctionBegin; 192 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 193 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 194 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 195 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 196 ierr = MPI_Reduce(&nasm->n,&N,1,MPIU_INT,MPIU_SUM,0,comm);CHKERRQ(ierr); 197 if (iascii) { 198 ierr = PetscViewerASCIIPrintf(viewer, " Nonlinear Additive Schwarz: total subdomain blocks = %D\n",N);CHKERRQ(ierr); 199 ierr = PetscViewerASCIIPrintf(viewer," Nonlinear Additive Schwarz: restriction/interpolation type - %s\n",SNESNASMTypes[nasm->type]);CHKERRQ(ierr); 200 ierr = PetscViewerASCIIPrintf(viewer," Nonlinear Additive Schwarz: subSNES iterations: %D subKSP iterations: %D\n",SNESNASMTypes[nasm->type]);CHKERRQ(ierr); 201 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 202 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] number of local blocks = %D\n",(int)rank,nasm->n);CHKERRQ(ierr); 203 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 204 ierr = PetscViewerASCIIPrintf(viewer," Local SNES objects:\n");CHKERRQ(ierr); 205 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 206 if (!rank) { 207 for (i=0; i<nasm->n; i++) { 208 ierr = PetscViewerGetSingleton(viewer,&sviewer);CHKERRQ(ierr); 209 ierr = SNESView(nasm->subsnes[i],sviewer);CHKERRQ(ierr); 210 ierr = PetscViewerRestoreSingleton(viewer,&sviewer);CHKERRQ(ierr); 211 } 212 } 213 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 214 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 215 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 216 } else if (isstring) { 217 ierr = PetscViewerStringSPrintf(viewer," blocks=%D,type=%s",N,SNESNASMTypes[nasm->type]);CHKERRQ(ierr); 218 ierr = PetscViewerGetSingleton(viewer,&sviewer);CHKERRQ(ierr); 219 if (nasm->subsnes && !rank) {ierr = SNESView(nasm->subsnes[0],sviewer);CHKERRQ(ierr);} 220 ierr = PetscViewerRestoreSingleton(viewer,&sviewer);CHKERRQ(ierr); 221 } 222 PetscFunctionReturn(0); 223 } 224 225 #undef __FUNCT__ 226 #define __FUNCT__ "SNESNASMSetSubdomains" 227 /*@ 228 SNESNASMSetSubdomains - Manually Set the context required to restrict and solve subdomain problems. 229 230 Not Collective 231 232 Input Parameters: 233 + SNES - the SNES context 234 . n - the number of local subdomains 235 . subsnes - solvers defined on the local subdomains 236 . iscatter - scatters into the nonoverlapping portions of the local subdomains 237 . oscatter - scatters into the overlapping portions of the local subdomains 238 - gscatter - scatters into the (ghosted) local vector of the local subdomain 239 240 Level: intermediate 241 242 .keywords: SNES, NASM 243 244 .seealso: SNESNASM, SNESNASMGetSubdomains() 245 @*/ 246 PetscErrorCode SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[]) 247 { 248 PetscErrorCode ierr; 249 PetscErrorCode (*f)(SNES,PetscInt,SNES*,VecScatter*,VecScatter*,VecScatter*); 250 251 PetscFunctionBegin; 252 ierr = PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",(void (**)(void))&f);CHKERRQ(ierr); 253 ierr = (f)(snes,n,subsnes,iscatter,oscatter,gscatter);CHKERRQ(ierr); 254 PetscFunctionReturn(0); 255 } 256 257 EXTERN_C_BEGIN 258 #undef __FUNCT__ 259 #define __FUNCT__ "SNESNASMSetSubdomains_NASM" 260 PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[]) 261 { 262 PetscInt i; 263 PetscErrorCode ierr; 264 SNES_NASM *nasm = (SNES_NASM*)snes->data; 265 266 PetscFunctionBegin; 267 if (snes->setupcalled) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"SNESNASMSetSubdomains() should be called before calling SNESSetUp()."); 268 269 /* tear down the previously set things */ 270 ierr = SNESReset(snes);CHKERRQ(ierr); 271 272 nasm->n = n; 273 if (oscatter) { 274 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)oscatter[i]);CHKERRQ(ierr);} 275 } 276 if (iscatter) { 277 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)iscatter[i]);CHKERRQ(ierr);} 278 } 279 if (gscatter) { 280 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)gscatter[i]);CHKERRQ(ierr);} 281 } 282 if (oscatter) { 283 ierr = PetscMalloc(n*sizeof(IS),&nasm->oscatter);CHKERRQ(ierr); 284 for (i=0; i<n; i++) { 285 nasm->oscatter[i] = oscatter[i]; 286 } 287 } 288 if (iscatter) { 289 ierr = PetscMalloc(n*sizeof(IS),&nasm->iscatter);CHKERRQ(ierr); 290 for (i=0; i<n; i++) { 291 nasm->iscatter[i] = iscatter[i]; 292 } 293 } 294 if (gscatter) { 295 ierr = PetscMalloc(n*sizeof(IS),&nasm->gscatter);CHKERRQ(ierr); 296 for (i=0; i<n; i++) { 297 nasm->gscatter[i] = gscatter[i]; 298 } 299 } 300 301 if (subsnes) { 302 ierr = PetscMalloc(n*sizeof(SNES),&nasm->subsnes);CHKERRQ(ierr); 303 for (i=0; i<n; i++) { 304 nasm->subsnes[i] = subsnes[i]; 305 } 306 } 307 PetscFunctionReturn(0); 308 } 309 EXTERN_C_END 310 311 #undef __FUNCT__ 312 #define __FUNCT__ "SNESNASMGetSubdomains" 313 /*@ 314 SNESNASMGetSubdomains - Get the local subdomain context. 315 316 Not Collective 317 318 Input Parameters: 319 . SNES - the SNES context 320 321 Output Parameters: 322 + n - the number of local subdomains 323 . subsnes - solvers defined on the local subdomains 324 . iscatter - scatters into the nonoverlapping portions of the local subdomains 325 . oscatter - scatters into the overlapping portions of the local subdomains 326 - gscatter - scatters into the (ghosted) local vector of the local subdomain 327 328 Level: intermediate 329 330 .keywords: SNES, NASM 331 332 .seealso: SNESNASM, SNESNASMSetSubdomains() 333 @*/ 334 PetscErrorCode SNESNASMGetSubdomains(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[]) 335 { 336 PetscErrorCode ierr; 337 PetscErrorCode (*f)(SNES,PetscInt*,SNES**,VecScatter**,VecScatter**,VecScatter**); 338 339 PetscFunctionBegin; 340 ierr = PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",(void (**)(void))&f);CHKERRQ(ierr); 341 ierr = (f)(snes,n,subsnes,iscatter,oscatter,gscatter);CHKERRQ(ierr); 342 PetscFunctionReturn(0); 343 } 344 345 EXTERN_C_BEGIN 346 #undef __FUNCT__ 347 #define __FUNCT__ "SNESNASMGetSubdomains_NASM" 348 PetscErrorCode SNESNASMGetSubdomains_NASM(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[]) 349 { 350 SNES_NASM *nasm = (SNES_NASM*)snes->data; 351 352 PetscFunctionBegin; 353 if (n) *n = nasm->n; 354 if (oscatter) *oscatter = nasm->oscatter; 355 if (iscatter) *iscatter = nasm->iscatter; 356 if (gscatter) *gscatter = nasm->gscatter; 357 if (subsnes) *subsnes = nasm->subsnes; 358 PetscFunctionReturn(0); 359 } 360 EXTERN_C_END 361 362 #undef __FUNCT__ 363 #define __FUNCT__ "SNESNASMGetSubdomainVecs" 364 /*@ 365 SNESNASMGetSubdomainVecs - Get the processor-local subdomain vectors 366 367 Not Collective 368 369 Input Parameters: 370 . SNES - the SNES context 371 372 Output Parameters: 373 + n - the number of local subdomains 374 . x - The subdomain solution vector 375 . y - The subdomain step vector 376 . b - The subdomain RHS vector 377 - xl - The subdomain local vectors (ghosted) 378 379 Level: developer 380 381 .keywords: SNES, NASM 382 383 .seealso: SNESNASM, SNESNASMGetSubdomains() 384 @*/ 385 PetscErrorCode SNESNASMGetSubdomainVecs(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b, Vec **xl) 386 { 387 PetscErrorCode ierr; 388 PetscErrorCode (*f)(SNES,PetscInt*,Vec**,Vec**,Vec**,Vec**); 389 390 PetscFunctionBegin; 391 ierr = PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",(void (**)(void))&f);CHKERRQ(ierr); 392 ierr = (f)(snes,n,x,y,b,xl);CHKERRQ(ierr); 393 PetscFunctionReturn(0); 394 } 395 396 EXTERN_C_BEGIN 397 #undef __FUNCT__ 398 #define __FUNCT__ "SNESNASMGetSubdomainVecs_NASM" 399 PetscErrorCode SNESNASMGetSubdomainVecs_NASM(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b,Vec **xl) 400 { 401 SNES_NASM *nasm = (SNES_NASM*)snes->data; 402 403 PetscFunctionBegin; 404 if (n) *n = nasm->n; 405 if (x) *x = nasm->x; 406 if (y) *y = nasm->y; 407 if (b) *b = nasm->b; 408 if (xl) *xl = nasm->xl; 409 PetscFunctionReturn(0); 410 } 411 EXTERN_C_END 412 413 #undef __FUNCT__ 414 #define __FUNCT__ "SNESNASMSetComputeFinalJacobian" 415 /*@ 416 SNESNASMSetComputeFinalJacobian - Schedules the computation of the global and subdomain jacobians upon convergence 417 418 Collective on SNES 419 420 Input Parameters: 421 + SNES - the SNES context 422 - flg - indication of whether to compute the jacobians or not 423 424 Level: developer 425 426 Notes: This is used almost exclusively in the implementation of ASPIN, where the converged subdomain and global jacobian 427 is needed at each linear iteration. 428 429 .keywords: SNES, NASM, ASPIN 430 431 .seealso: SNESNASM, SNESNASMGetSubdomains() 432 @*/ 433 PetscErrorCode SNESNASMSetComputeFinalJacobian(SNES snes,PetscBool flg) 434 { 435 PetscErrorCode (*f)(SNES,PetscBool); 436 PetscErrorCode ierr; 437 438 PetscFunctionBegin; 439 ierr = PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",(void (**)(void))&f);CHKERRQ(ierr); 440 ierr = (f)(snes,flg);CHKERRQ(ierr); 441 PetscFunctionReturn(0); 442 } 443 444 EXTERN_C_BEGIN 445 #undef __FUNCT__ 446 #define __FUNCT__ "SNESNASMSetComputeFinalJacobian_NASM" 447 PetscErrorCode SNESNASMSetComputeFinalJacobian_NASM(SNES snes,PetscBool flg) 448 { 449 SNES_NASM *nasm = (SNES_NASM*)snes->data; 450 451 PetscFunctionBegin; 452 nasm->finaljacobian = flg; 453 PetscFunctionReturn(0); 454 } 455 EXTERN_C_END 456 457 458 #undef __FUNCT__ 459 #define __FUNCT__ "SNESNASMSolveLocal_Private" 460 PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X) 461 { 462 SNES_NASM *nasm = (SNES_NASM*)snes->data; 463 SNES subsnes; 464 PetscInt i; 465 PetscErrorCode ierr; 466 Vec Xlloc,Xl,Bl,Yl; 467 VecScatter iscat,oscat,gscat; 468 DM dm,subdm; 469 470 PetscFunctionBegin; 471 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 472 ierr = VecSet(Y,0);CHKERRQ(ierr); 473 474 if (nasm->eventrestrictinterp) {ierr = PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);CHKERRQ(ierr);} 475 for (i=0; i<nasm->n; i++) { 476 /* scatter the solution to the local solution */ 477 Xlloc = nasm->xl[i]; 478 gscat = nasm->gscatter[i]; 479 oscat = nasm->oscatter[i]; 480 ierr = VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 481 if (B) { 482 /* scatter the RHS to the local RHS */ 483 Bl = nasm->b[i]; 484 ierr = VecScatterBegin(oscat,B,Bl,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 485 } 486 } 487 if (nasm->eventrestrictinterp) {ierr = PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);CHKERRQ(ierr);} 488 489 for (i=0; i<nasm->n; i++) { 490 Xlloc = nasm->xl[i]; 491 gscat = nasm->gscatter[i]; 492 oscat = nasm->oscatter[i]; 493 ierr = VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 494 if (B) { 495 Bl = nasm->b[i]; 496 ierr = VecScatterEnd(oscat,B,Bl,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 497 } 498 } 499 500 if (nasm->eventsubsolve) {ierr = PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);CHKERRQ(ierr);} 501 for (i=0; i<nasm->n; i++) { 502 Xl = nasm->x[i]; 503 Xlloc = nasm->xl[i]; 504 Yl = nasm->y[i]; 505 subsnes = nasm->subsnes[i]; 506 ierr = SNESGetDM(subsnes,&subdm);CHKERRQ(ierr); 507 iscat = nasm->iscatter[i]; 508 oscat = nasm->oscatter[i]; 509 gscat = nasm->gscatter[i]; 510 ierr = DMSubDomainRestrict(dm,oscat,gscat,subdm);CHKERRQ(ierr); 511 if (B) { 512 Bl = nasm->b[i]; 513 } else { 514 Bl = NULL; 515 } 516 ierr = DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);CHKERRQ(ierr); 517 ierr = DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);CHKERRQ(ierr); 518 ierr = VecCopy(Xl,Yl);CHKERRQ(ierr); 519 ierr = SNESSolve(subsnes,Bl,Yl);CHKERRQ(ierr); 520 ierr = VecAXPY(Yl,-1.0,Xl);CHKERRQ(ierr); 521 } 522 if (nasm->eventsubsolve) {ierr = PetscLogEventEnd(nasm->eventsubsolve,snes,0,0,0);CHKERRQ(ierr);} 523 524 ierr = MPI_Barrier(PetscObjectComm((PetscObject)snes));CHKERRQ(ierr); 525 526 if (nasm->eventrestrictinterp) {ierr = PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);CHKERRQ(ierr);} 527 for (i=0; i<nasm->n; i++) { 528 Yl = nasm->y[i]; 529 iscat = nasm->iscatter[i]; 530 oscat = nasm->oscatter[i]; 531 if (nasm->type == PC_ASM_BASIC) { 532 ierr = VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 533 } else if (nasm->type == PC_ASM_RESTRICT) { 534 ierr = VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 535 } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM"); 536 } 537 538 for (i=0; i<nasm->n; i++) { 539 Yl = nasm->y[i]; 540 iscat = nasm->iscatter[i]; 541 oscat = nasm->oscatter[i]; 542 if (nasm->type == PC_ASM_BASIC) { 543 ierr = VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 544 } else if (nasm->type == PC_ASM_RESTRICT) { 545 ierr = VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 546 } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM"); 547 } 548 if (nasm->eventrestrictinterp) {ierr = PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);CHKERRQ(ierr);} 549 550 ierr = MPI_Barrier(PetscObjectComm((PetscObject)snes));CHKERRQ(ierr); 551 552 ierr = VecAXPY(X,1.0,Y);CHKERRQ(ierr); 553 PetscFunctionReturn(0); 554 } 555 556 #undef __FUNCT__ 557 #define __FUNCT__ "SNESNASMComputeFinalJacobian_Private" 558 PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec X) 559 { 560 SNES_NASM *nasm = (SNES_NASM*)snes->data; 561 SNES subsnes; 562 PetscInt i; 563 PetscErrorCode ierr; 564 Vec Xlloc,Xl; 565 VecScatter oscat,gscat; 566 DM dm,subdm; 567 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 568 569 PetscFunctionBegin; 570 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 571 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 572 if (nasm->eventrestrictinterp) {ierr = PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);CHKERRQ(ierr);} 573 for (i=0; i<nasm->n; i++) { 574 /* scatter the solution to the local solution */ 575 Xlloc = nasm->xl[i]; 576 gscat = nasm->gscatter[i]; 577 oscat = nasm->oscatter[i]; 578 ierr = VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 579 } 580 if (nasm->eventrestrictinterp) {ierr = PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);CHKERRQ(ierr);} 581 582 for (i=0; i<nasm->n; i++) { 583 Xlloc = nasm->xl[i]; 584 gscat = nasm->gscatter[i]; 585 oscat = nasm->oscatter[i]; 586 ierr = VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 587 } 588 589 for (i=0; i<nasm->n; i++) { 590 Xl = nasm->x[i]; 591 Xlloc = nasm->xl[i]; 592 subsnes = nasm->subsnes[i]; 593 ierr = SNESGetDM(subsnes,&subdm);CHKERRQ(ierr); 594 oscat = nasm->oscatter[i]; 595 gscat = nasm->gscatter[i]; 596 ierr = DMSubDomainRestrict(dm,oscat,gscat,subdm);CHKERRQ(ierr); 597 ierr = DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);CHKERRQ(ierr); 598 ierr = DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);CHKERRQ(ierr); 599 ierr = SNESComputeJacobian(subsnes,Xl,&snes->jacobian,&subsnes->jacobian_pre,&flg);CHKERRQ(ierr); 600 ierr = KSPSetOperators(subsnes->ksp,subsnes->jacobian,subsnes->jacobian_pre,flg);CHKERRQ(ierr); 601 } 602 603 ierr = MPI_Barrier(PetscObjectComm((PetscObject)snes));CHKERRQ(ierr); 604 PetscFunctionReturn(0); 605 } 606 607 #undef __FUNCT__ 608 #define __FUNCT__ "SNESSolve_NASM" 609 PetscErrorCode SNESSolve_NASM(SNES snes) 610 { 611 Vec F; 612 Vec X; 613 Vec B; 614 Vec Y; 615 PetscInt i; 616 PetscReal fnorm; 617 PetscErrorCode ierr; 618 SNESNormType normtype; 619 SNES_NASM *nasm = (SNES_NASM*)snes->data; 620 621 PetscFunctionBegin; 622 X = snes->vec_sol; 623 Y = snes->vec_sol_update; 624 F = snes->vec_func; 625 B = snes->vec_rhs; 626 627 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 628 snes->iter = 0; 629 snes->norm = 0.; 630 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 631 snes->reason = SNES_CONVERGED_ITERATING; 632 ierr = SNESGetNormType(snes, &normtype);CHKERRQ(ierr); 633 if (normtype == SNES_NORM_FUNCTION || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) { 634 /* compute the initial function and preconditioned update delX */ 635 if (!snes->vec_func_init_set) { 636 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 637 if (snes->domainerror) { 638 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 639 PetscFunctionReturn(0); 640 } 641 } else snes->vec_func_init_set = PETSC_FALSE; 642 643 /* convergence test */ 644 if (!snes->norm_init_set) { 645 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 646 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 647 } else { 648 fnorm = snes->norm_init; 649 snes->norm_init_set = PETSC_FALSE; 650 } 651 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 652 snes->iter = 0; 653 snes->norm = fnorm; 654 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 655 ierr = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr); 656 ierr = SNESMonitor(snes,0,snes->norm);CHKERRQ(ierr); 657 658 /* set parameter for default relative tolerance convergence test */ 659 snes->ttol = fnorm*snes->rtol; 660 661 /* test convergence */ 662 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 663 if (snes->reason) PetscFunctionReturn(0); 664 } else { 665 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 666 ierr = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr); 667 ierr = SNESMonitor(snes,0,snes->norm);CHKERRQ(ierr); 668 } 669 670 /* Call general purpose update function */ 671 if (snes->ops->update) { 672 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 673 } 674 675 for (i = 0; i < snes->max_its; i++) { 676 ierr = SNESNASMSolveLocal_Private(snes,B,Y,X);CHKERRQ(ierr); 677 if (normtype == SNES_NORM_FUNCTION || ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY))) { 678 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 679 if (snes->domainerror) { 680 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 681 break; 682 } 683 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 684 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 685 } 686 /* Monitor convergence */ 687 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 688 snes->iter = i+1; 689 snes->norm = fnorm; 690 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 691 ierr = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr); 692 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 693 /* Test for convergence */ 694 if (normtype == SNES_NORM_FUNCTION) ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 695 if (snes->reason) break; 696 /* Call general purpose update function */ 697 if (snes->ops->update) { 698 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 699 } 700 } 701 if (nasm->finaljacobian) {ierr = SNESNASMComputeFinalJacobian_Private(snes,X);CHKERRQ(ierr);} 702 if (normtype == SNES_NORM_FUNCTION) { 703 if (i == snes->max_its) { 704 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);CHKERRQ(ierr); 705 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 706 } 707 } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */ 708 PetscFunctionReturn(0); 709 } 710 711 /*MC 712 SNESNASM - Nonlinear Additive Schwartz 713 714 Options Database: 715 + -snes_nasm_log - enable logging events for the communication and solve stages 716 . -snes_nasm_type <basic,restrict> - type of subdomain update used 717 . -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate 718 . -sub_snes_ - options prefix of the subdomain nonlinear solves 719 . -sub_ksp_ - options prefix of the subdomain Krylov solver 720 - -sub_pc_ - options prefix of the subdomain preconditioner 721 722 Level: advanced 723 724 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 725 M*/ 726 727 EXTERN_C_BEGIN 728 #undef __FUNCT__ 729 #define __FUNCT__ "SNESCreate_NASM" 730 PetscErrorCode SNESCreate_NASM(SNES snes) 731 { 732 SNES_NASM *nasm; 733 PetscErrorCode ierr; 734 735 PetscFunctionBegin; 736 ierr = PetscNewLog(snes, SNES_NASM, &nasm);CHKERRQ(ierr); 737 snes->data = (void*)nasm; 738 739 nasm->n = PETSC_DECIDE; 740 nasm->subsnes = 0; 741 nasm->x = 0; 742 nasm->xl = 0; 743 nasm->y = 0; 744 nasm->b = 0; 745 nasm->oscatter = 0; 746 nasm->iscatter = 0; 747 nasm->gscatter = 0; 748 749 nasm->type = PC_ASM_BASIC; 750 nasm->finaljacobian = PETSC_FALSE; 751 752 snes->ops->destroy = SNESDestroy_NASM; 753 snes->ops->setup = SNESSetUp_NASM; 754 snes->ops->setfromoptions = SNESSetFromOptions_NASM; 755 snes->ops->view = SNESView_NASM; 756 snes->ops->solve = SNESSolve_NASM; 757 snes->ops->reset = SNESReset_NASM; 758 759 snes->usesksp = PETSC_FALSE; 760 snes->usespc = PETSC_FALSE; 761 762 nasm->eventrestrictinterp = 0; 763 nasm->eventsubsolve = 0; 764 765 if (!snes->tolerancesset) { 766 snes->max_its = 10000; 767 snes->max_funcs = 10000; 768 } 769 770 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESNASMSetSubdomains_C","SNESNASMSetSubdomains_NASM", 771 SNESNASMSetSubdomains_NASM);CHKERRQ(ierr); 772 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESNASMGetSubdomains_C","SNESNASMGetSubdomains_NASM", 773 SNESNASMGetSubdomains_NASM);CHKERRQ(ierr); 774 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESNASMGetSubdomainVecs_C","SNESNASMGetSubdomainVecs_NASM", 775 SNESNASMGetSubdomainVecs_NASM);CHKERRQ(ierr); 776 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C","SNESNASMSetComputeFinalJacobian_NASM", 777 SNESNASMSetComputeFinalJacobian_NASM);CHKERRQ(ierr); 778 PetscFunctionReturn(0); 779 } 780 EXTERN_C_END 781