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,ddm; 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) { 103 ierr = DMCreateDomainDecompositionDM(dm,"default",&ddm);CHKERRQ(ierr); 104 if (!ddm) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM has no default decomposition defined. Set subsolves manually with SNESNASMSetSubdomains()."); 105 ierr = SNESSetDM(snes,ddm);CHKERRQ(ierr); 106 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 107 ierr = DMCreateDomainDecomposition(dm,&nasm->n,NULL,NULL,NULL,&subdms);CHKERRQ(ierr); 108 } 109 if (!subdms) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM has no default decomposition defined. Set subsolves manually with SNESNASMSetSubdomains()."); 110 ierr = DMCreateDomainDecompositionScatters(dm,nasm->n,subdms,&nasm->iscatter,&nasm->oscatter,&nasm->gscatter);CHKERRQ(ierr); 111 112 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 113 ierr = PetscMalloc(nasm->n*sizeof(SNES),&nasm->subsnes);CHKERRQ(ierr); 114 115 for (i=0; i<nasm->n; i++) { 116 ierr = SNESCreate(PETSC_COMM_SELF,&nasm->subsnes[i]);CHKERRQ(ierr); 117 ierr = SNESAppendOptionsPrefix(nasm->subsnes[i],optionsprefix);CHKERRQ(ierr); 118 ierr = SNESAppendOptionsPrefix(nasm->subsnes[i],"sub_");CHKERRQ(ierr); 119 ierr = SNESSetDM(nasm->subsnes[i],subdms[i]);CHKERRQ(ierr); 120 ierr = SNESSetFromOptions(nasm->subsnes[i]);CHKERRQ(ierr); 121 ierr = DMDestroy(&subdms[i]);CHKERRQ(ierr); 122 } 123 ierr = PetscFree(subdms);CHKERRQ(ierr); 124 } else SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot construct local problems automatically without a DM!"); 125 } else SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set subproblems manually if there is no DM!"); 126 /* allocate the global vectors */ 127 ierr = PetscMalloc(nasm->n*sizeof(Vec),&nasm->x);CHKERRQ(ierr); 128 ierr = PetscMalloc(nasm->n*sizeof(Vec),&nasm->xl);CHKERRQ(ierr); 129 ierr = PetscMalloc(nasm->n*sizeof(Vec),&nasm->y);CHKERRQ(ierr); 130 ierr = PetscMalloc(nasm->n*sizeof(Vec),&nasm->b);CHKERRQ(ierr); 131 132 for (i=0; i<nasm->n; i++) { 133 DM subdm; 134 ierr = SNESGetFunction(nasm->subsnes[i],&F,NULL,NULL);CHKERRQ(ierr); 135 ierr = VecDuplicate(F,&nasm->x[i]);CHKERRQ(ierr); 136 ierr = VecDuplicate(F,&nasm->y[i]);CHKERRQ(ierr); 137 ierr = VecDuplicate(F,&nasm->b[i]);CHKERRQ(ierr); 138 ierr = SNESGetDM(nasm->subsnes[i],&subdm);CHKERRQ(ierr); 139 ierr = DMCreateLocalVector(subdm,&nasm->xl[i]);CHKERRQ(ierr); 140 ierr = DMGlobalToLocalHookAdd(subdm,DMGlobalToLocalSubDomainDirichletHook_Private,NULL,nasm->xl[i]);CHKERRQ(ierr); 141 } 142 PetscFunctionReturn(0); 143 } 144 145 #undef __FUNCT__ 146 #define __FUNCT__ "SNESSetFromOptions_NASM" 147 PetscErrorCode SNESSetFromOptions_NASM(SNES snes) 148 { 149 PetscErrorCode ierr; 150 DM dm,ddm; 151 char ddm_name[1024]; 152 PCASMType asmtype; 153 PetscBool flg,monflg; 154 SNES_NASM *nasm = (SNES_NASM*)snes->data; 155 156 PetscFunctionBegin; 157 ierr = PetscOptionsHead("Nonlinear Additive Schwartz options");CHKERRQ(ierr); 158 ierr = PetscOptionsEnum("-snes_nasm_type","Type of restriction/extension","",SNESNASMTypes,(PetscEnum)nasm->type,(PetscEnum*)&asmtype,&flg);CHKERRQ(ierr); 159 if (flg) nasm->type = asmtype; 160 flg = PETSC_FALSE; 161 monflg = PETSC_TRUE; 162 ierr = PetscOptionsBool("-snes_nasm_log","Log times for subSNES solves and restriction","",monflg,&monflg,&flg);CHKERRQ(ierr); 163 if (flg) { 164 ierr = PetscLogEventRegister("SNESNASMSubSolve",((PetscObject)snes)->classid,&nasm->eventsubsolve);CHKERRQ(ierr); 165 ierr = PetscLogEventRegister("SNESNASMRestrict",((PetscObject)snes)->classid,&nasm->eventrestrictinterp);CHKERRQ(ierr); 166 } 167 ierr = PetscOptionsString("-snes_nasm_decomposition", "Name of the DM defining the composition", "SNESSetDM", ddm_name, ddm_name,1024,&flg);CHKERRQ(ierr); 168 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 169 if (flg) { 170 if (dm) { 171 ierr = DMCreateDomainDecompositionDM(dm, ddm_name, &ddm);CHKERRQ(ierr); 172 if (!ddm) SETERRQ1(((PetscObject)snes)->comm, PETSC_ERR_ARG_WRONGSTATE, "Unknown DM decomposition name %s", ddm_name); 173 ierr = PetscInfo(snes,"Using domain decomposition DM defined using options database\n");CHKERRQ(ierr); 174 ierr = SNESSetDM(snes,ddm);CHKERRQ(ierr); 175 } else SETERRQ(((PetscObject)snes)->comm, PETSC_ERR_ARG_WRONGSTATE, "No DM to decompose"); 176 } 177 ierr = PetscOptionsTail();CHKERRQ(ierr); 178 PetscFunctionReturn(0); 179 } 180 181 #undef __FUNCT__ 182 #define __FUNCT__ "SNESView_NASM" 183 PetscErrorCode SNESView_NASM(SNES snes, PetscViewer viewer) 184 { 185 PetscFunctionBegin; 186 SNES_NASM *nasm = (SNES_NASM*)snes->data; 187 PetscErrorCode ierr; 188 PetscMPIInt rank; 189 PetscInt i,N; 190 PetscBool iascii,isstring; 191 PetscViewer sviewer; 192 MPI_Comm comm = ((PetscObject)snes)->comm; 193 194 PetscFunctionBegin; 195 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 196 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 197 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 198 ierr = MPI_Reduce(&nasm->n,&N,1,MPIU_INT,MPIU_SUM,0,comm);CHKERRQ(ierr); 199 if (iascii) { 200 ierr = PetscViewerASCIIPrintf(viewer, " Nonlinear Additive Schwarz: total subdomain blocks = %D\n",N);CHKERRQ(ierr); 201 ierr = PetscViewerASCIIPrintf(viewer," Nonlinear Additive Schwarz: restriction/interpolation type - %s\n",SNESNASMTypes[nasm->type]);CHKERRQ(ierr); 202 ierr = PetscViewerASCIIPrintf(viewer," Nonlinear Additive Schwarz: subSNES iterations: %D subKSP iterations: %D\n",SNESNASMTypes[nasm->type]);CHKERRQ(ierr); 203 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 204 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] number of local blocks = %D\n",(int)rank,nasm->n);CHKERRQ(ierr); 205 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 206 ierr = PetscViewerASCIIPrintf(viewer," Local SNES objects:\n");CHKERRQ(ierr); 207 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 208 if (!rank) { 209 for (i=0; i<nasm->n; i++) { 210 ierr = PetscViewerGetSingleton(viewer,&sviewer);CHKERRQ(ierr); 211 ierr = SNESView(nasm->subsnes[i],sviewer);CHKERRQ(ierr); 212 ierr = PetscViewerRestoreSingleton(viewer,&sviewer);CHKERRQ(ierr); 213 } 214 } 215 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 216 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 217 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 218 } else if (isstring) { 219 ierr = PetscViewerStringSPrintf(viewer," blocks=%D,type=%s",N,SNESNASMTypes[nasm->type]);CHKERRQ(ierr); 220 ierr = PetscViewerGetSingleton(viewer,&sviewer);CHKERRQ(ierr); 221 if (nasm->subsnes && !rank) {ierr = SNESView(nasm->subsnes[0],sviewer);CHKERRQ(ierr);} 222 ierr = PetscViewerRestoreSingleton(viewer,&sviewer);CHKERRQ(ierr); 223 } 224 PetscFunctionReturn(0); 225 } 226 227 #undef __FUNCT__ 228 #define __FUNCT__ "SNESNASMSetSubdomains" 229 PetscErrorCode SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[]) 230 { 231 PetscErrorCode ierr; 232 PetscErrorCode (*f)(SNES,PetscInt,SNES*,VecScatter*,VecScatter*,VecScatter*); 233 234 PetscFunctionBegin; 235 ierr = PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",(void (**)(void))&f);CHKERRQ(ierr); 236 ierr = (f)(snes,n,subsnes,iscatter,oscatter,gscatter);CHKERRQ(ierr); 237 PetscFunctionReturn(0); 238 } 239 240 EXTERN_C_BEGIN 241 #undef __FUNCT__ 242 #define __FUNCT__ "SNESNASMSetSubdomains_NASM" 243 PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[]) 244 { 245 PetscInt i; 246 PetscErrorCode ierr; 247 SNES_NASM *nasm = (SNES_NASM*)snes->data; 248 249 PetscFunctionBegin; 250 if (snes->setupcalled) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"SNESNASMSetSubdomains() should be called before calling SNESSetUp()."); 251 252 /* tear down the previously set things */ 253 ierr = SNESReset(snes);CHKERRQ(ierr); 254 255 nasm->n = n; 256 if (oscatter) { 257 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)oscatter[i]);CHKERRQ(ierr);} 258 } 259 if (iscatter) { 260 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)iscatter[i]);CHKERRQ(ierr);} 261 } 262 if (gscatter) { 263 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)gscatter[i]);CHKERRQ(ierr);} 264 } 265 if (oscatter) { 266 ierr = PetscMalloc(n*sizeof(IS),&nasm->oscatter);CHKERRQ(ierr); 267 for (i=0; i<n; i++) { 268 nasm->oscatter[i] = oscatter[i]; 269 } 270 } 271 if (iscatter) { 272 ierr = PetscMalloc(n*sizeof(IS),&nasm->iscatter);CHKERRQ(ierr); 273 for (i=0; i<n; i++) { 274 nasm->iscatter[i] = iscatter[i]; 275 } 276 } 277 if (gscatter) { 278 ierr = PetscMalloc(n*sizeof(IS),&nasm->gscatter);CHKERRQ(ierr); 279 for (i=0; i<n; i++) { 280 nasm->gscatter[i] = gscatter[i]; 281 } 282 } 283 284 if (subsnes) { 285 ierr = PetscMalloc(n*sizeof(SNES),&nasm->subsnes);CHKERRQ(ierr); 286 for (i=0; i<n; i++) { 287 nasm->subsnes[i] = subsnes[i]; 288 } 289 } 290 PetscFunctionReturn(0); 291 } 292 EXTERN_C_END 293 294 #undef __FUNCT__ 295 #define __FUNCT__ "SNESNASMSolveLocal_Private" 296 PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X) 297 { 298 SNES_NASM *nasm = (SNES_NASM*)snes->data; 299 SNES subsnes; 300 PetscInt i; 301 PetscErrorCode ierr; 302 Vec Xlloc,Xl,Bl,Yl; 303 VecScatter iscat,oscat,gscat; 304 DM dm,subdm; 305 306 PetscFunctionBegin; 307 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 308 ierr = VecSet(Y,0);CHKERRQ(ierr); 309 310 if (nasm->eventrestrictinterp) {ierr = PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);CHKERRQ(ierr);} 311 for (i=0; i<nasm->n; i++) { 312 /* scatter the solution to the local solution */ 313 Xlloc = nasm->xl[i]; 314 gscat = nasm->gscatter[i]; 315 oscat = nasm->oscatter[i]; 316 ierr = VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 317 if (B) { 318 /* scatter the RHS to the local RHS */ 319 Bl = nasm->b[i]; 320 ierr = VecScatterBegin(oscat,B,Bl,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 321 } else { 322 Bl = NULL; 323 } 324 } 325 if (nasm->eventrestrictinterp) {ierr = PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);CHKERRQ(ierr);} 326 327 for (i=0; i<nasm->n; i++) { 328 Xlloc = nasm->xl[i]; 329 ierr = VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 330 if (B) { 331 ierr = VecScatterEnd(oscat,B,Bl,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 332 } 333 } 334 335 if (nasm->eventsubsolve) {ierr = PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);CHKERRQ(ierr);} 336 for (i=0; i<nasm->n; i++) { 337 Xl = nasm->x[i]; 338 Xlloc = nasm->xl[i]; 339 Yl = nasm->y[i]; 340 subsnes = nasm->subsnes[i]; 341 ierr = SNESGetDM(subsnes,&subdm);CHKERRQ(ierr); 342 iscat = nasm->iscatter[i]; 343 oscat = nasm->oscatter[i]; 344 gscat = nasm->gscatter[i]; 345 ierr = DMSubDomainRestrict(dm,oscat,gscat,subdm);CHKERRQ(ierr); 346 347 ierr = DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);CHKERRQ(ierr); 348 ierr = DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);CHKERRQ(ierr); 349 ierr = VecCopy(Xl,Yl);CHKERRQ(ierr); 350 ierr = SNESSolve(subsnes,Bl,Yl);CHKERRQ(ierr); 351 ierr = VecAXPY(Yl,-1.0,Xl);CHKERRQ(ierr); 352 } 353 if (nasm->eventsubsolve) {ierr = PetscLogEventEnd(nasm->eventsubsolve,snes,0,0,0);CHKERRQ(ierr);} 354 355 ierr = MPI_Barrier(((PetscObject)snes)->comm);CHKERRQ(ierr); 356 357 if (nasm->eventrestrictinterp) {ierr = PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);CHKERRQ(ierr);} 358 for (i=0; i<nasm->n; i++) { 359 Yl = nasm->y[i]; 360 iscat = nasm->iscatter[i]; 361 oscat = nasm->oscatter[i]; 362 if (nasm->type == PC_ASM_BASIC) { 363 ierr = VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 364 } else if (nasm->type == PC_ASM_RESTRICT) { 365 ierr = VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 366 } else SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM"); 367 } 368 369 for (i=0; i<nasm->n; i++) { 370 Yl = nasm->y[i]; 371 iscat = nasm->iscatter[i]; 372 oscat = nasm->oscatter[i]; 373 if (nasm->type == PC_ASM_BASIC) { 374 ierr = VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 375 } else if (nasm->type == PC_ASM_RESTRICT) { 376 ierr = VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 377 } else SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM"); 378 } 379 if (nasm->eventrestrictinterp) {ierr = PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);CHKERRQ(ierr);} 380 381 ierr = MPI_Barrier(((PetscObject)snes)->comm);CHKERRQ(ierr); 382 383 ierr = VecAXPY(X,1.0,Y);CHKERRQ(ierr); 384 PetscFunctionReturn(0); 385 } 386 387 #undef __FUNCT__ 388 #define __FUNCT__ "SNESSolve_NASM" 389 PetscErrorCode SNESSolve_NASM(SNES snes) 390 { 391 Vec F; 392 Vec X; 393 Vec B; 394 Vec Y; 395 PetscInt i; 396 PetscReal fnorm; 397 PetscErrorCode ierr; 398 SNESNormType normtype; 399 400 PetscFunctionBegin; 401 X = snes->vec_sol; 402 Y = snes->vec_sol_update; 403 F = snes->vec_func; 404 B = snes->vec_rhs; 405 406 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 407 snes->iter = 0; 408 snes->norm = 0.; 409 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 410 snes->reason = SNES_CONVERGED_ITERATING; 411 ierr = SNESGetNormType(snes, &normtype);CHKERRQ(ierr); 412 if (normtype == SNES_NORM_FUNCTION || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) { 413 /* compute the initial function and preconditioned update delX */ 414 if (!snes->vec_func_init_set) { 415 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 416 if (snes->domainerror) { 417 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 418 PetscFunctionReturn(0); 419 } 420 } else snes->vec_func_init_set = PETSC_FALSE; 421 422 /* convergence test */ 423 if (!snes->norm_init_set) { 424 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 425 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 426 } else { 427 fnorm = snes->norm_init; 428 snes->norm_init_set = PETSC_FALSE; 429 } 430 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 431 snes->iter = 0; 432 snes->norm = fnorm; 433 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 434 ierr = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr); 435 ierr = SNESMonitor(snes,0,snes->norm);CHKERRQ(ierr); 436 437 /* set parameter for default relative tolerance convergence test */ 438 snes->ttol = fnorm*snes->rtol; 439 440 /* test convergence */ 441 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 442 if (snes->reason) PetscFunctionReturn(0); 443 } else { 444 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 445 ierr = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr); 446 ierr = SNESMonitor(snes,0,snes->norm);CHKERRQ(ierr); 447 } 448 449 /* Call general purpose update function */ 450 if (snes->ops->update) { 451 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 452 } 453 454 for (i = 0; i < snes->max_its; i++) { 455 ierr = SNESNASMSolveLocal_Private(snes,B,Y,X);CHKERRQ(ierr); 456 if (normtype == SNES_NORM_FUNCTION || ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY))) { 457 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 458 if (snes->domainerror) { 459 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 460 PetscFunctionReturn(0); 461 } 462 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 463 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 464 } 465 /* Monitor convergence */ 466 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 467 snes->iter = i+1; 468 snes->norm = fnorm; 469 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 470 ierr = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr); 471 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 472 /* Test for convergence */ 473 if (normtype == SNES_NORM_FUNCTION) ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 474 if (snes->reason) PetscFunctionReturn(0); 475 /* Call general purpose update function */ 476 if (snes->ops->update) { 477 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 478 } 479 } 480 if (normtype == SNES_NORM_FUNCTION) { 481 if (i == snes->max_its) { 482 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);CHKERRQ(ierr); 483 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 484 } 485 } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */ 486 PetscFunctionReturn(0); 487 } 488 489 /*MC 490 SNESNASM - Nonlinear Additive Schwartz 491 492 Level: advanced 493 494 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 495 M*/ 496 497 EXTERN_C_BEGIN 498 #undef __FUNCT__ 499 #define __FUNCT__ "SNESCreate_NASM" 500 PetscErrorCode SNESCreate_NASM(SNES snes) 501 { 502 SNES_NASM *nasm; 503 PetscErrorCode ierr; 504 505 PetscFunctionBegin; 506 ierr = PetscNewLog(snes, SNES_NASM, &nasm);CHKERRQ(ierr); 507 snes->data = (void*)nasm; 508 509 nasm->n = PETSC_DECIDE; 510 nasm->subsnes = 0; 511 nasm->x = 0; 512 nasm->xl = 0; 513 nasm->y = 0; 514 nasm->b = 0; 515 nasm->oscatter = 0; 516 nasm->iscatter = 0; 517 nasm->gscatter = 0; 518 519 nasm->type = PC_ASM_BASIC; 520 521 snes->ops->destroy = SNESDestroy_NASM; 522 snes->ops->setup = SNESSetUp_NASM; 523 snes->ops->setfromoptions = SNESSetFromOptions_NASM; 524 snes->ops->view = SNESView_NASM; 525 snes->ops->solve = SNESSolve_NASM; 526 snes->ops->reset = SNESReset_NASM; 527 528 snes->usesksp = PETSC_FALSE; 529 snes->usespc = PETSC_FALSE; 530 531 nasm->eventrestrictinterp = 0; 532 nasm->eventsubsolve = 0; 533 534 if (!snes->tolerancesset) { 535 snes->max_its = 10000; 536 snes->max_funcs = 10000; 537 } 538 539 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESNASMSetSubdomains_C","SNESNASMSetSubdomains_NASM", 540 SNESNASMSetSubdomains_NASM);CHKERRQ(ierr); 541 PetscFunctionReturn(0); 542 } 543 EXTERN_C_END 544