Actual source code: nasm.c

petsc-3.9.4 2018-09-11
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 Schwartz 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: .keywords: SNES, NASM

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

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

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

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

311: /*@
312:    SNESNASMGetType - Get the type of subdomain update used

314:    Logically Collective on SNES

316:    Input Parameters:
317: .  SNES - the SNES context

319:    Output Parameters:
320: .  type - the type of update

322:    Level: intermediate

324: .keywords: SNES, NASM

326: .seealso: SNESNASM, SNESNASMSetType(), PCASMGetType()
327: @*/
328: PetscErrorCode SNESNASMGetType(SNES snes,PCASMType *type)
329: {

333:   PetscUseMethod(snes,"SNESNASMGetType_C",(SNES,PCASMType*),(snes,type));
334:   return(0);
335: }

337: static PetscErrorCode SNESNASMGetType_NASM(SNES snes,PCASMType *type)
338: {
339:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

342:   *type = nasm->type;
343:   return(0);
344: }

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

349:    Not Collective

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

359:    Level: intermediate

361: .keywords: SNES, NASM

363: .seealso: SNESNASM, SNESNASMGetSubdomains()
364: @*/
365: PetscErrorCode SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
366: {
368:   PetscErrorCode (*f)(SNES,PetscInt,SNES*,VecScatter*,VecScatter*,VecScatter*);

371:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",&f);
372:   if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
373:   return(0);
374: }

376: static PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
377: {
378:   PetscInt       i;
380:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

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

385:   /* tear down the previously set things */
386:   SNESReset(snes);

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

419:   if (subsnes) {
420:     PetscMalloc1(n,&nasm->subsnes);
421:     for (i=0; i<n; i++) {
422:       nasm->subsnes[i] = subsnes[i];
423:     }
424:     nasm->same_local_solves = PETSC_FALSE;
425:   }
426:   return(0);
427: }

429: /*@
430:    SNESNASMGetSubdomains - Get the local subdomain context.

432:    Not Collective

434:    Input Parameters:
435: .  SNES - the SNES context

437:    Output Parameters:
438: +  n - the number of local subdomains
439: .  subsnes - solvers defined on the local subdomains
440: .  iscatter - scatters into the nonoverlapping portions of the local subdomains
441: .  oscatter - scatters into the overlapping portions of the local subdomains
442: -  gscatter - scatters into the (ghosted) local vector of the local subdomain

444:    Level: intermediate

446: .keywords: SNES, NASM

448: .seealso: SNESNASM, SNESNASMSetSubdomains()
449: @*/
450: PetscErrorCode SNESNASMGetSubdomains(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
451: {
453:   PetscErrorCode (*f)(SNES,PetscInt*,SNES**,VecScatter**,VecScatter**,VecScatter**);

456:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",&f);
457:   if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
458:   return(0);
459: }

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

466:   if (n) *n = nasm->n;
467:   if (oscatter) *oscatter = nasm->oscatter;
468:   if (iscatter) *iscatter = nasm->iscatter;
469:   if (gscatter) *gscatter = nasm->gscatter;
470:   if (subsnes)  {
471:     *subsnes  = nasm->subsnes;
472:     nasm->same_local_solves = PETSC_FALSE;
473:   }
474:   return(0);
475: }

477: /*@
478:    SNESNASMGetSubdomainVecs - Get the processor-local subdomain vectors

480:    Not Collective

482:    Input Parameters:
483: .  SNES - the SNES context

485:    Output Parameters:
486: +  n - the number of local subdomains
487: .  x - The subdomain solution vector
488: .  y - The subdomain step vector
489: .  b - The subdomain RHS vector
490: -  xl - The subdomain local vectors (ghosted)

492:    Level: developer

494: .keywords: SNES, NASM

496: .seealso: SNESNASM, SNESNASMGetSubdomains()
497: @*/
498: PetscErrorCode SNESNASMGetSubdomainVecs(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b, Vec **xl)
499: {
501:   PetscErrorCode (*f)(SNES,PetscInt*,Vec**,Vec**,Vec**,Vec**);

504:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",&f);
505:   if (f) {(f)(snes,n,x,y,b,xl);}
506:   return(0);
507: }

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

