Actual source code: nasm.c

  1: #include <petsc/private/snesimpl.h>
  2: #include <petscdm.h>

  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:   Vec        weight;              /* weighting for adding updates on overlaps, in global space */
 12:   VecScatter *oscatter;           /* scatter from global space to the subdomain global space */
 13:   VecScatter *oscatter_copy;      /* copy of the above */
 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:   PCASMType  type;                /* ASM type */
 17:   PetscBool  usesdm;              /* use the DM for setting up the subproblems */
 18:   PetscBool  finaljacobian;       /* compute the jacobian of the converged solution */
 19:   PetscReal  damping;             /* damping parameter for updates from the blocks */
 20:   PetscBool  weight_set;          /* use a weight in the overlap updates */

 22:   /* logging events */
 23:   PetscLogEvent eventrestrictinterp;
 24:   PetscLogEvent eventsubsolve;

 26:   PetscInt      fjtype;            /* type of computed jacobian */
 27:   Vec           xinit;             /* initial solution in case the final jacobian type is computed as first */
 28: } SNES_NASM;

 30: const char *const SNESNASMTypes[] = {"NONE","RESTRICT","INTERPOLATE","BASIC","PCASMType","PC_ASM_",NULL};
 31: const char *const SNESNASMFJTypes[] = {"FINALOUTER","FINALINNER","INITIAL"};

 33: static PetscErrorCode SNESReset_NASM(SNES snes)
 34: {
 35:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
 37:   PetscInt       i;

 40:   for (i=0; i<nasm->n; i++) {
 41:     if (nasm->xl) { VecDestroy(&nasm->xl[i]); }
 42:     if (nasm->x) { VecDestroy(&nasm->x[i]); }
 43:     if (nasm->y) { VecDestroy(&nasm->y[i]); }
 44:     if (nasm->b) { VecDestroy(&nasm->b[i]); }

 46:     if (nasm->subsnes) { SNESDestroy(&nasm->subsnes[i]); }
 47:     if (nasm->oscatter) { VecScatterDestroy(&nasm->oscatter[i]); }
 48:     if (nasm->oscatter_copy) { VecScatterDestroy(&nasm->oscatter_copy[i]); }
 49:     if (nasm->iscatter) { VecScatterDestroy(&nasm->iscatter[i]); }
 50:     if (nasm->gscatter) { VecScatterDestroy(&nasm->gscatter[i]); }
 51:   }

 53:   PetscFree(nasm->x);
 54:   PetscFree(nasm->xl);
 55:   PetscFree(nasm->y);
 56:   PetscFree(nasm->b);

 58:   if (nasm->xinit) {VecDestroy(&nasm->xinit);}

 60:   PetscFree(nasm->subsnes);
 61:   PetscFree(nasm->oscatter);
 62:   PetscFree(nasm->oscatter_copy);
 63:   PetscFree(nasm->iscatter);
 64:   PetscFree(nasm->gscatter);

 66:   if (nasm->weight_set) {
 67:     VecDestroy(&nasm->weight);
 68:   }

 70:   nasm->eventrestrictinterp = 0;
 71:   nasm->eventsubsolve = 0;
 72:   return(0);
 73: }

 75: static PetscErrorCode SNESDestroy_NASM(SNES snes)
 76: {

 80:   SNESReset_NASM(snes);
 81:   PetscFree(snes->data);
 82:   return(0);
 83: }

 85: static PetscErrorCode DMGlobalToLocalSubDomainDirichletHook_Private(DM dm,Vec g,InsertMode mode,Vec l,void *ctx)
 86: {
 88:   Vec            bcs = (Vec)ctx;

 91:   VecCopy(bcs,l);
 92:   return(0);
 93: }

 95: static PetscErrorCode SNESSetUp_NASM(SNES snes)
 96: {
 97:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
 99:   DM             dm,subdm;
100:   DM             *subdms;
101:   PetscInt       i;
102:   const char     *optionsprefix;
103:   Vec            F;
104:   PetscMPIInt    size;
105:   KSP            ksp;
106:   PC             pc;

109:   if (!nasm->subsnes) {
110:     SNESGetDM(snes,&dm);
111:     if (dm) {
112:       nasm->usesdm = PETSC_TRUE;
113:       DMCreateDomainDecomposition(dm,&nasm->n,NULL,NULL,NULL,&subdms);
114:       if (!subdms) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM has no default decomposition defined.  Set subsolves manually with SNESNASMSetSubdomains().");
115:       DMCreateDomainDecompositionScatters(dm,nasm->n,subdms,&nasm->iscatter,&nasm->oscatter,&nasm->gscatter);
116:       PetscMalloc1(nasm->n, &nasm->oscatter_copy);
117:       for (i=0; i<nasm->n; i++) {
118:         VecScatterCopy(nasm->oscatter[i], &nasm->oscatter_copy[i]);
119:       }

121:       SNESGetOptionsPrefix(snes, &optionsprefix);
122:       PetscMalloc1(nasm->n,&nasm->subsnes);
123:       for (i=0; i<nasm->n; i++) {
124:         SNESCreate(PETSC_COMM_SELF,&nasm->subsnes[i]);
125:         PetscObjectIncrementTabLevel((PetscObject)nasm->subsnes[i], (PetscObject)snes, 1);
126:         SNESAppendOptionsPrefix(nasm->subsnes[i],optionsprefix);
127:         SNESAppendOptionsPrefix(nasm->subsnes[i],"sub_");
128:         SNESSetDM(nasm->subsnes[i],subdms[i]);
129:         MPI_Comm_size(PetscObjectComm((PetscObject)nasm->subsnes[i]),&size);
130:         if (size == 1) {
131:           SNESGetKSP(nasm->subsnes[i],&ksp);
132:           KSPGetPC(ksp,&pc);
133:           KSPSetType(ksp,KSPPREONLY);
134:           PCSetType(pc,PCLU);
135:         }
136:         SNESSetFromOptions(nasm->subsnes[i]);
137:         DMDestroy(&subdms[i]);
138:       }
139:       PetscFree(subdms);
140:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Cannot construct local problems automatically without a DM!");
141:   } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set subproblems manually if there is no DM!");
142:   /* allocate the global vectors */
143:   if (!nasm->x) {
144:     PetscCalloc1(nasm->n,&nasm->x);
145:   }
146:   if (!nasm->xl) {
147:     PetscCalloc1(nasm->n,&nasm->xl);
148:   }
149:   if (!nasm->y) {
150:     PetscCalloc1(nasm->n,&nasm->y);
151:   }
152:   if (!nasm->b) {
153:     PetscCalloc1(nasm->n,&nasm->b);
154:   }

156:   for (i=0; i<nasm->n; i++) {
157:     SNESGetFunction(nasm->subsnes[i],&F,NULL,NULL);
158:     if (!nasm->x[i]) {VecDuplicate(F,&nasm->x[i]);}
159:     if (!nasm->y[i]) {VecDuplicate(F,&nasm->y[i]);}
160:     if (!nasm->b[i]) {VecDuplicate(F,&nasm->b[i]);}
161:     if (!nasm->xl[i]) {
162:       SNESGetDM(nasm->subsnes[i],&subdm);
163:       DMCreateLocalVector(subdm,&nasm->xl[i]);
164:       DMGlobalToLocalHookAdd(subdm,DMGlobalToLocalSubDomainDirichletHook_Private,NULL,nasm->xl[i]);
165:     }
166:   }
167:   if (nasm->finaljacobian) {
168:     SNESSetUpMatrices(snes);
169:     if (nasm->fjtype == 2) {
170:       VecDuplicate(snes->vec_sol,&nasm->xinit);
171:     }
172:     for (i=0; i<nasm->n;i++) {
173:       SNESSetUpMatrices(nasm->subsnes[i]);
174:     }
175:   }
176:   return(0);
177: }

179: static PetscErrorCode SNESSetFromOptions_NASM(PetscOptionItems *PetscOptionsObject,SNES snes)
180: {
181:   PetscErrorCode    ierr;
182:   PCASMType         asmtype;
183:   PetscBool         flg,monflg;
184:   SNES_NASM         *nasm = (SNES_NASM*)snes->data;

187:   PetscOptionsHead(PetscOptionsObject,"Nonlinear Additive Schwarz options");
188:   PetscOptionsEnum("-snes_nasm_type","Type of restriction/extension","",SNESNASMTypes,(PetscEnum)nasm->type,(PetscEnum*)&asmtype,&flg);
189:   if (flg) {SNESNASMSetType(snes,asmtype);}
190:   flg    = PETSC_FALSE;
191:   monflg = PETSC_TRUE;
192:   PetscOptionsReal("-snes_nasm_damping","The new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)","SNESNASMSetDamping",nasm->damping,&nasm->damping,&flg);
193:   if (flg) {SNESNASMSetDamping(snes,nasm->damping);}
194:   PetscOptionsDeprecated("-snes_nasm_sub_view",NULL,"3.15","Use -snes_view ::ascii_info_detail");
195:   PetscOptionsBool("-snes_nasm_finaljacobian","Compute the global jacobian of the final iterate (for ASPIN)","",nasm->finaljacobian,&nasm->finaljacobian,NULL);
196:   PetscOptionsEList("-snes_nasm_finaljacobian_type","The type of the final jacobian computed.","",SNESNASMFJTypes,3,SNESNASMFJTypes[0],&nasm->fjtype,NULL);
197:   PetscOptionsBool("-snes_nasm_log","Log times for subSNES solves and restriction","",monflg,&monflg,&flg);
198:   if (flg) {
199:     PetscLogEventRegister("SNESNASMSubSolve",((PetscObject)snes)->classid,&nasm->eventsubsolve);
200:     PetscLogEventRegister("SNESNASMRestrict",((PetscObject)snes)->classid,&nasm->eventrestrictinterp);
201:   }
202:   PetscOptionsTail();
203:   return(0);
204: }

206: static PetscErrorCode SNESView_NASM(SNES snes, PetscViewer viewer)
207: {
208:   SNES_NASM         *nasm = (SNES_NASM*)snes->data;
209:   PetscErrorCode    ierr;
210:   PetscMPIInt       rank,size;
211:   PetscInt          i,N,bsz;
212:   PetscBool         iascii,isstring;
213:   PetscViewer       sviewer;
214:   MPI_Comm          comm;
215:   PetscViewerFormat format;
216:   const char        *prefix;

219:   PetscObjectGetComm((PetscObject)snes,&comm);
220:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
221:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
222:   MPI_Comm_rank(comm,&rank);
223:   MPI_Comm_size(comm,&size);
224:   MPIU_Allreduce(&nasm->n,&N,1,MPIU_INT,MPI_SUM,comm);
225:   if (iascii) {
226:     PetscViewerASCIIPrintf(viewer, "  total subdomain blocks = %D\n",N);
227:     PetscViewerGetFormat(viewer,&format);
228:     if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
229:       if (nasm->subsnes) {
230:         PetscViewerASCIIPrintf(viewer,"  Local solver information for first block on rank 0:\n");
231:         SNESGetOptionsPrefix(snes,&prefix);
232:         PetscViewerASCIIPrintf(viewer,"  Use -%ssnes_view ::ascii_info_detail to display information for all blocks\n",prefix?prefix:"");
233:         PetscViewerASCIIPushTab(viewer);
234:         PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
235:         if (!rank) {
236:           PetscViewerASCIIPushTab(viewer);
237:           SNESView(nasm->subsnes[0],sviewer);
238:           PetscViewerASCIIPopTab(viewer);
239:         }
240:         PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
241:         PetscViewerASCIIPopTab(viewer);
242:       }
243:     } else {
244:       /* print the solver on each block */
245:       PetscViewerASCIIPushSynchronized(viewer);
246:       PetscViewerASCIISynchronizedPrintf(viewer,"  [%d] number of local blocks = %D\n",(int)rank,nasm->n);
247:       PetscViewerFlush(viewer);
248:       PetscViewerASCIIPopSynchronized(viewer);
249:       PetscViewerASCIIPrintf(viewer,"  Local solver information for each block is in the following SNES objects:\n");
250:       PetscViewerASCIIPushTab(viewer);
251:       PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
252:       PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
253:       for (i=0; i<nasm->n; i++) {
254:         VecGetLocalSize(nasm->x[i],&bsz);
255:         PetscViewerASCIIPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);
256:         SNESView(nasm->subsnes[i],sviewer);
257:         PetscViewerASCIIPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
258:       }
259:       PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
260:       PetscViewerFlush(viewer);
261:       PetscViewerASCIIPopTab(viewer);
262:     }
263:   } else if (isstring) {
264:     PetscViewerStringSPrintf(viewer," blocks=%D,type=%s",N,SNESNASMTypes[nasm->type]);
265:     PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
266:     if (nasm->subsnes && !rank) {SNESView(nasm->subsnes[0],sviewer);}
267:     PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
268:   }
269:   return(0);
270: }

