Actual source code: nasm.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  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  same_local_solves;   /* flag to determine if the solvers have been individually modified */
 21:   PetscBool  weight_set;          /* use a weight in the overlap updates */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

188:   PetscOptionsHead(PetscOptionsObject,"Nonlinear Additive Schwarz options");
189:   PetscOptionsEnum("-snes_nasm_type","Type of restriction/extension","",SNESNASMTypes,(PetscEnum)nasm->type,(PetscEnum*)&asmtype,&flg);
190:   if (flg) {SNESNASMSetType(snes,asmtype);}
191:   flg    = PETSC_FALSE;
192:   monflg = PETSC_TRUE;
193:   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);
194:   if (flg) {SNESNASMSetDamping(snes,nasm->damping);}
195:   subviewflg = PETSC_FALSE;
196:   PetscOptionsBool("-snes_nasm_sub_view","Print detailed information for every processor when using -snes_view","",subviewflg,&subviewflg,&flg);
197:   if (flg) {
198:     nasm->same_local_solves = PETSC_FALSE;
199:     if (!subviewflg) {
200:       nasm->same_local_solves = PETSC_TRUE;
201:     }
202:   }
203:   PetscOptionsBool("-snes_nasm_finaljacobian","Compute the global jacobian of the final iterate (for ASPIN)","",nasm->finaljacobian,&nasm->finaljacobian,NULL);
204:   PetscOptionsEList("-snes_nasm_finaljacobian_type","The type of the final jacobian computed.","",SNESNASMFJTypes,3,SNESNASMFJTypes[0],&nasm->fjtype,NULL);
205:   PetscOptionsBool("-snes_nasm_log","Log times for subSNES solves and restriction","",monflg,&monflg,&flg);
206:   if (flg) {
207:     PetscLogEventRegister("SNESNASMSubSolve",((PetscObject)snes)->classid,&nasm->eventsubsolve);
208:     PetscLogEventRegister("SNESNASMRestrict",((PetscObject)snes)->classid,&nasm->eventrestrictinterp);
209:   }
210:   PetscOptionsTail();
211:   return(0);
212: }

214: static PetscErrorCode SNESView_NASM(SNES snes, PetscViewer viewer)
215: {
216:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
218:   PetscMPIInt    rank,size;
219:   PetscInt       i,N,bsz;
220:   PetscBool      iascii,isstring;
221:   PetscViewer    sviewer;
222:   MPI_Comm       comm;

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

275: /*@
276:    SNESNASMSetType - Set the type of subdomain update used

278:    Logically Collective on SNES

280:    Input Parameters:
281: +  SNES - the SNES context
282: -  type - the type of update, PC_ASM_BASIC or PC_ASM_RESTRICT

284:    Level: intermediate

286: .seealso: SNESNASM, SNESNASMGetType(), PCASMSetType()
287: @*/
288: PetscErrorCode SNESNASMSetType(SNES snes,PCASMType type)
289: {
291:   PetscErrorCode (*f)(SNES,PCASMType);

294:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetType_C",&f);
295:   if (f) {(f)(snes,type);}
296:   return(0);
297: }

299: static PetscErrorCode SNESNASMSetType_NASM(SNES snes,PCASMType type)
300: {
301:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

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

309: /*@
310:    SNESNASMGetType - Get the type of subdomain update used

312:    Logically Collective on SNES

314:    Input Parameters:
315: .  SNES - the SNES context

317:    Output Parameters:
318: .  type - the type of update

320:    Level: intermediate

322: .seealso: SNESNASM, SNESNASMSetType(), PCASMGetType()
323: @*/
324: PetscErrorCode SNESNASMGetType(SNES snes,PCASMType *type)
325: {

329:   PetscUseMethod(snes,"SNESNASMGetType_C",(SNES,PCASMType*),(snes,type));
330:   return(0);
331: }

333: static PetscErrorCode SNESNASMGetType_NASM(SNES snes,PCASMType *type)
334: {
335:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

338:   *type = nasm->type;
339:   return(0);
340: }

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

345:    Not Collective

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

355:    Level: intermediate

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

365:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",&f);
366:   if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
367:   return(0);
368: }

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

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

379:   /* tear down the previously set things */
380:   SNESReset(snes);

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

413:   if (subsnes) {
414:     PetscMalloc1(n,&nasm->subsnes);
415:     for (i=0; i<n; i++) {
416:       nasm->subsnes[i] = subsnes[i];
417:     }
418:     nasm->same_local_solves = PETSC_FALSE;
419:   }
420:   return(0);
421: }

