Actual source code: nasm.c

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

  4: typedef struct {
  5:   PetscInt   n;                   /* local subdomains */
  6:   SNES       *subsnes;            /* nonlinear solvers for each subdomain */
  7:   Vec        *x;                  /* solution vectors */
  8:   Vec        *xl;                 /* solution local vectors */
  9:   Vec        *y;                  /* step vectors */
 10:   Vec        *b;                  /* rhs vectors */
 11:   Vec        weight;              /* weighting for adding updates on overlaps, in global space */
 12:   VecScatter *oscatter;           /* scatter from global space to the subdomain global space */
 13:   VecScatter *oscatter_copy;      /* copy of the above */
 14:   VecScatter *iscatter;           /* scatter from global space to the nonoverlapping subdomain space */
 15:   VecScatter *gscatter;           /* scatter from global space to the subdomain local space */
 16:   PCASMType  type;                /* ASM type */
 17:   PetscBool  usesdm;              /* use the DM for setting up the subproblems */
 18:   PetscBool  finaljacobian;       /* compute the jacobian of the converged solution */
 19:   PetscReal  damping;             /* damping parameter for updates from the blocks */
 20:   PetscBool  weight_set;          /* use a weight in the overlap updates */

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

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

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

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

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

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

 51:   PetscFree(nasm->x);
 52:   PetscFree(nasm->xl);
 53:   PetscFree(nasm->y);
 54:   PetscFree(nasm->b);

 56:   if (nasm->xinit) VecDestroy(&nasm->xinit);

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

 64:   if (nasm->weight_set) {
 65:     VecDestroy(&nasm->weight);
 66:   }

 68:   nasm->eventrestrictinterp = 0;
 69:   nasm->eventsubsolve = 0;
 70:   return 0;
 71: }

 73: static PetscErrorCode SNESDestroy_NASM(SNES snes)
 74: {
 75:   SNESReset_NASM(snes);
 76:   PetscFree(snes->data);
 77:   return 0;
 78: }

 80: static PetscErrorCode DMGlobalToLocalSubDomainDirichletHook_Private(DM dm,Vec g,InsertMode mode,Vec l,void *ctx)
 81: {
 82:   Vec            bcs = (Vec)ctx;

 84:   VecCopy(bcs,l);
 85:   return 0;
 86: }

 88: static PetscErrorCode SNESSetUp_NASM(SNES snes)
 89: {
 90:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
 91:   DM             dm,subdm;
 92:   DM             *subdms;
 93:   PetscInt       i;
 94:   const char     *optionsprefix;
 95:   Vec            F;
 96:   PetscMPIInt    size;
 97:   KSP            ksp;
 98:   PC             pc;

100:   if (!nasm->subsnes) {
101:     SNESGetDM(snes,&dm);
102:     if (dm) {
103:       nasm->usesdm = PETSC_TRUE;
104:       DMCreateDomainDecomposition(dm,&nasm->n,NULL,NULL,NULL,&subdms);
106:       DMCreateDomainDecompositionScatters(dm,nasm->n,subdms,&nasm->iscatter,&nasm->oscatter,&nasm->gscatter);
107:       PetscMalloc1(nasm->n, &nasm->oscatter_copy);
108:       for (i=0; i<nasm->n; i++) {
109:         VecScatterCopy(nasm->oscatter[i], &nasm->oscatter_copy[i]);
110:       }

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

147:   for (i=0; i<nasm->n; i++) {
148:     SNESGetFunction(nasm->subsnes[i],&F,NULL,NULL);
149:     if (!nasm->x[i]) VecDuplicate(F,&nasm->x[i]);
150:     if (!nasm->y[i]) VecDuplicate(F,&nasm->y[i]);
151:     if (!nasm->b[i]) VecDuplicate(F,&nasm->b[i]);
152:     if (!nasm->xl[i]) {
153:       SNESGetDM(nasm->subsnes[i],&subdm);
154:       DMCreateLocalVector(subdm,&nasm->xl[i]);
155:       DMGlobalToLocalHookAdd(subdm,DMGlobalToLocalSubDomainDirichletHook_Private,NULL,nasm->xl[i]);
156:     }
157:   }
158:   if (nasm->finaljacobian) {
159:     SNESSetUpMatrices(snes);
160:     if (nasm->fjtype == 2) {
161:       VecDuplicate(snes->vec_sol,&nasm->xinit);
162:     }
163:     for (i=0; i<nasm->n;i++) {
164:       SNESSetUpMatrices(nasm->subsnes[i]);
165:     }
166:   }
167:   return 0;
168: }

170: static PetscErrorCode SNESSetFromOptions_NASM(PetscOptionItems *PetscOptionsObject,SNES snes)
171: {
172:   PCASMType         asmtype;
173:   PetscBool         flg,monflg;
174:   SNES_NASM         *nasm = (SNES_NASM*)snes->data;

176:   PetscOptionsHead(PetscOptionsObject,"Nonlinear Additive Schwarz options");
177:   PetscOptionsEnum("-snes_nasm_type","Type of restriction/extension","",SNESNASMTypes,(PetscEnum)nasm->type,(PetscEnum*)&asmtype,&flg);
178:   if (flg) SNESNASMSetType(snes,asmtype);
179:   flg    = PETSC_FALSE;
180:   monflg = PETSC_TRUE;
181:   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);
182:   if (flg) SNESNASMSetDamping(snes,nasm->damping);
183:   PetscOptionsDeprecated("-snes_nasm_sub_view",NULL,"3.15","Use -snes_view ::ascii_info_detail");
184:   PetscOptionsBool("-snes_nasm_finaljacobian","Compute the global jacobian of the final iterate (for ASPIN)","",nasm->finaljacobian,&nasm->finaljacobian,NULL);
185:   PetscOptionsEList("-snes_nasm_finaljacobian_type","The type of the final jacobian computed.","",SNESNASMFJTypes,3,SNESNASMFJTypes[0],&nasm->fjtype,NULL);
186:   PetscOptionsBool("-snes_nasm_log","Log times for subSNES solves and restriction","",monflg,&monflg,&flg);
187:   if (flg) {
188:     PetscLogEventRegister("SNESNASMSubSolve",((PetscObject)snes)->classid,&nasm->eventsubsolve);
189:     PetscLogEventRegister("SNESNASMRestrict",((PetscObject)snes)->classid,&nasm->eventrestrictinterp);
190:   }
191:   PetscOptionsTail();
192:   return 0;
193: }

195: static PetscErrorCode SNESView_NASM(SNES snes, PetscViewer viewer)
196: {
197:   SNES_NASM         *nasm = (SNES_NASM*)snes->data;
198:   PetscMPIInt       rank,size;
199:   PetscInt          i,N,bsz;
200:   PetscBool         iascii,isstring;
201:   PetscViewer       sviewer;
202:   MPI_Comm          comm;
203:   PetscViewerFormat format;
204:   const char        *prefix;

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

259: /*@
260:    SNESNASMSetType - Set the type of subdomain update used

262:    Logically Collective on SNES

264:    Input Parameters:
265: +  SNES - the SNES context
266: -  type - the type of update, PC_ASM_BASIC or PC_ASM_RESTRICT

268:    Level: intermediate

270: .seealso: SNESNASM, SNESNASMGetType(), PCASMSetType()
271: @*/
272: PetscErrorCode SNESNASMSetType(SNES snes,PCASMType type)
273: {
274:   PetscErrorCode (*f)(SNES,PCASMType);

276:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetType_C",&f);
277:   if (f) (f)(snes,type);
278:   return 0;
279: }

281: static PetscErrorCode SNESNASMSetType_NASM(SNES snes,PCASMType type)
282: {
283:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

286:   nasm->type = type;
287:   return 0;
288: }

290: /*@
291:    SNESNASMGetType - Get the type of subdomain update used

293:    Logically Collective on SNES

295:    Input Parameters:
296: .  SNES - the SNES context

298:    Output Parameters:
299: .  type - the type of update

301:    Level: intermediate

303: .seealso: SNESNASM, SNESNASMSetType(), PCASMGetType()
304: @*/
305: PetscErrorCode SNESNASMGetType(SNES snes,PCASMType *type)
306: {
307:   PetscUseMethod(snes,"SNESNASMGetType_C",(SNES,PCASMType*),(snes,type));
308:   return 0;
309: }

311: static PetscErrorCode SNESNASMGetType_NASM(SNES snes,PCASMType *type)
312: {
313:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

315:   *type = nasm->type;
316:   return 0;
317: }

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

322:    Not Collective

324:    Input Parameters:
325: +  SNES - the SNES context
326: .  n - the number of local subdomains
327: .  subsnes - solvers defined on the local subdomains
328: .  iscatter - scatters into the nonoverlapping portions of the local subdomains
329: .  oscatter - scatters into the overlapping portions of the local subdomains
330: -  gscatter - scatters into the (ghosted) local vector of the local subdomain

332:    Level: intermediate

334: .seealso: SNESNASM, SNESNASMGetSubdomains()
335: @*/
336: PetscErrorCode SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
337: {
338:   PetscErrorCode (*f)(SNES,PetscInt,SNES*,VecScatter*,VecScatter*,VecScatter*);

340:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",&f);
341:   if (f) (f)(snes,n,subsnes,iscatter,oscatter,gscatter);
342:   return 0;
343: }

345: static PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
346: {
347:   PetscInt       i;
348:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;


352:   /* tear down the previously set things */
353:   SNESReset(snes);

355:   nasm->n = n;
356:   if (oscatter) {
357:     for (i=0; i<n; i++) PetscObjectReference((PetscObject)oscatter[i]);
358:   }
359:   if (iscatter) {
360:     for (i=0; i<n; i++) PetscObjectReference((PetscObject)iscatter[i]);
361:   }
362:   if (gscatter) {
363:     for (i=0; i<n; i++) PetscObjectReference((PetscObject)gscatter[i]);
364:   }
365:   if (oscatter) {
366:     PetscMalloc1(n,&nasm->oscatter);
367:     PetscMalloc1(n,&nasm->oscatter_copy);
368:     for (i=0; i<n; i++) {
369:       nasm->oscatter[i] = oscatter[i];
370:       VecScatterCopy(oscatter[i], &nasm->oscatter_copy[i]);
371:     }
372:   }
373:   if (iscatter) {
374:     PetscMalloc1(n,&nasm->iscatter);
375:     for (i=0; i<n; i++) {
376:       nasm->iscatter[i] = iscatter[i];
377:     }
378:   }
379:   if (gscatter) {
380:     PetscMalloc1(n,&nasm->gscatter);
381:     for (i=0; i<n; i++) {
382:       nasm->gscatter[i] = gscatter[i];
383:     }
384:   }

386:   if (subsnes) {
387:     PetscMalloc1(n,&nasm->subsnes);
388:     for (i=0; i<n; i++) {
389:       nasm->subsnes[i] = subsnes[i];
390:     }
391:   }
392:   return 0;
393: }

395: /*@
396:    SNESNASMGetSubdomains - Get the local subdomain context.

398:    Not Collective

400:    Input Parameter:
401: .  SNES - the SNES context

403:    Output Parameters:
404: +  n - the number of local subdomains
405: .  subsnes - solvers defined on the local subdomains
406: .  iscatter - scatters into the nonoverlapping portions of the local subdomains
407: .  oscatter - scatters into the overlapping portions of the local subdomains
408: -  gscatter - scatters into the (ghosted) local vector of the local subdomain

410:    Level: intermediate

412: .seealso: SNESNASM, SNESNASMSetSubdomains()
413: @*/
414: PetscErrorCode SNESNASMGetSubdomains(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
415: {
416:   PetscErrorCode (*f)(SNES,PetscInt*,SNES**,VecScatter**,VecScatter**,VecScatter**);

418:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",&f);
419:   if (f) (f)(snes,n,subsnes,iscatter,oscatter,gscatter);
420:   return 0;
421: }

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

427:   if (n) *n = nasm->n;
428:   if (oscatter) *oscatter = nasm->oscatter;
429:   if (iscatter) *iscatter = nasm->iscatter;
430:   if (gscatter) *gscatter = nasm->gscatter;
431:   if (subsnes)  *subsnes  = nasm->subsnes;
432:   return 0;
433: }

435: /*@
436:    SNESNASMGetSubdomainVecs - Get the processor-local subdomain vectors

438:    Not Collective

440:    Input Parameter:
441: .  SNES - the SNES context

443:    Output Parameters:
444: +  n - the number of local subdomains
445: .  x - The subdomain solution vector
446: .  y - The subdomain step vector
447: .  b - The subdomain RHS vector
448: -  xl - The subdomain local vectors (ghosted)

450:    Level: developer

452: .seealso: SNESNASM, SNESNASMGetSubdomains()
453: @*/
454: PetscErrorCode SNESNASMGetSubdomainVecs(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b, Vec **xl)
455: {
456:   PetscErrorCode (*f)(SNES,PetscInt*,Vec**,Vec**,Vec**,Vec**);

458:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",&f);
459:   if (f) (f)(snes,n,x,y,b,xl);
460:   return 0;
461: }

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

467:   if (n)  *n  = nasm->n;
468:   if (x)  *x  = nasm->x;
469:   if (y)  *y  = nasm->y;
470:   if (b)  *b  = nasm->b;
471:   if (xl) *xl = nasm->xl;
472:   return 0;
473: }

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

478:    Collective on SNES

480:    Input Parameters:
481: +  SNES - the SNES context
482: -  flg - indication of whether to compute the Jacobians or not

484:    Level: developer

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

490: .seealso: SNESNASM, SNESNASMGetSubdomains()
491: @*/
492: PetscErrorCode SNESNASMSetComputeFinalJacobian(SNES snes,PetscBool flg)
493: {
494:   PetscErrorCode (*f)(SNES,PetscBool);

496:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",&f);
497:   if (f) (f)(snes,flg);
498:   return 0;
499: }

501: static PetscErrorCode SNESNASMSetComputeFinalJacobian_NASM(SNES snes,PetscBool flg)
502: {
503:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

505:   nasm->finaljacobian = flg;
506:   return 0;
507: }

509: /*@
510:    SNESNASMSetDamping - Sets the update damping for NASM

512:    Logically collective on SNES

514:    Input Parameters:
515: +  SNES - the SNES context
516: -  dmp - damping

518:    Level: intermediate

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

523: .seealso: SNESNASM, SNESNASMGetDamping()
524: @*/
525: PetscErrorCode SNESNASMSetDamping(SNES snes,PetscReal dmp)
526: {
527:   PetscErrorCode (*f)(SNES,PetscReal);

529:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetDamping_C",(void (**)(void))&f);
530:   if (f) (f)(snes,dmp);
531:   return 0;
532: }

534: static PetscErrorCode SNESNASMSetDamping_NASM(SNES snes,PetscReal dmp)
535: {
536:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

538:   nasm->damping = dmp;
539:   return 0;
540: }

542: /*@
543:    SNESNASMGetDamping - Gets the update damping for NASM

545:    Not Collective

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

551:    Level: intermediate

553: .seealso: SNESNASM, SNESNASMSetDamping()
554: @*/
555: PetscErrorCode SNESNASMGetDamping(SNES snes,PetscReal *dmp)
556: {
557:   PetscUseMethod(snes,"SNESNASMGetDamping_C",(SNES,PetscReal*),(snes,dmp));
558:   return 0;
559: }

561: static PetscErrorCode SNESNASMGetDamping_NASM(SNES snes,PetscReal *dmp)
562: {
563:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

565:   *dmp = nasm->damping;
566:   return 0;
567: }

569: /*
570:   Input Parameters:
571: + snes - The solver
572: . B - The RHS vector
573: - X - The initial guess

575:   Output Parameters:
576: . Y - The solution update

578:   TODO: All scatters should be packed into one
579: */
580: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X)
581: {
582:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
583:   SNES           subsnes;
584:   PetscInt       i;
585:   PetscReal      dmp;
586:   Vec            Xl,Bl,Yl,Xlloc;
587:   VecScatter     iscat,oscat,gscat,oscat_copy;
588:   DM             dm,subdm;
589:   PCASMType      type;

591:   SNESNASMGetType(snes,&type);
592:   SNESGetDM(snes,&dm);
593:   VecSet(Y,0);
594:   if (nasm->eventrestrictinterp) PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);
595:   for (i=0; i<nasm->n; i++) {
596:     /* scatter the solution to the global solution and the local solution */
597:     Xl      = nasm->x[i];
598:     Xlloc   = nasm->xl[i];
599:     oscat   = nasm->oscatter[i];
600:     oscat_copy = nasm->oscatter_copy[i];
601:     gscat   = nasm->gscatter[i];
602:     VecScatterBegin(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
603:     VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
604:     if (B) {
605:       /* scatter the RHS to the local RHS */
606:       Bl   = nasm->b[i];
607:       VecScatterBegin(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
608:     }
609:   }
610:   if (nasm->eventrestrictinterp) PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);

612:   if (nasm->eventsubsolve) PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);
613:   for (i=0; i<nasm->n; i++) {
614:     Xl    = nasm->x[i];
615:     Xlloc = nasm->xl[i];
616:     Yl    = nasm->y[i];
617:     subsnes = nasm->subsnes[i];
618:     SNESGetDM(subsnes,&subdm);
619:     iscat   = nasm->iscatter[i];
620:     oscat   = nasm->oscatter[i];
621:     oscat_copy = nasm->oscatter_copy[i];
622:     gscat   = nasm->gscatter[i];
623:     VecScatterEnd(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
624:     VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
625:     if (B) {
626:       Bl   = nasm->b[i];
627:       VecScatterEnd(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
628:     } else Bl = NULL;

630:     DMSubDomainRestrict(dm,oscat,gscat,subdm);
631:     VecCopy(Xl,Yl);
632:     SNESSolve(subsnes,Bl,Xl);
633:     VecAYPX(Yl,-1.0,Xl);
634:     VecScale(Yl, nasm->damping);
635:     if (type == PC_ASM_BASIC) {
636:       VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
637:       VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
638:     } else if (type == PC_ASM_RESTRICT) {
639:       VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
640:       VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
641:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
642:   }
643:   if (nasm->eventsubsolve) PetscLogEventEnd(nasm->eventsubsolve,snes,0,0,0);
644:   if (nasm->eventrestrictinterp) PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);
645:   if (nasm->weight_set) {
646:     VecPointwiseMult(Y,Y,nasm->weight);
647:   }
648:   if (nasm->eventrestrictinterp) PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);
649:   SNESNASMGetDamping(snes,&dmp);
650:   VecAXPY(X,dmp,Y);
651:   return 0;
652: }

