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