423: /*@
424:    SNESNASMGetSubdomains - Get the local subdomain context.

426:    Not Collective

428:    Input Parameters:
429: .  SNES - the SNES context

431:    Output Parameters:
432: +  n - the number of local subdomains
433: .  subsnes - solvers defined on the local subdomains
434: .  iscatter - scatters into the nonoverlapping portions of the local subdomains
435: .  oscatter - scatters into the overlapping portions of the local subdomains
436: -  gscatter - scatters into the (ghosted) local vector of the local subdomain

438:    Level: intermediate

440: .seealso: SNESNASM, SNESNASMSetSubdomains()
441: @*/
442: PetscErrorCode SNESNASMGetSubdomains(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
443: {
445:   PetscErrorCode (*f)(SNES,PetscInt*,SNES**,VecScatter**,VecScatter**,VecScatter**);

448:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",&f);
449:   if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
450:   return(0);
451: }

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

458:   if (n) *n = nasm->n;
459:   if (oscatter) *oscatter = nasm->oscatter;
460:   if (iscatter) *iscatter = nasm->iscatter;
461:   if (gscatter) *gscatter = nasm->gscatter;
462:   if (subsnes)  {
463:     *subsnes  = nasm->subsnes;
464:     nasm->same_local_solves = PETSC_FALSE;
465:   }
466:   return(0);
467: }

469: /*@
470:    SNESNASMGetSubdomainVecs - Get the processor-local subdomain vectors

472:    Not Collective

474:    Input Parameters:
475: .  SNES - the SNES context

477:    Output Parameters:
478: +  n - the number of local subdomains
479: .  x - The subdomain solution vector
480: .  y - The subdomain step vector
481: .  b - The subdomain RHS vector
482: -  xl - The subdomain local vectors (ghosted)

484:    Level: developer

486: .seealso: SNESNASM, SNESNASMGetSubdomains()
487: @*/
488: PetscErrorCode SNESNASMGetSubdomainVecs(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b, Vec **xl)
489: {
491:   PetscErrorCode (*f)(SNES,PetscInt*,Vec**,Vec**,Vec**,Vec**);

494:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",&f);
495:   if (f) {(f)(snes,n,x,y,b,xl);}
496:   return(0);
497: }

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

504:   if (n)  *n  = nasm->n;
505:   if (x)  *x  = nasm->x;
506:   if (y)  *y  = nasm->y;
507:   if (b)  *b  = nasm->b;
508:   if (xl) *xl = nasm->xl;
509:   return(0);
510: }

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

515:    Collective on SNES

517:    Input Parameters:
518: +  SNES - the SNES context
519: -  flg - indication of whether to compute the Jacobians or not

521:    Level: developer

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

527: .seealso: SNESNASM, SNESNASMGetSubdomains()
528: @*/
529: PetscErrorCode SNESNASMSetComputeFinalJacobian(SNES snes,PetscBool flg)
530: {
531:   PetscErrorCode (*f)(SNES,PetscBool);

535:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",&f);
536:   if (f) {(f)(snes,flg);}
537:   return(0);
538: }

540: static PetscErrorCode SNESNASMSetComputeFinalJacobian_NASM(SNES snes,PetscBool flg)
541: {
542:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

545:   nasm->finaljacobian = flg;
546:   return(0);
547: }

549: /*@
550:    SNESNASMSetDamping - Sets the update damping for NASM

552:    Logically collective on SNES

554:    Input Parameters:
555: +  SNES - the SNES context
556: -  dmp - damping

558:    Level: intermediate

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

563: .seealso: SNESNASM, SNESNASMGetDamping()
564: @*/
565: PetscErrorCode SNESNASMSetDamping(SNES snes,PetscReal dmp)
566: {
567:   PetscErrorCode (*f)(SNES,PetscReal);

571:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetDamping_C",(void (**)(void))&f);
572:   if (f) {(f)(snes,dmp);}
573:   return(0);
574: }

576: static PetscErrorCode SNESNASMSetDamping_NASM(SNES snes,PetscReal dmp)
577: {
578:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

581:   nasm->damping = dmp;
582:   return(0);
583: }

585: /*@
586:    SNESNASMGetDamping - Gets the update damping for NASM

588:    Not Collective

590:    Input Parameters:
591: +  SNES - the SNES context
592: -  dmp - damping

594:    Level: intermediate

596: .seealso: SNESNASM, SNESNASMSetDamping()
597: @*/
598: PetscErrorCode SNESNASMGetDamping(SNES snes,PetscReal *dmp)
599: {

603:   PetscUseMethod(snes,"SNESNASMGetDamping_C",(SNES,PetscReal*),(snes,dmp));
604:   return(0);
605: }

