Actual source code: nasm.c

petsc-3.11.2 2019-05-18
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:   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:
573:     The new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)

575: .keywords: SNES, NASM, damping

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

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

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

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

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

602:    Not Collective

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

608:    Level: intermediate

610: .keywords: SNES, NASM, damping

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

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

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

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


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

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

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

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


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

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

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

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

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


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

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

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

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

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

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

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

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

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

883:    Level: advanced

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

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

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

899: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
900: {
901:   SNES_NASM      *nasm;

905:   PetscNewLog(snes,&nasm);
906:   snes->data = (void*)nasm;

908:   nasm->n        = PETSC_DECIDE;
909:   nasm->subsnes  = 0;
910:   nasm->x        = 0;
911:   nasm->xl       = 0;
912:   nasm->y        = 0;
913:   nasm->b        = 0;
914:   nasm->oscatter = 0;
915:   nasm->oscatter_copy = 0;
916:   nasm->iscatter = 0;
917:   nasm->gscatter = 0;
918:   nasm->damping  = 1.;

920:   nasm->type              = PC_ASM_BASIC;
921:   nasm->finaljacobian     = PETSC_FALSE;
922:   nasm->same_local_solves = PETSC_TRUE;
923:   nasm->weight_set        = PETSC_FALSE;

925:   snes->ops->destroy        = SNESDestroy_NASM;
926:   snes->ops->setup          = SNESSetUp_NASM;
927:   snes->ops->setfromoptions = SNESSetFromOptions_NASM;
928:   snes->ops->view           = SNESView_NASM;
929:   snes->ops->solve          = SNESSolve_NASM;
930:   snes->ops->reset          = SNESReset_NASM;

932:   snes->usesksp = PETSC_FALSE;
933:   snes->usesnpc = PETSC_FALSE;

935:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

937:   nasm->fjtype              = 0;
938:   nasm->xinit               = NULL;
939:   nasm->eventrestrictinterp = 0;
940:   nasm->eventsubsolve       = 0;

942:   if (!snes->tolerancesset) {
943:     snes->max_its   = 10000;
944:     snes->max_funcs = 10000;
945:   }

947:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetType_C",SNESNASMSetType_NASM);
948:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetType_C",SNESNASMGetType_NASM);
949:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
950:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
951:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
952:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
953:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
954:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
955:   return(0);
956: }

958: /*@
959:    SNESNASMGetSNES - Gets a subsolver

961:    Not collective

963:    Input Parameters:
964: +  snes - the SNES context
965: -  i - the number of the subsnes to get

967:    Output Parameters:
968: .  subsnes - the subsolver context

970:    Level: intermediate

972: .keywords: SNES, NASM

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

981:   if (i < 0 || i >= nasm->n) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"No such subsolver");
982:   *subsnes = nasm->subsnes[i];
983:   return(0);
984: }

986: /*@
987:    SNESNASMGetNumber - Gets number of subsolvers

989:    Not collective

991:    Input Parameters:
992: .  snes - the SNES context

994:    Output Parameters:
995: .  n - the number of subsolvers

997:    Level: intermediate

999: .keywords: SNES, NASM

1001: .seealso: SNESNASM, SNESNASMGetSNES()
1002: @*/
1003: PetscErrorCode SNESNASMGetNumber(SNES snes,PetscInt *n)
1004: {
1005:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

1008:   *n = nasm->n;
1009:   return(0);
1010: }

1012: /*@
1013:    SNESNASMSetWeight - Sets weight to use when adding overlapping updates

1015:    Collective

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

1021:    Level: intermediate

1023: .keywords: SNES, NASM

1025: .seealso: SNESNASM
1026: @*/
1027: PetscErrorCode SNESNASMSetWeight(SNES snes,Vec weight)
1028: {
1029:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;


1034:   VecDestroy(&nasm->weight);
1035:   nasm->weight_set = PETSC_TRUE;
1036:   nasm->weight     = weight;
1037:   PetscObjectReference((PetscObject)nasm->weight);

1039:   return(0);
1040: }