272: /*@
273:    SNESNASMSetType - Set the type of subdomain update used

275:    Logically Collective on SNES

277:    Input Parameters:
278: +  SNES - the SNES context
279: -  type - the type of update, PC_ASM_BASIC or PC_ASM_RESTRICT

281:    Level: intermediate

283: .seealso: SNESNASM, SNESNASMGetType(), PCASMSetType()
284: @*/
285: PetscErrorCode SNESNASMSetType(SNES snes,PCASMType type)
286: {
288:   PetscErrorCode (*f)(SNES,PCASMType);

291:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetType_C",&f);
292:   if (f) {(f)(snes,type);}
293:   return(0);
294: }

296: static PetscErrorCode SNESNASMSetType_NASM(SNES snes,PCASMType type)
297: {
298:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

301:   if (type != PC_ASM_BASIC && type != PC_ASM_RESTRICT) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"SNESNASM only supports basic and restrict types");
302:   nasm->type = type;
303:   return(0);
304: }

306: /*@
307:    SNESNASMGetType - Get the type of subdomain update used

309:    Logically Collective on SNES

311:    Input Parameters:
312: .  SNES - the SNES context

314:    Output Parameters:
315: .  type - the type of update

317:    Level: intermediate

319: .seealso: SNESNASM, SNESNASMSetType(), PCASMGetType()
320: @*/
321: PetscErrorCode SNESNASMGetType(SNES snes,PCASMType *type)
322: {

326:   PetscUseMethod(snes,"SNESNASMGetType_C",(SNES,PCASMType*),(snes,type));
327:   return(0);
328: }

