Actual source code: nasm.c

petsc-3.11.1 2019-04-17
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:
534:     This is used almost exclusively in the implementation of ASPIN, where the converged subdomain and global jacobian
535:    is needed at each linear iteration.

537: .keywords: SNES, NASM, ASPIN

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

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

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

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

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

565:    Logically collective on SNES

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

571:    Level: intermediate

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

576: .keywords: SNES, NASM, damping

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

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

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

596:   nasm->damping = dmp;
597:   return(0);
598: }

600: /*@
601:    SNESNASMGetDamping - Gets the update damping for NASM

603:    Not Collective

605:    Input Parameters:
606: +  SNES - the SNES context
607: -  dmp - damping

609:    Level: intermediate

611: .keywords: SNES, NASM, damping

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

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

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

629:   *dmp = nasm->damping;
630:   return(0);
631: }


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

640:   Output Parameters:
641: . Y - The solution update

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

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


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

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

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

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

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


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

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

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

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

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

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

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

871: /*MC
872:   SNESNASM - Nonlinear Additive Schwartz

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

884:    Level: advanced

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

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

893: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
894: {
895:   SNES_NASM      *nasm;

899:   PetscNewLog(snes,&nasm);
900:   snes->data = (void*)nasm;

902:   nasm->n        = PETSC_DECIDE;
903:   nasm->subsnes  = 0;
904:   nasm->x        = 0;
905:   nasm->xl       = 0;
906:   nasm->y        = 0;
907:   nasm->b        = 0;
908:   nasm->oscatter = 0;
909:   nasm->oscatter_copy = 0;
910:   nasm->iscatter = 0;
911:   nasm->gscatter = 0;
912:   nasm->damping  = 1.;

914:   nasm->type              = PC_ASM_BASIC;
915:   nasm->finaljacobian     = PETSC_FALSE;
916:   nasm->same_local_solves = PETSC_TRUE;
917:   nasm->weight_set        = PETSC_FALSE;

919:   snes->ops->destroy        = SNESDestroy_NASM;
920:   snes->ops->setup          = SNESSetUp_NASM;
921:   snes->ops->setfromoptions = SNESSetFromOptions_NASM;
922:   snes->ops->view           = SNESView_NASM;
923:   snes->ops->solve          = SNESSolve_NASM;
924:   snes->ops->reset          = SNESReset_NASM;

926:   snes->usesksp = PETSC_FALSE;
927:   snes->usesnpc = PETSC_FALSE;

929:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

931:   nasm->fjtype              = 0;
932:   nasm->xinit               = NULL;
933:   nasm->eventrestrictinterp = 0;
934:   nasm->eventsubsolve       = 0;

936:   if (!snes->tolerancesset) {
937:     snes->max_its   = 10000;
938:     snes->max_funcs = 10000;
939:   }

941:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetType_C",SNESNASMSetType_NASM);
942:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetType_C",SNESNASMGetType_NASM);
943:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
944:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
945:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
946:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
947:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
948:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
949:   return(0);
950: }

952: /*@
953:    SNESNASMGetSNES - Gets a subsolver

955:    Not collective

957:    Input Parameters:
958: +  snes - the SNES context
959: -  i - the number of the subsnes to get

961:    Output Parameters:
962: .  subsnes - the subsolver context

964:    Level: intermediate

966: .keywords: SNES, NASM

968: .seealso: SNESNASM, SNESNASMGetNumber()
969: @*/
970: PetscErrorCode SNESNASMGetSNES(SNES snes,PetscInt i,SNES *subsnes)
971: {
972:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

975:   if (i < 0 || i >= nasm->n) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"No such subsolver");
976:   *subsnes = nasm->subsnes[i];
977:   return(0);
978: }

980: /*@
981:    SNESNASMGetNumber - Gets number of subsolvers

983:    Not collective

985:    Input Parameters:
986: .  snes - the SNES context

988:    Output Parameters:
989: .  n - the number of subsolvers

991:    Level: intermediate

993: .keywords: SNES, NASM

995: .seealso: SNESNASM, SNESNASMGetSNES()
996: @*/
997: PetscErrorCode SNESNASMGetNumber(SNES snes,PetscInt *n)
998: {
999:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

1002:   *n = nasm->n;
1003:   return(0);
1004: }

1006: /*@
1007:    SNESNASMSetWeight - Sets weight to use when adding overlapping updates

1009:    Collective

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

1015:    Level: intermediate

1017: .keywords: SNES, NASM

1019: .seealso: SNESNASM
1020: @*/
1021: PetscErrorCode SNESNASMSetWeight(SNES snes,Vec weight)
1022: {
1023:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;


1028:   VecDestroy(&nasm->weight);
1029:   nasm->weight_set = PETSC_TRUE;
1030:   nasm->weight     = weight;
1031:   PetscObjectReference((PetscObject)nasm->weight);

1033:   return(0);
1034: }