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 8 Vec *x; /* solution vectors */ 9 Vec *xl; /* solution local vectors */ 10 Vec *y; /* step vectors */ 11 Vec *b; /* rhs vectors */ 12 13 VecScatter *oscatter; /* scatter from global space to the subdomain global space */ 14 VecScatter *iscatter; /* scatter from global space to the nonoverlapping subdomain space */ 15 VecScatter *gscatter; /* scatter from global space to the subdomain local space */ 16 17 PCASMType type; /* ASM type */ 18 19 PetscBool usesdm; /* use the DM for setting up the subproblems */ 20 } SNES_NASM; 21 22 #undef __FUNCT__ 23 #define __FUNCT__ "SNESReset_NASM" 24 PetscErrorCode SNESReset_NASM(SNES snes) 25 { 26 SNES_NASM *nasm = (SNES_NASM *)snes->data; 27 PetscErrorCode ierr; 28 PetscInt i; 29 PetscFunctionBegin; 30 for (i=0;i<nasm->n;i++) { 31 if (nasm->x) { ierr = VecDestroy(&nasm->x[i]);CHKERRQ(ierr); } 32 if (nasm->xl){ ierr = VecDestroy(&nasm->xl[i]);CHKERRQ(ierr); } 33 if (nasm->y) { ierr = VecDestroy(&nasm->y[i]);CHKERRQ(ierr); } 34 if (nasm->b) { ierr = VecDestroy(&nasm->b[i]);CHKERRQ(ierr); } 35 36 if (nasm->subsnes) { ierr = SNESDestroy(&nasm->subsnes[i]);CHKERRQ(ierr); } 37 if (nasm->oscatter) { ierr = VecScatterDestroy(&nasm->oscatter[i]);CHKERRQ(ierr); } 38 if (nasm->iscatter) { ierr = VecScatterDestroy(&nasm->iscatter[i]);CHKERRQ(ierr); } 39 if (nasm->gscatter) { ierr = VecScatterDestroy(&nasm->gscatter[i]);CHKERRQ(ierr); } 40 } 41 42 if (nasm->x) {ierr = PetscFree(nasm->x);CHKERRQ(ierr);} 43 if (nasm->xl) {ierr = PetscFree(nasm->xl);CHKERRQ(ierr);} 44 if (nasm->y) {ierr = PetscFree(nasm->y);CHKERRQ(ierr);} 45 if (nasm->b) {ierr = PetscFree(nasm->b);CHKERRQ(ierr);} 46 47 if (nasm->subsnes) {ierr = PetscFree(nasm->subsnes);CHKERRQ(ierr);} 48 if (nasm->oscatter) {ierr = PetscFree(nasm->oscatter);CHKERRQ(ierr);} 49 if (nasm->iscatter) {ierr = PetscFree(nasm->iscatter);CHKERRQ(ierr);} 50 if (nasm->gscatter) {ierr = PetscFree(nasm->gscatter);CHKERRQ(ierr);} 51 52 PetscFunctionReturn(0); 53 } 54 55 #undef __FUNCT__ 56 #define __FUNCT__ "SNESDestroy_NASM" 57 PetscErrorCode SNESDestroy_NASM(SNES snes) 58 { 59 PetscErrorCode ierr; 60 PetscFunctionBegin; 61 ierr = SNESReset_NASM(snes);CHKERRQ(ierr); 62 ierr = PetscFree(snes->data);CHKERRQ(ierr); 63 PetscFunctionReturn(0); 64 } 65 66 #undef __FUNCT__ 67 #define __FUNCT__ "DMGlobalToLocalSubDomainDirichletHook_Private" 68 PetscErrorCode DMGlobalToLocalSubDomainDirichletHook_Private(DM dm,Vec g,InsertMode mode,Vec l,void *ctx) { 69 PetscErrorCode ierr; 70 Vec bcs = (Vec)ctx; 71 PetscFunctionBegin; 72 ierr = VecCopy(bcs,l);CHKERRQ(ierr); 73 PetscFunctionReturn(0); 74 } 75 76 #undef __FUNCT__ 77 #define __FUNCT__ "SNESSetUp_NASM" 78 PetscErrorCode SNESSetUp_NASM(SNES snes) 79 { 80 SNES_NASM *nasm = (SNES_NASM *)snes->data; 81 PetscErrorCode ierr; 82 DM dm; 83 DM *subdms; 84 PetscInt i; 85 const char *optionsprefix; 86 Vec F; 87 88 PetscFunctionBegin; 89 90 if (!nasm->subsnes) { 91 if (snes->dm) { 92 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 93 nasm->usesdm = PETSC_TRUE; 94 /* create the subdomain dms on this proc */ 95 ierr = DMCreateDomainDecomposition(dm,&nasm->n,PETSC_NULL,PETSC_NULL,PETSC_NULL,&subdms);CHKERRQ(ierr); 96 ierr = DMCreateDomainDecompositionScatters(dm,nasm->n,subdms,&nasm->iscatter,&nasm->oscatter,&nasm->gscatter);CHKERRQ(ierr); 97 98 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 99 ierr = PetscMalloc(nasm->n*sizeof(SNES),&nasm->subsnes);CHKERRQ(ierr); 100 101 for (i=0;i<nasm->n;i++) { 102 ierr = SNESCreate(PETSC_COMM_SELF,&nasm->subsnes[i]);CHKERRQ(ierr); 103 ierr = SNESAppendOptionsPrefix(nasm->subsnes[i],optionsprefix);CHKERRQ(ierr); 104 ierr = SNESAppendOptionsPrefix(nasm->subsnes[i],"sub_");CHKERRQ(ierr); 105 ierr = SNESSetDM(nasm->subsnes[i],subdms[i]);CHKERRQ(ierr); 106 ierr = SNESSetFromOptions(nasm->subsnes[i]);CHKERRQ(ierr); 107 ierr = DMDestroy(&subdms[i]);CHKERRQ(ierr); 108 } 109 ierr = PetscFree(subdms);CHKERRQ(ierr); 110 } else { 111 SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot construct local problems automatically without a DM!"); 112 } 113 } else { 114 SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set subproblems manually if there is no DM!"); 115 } 116 /* allocate the global vectors */ 117 ierr = PetscMalloc(nasm->n*sizeof(Vec),&nasm->x);CHKERRQ(ierr); 118 ierr = PetscMalloc(nasm->n*sizeof(Vec),&nasm->xl);CHKERRQ(ierr); 119 ierr = PetscMalloc(nasm->n*sizeof(Vec),&nasm->y);CHKERRQ(ierr); 120 ierr = PetscMalloc(nasm->n*sizeof(Vec),&nasm->b);CHKERRQ(ierr); 121 122 for (i=0;i<nasm->n;i++) { 123 DM subdm; 124 ierr = SNESGetFunction(nasm->subsnes[i],&F,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 125 ierr = VecDuplicate(F,&nasm->x[i]);CHKERRQ(ierr); 126 ierr = VecDuplicate(F,&nasm->y[i]);CHKERRQ(ierr); 127 ierr = VecDuplicate(F,&nasm->b[i]);CHKERRQ(ierr); 128 ierr = SNESGetDM(nasm->subsnes[i],&subdm);CHKERRQ(ierr); 129 ierr = DMCreateLocalVector(subdm,&nasm->xl[i]);CHKERRQ(ierr); 130 ierr = DMGlobalToLocalHookAdd(subdm,DMGlobalToLocalSubDomainDirichletHook_Private,PETSC_NULL,nasm->xl[i]);CHKERRQ(ierr); 131 } 132 133 PetscFunctionReturn(0); 134 } 135 136 #undef __FUNCT__ 137 #define __FUNCT__ "SNESSetFromOptions_NASM" 138 PetscErrorCode SNESSetFromOptions_NASM(SNES snes) 139 { 140 PetscErrorCode ierr; 141 PCASMType asmtype; 142 const char *const SNESNASMTypes[] = {"NONE","RESTRICT","INTERPOLATE","BASIC","PCASMType","PC_ASM_",0}; 143 PetscBool flg; 144 SNES_NASM *nasm = (SNES_NASM *)snes->data; 145 PetscFunctionBegin; 146 ierr = PetscOptionsHead("Nonlinear Additive Schwartz options");CHKERRQ(ierr); 147 ierr = PetscOptionsEnum("-snes_nasm_type","Type of restriction/extension","",SNESNASMTypes,(PetscEnum)nasm->type,(PetscEnum*)&asmtype,&flg);CHKERRQ(ierr); 148 if (flg) {nasm->type = asmtype;} 149 ierr = PetscOptionsTail();CHKERRQ(ierr); 150 PetscFunctionReturn(0); 151 } 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "SNESView_NASM" 155 PetscErrorCode SNESView_NASM(SNES snes, PetscViewer viewer) 156 { 157 PetscFunctionBegin; 158 PetscFunctionReturn(0); 159 } 160 161 #undef __FUNCT__ 162 #define __FUNCT__ "SNESNASMSetSubdomains" 163 PetscErrorCode SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[]) { 164 PetscErrorCode ierr; 165 PetscErrorCode (*f)(SNES,PetscInt,SNES*,VecScatter*,VecScatter*,VecScatter*); 166 PetscFunctionBegin; 167 ierr = PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",(void (**)(void))&f);CHKERRQ(ierr); 168 ierr = (f)(snes,n,subsnes,iscatter,oscatter,gscatter);CHKERRQ(ierr); 169 PetscFunctionReturn(0); 170 } 171 172 EXTERN_C_BEGIN 173 #undef __FUNCT__ 174 #define __FUNCT__ "SNESNASMSetSubdomains_NASM" 175 PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[]) { 176 PetscInt i; 177 PetscErrorCode ierr; 178 SNES_NASM *nasm = (SNES_NASM *)snes->data; 179 PetscFunctionBegin; 180 if (snes->setupcalled)SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"SNESNASMSetSubdomains() should be called before calling SNESSetUp()."); 181 182 /* tear down the previously set things */ 183 ierr = SNESReset(snes);CHKERRQ(ierr); 184 185 nasm->n = n; 186 if (oscatter) { 187 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)oscatter[i]);CHKERRQ(ierr);} 188 } 189 if (iscatter) { 190 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)iscatter[i]);CHKERRQ(ierr);} 191 } 192 if (gscatter) { 193 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)gscatter[i]);CHKERRQ(ierr);} 194 } 195 if (oscatter) { 196 ierr = PetscMalloc(n*sizeof(IS),&nasm->oscatter);CHKERRQ(ierr); 197 for (i=0; i<n; i++) { 198 nasm->oscatter[i] = oscatter[i]; 199 } 200 } 201 if (iscatter) { 202 ierr = PetscMalloc(n*sizeof(IS),&nasm->iscatter);CHKERRQ(ierr); 203 for (i=0; i<n; i++) { 204 nasm->iscatter[i] = iscatter[i]; 205 } 206 } 207 if (gscatter) { 208 ierr = PetscMalloc(n*sizeof(IS),&nasm->gscatter);CHKERRQ(ierr); 209 for (i=0; i<n; i++) { 210 nasm->gscatter[i] = gscatter[i]; 211 } 212 } 213 214 if (subsnes) { 215 ierr = PetscMalloc(n*sizeof(SNES),&nasm->subsnes);CHKERRQ(ierr); 216 for (i=0; i<n; i++) { 217 nasm->subsnes[i] = subsnes[i]; 218 } 219 } 220 PetscFunctionReturn(0); 221 } 222 EXTERN_C_END 223 224 #undef __FUNCT__ 225 #define __FUNCT__ "SNESNASMSolveLocal_Private" 226 PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X) { 227 SNES_NASM *nasm = (SNES_NASM *)snes->data; 228 SNES subsnes; 229 PetscInt i; 230 PetscErrorCode ierr; 231 Vec Xlloc,Xl,Bl,Yl; 232 VecScatter iscat,oscat,gscat; 233 DM dm,subdm; 234 PetscFunctionBegin; 235 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 236 ierr = VecSet(Y,0);CHKERRQ(ierr); 237 for(i=0;i<nasm->n;i++) { 238 subsnes = nasm->subsnes[i]; 239 ierr = SNESGetDM(subsnes,&subdm);CHKERRQ(ierr); 240 iscat = nasm->iscatter[i]; 241 oscat = nasm->oscatter[i]; 242 gscat = nasm->gscatter[i]; 243 ierr = DMSubDomainRestrict(dm,oscat,gscat,subdm);CHKERRQ(ierr); 244 245 /* scatter the solution to the local solution */ 246 Xl = nasm->x[i]; 247 Xlloc = nasm->xl[i]; 248 Yl = nasm->y[i]; 249 250 ierr = VecScatterBegin(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 251 ierr = VecScatterEnd(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 252 253 ierr = VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 254 ierr = VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 255 256 ierr = VecCopy(Xl,Yl);CHKERRQ(ierr); 257 258 if (B) { 259 /* scatter the RHS to the local RHS */ 260 Bl = nasm->b[i]; 261 ierr = VecScatterBegin(oscat,B,Bl,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 262 ierr = VecScatterEnd(oscat,B,Bl,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 263 } else { 264 Bl = PETSC_NULL; 265 } 266 ierr = SNESSolve(subsnes,Bl,Yl);CHKERRQ(ierr); 267 268 ierr = VecAXPY(Yl,-1.0,Xl);CHKERRQ(ierr); 269 270 if (nasm->type == PC_ASM_BASIC) { 271 ierr = VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 272 ierr = VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 273 } else if (nasm->type == PC_ASM_RESTRICT) { 274 ierr = VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 275 ierr = VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 276 } else { 277 SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Only basic and interpolate types are supported for SNESNASM"); 278 } 279 } 280 281 ierr = VecAssemblyBegin(Y);CHKERRQ(ierr); 282 ierr = VecAssemblyEnd(Y);CHKERRQ(ierr); 283 284 ierr = VecAXPY(X,1.0,Y);CHKERRQ(ierr); 285 286 PetscFunctionReturn(0); 287 } 288 289 #undef __FUNCT__ 290 #define __FUNCT__ "SNESSolve_NASM" 291 PetscErrorCode SNESSolve_NASM(SNES snes) 292 { 293 Vec F; 294 Vec X; 295 Vec B; 296 Vec Y; 297 PetscInt i; 298 PetscReal fnorm; 299 PetscErrorCode ierr; 300 SNESNormType normtype; 301 302 PetscFunctionBegin; 303 304 X = snes->vec_sol; 305 Y = snes->vec_sol_update; 306 F = snes->vec_func; 307 B = snes->vec_rhs; 308 309 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 310 snes->iter = 0; 311 snes->norm = 0.; 312 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 313 snes->reason = SNES_CONVERGED_ITERATING; 314 ierr = SNESGetNormType(snes, &normtype);CHKERRQ(ierr); 315 if (normtype == SNES_NORM_FUNCTION || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) { 316 /* compute the initial function and preconditioned update delX */ 317 if (!snes->vec_func_init_set) { 318 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 319 if (snes->domainerror) { 320 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 321 PetscFunctionReturn(0); 322 } 323 } else { 324 snes->vec_func_init_set = PETSC_FALSE; 325 } 326 327 /* convergence test */ 328 if (!snes->norm_init_set) { 329 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 330 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 331 } else { 332 fnorm = snes->norm_init; 333 snes->norm_init_set = PETSC_FALSE; 334 } 335 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 336 snes->iter = 0; 337 snes->norm = fnorm; 338 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 339 SNESLogConvHistory(snes,snes->norm,0); 340 ierr = SNESMonitor(snes,0,snes->norm);CHKERRQ(ierr); 341 342 /* set parameter for default relative tolerance convergence test */ 343 snes->ttol = fnorm*snes->rtol; 344 345 /* test convergence */ 346 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 347 if (snes->reason) PetscFunctionReturn(0); 348 } else { 349 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 350 SNESLogConvHistory(snes,snes->norm,0); 351 ierr = SNESMonitor(snes,0,snes->norm);CHKERRQ(ierr); 352 } 353 354 /* Call general purpose update function */ 355 if (snes->ops->update) { 356 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 357 } 358 359 for (i = 0; i < snes->max_its; i++) { 360 ierr = SNESNASMSolveLocal_Private(snes,B,Y,X);CHKERRQ(ierr); 361 if (normtype == SNES_NORM_FUNCTION || ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY))) { 362 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 363 if (snes->domainerror) { 364 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 365 PetscFunctionReturn(0); 366 } 367 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 368 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 369 } 370 /* Monitor convergence */ 371 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 372 snes->iter = i+1; 373 snes->norm = fnorm; 374 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 375 SNESLogConvHistory(snes,snes->norm,0); 376 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 377 /* Test for convergence */ 378 if (normtype == SNES_NORM_FUNCTION) ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 379 if (snes->reason) PetscFunctionReturn(0); 380 /* Call general purpose update function */ 381 if (snes->ops->update) { 382 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 383 } 384 } 385 if (normtype == SNES_NORM_FUNCTION) { 386 if (i == snes->max_its) { 387 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);CHKERRQ(ierr); 388 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 389 } 390 } else { 391 if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */ 392 } 393 PetscFunctionReturn(0); 394 } 395 396 /*MC 397 SNESNASM - Nonlinear Additive Schwartz 398 399 Level: advanced 400 401 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 402 M*/ 403 404 EXTERN_C_BEGIN 405 #undef __FUNCT__ 406 #define __FUNCT__ "SNESCreate_NASM" 407 PetscErrorCode SNESCreate_NASM(SNES snes) 408 { 409 SNES_NASM *nasm; 410 PetscErrorCode ierr; 411 412 PetscFunctionBegin; 413 414 415 ierr = PetscNewLog(snes, SNES_NASM, &nasm);CHKERRQ(ierr); 416 snes->data = (void*)nasm; 417 418 nasm->n = PETSC_DECIDE; 419 nasm->subsnes = 0; 420 nasm->x = 0; 421 nasm->xl = 0; 422 nasm->y = 0; 423 nasm->b = 0; 424 nasm->oscatter = 0; 425 nasm->iscatter = 0; 426 nasm->gscatter = 0; 427 428 nasm->type = PC_ASM_BASIC; 429 430 snes->ops->destroy = SNESDestroy_NASM; 431 snes->ops->setup = SNESSetUp_NASM; 432 snes->ops->setfromoptions = SNESSetFromOptions_NASM; 433 snes->ops->view = SNESView_NASM; 434 snes->ops->solve = SNESSolve_NASM; 435 snes->ops->reset = SNESReset_NASM; 436 437 snes->usesksp = PETSC_FALSE; 438 snes->usespc = PETSC_FALSE; 439 440 if (!snes->tolerancesset) { 441 snes->max_its = 10000; 442 snes->max_funcs = 10000; 443 } 444 445 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESNASMSetSubdomains_C","SNESNASMSetSubdomains_NASM", 446 SNESNASMSetSubdomains_NASM);CHKERRQ(ierr); 447 448 PetscFunctionReturn(0); 449 } 450 EXTERN_C_END 451