330: static PetscErrorCode SNESNASMGetType_NASM(SNES snes,PCASMType *type)
331: {
332:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

335:   *type = nasm->type;
336:   return(0);
337: }

339: /*@
340:    SNESNASMSetSubdomains - Manually Set the context required to restrict and solve subdomain problems.

342:    Not Collective

344:    Input Parameters:
345: +  SNES - the SNES context
346: .  n - the number of local subdomains
347: .  subsnes - solvers defined on the local subdomains
348: .  iscatter - scatters into the nonoverlapping portions of the local subdomains
349: .  oscatter - scatters into the overlapping portions of the local subdomains
350: -  gscatter - scatters into the (ghosted) local vector of the local subdomain

352:    Level: intermediate

354: .seealso: SNESNASM, SNESNASMGetSubdomains()
355: @*/
356: PetscErrorCode SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
357: {
359:   PetscErrorCode (*f)(SNES,PetscInt,SNES*,VecScatter*,VecScatter*,VecScatter*);

362:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",&f);
363:   if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
364:   return(0);
365: }

367: static PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
368: {
369:   PetscInt       i;
371:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

374:   if (snes->setupcalled) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"SNESNASMSetSubdomains() should be called before calling SNESSetUp().");

376:   /* tear down the previously set things */
377:   SNESReset(snes);