607: static PetscErrorCode SNESNASMGetDamping_NASM(SNES snes,PetscReal *dmp)
608: {
609:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

612:   *dmp = nasm->damping;
613:   return(0);
614: }


617: /*
618:   Input Parameters:
619: + snes - The solver
620: . B - The RHS vector
621: - X - The initial guess

623:   Output Parameters:
624: . Y - The solution update

626:   TODO: All scatters should be packed into one
627: */
628: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X)
629: {
630:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
631:   SNES           subsnes;
632:   PetscInt       i;
633:   PetscReal      dmp;
635:   Vec            Xl,Bl,Yl,Xlloc;
636:   VecScatter     iscat,oscat,gscat,oscat_copy;
637:   DM             dm,subdm;
638:   PCASMType      type;

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


663:   if (nasm->eventsubsolve) {PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);}
664:   for (i=0; i<nasm->n; i++) {
665:     Xl    = nasm->x[i];
666:     Xlloc = nasm->xl[i];
667:     Yl    = nasm->y[i];
668:     subsnes = nasm->subsnes[i];
669:     SNESGetDM(subsnes,&subdm);
670:     iscat   = nasm->iscatter[i];
671:     oscat   = nasm->oscatter[i];
672:     oscat_copy = nasm->oscatter_copy[i];
673:     gscat   = nasm->gscatter[i];
674:     VecScatterEnd(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
675:     VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
676:     if (B) {
677:       Bl   = nasm->b[i];
678:       VecScatterEnd(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
679:     } else Bl = NULL;

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

705: static PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
706: {
707:   Vec            X = Xfinal;
708:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
709:   SNES           subsnes;
710:   PetscInt       i,lag = 1;
712:   Vec            Xlloc,Xl,Fl,F;
713:   VecScatter     oscat,gscat;
714:   DM             dm,subdm;

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

754: static PetscErrorCode SNESSolve_NASM(SNES snes)
755: {
756:   Vec              F;
757:   Vec              X;
758:   Vec              B;
759:   Vec              Y;
760:   PetscInt         i;
761:   PetscReal        fnorm = 0.0;
762:   PetscErrorCode   ierr;
763:   SNESNormSchedule normschedule;
764:   SNES_NASM        *nasm = (SNES_NASM*)snes->data;


768:   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);

770:   PetscCitationsRegister(SNESCitation,&SNEScite);
771:   X = snes->vec_sol;
772:   Y = snes->vec_sol_update;
773:   F = snes->vec_func;
774:   B = snes->vec_rhs;

776:   PetscObjectSAWsTakeAccess((PetscObject)snes);
777:   snes->iter   = 0;
778:   snes->norm   = 0.;
779:   PetscObjectSAWsGrantAccess((PetscObject)snes);
780:   snes->reason = SNES_CONVERGED_ITERATING;
781:   SNESGetNormSchedule(snes, &normschedule);
782:   if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
783:     /* compute the initial function and preconditioned update delX */
784:     if (!snes->vec_func_init_set) {
785:       SNESComputeFunction(snes,X,F);
786:     } else snes->vec_func_init_set = PETSC_FALSE;

788:     VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
789:     SNESCheckFunctionNorm(snes,fnorm);
790:     PetscObjectSAWsTakeAccess((PetscObject)snes);
791:     snes->iter = 0;
792:     snes->norm = fnorm;
793:     PetscObjectSAWsGrantAccess((PetscObject)snes);
794:     SNESLogConvergenceHistory(snes,snes->norm,0);
795:     SNESMonitor(snes,0,snes->norm);

797:     /* test convergence */
798:     (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
799:     if (snes->reason) return(0);
800:   } else {
801:     PetscObjectSAWsGrantAccess((PetscObject)snes);
802:     SNESLogConvergenceHistory(snes,snes->norm,0);
803:     SNESMonitor(snes,0,snes->norm);
804:   }

806:   /* Call general purpose update function */
807:   if (snes->ops->update) {
808:     (*snes->ops->update)(snes, snes->iter);
809:   }
810:   /* copy the initial solution over for later */
811:   if (nasm->fjtype == 2) {VecCopy(X,nasm->xinit);}

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

846: /*MC
847:   SNESNASM - Nonlinear Additive Schwarz

849:    Options Database:
850: +  -snes_nasm_log - enable logging events for the communication and solve stages
851: .  -snes_nasm_type <basic,restrict> - type of subdomain update used
852: .  -snes_asm_damping <dmp> - the new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
853: .  -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
854: .  -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> - pick state the jacobian is calculated at
855: .  -sub_snes_ - options prefix of the subdomain nonlinear solves
856: .  -sub_ksp_ - options prefix of the subdomain Krylov solver
857: -  -sub_pc_ - options prefix of the subdomain preconditioner

859:    Level: advanced

861:    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
862:        false and SNESView() and -snes_view do not display a KSP object. However, if the flag nasm->finaljacobian is set (for example, if
863:        NASM is used as a nonlinear preconditioner for  KSPASPIN) then SNESSetUpMatrices() is called to generate the Jacobian (needed by KSPASPIN)
864:        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
865:        used by SNESASPIN they share the same Jacobian matrices because SNESSetUp() (called on the outer SNES KSPASPIN) causes the inner SNES
866:        object (in this case SNESNASM) to inherit the outer Jacobian matrices.

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

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

875: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
876: {
877:   SNES_NASM      *nasm;

881:   PetscNewLog(snes,&nasm);
882:   snes->data = (void*)nasm;

884:   nasm->n        = PETSC_DECIDE;
885:   nasm->subsnes  = 0;
886:   nasm->x        = 0;
887:   nasm->xl       = 0;
888:   nasm->y        = 0;
889:   nasm->b        = 0;
890:   nasm->oscatter = 0;
891:   nasm->oscatter_copy = 0;
892:   nasm->iscatter = 0;
893:   nasm->gscatter = 0;
894:   nasm->damping  = 1.;

896:   nasm->type              = PC_ASM_BASIC;
897:   nasm->finaljacobian     = PETSC_FALSE;
898:   nasm->same_local_solves = PETSC_TRUE;
899:   nasm->weight_set        = PETSC_FALSE;

901:   snes->ops->destroy        = SNESDestroy_NASM;
902:   snes->ops->setup          = SNESSetUp_NASM;
903:   snes->ops->setfromoptions = SNESSetFromOptions_NASM;
904:   snes->ops->view           = SNESView_NASM;
905:   snes->ops->solve          = SNESSolve_NASM;
906:   snes->ops->reset          = SNESReset_NASM;

908:   snes->usesksp = PETSC_FALSE;
909:   snes->usesnpc = PETSC_FALSE;

911:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

913:   nasm->fjtype              = 0;
914:   nasm->xinit               = NULL;
915:   nasm->eventrestrictinterp = 0;
916:   nasm->eventsubsolve       = 0;

918:   if (!snes->tolerancesset) {
919:     snes->max_its   = 10000;
920:     snes->max_funcs = 10000;
921:   }

923:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetType_C",SNESNASMSetType_NASM);
924:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetType_C",SNESNASMGetType_NASM);
925:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
926:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
927:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
928:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
929:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
930:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
931:   return(0);
932: }

934: /*@
935:    SNESNASMGetSNES - Gets a subsolver

937:    Not collective

939:    Input Parameters:
940: +  snes - the SNES context
941: -  i - the number of the subsnes to get

943:    Output Parameters:
944: .  subsnes - the subsolver context

946:    Level: intermediate

948: .seealso: SNESNASM, SNESNASMGetNumber()
949: @*/
950: PetscErrorCode SNESNASMGetSNES(SNES snes,PetscInt i,SNES *subsnes)
951: {
952:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

955:   if (i < 0 || i >= nasm->n) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"No such subsolver");
956:   *subsnes = nasm->subsnes[i];
957:   return(0);
958: }

960: /*@
961:    SNESNASMGetNumber - Gets number of subsolvers

963:    Not collective

965:    Input Parameters:
966: .  snes - the SNES context

968:    Output Parameters:
969: .  n - the number of subsolvers

971:    Level: intermediate

973: .seealso: SNESNASM, SNESNASMGetSNES()
974: @*/
975: PetscErrorCode SNESNASMGetNumber(SNES snes,PetscInt *n)
976: {
977:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

980:   *n = nasm->n;
981:   return(0);
982: }

984: /*@
985:    SNESNASMSetWeight - Sets weight to use when adding overlapping updates

987:    Collective

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

993:    Level: intermediate

995: .seealso: SNESNASM
996: @*/
997: PetscErrorCode SNESNASMSetWeight(SNES snes,Vec weight)
998: {
999:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;


1004:   VecDestroy(&nasm->weight);
1005:   nasm->weight_set = PETSC_TRUE;
1006:   nasm->weight     = weight;
1007:   PetscObjectReference((PetscObject)nasm->weight);

1009:   return(0);
1010: }