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