379:   nasm->n = n;
380:   if (oscatter) {
381:     for (i=0; i<n; i++) {PetscObjectReference((PetscObject)oscatter[i]);}
382:   }
383:   if (iscatter) {
384:     for (i=0; i<n; i++) {PetscObjectReference((PetscObject)iscatter[i]);}
385:   }
386:   if (gscatter) {
387:     for (i=0; i<n; i++) {PetscObjectReference((PetscObject)gscatter[i]);}
388:   }
389:   if (oscatter) {
390:     PetscMalloc1(n,&nasm->oscatter);
391:     PetscMalloc1(n,&nasm->oscatter_copy);
392:     for (i=0; i<n; i++) {
393:       nasm->oscatter[i] = oscatter[i];
394:       VecScatterCopy(oscatter[i], &nasm->oscatter_copy[i]);
395:     }
396:   }
397:   if (iscatter) {
398:     PetscMalloc1(n,&nasm->iscatter);
399:     for (i=0; i<n; i++) {
400:       nasm->iscatter[i] = iscatter[i];
401:     }
402:   }
403:   if (gscatter) {
404:     PetscMalloc1(n,&nasm->gscatter);
405:     for (i=0; i<n; i++) {
406:       nasm->gscatter[i] = gscatter[i];
407:     }
408:   }

410:   if (subsnes) {
411:     PetscMalloc1(n,&nasm->subsnes);
412:     for (i=0; i<n; i++) {
413:       nasm->subsnes[i] = subsnes[i];
414:     }
415:   }
416:   return(0);
417: }