514:   if (n)  *n  = nasm->n;
515:   if (x)  *x  = nasm->x;
516:   if (y)  *y  = nasm->y;
517:   if (b)  *b  = nasm->b;
518:   if (xl) *xl = nasm->xl;
519:   return(0);
520: }

522: /*@
523:    SNESNASMSetComputeFinalJacobian - Schedules the computation of the global and subdomain jacobians upon convergence

525:    Collective on SNES

527:    Input Parameters:
528: +  SNES - the SNES context
529: -  flg - indication of whether to compute the jacobians or not

531:    Level: developer

533:    Notes: This is used almost exclusively in the implementation of ASPIN, where the converged subdomain and global jacobian
534:    is needed at each linear iteration.

536: .keywords: SNES, NASM, ASPIN

538: .seealso: SNESNASM, SNESNASMGetSubdomains()
539: @*/
540: PetscErrorCode SNESNASMSetComputeFinalJacobian(SNES snes,PetscBool flg)
541: {
542:   PetscErrorCode (*f)(SNES,PetscBool);

546:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",&f);
547:   if (f) {(f)(snes,flg);}
548:   return(0);
549: }

551: static PetscErrorCode SNESNASMSetComputeFinalJacobian_NASM(SNES snes,PetscBool flg)
552: {
553:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

556:   nasm->finaljacobian = flg;
557:   if (flg) snes->usesksp = PETSC_TRUE;
558:   return(0);
559: }

561: /*@
562:    SNESNASMSetDamping - Sets the update damping for NASM

564:    Logically collective on SNES

566:    Input Parameters:
567: +  SNES - the SNES context
568: -  dmp - damping

570:    Level: intermediate

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

574: .keywords: SNES, NASM, damping

576: .seealso: SNESNASM, SNESNASMGetDamping()
577: @*/
578: PetscErrorCode SNESNASMSetDamping(SNES snes,PetscReal dmp)
579: {
580:   PetscErrorCode (*f)(SNES,PetscReal);

584:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetDamping_C",(void (**)(void))&f);
585:   if (f) {(f)(snes,dmp);}
586:   return(0);
587: }

589: static PetscErrorCode SNESNASMSetDamping_NASM(SNES snes,PetscReal dmp)
590: {
591:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

594:   nasm->damping = dmp;
595:   return(0);
596: }

598: /*@
599:    SNESNASMGetDamping - Gets the update damping for NASM

601:    Not Collective

603:    Input Parameters:
604: +  SNES - the SNES context
605: -  dmp - damping

607:    Level: intermediate

609: .keywords: SNES, NASM, damping

611: .seealso: SNESNASM, SNESNASMSetDamping()
612: @*/
613: PetscErrorCode SNESNASMGetDamping(SNES snes,PetscReal *dmp)
614: {

618:   PetscUseMethod(snes,"SNESNASMGetDamping_C",(SNES,PetscReal*),(snes,dmp));
619:   return(0);
620: }

622: static PetscErrorCode SNESNASMGetDamping_NASM(SNES snes,PetscReal *dmp)
623: {
624:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

627:   *dmp = nasm->damping;
628:   return(0);
629: }


