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