419: /*@
420:    SNESNASMGetSubdomains - Get the local subdomain context.

422:    Not Collective

424:    Input Parameters:
425: .  SNES - the SNES context

427:    Output Parameters:
428: +  n - the number of local subdomains
429: .  subsnes - solvers defined on the local subdomains
430: .  iscatter - scatters into the nonoverlapping portions of the local subdomains
431: .  oscatter - scatters into the overlapping portions of the local subdomains
432: -  gscatter - scatters into the (ghosted) local vector of the local subdomain

434:    Level: intermediate

436: .seealso: SNESNASM, SNESNASMSetSubdomains()
437: @*/
438: PetscErrorCode SNESNASMGetSubdomains(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
439: {
441:   PetscErrorCode (*f)(SNES,PetscInt*,SNES**,VecScatter**,VecScatter**,VecScatter**);

444:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",&f);
445:   if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
446:   return(0);
447: }

449: static PetscErrorCode SNESNASMGetSubdomains_NASM(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
450: {
451:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

454:   if (n) *n = nasm->n;
455:   if (oscatter) *oscatter = nasm->oscatter;
456:   if (iscatter) *iscatter = nasm->iscatter;
457:   if (gscatter) *gscatter = nasm->gscatter;
458:   if (subsnes)  *subsnes  = nasm->subsnes;
459:   return(0);
460: }

462: /*@
463:    SNESNASMGetSubdomainVecs - Get the processor-local subdomain vectors

465:    Not Collective

467:    Input Parameters:
468: .  SNES - the SNES context

470:    Output Parameters:
471: +  n - the number of local subdomains
472: .  x - The subdomain solution vector
473: .  y - The subdomain step vector
474: .  b - The subdomain RHS vector
475: -  xl - The subdomain local vectors (ghosted)

477:    Level: developer

479: .seealso: SNESNASM, SNESNASMGetSubdomains()
480: @*/
481: PetscErrorCode SNESNASMGetSubdomainVecs(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b, Vec **xl)
482: {
484:   PetscErrorCode (*f)(SNES,PetscInt*,Vec**,Vec**,Vec**,Vec**);

487:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",&f);
488:   if (f) {(f)(snes,n,x,y,b,xl);}
489:   return(0);
490: }

492: static PetscErrorCode SNESNASMGetSubdomainVecs_NASM(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b,Vec **xl)
493: {
494:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

497:   if (n)  *n  = nasm->n;
498:   if (x)  *x  = nasm->x;
499:   if (y)  *y  = nasm->y;
500:   if (b)  *b  = nasm->b;
501:   if (xl) *xl = nasm->xl;
502:   return(0);
503: }

505: /*@
506:    SNESNASMSetComputeFinalJacobian - Schedules the computation of the global and subdomain Jacobians upon convergence

508:    Collective on SNES

510:    Input Parameters:
511: +  SNES - the SNES context
512: -  flg - indication of whether to compute the Jacobians or not

514:    Level: developer

516:    Notes:
517:    This is used almost exclusively in the implementation of ASPIN, where the converged subdomain and global Jacobian
518:    is needed at each linear iteration.

520: .seealso: SNESNASM, SNESNASMGetSubdomains()
521: @*/
522: PetscErrorCode SNESNASMSetComputeFinalJacobian(SNES snes,PetscBool flg)
523: {
524:   PetscErrorCode (*f)(SNES,PetscBool);

528:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",&f);
529:   if (f) {(f)(snes,flg);}
530:   return(0);
531: }

533: static PetscErrorCode SNESNASMSetComputeFinalJacobian_NASM(SNES snes,PetscBool flg)
534: {
535:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

538:   nasm->finaljacobian = flg;
539:   return(0);
540: }

542: /*@
543:    SNESNASMSetDamping - Sets the update damping for NASM

545:    Logically collective on SNES

547:    Input Parameters:
548: +  SNES - the SNES context
549: -  dmp - damping

551:    Level: intermediate

553:    Notes:
554:     The new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)

556: .seealso: SNESNASM, SNESNASMGetDamping()
557: @*/
558: PetscErrorCode SNESNASMSetDamping(SNES snes,PetscReal dmp)
559: {
560:   PetscErrorCode (*f)(SNES,PetscReal);

564:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetDamping_C",(void (**)(void))&f);
565:   if (f) {(f)(snes,dmp);}
566:   return(0);
567: }

569: static PetscErrorCode SNESNASMSetDamping_NASM(SNES snes,PetscReal dmp)
570: {
571:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

574:   nasm->damping = dmp;
575:   return(0);
576: }

578: /*@
579:    SNESNASMGetDamping - Gets the update damping for NASM

581:    Not Collective

583:    Input Parameters:
584: +  SNES - the SNES context
585: -  dmp - damping

587:    Level: intermediate

589: .seealso: SNESNASM, SNESNASMSetDamping()
590: @*/
591: PetscErrorCode SNESNASMGetDamping(SNES snes,PetscReal *dmp)
592: {

596:   PetscUseMethod(snes,"SNESNASMGetDamping_C",(SNES,PetscReal*),(snes,dmp));
597:   return(0);
598: }

600: static PetscErrorCode SNESNASMGetDamping_NASM(SNES snes,PetscReal *dmp)
601: {
602:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

605:   *dmp = nasm->damping;
606:   return(0);
607: }

609: /*
610:   Input Parameters:
611: + snes - The solver
612: . B - The RHS vector
613: - X - The initial guess

615:   Output Parameters:
616: . Y - The solution update

618:   TODO: All scatters should be packed into one
619: */
620: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X)
621: {
622:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
623:   SNES           subsnes;
624:   PetscInt       i;
625:   PetscReal      dmp;
627:   Vec            Xl,Bl,Yl,Xlloc;
628:   VecScatter     iscat,oscat,gscat,oscat_copy;
629:   DM             dm,subdm;
630:   PCASMType      type;

633:   SNESNASMGetType(snes,&type);
634:   SNESGetDM(snes,&dm);
635:   VecSet(Y,0);
636:   if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
637:   for (i=0; i<nasm->n; i++) {
638:     /* scatter the solution to the global solution and the local solution */
639:     Xl      = nasm->x[i];
640:     Xlloc   = nasm->xl[i];
641:     oscat   = nasm->oscatter[i];
642:     oscat_copy = nasm->oscatter_copy[i];
643:     gscat   = nasm->gscatter[i];
644:     VecScatterBegin(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
645:     VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
646:     if (B) {
647:       /* scatter the RHS to the local RHS */
648:       Bl   = nasm->b[i];
649:       VecScatterBegin(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
650:     }
651:   }
652:   if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}

654:   if (nasm->eventsubsolve) {PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);}
655:   for (i=0; i<nasm->n; i++) {
656:     Xl    = nasm->x[i];
657:     Xlloc = nasm->xl[i];
658:     Yl    = nasm->y[i];
659:     subsnes = nasm->subsnes[i];
660:     SNESGetDM(subsnes,&subdm);
661:     iscat   = nasm->iscatter[i];
662:     oscat   = nasm->oscatter[i];
663:     oscat_copy = nasm->oscatter_copy[i];
664:     gscat   = nasm->gscatter[i];
665:     VecScatterEnd(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
666:     VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
667:     if (B) {
668:       Bl   = nasm->b[i];
669:       VecScatterEnd(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
670:     } else Bl = NULL;

672:     DMSubDomainRestrict(dm,oscat,gscat,subdm);
673:     VecCopy(Xl,Yl);
674:     SNESSolve(subsnes,Bl,Xl);
675:     VecAYPX(Yl,-1.0,Xl);
676:     VecScale(Yl, nasm->damping);
677:     if (type == PC_ASM_BASIC) {
678:       VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
679:       VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
680:     } else if (type == PC_ASM_RESTRICT) {
681:       VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
682:       VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
683:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
684:   }
685:   if (nasm->eventsubsolve) {PetscLogEventEnd(nasm->eventsubsolve,snes,0,0,0);}
686:   if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
687:   if (nasm->weight_set) {
688:     VecPointwiseMult(Y,Y,nasm->weight);
689:   }
690:   if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
691:   SNESNASMGetDamping(snes,&dmp);
692:   VecAXPY(X,dmp,Y);
693:   return(0);
694: }

696: static PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
697: {
698:   Vec            X = Xfinal;
699:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
700:   SNES           subsnes;
701:   PetscInt       i,lag = 1;
703:   Vec            Xlloc,Xl,Fl,F;
704:   VecScatter     oscat,gscat;
705:   DM             dm,subdm;

708:   if (nasm->fjtype == 2) X = nasm->xinit;
709:   F = snes->vec_func;
710:   if (snes->normschedule == SNES_NORM_NONE) {SNESComputeFunction(snes,X,F);}
711:   SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
712:   SNESGetDM(snes,&dm);
713:   if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
714:   if (nasm->fjtype != 1) {
715:     for (i=0; i<nasm->n; i++) {
716:       Xlloc = nasm->xl[i];
717:       gscat = nasm->gscatter[i];
718:       VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
719:     }
720:   }
721:   if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
722:   for (i=0; i<nasm->n; i++) {
723:     Fl      = nasm->subsnes[i]->vec_func;
724:     Xl      = nasm->x[i];
725:     Xlloc   = nasm->xl[i];
726:     subsnes = nasm->subsnes[i];
727:     oscat   = nasm->oscatter[i];
728:     gscat   = nasm->gscatter[i];
729:     if (nasm->fjtype != 1) {VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);}
730:     SNESGetDM(subsnes,&subdm);
731:     DMSubDomainRestrict(dm,oscat,gscat,subdm);
732:     if (nasm->fjtype != 1) {
733:       DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
734:       DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
735:     }
736:     if (subsnes->lagjacobian == -1)    subsnes->lagjacobian = -2;
737:     else if (subsnes->lagjacobian > 1) lag = subsnes->lagjacobian;
738:     SNESComputeFunction(subsnes,Xl,Fl);
739:     SNESComputeJacobian(subsnes,Xl,subsnes->jacobian,subsnes->jacobian_pre);
740:     if (lag > 1) subsnes->lagjacobian = lag;
741:   }
742:   return(0);
743: }

745: static PetscErrorCode SNESSolve_NASM(SNES snes)
746: {
747:   Vec              F;
748:   Vec              X;
749:   Vec              B;
750:   Vec              Y;
751:   PetscInt         i;
752:   PetscReal        fnorm = 0.0;
753:   PetscErrorCode   ierr;
754:   SNESNormSchedule normschedule;
755:   SNES_NASM        *nasm = (SNES_NASM*)snes->data;


759:   if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);

761:   PetscCitationsRegister(SNESCitation,&SNEScite);
762:   X = snes->vec_sol;
763:   Y = snes->vec_sol_update;
764:   F = snes->vec_func;
765:   B = snes->vec_rhs;

767:   PetscObjectSAWsTakeAccess((PetscObject)snes);
768:   snes->iter   = 0;
769:   snes->norm   = 0.;
770:   PetscObjectSAWsGrantAccess((PetscObject)snes);
771:   snes->reason = SNES_CONVERGED_ITERATING;
772:   SNESGetNormSchedule(snes, &normschedule);
773:   if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
774:     /* compute the initial function and preconditioned update delX */
775:     if (!snes->vec_func_init_set) {
776:       SNESComputeFunction(snes,X,F);
777:     } else snes->vec_func_init_set = PETSC_FALSE;

779:     VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
780:     SNESCheckFunctionNorm(snes,fnorm);
781:     PetscObjectSAWsTakeAccess((PetscObject)snes);
782:     snes->iter = 0;
783:     snes->norm = fnorm;
784:     PetscObjectSAWsGrantAccess((PetscObject)snes);
785:     SNESLogConvergenceHistory(snes,snes->norm,0);
786:     SNESMonitor(snes,0,snes->norm);

788:     /* test convergence */
789:     (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
790:     if (snes->reason) return(0);
791:   } else {
792:     PetscObjectSAWsGrantAccess((PetscObject)snes);
793:     SNESLogConvergenceHistory(snes,snes->norm,0);
794:     SNESMonitor(snes,0,snes->norm);
795:   }

797:   /* Call general purpose update function */
798:   if (snes->ops->update) {
799:     (*snes->ops->update)(snes, snes->iter);
800:   }
801:   /* copy the initial solution over for later */
802:   if (nasm->fjtype == 2) {VecCopy(X,nasm->xinit);}

804:   for (i=0; i < snes->max_its; i++) {
805:     SNESNASMSolveLocal_Private(snes,B,Y,X);
806:     if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
807:       SNESComputeFunction(snes,X,F);
808:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
809:       SNESCheckFunctionNorm(snes,fnorm);
810:     }
811:     /* Monitor convergence */
812:     PetscObjectSAWsTakeAccess((PetscObject)snes);
813:     snes->iter = i+1;
814:     snes->norm = fnorm;
815:     PetscObjectSAWsGrantAccess((PetscObject)snes);
816:     SNESLogConvergenceHistory(snes,snes->norm,0);
817:     SNESMonitor(snes,snes->iter,snes->norm);
818:     /* Test for convergence */
819:     if (normschedule == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);}
820:     if (snes->reason) break;
821:     /* Call general purpose update function */
822:     if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
823:   }
824:   if (nasm->finaljacobian) {
825:     SNESNASMComputeFinalJacobian_Private(snes,X);
826:     SNESCheckJacobianDomainerror(snes);
827:   }
828:   if (normschedule == SNES_NORM_ALWAYS) {
829:     if (i == snes->max_its) {
830:       PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
831:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
832:     }
833:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */
834:   return(0);
835: }

837: /*MC
838:   SNESNASM - Nonlinear Additive Schwarz

840:    Options Database:
841: +  -snes_nasm_log - enable logging events for the communication and solve stages
842: .  -snes_nasm_type <basic,restrict> - type of subdomain update used
843: .  -snes_asm_damping <dmp> - the new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
844: .  -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
845: .  -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> - pick state the jacobian is calculated at
846: .  -sub_snes_ - options prefix of the subdomain nonlinear solves
847: .  -sub_ksp_ - options prefix of the subdomain Krylov solver
848: -  -sub_pc_ - options prefix of the subdomain preconditioner

850:    Level: advanced

852:    Developer Note: This is a non-Newton based nonlinear solver that does not directly require a Jacobian; hence the flag snes->usesksp is set to
853:        false and SNESView() and -snes_view do not display a KSP object. However, if the flag nasm->finaljacobian is set (for example, if
854:        NASM is used as a nonlinear preconditioner for  KSPASPIN) then SNESSetUpMatrices() is called to generate the Jacobian (needed by KSPASPIN)
855:        and this utilizes the KSP for storing the matrices, but the KSP is never used for solving a linear system. Note that when SNESNASM is
856:        used by SNESASPIN they share the same Jacobian matrices because SNESSetUp() (called on the outer SNES KSPASPIN) causes the inner SNES
857:        object (in this case SNESNASM) to inherit the outer Jacobian matrices.

859:    References:
860: .  1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
861:    SIAM Review, 57(4), 2015

863: .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types), SNESNASMSetType(), SNESNASMGetType(), SNESNASMSetSubdomains(), SNESNASMGetSubdomains(), SNESNASMGetSubdomainVecs(), SNESNASMSetComputeFinalJacobian(), SNESNASMSetDamping(), SNESNASMGetDamping()
864: M*/

866: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
867: {
868:   SNES_NASM      *nasm;

872:   PetscNewLog(snes,&nasm);
873:   snes->data = (void*)nasm;

875:   nasm->n        = PETSC_DECIDE;
876:   nasm->subsnes  = NULL;
877:   nasm->x        = NULL;
878:   nasm->xl       = NULL;
879:   nasm->y        = NULL;
880:   nasm->b        = NULL;
881:   nasm->oscatter = NULL;
882:   nasm->oscatter_copy = NULL;
883:   nasm->iscatter = NULL;
884:   nasm->gscatter = NULL;
885:   nasm->damping  = 1.;

887:   nasm->type              = PC_ASM_BASIC;
888:   nasm->finaljacobian     = PETSC_FALSE;
889:   nasm->weight_set        = PETSC_FALSE;

891:   snes->ops->destroy        = SNESDestroy_NASM;
892:   snes->ops->setup          = SNESSetUp_NASM;
893:   snes->ops->setfromoptions = SNESSetFromOptions_NASM;
894:   snes->ops->view           = SNESView_NASM;
895:   snes->ops->solve          = SNESSolve_NASM;
896:   snes->ops->reset          = SNESReset_NASM;

898:   snes->usesksp = PETSC_FALSE;
899:   snes->usesnpc = PETSC_FALSE;

901:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

903:   nasm->fjtype              = 0;
904:   nasm->xinit               = NULL;
905:   nasm->eventrestrictinterp = 0;
906:   nasm->eventsubsolve       = 0;

908:   if (!snes->tolerancesset) {
909:     snes->max_its   = 10000;
910:     snes->max_funcs = 10000;
911:   }

913:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetType_C",SNESNASMSetType_NASM);
914:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetType_C",SNESNASMGetType_NASM);
915:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
916:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
917:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
918:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
919:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
920:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
921:   return(0);
922: }

924: /*@
925:    SNESNASMGetSNES - Gets a subsolver

927:    Not collective

929:    Input Parameters:
930: +  snes - the SNES context
931: -  i - the number of the subsnes to get

933:    Output Parameters:
934: .  subsnes - the subsolver context

936:    Level: intermediate

938: .seealso: SNESNASM, SNESNASMGetNumber()
939: @*/
940: PetscErrorCode SNESNASMGetSNES(SNES snes,PetscInt i,SNES *subsnes)
941: {
942:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

945:   if (i < 0 || i >= nasm->n) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"No such subsolver");
946:   *subsnes = nasm->subsnes[i];
947:   return(0);
948: }

950: /*@
951:    SNESNASMGetNumber - Gets number of subsolvers

953:    Not collective

955:    Input Parameters:
956: .  snes - the SNES context

958:    Output Parameters:
959: .  n - the number of subsolvers

961:    Level: intermediate

963: .seealso: SNESNASM, SNESNASMGetSNES()
964: @*/
965: PetscErrorCode SNESNASMGetNumber(SNES snes,PetscInt *n)
966: {
967:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

970:   *n = nasm->n;
971:   return(0);
972: }

974: /*@
975:    SNESNASMSetWeight - Sets weight to use when adding overlapping updates

977:    Collective

979:    Input Parameters:
980: +  snes - the SNES context
981: -  weight - the weights to use (typically 1/N for each dof, where N is the number of patches it appears in)

983:    Level: intermediate

985: .seealso: SNESNASM
986: @*/
987: PetscErrorCode SNESNASMSetWeight(SNES snes,Vec weight)
988: {
989:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;


994:   VecDestroy(&nasm->weight);
995:   nasm->weight_set = PETSC_TRUE;
996:   nasm->weight     = weight;
997:   PetscObjectReference((PetscObject)nasm->weight);

999:   return(0);
1000: }