654: static PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
655: {
656:   Vec            X = Xfinal;
657:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
658:   SNES           subsnes;
659:   PetscInt       i,lag = 1;
660:   Vec            Xlloc,Xl,Fl,F;
661:   VecScatter     oscat,gscat;
662:   DM             dm,subdm;

664:   if (nasm->fjtype == 2) X = nasm->xinit;
665:   F = snes->vec_func;
666:   if (snes->normschedule == SNES_NORM_NONE) SNESComputeFunction(snes,X,F);
667:   SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
668:   SNESGetDM(snes,&dm);
669:   if (nasm->eventrestrictinterp) PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);
670:   if (nasm->fjtype != 1) {
671:     for (i=0; i<nasm->n; i++) {
672:       Xlloc = nasm->xl[i];
673:       gscat = nasm->gscatter[i];
674:       VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
675:     }
676:   }
677:   if (nasm->eventrestrictinterp) PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);
678:   for (i=0; i<nasm->n; i++) {
679:     Fl      = nasm->subsnes[i]->vec_func;
680:     Xl      = nasm->x[i];
681:     Xlloc   = nasm->xl[i];
682:     subsnes = nasm->subsnes[i];
683:     oscat   = nasm->oscatter[i];
684:     gscat   = nasm->gscatter[i];
685:     if (nasm->fjtype != 1) VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
686:     SNESGetDM(subsnes,&subdm);
687:     DMSubDomainRestrict(dm,oscat,gscat,subdm);
688:     if (nasm->fjtype != 1) {
689:       DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
690:       DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
691:     }
692:     if (subsnes->lagjacobian == -1)    subsnes->lagjacobian = -2;
693:     else if (subsnes->lagjacobian > 1) lag = subsnes->lagjacobian;
694:     SNESComputeFunction(subsnes,Xl,Fl);
695:     SNESComputeJacobian(subsnes,Xl,subsnes->jacobian,subsnes->jacobian_pre);
696:     if (lag > 1) subsnes->lagjacobian = lag;
697:   }
698:   return 0;
699: }

