Actual source code: nasm.c

petsc-3.7.3 2016-07-24
Report Typos and Errors
  1: #include <petsc/private/snesimpl.h>             /*I   "petscsnes.h"   I*/
  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:   VecScatter *oscatter;           /* scatter from global space to the subdomain global space */
 12:   VecScatter *iscatter;           /* scatter from global space to the nonoverlapping subdomain space */
 13:   VecScatter *gscatter;           /* scatter from global space to the subdomain local space */
 14:   PCASMType  type;                /* ASM type */
 15:   PetscBool  usesdm;              /* use the DM for setting up the subproblems */
 16:   PetscBool  finaljacobian;       /* compute the jacobian of the converged solution */
 17:   PetscReal  damping;             /* damping parameter for updates from the blocks */
 18:   PetscBool  same_local_solves;   /* flag to determine if the solvers have been individually modified */

 20:   /* logging events */
 21:   PetscLogEvent eventrestrictinterp;
 22:   PetscLogEvent eventsubsolve;

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

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

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

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

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

 52:   if (nasm->x) {PetscFree(nasm->x);}
 53:   if (nasm->xl) {PetscFree(nasm->xl);}
 54:   if (nasm->y) {PetscFree(nasm->y);}
 55:   if (nasm->b) {PetscFree(nasm->b);}

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

 59:   if (nasm->subsnes) {PetscFree(nasm->subsnes);}
 60:   if (nasm->oscatter) {PetscFree(nasm->oscatter);}
 61:   if (nasm->iscatter) {PetscFree(nasm->iscatter);}
 62:   if (nasm->gscatter) {PetscFree(nasm->gscatter);}

 64:   nasm->eventrestrictinterp = 0;
 65:   nasm->eventsubsolve = 0;
 66:   return(0);
 67: }

 71: PetscErrorCode SNESDestroy_NASM(SNES snes)
 72: {

 76:   SNESReset_NASM(snes);
 77:   PetscFree(snes->data);
 78:   return(0);
 79: }

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

 89:   VecCopy(bcs,l);
 90:   return(0);
 91: }

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