632: /*
633:   Input Parameters:
634: + snes - The solver
635: . B - The RHS vector
636: - X - The initial guess

638:   Output Parameters:
639: . Y - The solution update

641:   TODO: All scatters should be packed into one
642: */
643: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X)
644: {
645:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
646:   SNES           subsnes;
647:   PetscInt       i;
648:   PetscReal      dmp;
650:   Vec            Xl,Bl,Yl,Xlloc;
651:   VecScatter     iscat,oscat,gscat,oscat_copy;
652:   DM             dm,subdm;
653:   PCASMType      type;

656:   SNESNASMGetType(snes,&type);
657:   SNESGetDM(snes,&dm);
658:   VecSet(Y,0);
659:   if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
660:   for (i=0; i<nasm->n; i++) {
661:     /* scatter the solution to the global solution and the local solution */
662:     Xl      = nasm->x[i];
663:     Xlloc   = nasm->xl[i];
664:     oscat   = nasm->oscatter[i];
665:     oscat_copy = nasm->oscatter_copy[i];
666:     gscat   = nasm->gscatter[i];
667:     VecScatterBegin(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
668:     VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
669:     if (B) {
670:       /* scatter the RHS to the local RHS */
671:       Bl   = nasm->b[i];
672:       VecScatterBegin(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
673:     }
674:   }
675:   if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}


678:   if (nasm->eventsubsolve) {PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);}
679:   for (i=0; i<nasm->n; i++) {
680:     Xl    = nasm->x[i];
681:     Xlloc = nasm->xl[i];
682:     Yl    = nasm->y[i];
683:     subsnes = nasm->subsnes[i];
684:     SNESGetDM(subsnes,&subdm);
685:     iscat   = nasm->iscatter[i];
686:     oscat   = nasm->oscatter[i];
687:     oscat_copy = nasm->oscatter_copy[i];
688:     gscat   = nasm->gscatter[i];
689:     VecScatterEnd(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
690:     VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
691:     if (B) {
692:       Bl   = nasm->b[i];
693:       VecScatterEnd(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
694:     } else Bl = NULL;

696:     DMSubDomainRestrict(dm,oscat,gscat,subdm);
697:     VecCopy(Xl,Yl);
698:     SNESSolve(subsnes,Bl,Xl);
699:     VecAYPX(Yl,-1.0,Xl);
700:     VecScale(Yl, nasm->damping);
701:     if (type == PC_ASM_BASIC) {
702:       VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
703:     } else if (type == PC_ASM_RESTRICT) {
704:       VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
705:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
706:   }
707:   if (nasm->eventsubsolve) {PetscLogEventEnd(nasm->eventsubsolve,snes,0,0,0);}
708:   if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
709:   for (i=0; i<nasm->n; i++) {
710:     Yl    = nasm->y[i];
711:     iscat   = nasm->iscatter[i];
712:     oscat   = nasm->oscatter[i];
713:     if (type == PC_ASM_BASIC) {
714:       VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
715:     } else if (type == PC_ASM_RESTRICT) {
716:       VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
717:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
718:   }
719:   if (nasm->weight_set) {
720:     VecPointwiseMult(Y,Y,nasm->weight);
721:   }
722:   if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
723:   SNESNASMGetDamping(snes,&dmp);
724:   VecAXPY(X,dmp,Y);
725:   return(0);
726: }

728: static PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
729: {
730:   Vec            X = Xfinal;
731:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
732:   SNES           subsnes;
733:   PetscInt       i,lag = 1;
735:   Vec            Xlloc,Xl,Fl,F;
736:   VecScatter     oscat,gscat;
737:   DM             dm,subdm;

740:   if (nasm->fjtype == 2) X = nasm->xinit;
741:   F = snes->vec_func;
742:   if (snes->normschedule == SNES_NORM_NONE) {SNESComputeFunction(snes,X,F);}
743:   SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
744:   SNESGetDM(snes,&dm);
745:   if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
746:   if (nasm->fjtype != 1) {
747:     for (i=0; i<nasm->n; i++) {
748:       Xlloc = nasm->xl[i];
749:       gscat = nasm->gscatter[i];
750:       VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
751:     }
752:   }
753:   if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
754:   for (i=0; i<nasm->n; i++) {
755:     Fl      = nasm->subsnes[i]->vec_func;
756:     Xl      = nasm->x[i];
757:     Xlloc   = nasm->xl[i];
758:     subsnes = nasm->subsnes[i];
759:     oscat   = nasm->oscatter[i];
760:     gscat   = nasm->gscatter[i];
761:     if (nasm->fjtype != 1) {VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);}
762:     SNESGetDM(subsnes,&subdm);
763:     DMSubDomainRestrict(dm,oscat,gscat,subdm);
764:     if (nasm->fjtype != 1) {
765:       DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
766:       DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
767:     }
768:     if (subsnes->lagjacobian == -1)    subsnes->lagjacobian = -2;
769:     else if (subsnes->lagjacobian > 1) lag = subsnes->lagjacobian;
770:     SNESComputeFunction(subsnes,Xl,Fl);
771:     SNESComputeJacobian(subsnes,Xl,subsnes->jacobian,subsnes->jacobian_pre);
772:     if (lag > 1) subsnes->lagjacobian = lag;
773:   }
774:   return(0);
775: }

777: static PetscErrorCode SNESSolve_NASM(SNES snes)
778: {
779:   Vec              F;
780:   Vec              X;
781:   Vec              B;
782:   Vec              Y;
783:   PetscInt         i;
784:   PetscReal        fnorm = 0.0;
785:   PetscErrorCode   ierr;
786:   SNESNormSchedule normschedule;
787:   SNES_NASM        *nasm = (SNES_NASM*)snes->data;


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

793:   PetscCitationsRegister(SNESCitation,&SNEScite);
794:   X = snes->vec_sol;
795:   Y = snes->vec_sol_update;
796:   F = snes->vec_func;
797:   B = snes->vec_rhs;

799:   PetscObjectSAWsTakeAccess((PetscObject)snes);
800:   snes->iter   = 0;
801:   snes->norm   = 0.;
802:   PetscObjectSAWsGrantAccess((PetscObject)snes);
803:   snes->reason = SNES_CONVERGED_ITERATING;
804:   SNESGetNormSchedule(snes, &normschedule);
805:   if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
806:     /* compute the initial function and preconditioned update delX */
807:     if (!snes->vec_func_init_set) {
808:       SNESComputeFunction(snes,X,F);
809:     } else snes->vec_func_init_set = PETSC_FALSE;

811:     VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
812:     SNESCheckFunctionNorm(snes,fnorm);
813:     PetscObjectSAWsTakeAccess((PetscObject)snes);
814:     snes->iter = 0;
815:     snes->norm = fnorm;
816:     PetscObjectSAWsGrantAccess((PetscObject)snes);
817:     SNESLogConvergenceHistory(snes,snes->norm,0);
818:     SNESMonitor(snes,0,snes->norm);

820:     /* test convergence */
821:     (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
822:     if (snes->reason) return(0);
823:   } else {
824:     PetscObjectSAWsGrantAccess((PetscObject)snes);
825:     SNESLogConvergenceHistory(snes,snes->norm,0);
826:     SNESMonitor(snes,0,snes->norm);
827:   }

829:   /* Call general purpose update function */
830:   if (snes->ops->update) {
831:     (*snes->ops->update)(snes, snes->iter);
832:   }
833:   /* copy the initial solution over for later */
834:   if (nasm->fjtype == 2) {VecCopy(X,nasm->xinit);}

836:   for (i=0; i < snes->max_its; i++) {
837:     SNESNASMSolveLocal_Private(snes,B,Y,X);
838:     if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
839:       SNESComputeFunction(snes,X,F);
840:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
841:       SNESCheckFunctionNorm(snes,fnorm);
842:     }
843:     /* Monitor convergence */
844:     PetscObjectSAWsTakeAccess((PetscObject)snes);
845:     snes->iter = i+1;
846:     snes->norm = fnorm;
847:     PetscObjectSAWsGrantAccess((PetscObject)snes);
848:     SNESLogConvergenceHistory(snes,snes->norm,0);
849:     SNESMonitor(snes,snes->iter,snes->norm);
850:     /* Test for convergence */
851:     if (normschedule == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);}
852:     if (snes->reason) break;
853:     /* Call general purpose update function */
854:     if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
855:   }
856:   if (nasm->finaljacobian) {SNESNASMComputeFinalJacobian_Private(snes,X);}
857:   if (normschedule == SNES_NORM_ALWAYS) {
858:     if (i == snes->max_its) {
859:       PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
860:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
861:     }
862:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */
863:   return(0);
864: }

866: /*MC
867:   SNESNASM - Nonlinear Additive Schwartz

869:    Options Database:
870: +  -snes_nasm_log - enable logging events for the communication and solve stages
871: .  -snes_nasm_type <basic,restrict> - type of subdomain update used
872: .  -snes_asm_damping <dmp> - the new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
873: .  -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
874: .  -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> - pick state the jacobian is calculated at
875: .  -sub_snes_ - options prefix of the subdomain nonlinear solves
876: .  -sub_ksp_ - options prefix of the subdomain Krylov solver
877: -  -sub_pc_ - options prefix of the subdomain preconditioner

879:    Level: advanced

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

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

888: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
889: {
890:   SNES_NASM      *nasm;

894:   PetscNewLog(snes,&nasm);
895:   snes->data = (void*)nasm;

897:   nasm->n        = PETSC_DECIDE;
898:   nasm->subsnes  = 0;
899:   nasm->x        = 0;
900:   nasm->xl       = 0;
901:   nasm->y        = 0;
902:   nasm->b        = 0;
903:   nasm->oscatter = 0;
904:   nasm->oscatter_copy = 0;
905:   nasm->iscatter = 0;
906:   nasm->gscatter = 0;
907:   nasm->damping  = 1.;

909:   nasm->type              = PC_ASM_BASIC;
910:   nasm->finaljacobian     = PETSC_FALSE;
911:   nasm->same_local_solves = PETSC_TRUE;
912:   nasm->weight_set        = PETSC_FALSE;

914:   snes->ops->destroy        = SNESDestroy_NASM;
915:   snes->ops->setup          = SNESSetUp_NASM;
916:   snes->ops->setfromoptions = SNESSetFromOptions_NASM;
917:   snes->ops->view           = SNESView_NASM;
918:   snes->ops->solve          = SNESSolve_NASM;
919:   snes->ops->reset          = SNESReset_NASM;

921:   snes->usesksp = PETSC_FALSE;
922:   snes->usesnpc = PETSC_FALSE;

924:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

926:   nasm->fjtype              = 0;
927:   nasm->xinit               = NULL;
928:   nasm->eventrestrictinterp = 0;
929:   nasm->eventsubsolve       = 0;

931:   if (!snes->tolerancesset) {
932:     snes->max_its   = 10000;
933:     snes->max_funcs = 10000;
934:   }

936:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetType_C",SNESNASMSetType_NASM);
937:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetType_C",SNESNASMGetType_NASM);
938:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
939:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
940:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
941:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
942:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
943:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
944:   return(0);
945: }

947: /*@
948:    SNESNASMGetSNES - Gets a subsolver

950:    Not collective

952:    Input Parameters:
953: +  snes - the SNES context
954: -  i - the number of the subsnes to get

956:    Output Parameters:
957: .  subsnes - the subsolver context

959:    Level: intermediate

961: .keywords: SNES, NASM

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

970:   if (i < 0 || i >= nasm->n) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"No such subsolver");
971:   *subsnes = nasm->subsnes[i];
972:   return(0);
973: }

975: /*@
976:    SNESNASMGetNumber - Gets number of subsolvers

978:    Not collective

980:    Input Parameters:
981: .  snes - the SNES context

983:    Output Parameters:
984: .  n - the number of subsolvers

986:    Level: intermediate

988: .keywords: SNES, NASM

990: .seealso: SNESNASM, SNESNASMGetSNES()
991: @*/
992: PetscErrorCode SNESNASMGetNumber(SNES snes,PetscInt *n)
993: {
994:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

997:   *n = nasm->n;
998:   return(0);
999: }

1001: /*@
1002:    SNESNASMSetWeight - Sets weight to use when adding overlapping updates

1004:    Collective

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

1010:    Level: intermediate

1012: .keywords: SNES, NASM

1014: .seealso: SNESNASM
1015: @*/
1016: PetscErrorCode SNESNASMSetWeight(SNES snes,Vec weight)
1017: {
1018:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;


1023:   VecDestroy(&nasm->weight);
1024:   nasm->weight_set = PETSC_TRUE;
1025:   nasm->weight     = weight;
1026:   PetscObjectReference((PetscObject)nasm->weight);

1028:   return(0);
1029: }