701: static PetscErrorCode SNESSolve_NASM(SNES snes)
702: {
703:   Vec              F;
704:   Vec              X;
705:   Vec              B;
706:   Vec              Y;
707:   PetscInt         i;
708:   PetscReal        fnorm = 0.0;
709:   SNESNormSchedule normschedule;
710:   SNES_NASM        *nasm = (SNES_NASM*)snes->data;



715:   PetscCitationsRegister(SNESCitation,&SNEScite);
716:   X = snes->vec_sol;
717:   Y = snes->vec_sol_update;
718:   F = snes->vec_func;
719:   B = snes->vec_rhs;

721:   PetscObjectSAWsTakeAccess((PetscObject)snes);
722:   snes->iter   = 0;
723:   snes->norm   = 0.;
724:   PetscObjectSAWsGrantAccess((PetscObject)snes);
725:   snes->reason = SNES_CONVERGED_ITERATING;
726:   SNESGetNormSchedule(snes, &normschedule);
727:   if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
728:     /* compute the initial function and preconditioned update delX */
729:     if (!snes->vec_func_init_set) {
730:       SNESComputeFunction(snes,X,F);
731:     } else snes->vec_func_init_set = PETSC_FALSE;

733:     VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
734:     SNESCheckFunctionNorm(snes,fnorm);
735:     PetscObjectSAWsTakeAccess((PetscObject)snes);
736:     snes->iter = 0;
737:     snes->norm = fnorm;
738:     PetscObjectSAWsGrantAccess((PetscObject)snes);
739:     SNESLogConvergenceHistory(snes,snes->norm,0);
740:     SNESMonitor(snes,0,snes->norm);

742:     /* test convergence */
743:     (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
744:     if (snes->reason) return 0;
745:   } else {
746:     PetscObjectSAWsGrantAccess((PetscObject)snes);
747:     SNESLogConvergenceHistory(snes,snes->norm,0);
748:     SNESMonitor(snes,0,snes->norm);
749:   }

751:   /* Call general purpose update function */
752:   if (snes->ops->update) {
753:     (*snes->ops->update)(snes, snes->iter);
754:   }
755:   /* copy the initial solution over for later */
756:   if (nasm->fjtype == 2) VecCopy(X,nasm->xinit);

758:   for (i=0; i < snes->max_its; i++) {
759:     SNESNASMSolveLocal_Private(snes,B,Y,X);
760:     if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
761:       SNESComputeFunction(snes,X,F);
762:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
763:       SNESCheckFunctionNorm(snes,fnorm);
764:     }
765:     /* Monitor convergence */
766:     PetscObjectSAWsTakeAccess((PetscObject)snes);
767:     snes->iter = i+1;
768:     snes->norm = fnorm;
769:     PetscObjectSAWsGrantAccess((PetscObject)snes);
770:     SNESLogConvergenceHistory(snes,snes->norm,0);
771:     SNESMonitor(snes,snes->iter,snes->norm);
772:     /* Test for convergence */
773:     if (normschedule == SNES_NORM_ALWAYS) (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
774:     if (snes->reason) break;
775:     /* Call general purpose update function */
776:     if (snes->ops->update) (*snes->ops->update)(snes, snes->iter);
777:   }
778:   if (nasm->finaljacobian) {
779:     SNESNASMComputeFinalJacobian_Private(snes,X);
780:     SNESCheckJacobianDomainerror(snes);
781:   }
782:   if (normschedule == SNES_NORM_ALWAYS) {
783:     if (i == snes->max_its) {
784:       PetscInfo(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
785:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
786:     }
787:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */
788:   return 0;
789: }

791: /*MC
792:   SNESNASM - Nonlinear Additive Schwarz

794:    Options Database:
795: +  -snes_nasm_log - enable logging events for the communication and solve stages
796: .  -snes_nasm_type <basic,restrict> - type of subdomain update used
797: .  -snes_asm_damping <dmp> - the new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
798: .  -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
799: .  -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> - pick state the jacobian is calculated at
800: .  -sub_snes_ - options prefix of the subdomain nonlinear solves
801: .  -sub_ksp_ - options prefix of the subdomain Krylov solver
802: -  -sub_pc_ - options prefix of the subdomain preconditioner

804:    Level: advanced

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

813:    References:
814: .  * - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
815:    SIAM Review, 57(4), 2015

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

820: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
821: {
822:   SNES_NASM      *nasm;

824:   PetscNewLog(snes,&nasm);
825:   snes->data = (void*)nasm;

827:   nasm->n        = PETSC_DECIDE;
828:   nasm->subsnes  = NULL;
829:   nasm->x        = NULL;
830:   nasm->xl       = NULL;
831:   nasm->y        = NULL;
832:   nasm->b        = NULL;
833:   nasm->oscatter = NULL;
834:   nasm->oscatter_copy = NULL;
835:   nasm->iscatter = NULL;
836:   nasm->gscatter = NULL;
837:   nasm->damping  = 1.;

839:   nasm->type              = PC_ASM_BASIC;
840:   nasm->finaljacobian     = PETSC_FALSE;
841:   nasm->weight_set        = PETSC_FALSE;

843:   snes->ops->destroy        = SNESDestroy_NASM;
844:   snes->ops->setup          = SNESSetUp_NASM;
845:   snes->ops->setfromoptions = SNESSetFromOptions_NASM;
846:   snes->ops->view           = SNESView_NASM;
847:   snes->ops->solve          = SNESSolve_NASM;
848:   snes->ops->reset          = SNESReset_NASM;

850:   snes->usesksp = PETSC_FALSE;
851:   snes->usesnpc = PETSC_FALSE;

853:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

855:   nasm->fjtype              = 0;
856:   nasm->xinit               = NULL;
857:   nasm->eventrestrictinterp = 0;
858:   nasm->eventsubsolve       = 0;

860:   if (!snes->tolerancesset) {
861:     snes->max_its   = 10000;
862:     snes->max_funcs = 10000;
863:   }

865:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetType_C",SNESNASMSetType_NASM);
866:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetType_C",SNESNASMGetType_NASM);
867:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
868:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
869:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
870:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
871:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
872:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
873:   return 0;
874: }

876: /*@
877:    SNESNASMGetSNES - Gets a subsolver

879:    Not collective

881:    Input Parameters:
882: +  snes - the SNES context
883: -  i - the number of the subsnes to get

885:    Output Parameters:
886: .  subsnes - the subsolver context

888:    Level: intermediate

890: .seealso: SNESNASM, SNESNASMGetNumber()
891: @*/
892: PetscErrorCode SNESNASMGetSNES(SNES snes,PetscInt i,SNES *subsnes)
893: {
894:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

897:   *subsnes = nasm->subsnes[i];
898:   return 0;
899: }

901: /*@
902:    SNESNASMGetNumber - Gets number of subsolvers

904:    Not collective

906:    Input Parameters:
907: .  snes - the SNES context

909:    Output Parameters:
910: .  n - the number of subsolvers

912:    Level: intermediate

914: .seealso: SNESNASM, SNESNASMGetSNES()
915: @*/
916: PetscErrorCode SNESNASMGetNumber(SNES snes,PetscInt *n)
917: {
918:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

920:   *n = nasm->n;
921:   return 0;
922: }

924: /*@
925:    SNESNASMSetWeight - Sets weight to use when adding overlapping updates

927:    Collective

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

933:    Level: intermediate

935: .seealso: SNESNASM
936: @*/
937: PetscErrorCode SNESNASMSetWeight(SNES snes,Vec weight)
938: {
939:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;


942:   VecDestroy(&nasm->weight);
943:   nasm->weight_set = PETSC_TRUE;
944:   nasm->weight     = weight;
945:   PetscObjectReference((PetscObject)nasm->weight);

947:   return 0;
948: }