109:   if (!nasm->subsnes) {
110:     SNESGetDM(snes,&dm);
111:     if (dm) {
112:       nasm->usesdm = PETSC_TRUE;
113:       DMCreateDomainDecomposition(dm,&nasm->n,NULL,NULL,NULL,&subdms);
114:       if (!subdms) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM has no default decomposition defined.  Set subsolves manually with SNESNASMSetSubdomains().");
115:       DMCreateDomainDecompositionScatters(dm,nasm->n,subdms,&nasm->iscatter,&nasm->oscatter,&nasm->gscatter);

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

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

176: PetscErrorCode SNESSetFromOptions_NASM(PetscOptionItems *PetscOptionsObject,SNES snes)
177: {
178:   PetscErrorCode    ierr;
179:   PCASMType         asmtype;
180:   PetscBool         flg,monflg,subviewflg;
181:   SNES_NASM         *nasm = (SNES_NASM*)snes->data;

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

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

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

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: }

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

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

315: /*@
316:    SNESNASMGetType - Get the type of subdomain update used

318:    Logically Collective on SNES

320:    Input Parameters:
321: .  SNES - the SNES context

323:    Output Parameters:
324: .  type - the type of update

326:    Level: intermediate

328: .keywords: SNES, NASM

330: .seealso: SNESNASM, SNESNASMSetType(), PCASMGetType()
331: @*/
332: PetscErrorCode SNESNASMGetType(SNES snes,PCASMType *type)
333: {
335:   PetscErrorCode (*f)(SNES,PCASMType*);

338:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetType_C",&f);
339:   if (f) {(f)(snes,type);}
340:   return(0);
341: }

345: PetscErrorCode SNESNASMGetType_NASM(SNES snes,PCASMType *type)
346: {
347:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

350:   *type = nasm->type;
351:   return(0);
352: }

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

359:    Not Collective

361:    Input Parameters:
362: +  SNES - the SNES context
363: .  n - the number of local subdomains
364: .  subsnes - solvers defined on the local subdomains
365: .  iscatter - scatters into the nonoverlapping portions of the local subdomains
366: .  oscatter - scatters into the overlapping portions of the local subdomains
367: -  gscatter - scatters into the (ghosted) local vector of the local subdomain

369:    Level: intermediate

371: .keywords: SNES, NASM

373: .seealso: SNESNASM, SNESNASMGetSubdomains()
374: @*/
375: PetscErrorCode SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
376: {
378:   PetscErrorCode (*f)(SNES,PetscInt,SNES*,VecScatter*,VecScatter*,VecScatter*);

381:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",&f);
382:   if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
383:   return(0);
384: }

388: PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
389: {
390:   PetscInt       i;
392:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

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

397:   /* tear down the previously set things */
398:   SNESReset(snes);

400:   nasm->n = n;
401:   if (oscatter) {
402:     for (i=0; i<n; i++) {PetscObjectReference((PetscObject)oscatter[i]);}
403:   }
404:   if (iscatter) {
405:     for (i=0; i<n; i++) {PetscObjectReference((PetscObject)iscatter[i]);}
406:   }
407:   if (gscatter) {
408:     for (i=0; i<n; i++) {PetscObjectReference((PetscObject)gscatter[i]);}
409:   }
410:   if (oscatter) {
411:     PetscMalloc1(n,&nasm->oscatter);
412:     for (i=0; i<n; i++) {
413:       nasm->oscatter[i] = oscatter[i];
414:     }
415:   }
416:   if (iscatter) {
417:     PetscMalloc1(n,&nasm->iscatter);
418:     for (i=0; i<n; i++) {
419:       nasm->iscatter[i] = iscatter[i];
420:     }
421:   }
422:   if (gscatter) {
423:     PetscMalloc1(n,&nasm->gscatter);
424:     for (i=0; i<n; i++) {
425:       nasm->gscatter[i] = gscatter[i];
426:     }
427:   }

429:   if (subsnes) {
430:     PetscMalloc1(n,&nasm->subsnes);
431:     for (i=0; i<n; i++) {
432:       nasm->subsnes[i] = subsnes[i];
433:     }
434:     nasm->same_local_solves = PETSC_FALSE;
435:   }
436:   return(0);
437: }

441: /*@
442:    SNESNASMGetSubdomains - Get the local subdomain context.

444:    Not Collective

446:    Input Parameters:
447: .  SNES - the SNES context

449:    Output Parameters:
450: +  n - the number of local subdomains
451: .  subsnes - solvers defined on the local subdomains
452: .  iscatter - scatters into the nonoverlapping portions of the local subdomains
453: .  oscatter - scatters into the overlapping portions of the local subdomains
454: -  gscatter - scatters into the (ghosted) local vector of the local subdomain

456:    Level: intermediate

458: .keywords: SNES, NASM

460: .seealso: SNESNASM, SNESNASMSetSubdomains()
461: @*/
462: PetscErrorCode SNESNASMGetSubdomains(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
463: {
465:   PetscErrorCode (*f)(SNES,PetscInt*,SNES**,VecScatter**,VecScatter**,VecScatter**);

468:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",&f);
469:   if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
470:   return(0);
471: }

475: PetscErrorCode SNESNASMGetSubdomains_NASM(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
476: {
477:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

480:   if (n) *n = nasm->n;
481:   if (oscatter) *oscatter = nasm->oscatter;
482:   if (iscatter) *iscatter = nasm->iscatter;
483:   if (gscatter) *gscatter = nasm->gscatter;
484:   if (subsnes)  {
485:     *subsnes  = nasm->subsnes;
486:     nasm->same_local_solves = PETSC_FALSE;
487:   }
488:   return(0);
489: }

493: /*@
494:    SNESNASMGetSubdomainVecs - Get the processor-local subdomain vectors

496:    Not Collective

498:    Input Parameters:
499: .  SNES - the SNES context

501:    Output Parameters:
502: +  n - the number of local subdomains
503: .  x - The subdomain solution vector
504: .  y - The subdomain step vector
505: .  b - The subdomain RHS vector
506: -  xl - The subdomain local vectors (ghosted)

508:    Level: developer

510: .keywords: SNES, NASM

512: .seealso: SNESNASM, SNESNASMGetSubdomains()
513: @*/
514: PetscErrorCode SNESNASMGetSubdomainVecs(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b, Vec **xl)
515: {
517:   PetscErrorCode (*f)(SNES,PetscInt*,Vec**,Vec**,Vec**,Vec**);

520:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",&f);
521:   if (f) {(f)(snes,n,x,y,b,xl);}
522:   return(0);
523: }

527: PetscErrorCode SNESNASMGetSubdomainVecs_NASM(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b,Vec **xl)
528: {
529:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

532:   if (n)  *n  = nasm->n;
533:   if (x)  *x  = nasm->x;
534:   if (y)  *y  = nasm->y;
535:   if (b)  *b  = nasm->b;
536:   if (xl) *xl = nasm->xl;
537:   return(0);
538: }

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

545:    Collective on SNES

547:    Input Parameters:
548: +  SNES - the SNES context
549: -  flg - indication of whether to compute the jacobians or not

551:    Level: developer

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

556: .keywords: SNES, NASM, ASPIN

558: .seealso: SNESNASM, SNESNASMGetSubdomains()
559: @*/
560: PetscErrorCode SNESNASMSetComputeFinalJacobian(SNES snes,PetscBool flg)
561: {
562:   PetscErrorCode (*f)(SNES,PetscBool);

566:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",&f);
567:   if (f) {(f)(snes,flg);}
568:   return(0);
569: }

573: PetscErrorCode SNESNASMSetComputeFinalJacobian_NASM(SNES snes,PetscBool flg)
574: {
575:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

578:   nasm->finaljacobian = flg;
579:   if (flg) snes->usesksp = PETSC_TRUE;
580:   return(0);
581: }

585: /*@
586:    SNESNASMSetDamping - Sets the update damping for NASM

588:    Logically collective on SNES

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

594:    Level: intermediate

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

598: .keywords: SNES, NASM, damping

600: .seealso: SNESNASM, SNESNASMGetDamping()
601: @*/
602: PetscErrorCode SNESNASMSetDamping(SNES snes,PetscReal dmp)
603: {
604:   PetscErrorCode (*f)(SNES,PetscReal);

608:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetDamping_C",(void (**)(void))&f);
609:   if (f) {(f)(snes,dmp);}
610:   return(0);
611: }

615: PetscErrorCode SNESNASMSetDamping_NASM(SNES snes,PetscReal dmp)
616: {
617:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

620:   nasm->damping = dmp;
621:   return(0);
622: }

626: /*@
627:    SNESNASMGetDamping - Gets the update damping for NASM

629:    Not Collective

631:    Input Parameters:
632: +  SNES - the SNES context
633: -  dmp - damping

635:    Level: intermediate

637: .keywords: SNES, NASM, damping

639: .seealso: SNESNASM, SNESNASMSetDamping()
640: @*/
641: PetscErrorCode SNESNASMGetDamping(SNES snes,PetscReal *dmp)
642: {
643:   PetscErrorCode (*f)(SNES,PetscReal*);

647:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetDamping_C",(void (**)(void))&f);
648:   if (f) {(f)(snes,dmp);}
649:   return(0);
650: }

654: PetscErrorCode SNESNASMGetDamping_NASM(SNES snes,PetscReal *dmp)
655: {
656:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

659:   *dmp = nasm->damping;
660:   return(0);
661: }


666: /*
667:   Input Parameters:
668: + snes - The solver
669: . B - The RHS vector
670: - X - The initial guess

672:   Output Parameters:
673: . Y - The solution update

675:   TODO: All scatters should be packed into one
676: */
677: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X)
678: {
679:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
680:   SNES           subsnes;
681:   PetscInt       i;
682:   PetscReal      dmp;
684:   Vec            Xlloc,Xl,Bl,Yl;
685:   VecScatter     iscat,oscat,gscat;
686:   DM             dm,subdm;
687:   PCASMType      type;

690:   SNESNASMGetType(snes,&type);
691:   SNESGetDM(snes,&dm);
692:   VecSet(Y,0);
693:   if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
694:   for (i=0; i<nasm->n; i++) {
695:     /* scatter the solution to the local solution */
696:     Xlloc = nasm->xl[i];
697:     gscat   = nasm->gscatter[i];
698:     oscat   = nasm->oscatter[i];
699:     VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
700:     if (B) {
701:       /* scatter the RHS to the local RHS */
702:       Bl   = nasm->b[i];
703:       VecScatterBegin(oscat,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
704:     }
705:   }
706:   if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}


709:   if (nasm->eventsubsolve) {PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);}
710:   for (i=0; i<nasm->n; i++) {
711:     Xl    = nasm->x[i];
712:     Xlloc = nasm->xl[i];
713:     Yl    = nasm->y[i];
714:     subsnes = nasm->subsnes[i];
715:     SNESGetDM(subsnes,&subdm);
716:     iscat   = nasm->iscatter[i];
717:     oscat   = nasm->oscatter[i];
718:     gscat   = nasm->gscatter[i];
719:     VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
720:     if (B) {
721:       Bl   = nasm->b[i];
722:       VecScatterEnd(oscat,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
723:     } else Bl = NULL;
724:     DMSubDomainRestrict(dm,oscat,gscat,subdm);
725:     /* Could scatter directly from X */
726:     DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
727:     DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
728:     VecCopy(Xl,Yl);
729:     SNESSolve(subsnes,Bl,Xl);
730:     VecAYPX(Yl,-1.0,Xl);
731:     if (type == PC_ASM_BASIC) {
732:       VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
733:     } else if (type == PC_ASM_RESTRICT) {
734:       VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
735:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
736:   }
737:   if (nasm->eventsubsolve) {PetscLogEventEnd(nasm->eventsubsolve,snes,0,0,0);}
738:   if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
739:   for (i=0; i<nasm->n; i++) {
740:     Yl    = nasm->y[i];
741:     iscat   = nasm->iscatter[i];
742:     oscat   = nasm->oscatter[i];
743:     if (type == PC_ASM_BASIC) {
744:       VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
745:     } else if (type == PC_ASM_RESTRICT) {
746:       VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
747:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
748:   }
749:   if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
750:   SNESNASMGetDamping(snes,&dmp);
751:   VecAXPY(X,dmp,Y);
752:   return(0);
753: }

757: PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
758: {
759:   Vec            X = Xfinal;
760:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
761:   SNES           subsnes;
762:   PetscInt       i,lag = 1;
764:   Vec            Xlloc,Xl,Fl,F;
765:   VecScatter     oscat,gscat;
766:   DM             dm,subdm;

769:   if (nasm->fjtype == 2) X = nasm->xinit;
770:   F = snes->vec_func;
771:   if (snes->normschedule == SNES_NORM_NONE) {SNESComputeFunction(snes,X,F);}
772:   SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
773:   SNESGetDM(snes,&dm);
774:   if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
775:   if (nasm->fjtype != 1) {
776:     for (i=0; i<nasm->n; i++) {
777:       Xlloc = nasm->xl[i];
778:       gscat = nasm->gscatter[i];
779:       VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
780:     }
781:   }
782:   if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
783:   for (i=0; i<nasm->n; i++) {
784:     Fl      = nasm->subsnes[i]->vec_func;
785:     Xl      = nasm->x[i];
786:     Xlloc   = nasm->xl[i];
787:     subsnes = nasm->subsnes[i];
788:     oscat   = nasm->oscatter[i];
789:     gscat   = nasm->gscatter[i];
790:     if (nasm->fjtype != 1) {VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);}
791:     SNESGetDM(subsnes,&subdm);
792:     DMSubDomainRestrict(dm,oscat,gscat,subdm);
793:     if (nasm->fjtype != 1) {
794:       DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
795:       DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
796:     }
797:     if (subsnes->lagjacobian == -1)    subsnes->lagjacobian = -2;
798:     else if (subsnes->lagjacobian > 1) lag = subsnes->lagjacobian;
799:     SNESComputeFunction(subsnes,Xl,Fl);
800:     SNESComputeJacobian(subsnes,Xl,subsnes->jacobian,subsnes->jacobian_pre);
801:     if (lag > 1) subsnes->lagjacobian = lag;
802:   }
803:   return(0);
804: }

808: PetscErrorCode SNESSolve_NASM(SNES snes)
809: {
810:   Vec              F;
811:   Vec              X;
812:   Vec              B;
813:   Vec              Y;
814:   PetscInt         i;
815:   PetscReal        fnorm = 0.0;
816:   PetscErrorCode   ierr;
817:   SNESNormSchedule normschedule;
818:   SNES_NASM        *nasm = (SNES_NASM*)snes->data;


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

824:   PetscCitationsRegister(SNESCitation,&SNEScite);
825:   X = snes->vec_sol;
826:   Y = snes->vec_sol_update;
827:   F = snes->vec_func;
828:   B = snes->vec_rhs;

830:   PetscObjectSAWsTakeAccess((PetscObject)snes);
831:   snes->iter   = 0;
832:   snes->norm   = 0.;
833:   PetscObjectSAWsGrantAccess((PetscObject)snes);
834:   snes->reason = SNES_CONVERGED_ITERATING;
835:   SNESGetNormSchedule(snes, &normschedule);
836:   if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
837:     /* compute the initial function and preconditioned update delX */
838:     if (!snes->vec_func_init_set) {
839:       SNESComputeFunction(snes,X,F);
840:     } else snes->vec_func_init_set = PETSC_FALSE;

842:     VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
843:     SNESCheckFunctionNorm(snes,fnorm);
844:     PetscObjectSAWsTakeAccess((PetscObject)snes);
845:     snes->iter = 0;
846:     snes->norm = fnorm;
847:     PetscObjectSAWsGrantAccess((PetscObject)snes);
848:     SNESLogConvergenceHistory(snes,snes->norm,0);
849:     SNESMonitor(snes,0,snes->norm);

851:     /* test convergence */
852:     (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
853:     if (snes->reason) return(0);
854:   } else {
855:     PetscObjectSAWsGrantAccess((PetscObject)snes);
856:     SNESLogConvergenceHistory(snes,snes->norm,0);
857:     SNESMonitor(snes,0,snes->norm);
858:   }

860:   /* Call general purpose update function */
861:   if (snes->ops->update) {
862:     (*snes->ops->update)(snes, snes->iter);
863:   }
864:   /* copy the initial solution over for later */
865:   if (nasm->fjtype == 2) {VecCopy(X,nasm->xinit);}

867:   for (i = 0; i < snes->max_its; i++) {
868:     SNESNASMSolveLocal_Private(snes,B,Y,X);
869:     if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
870:       SNESComputeFunction(snes,X,F);
871:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
872:       SNESCheckFunctionNorm(snes,fnorm);
873:     }
874:     /* Monitor convergence */
875:     PetscObjectSAWsTakeAccess((PetscObject)snes);
876:     snes->iter = i+1;
877:     snes->norm = fnorm;
878:     PetscObjectSAWsGrantAccess((PetscObject)snes);
879:     SNESLogConvergenceHistory(snes,snes->norm,0);
880:     SNESMonitor(snes,snes->iter,snes->norm);
881:     /* Test for convergence */
882:     if (normschedule == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);}
883:     if (snes->reason) break;
884:     /* Call general purpose update function */
885:     if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
886:   }
887:   if (nasm->finaljacobian) {SNESNASMComputeFinalJacobian_Private(snes,X);}
888:   if (normschedule == SNES_NORM_ALWAYS) {
889:     if (i == snes->max_its) {
890:       PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
891:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
892:     }
893:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */
894:   return(0);
895: }

897: /*MC
898:   SNESNASM - Nonlinear Additive Schwartz

900:    Options Database:
901: +  -snes_nasm_log - enable logging events for the communication and solve stages
902: .  -snes_nasm_type <basic,restrict> - type of subdomain update used
903: .  -snes_asm_damping <dmp> - the new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
904: .  -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
905: .  -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> - pick state the jacobian is calculated at
906: .  -sub_snes_ - options prefix of the subdomain nonlinear solves
907: .  -sub_ksp_ - options prefix of the subdomain Krylov solver
908: -  -sub_pc_ - options prefix of the subdomain preconditioner

910:    Level: advanced

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

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

921: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
922: {
923:   SNES_NASM      *nasm;

927:   PetscNewLog(snes,&nasm);
928:   snes->data = (void*)nasm;

930:   nasm->n        = PETSC_DECIDE;
931:   nasm->subsnes  = 0;
932:   nasm->x        = 0;
933:   nasm->xl       = 0;
934:   nasm->y        = 0;
935:   nasm->b        = 0;
936:   nasm->oscatter = 0;
937:   nasm->iscatter = 0;
938:   nasm->gscatter = 0;
939:   nasm->damping  = 1.;

941:   nasm->type = PC_ASM_BASIC;
942:   nasm->finaljacobian = PETSC_FALSE;
943:   nasm->same_local_solves = PETSC_TRUE;

945:   snes->ops->destroy        = SNESDestroy_NASM;
946:   snes->ops->setup          = SNESSetUp_NASM;
947:   snes->ops->setfromoptions = SNESSetFromOptions_NASM;
948:   snes->ops->view           = SNESView_NASM;
949:   snes->ops->solve          = SNESSolve_NASM;
950:   snes->ops->reset          = SNESReset_NASM;

952:   snes->usesksp = PETSC_FALSE;
953:   snes->usespc  = PETSC_FALSE;

955:   nasm->fjtype              = 0;
956:   nasm->xinit               = NULL;
957:   nasm->eventrestrictinterp = 0;
958:   nasm->eventsubsolve       = 0;

960:   if (!snes->tolerancesset) {
961:     snes->max_its   = 10000;
962:     snes->max_funcs = 10000;
963:   }

965:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetType_C",SNESNASMSetType_NASM);
966:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetType_C",SNESNASMGetType_NASM);
967:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
968:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
969:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
970:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
971:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
972:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
973:   return(0);
974: }