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