Actual source code: nasm.c
petsc-3.8.4 2018-03-24
1: #include <petsc/private/snesimpl.h>
2: #include <petscdm.h>
4: typedef struct {
5: PetscInt n; /* local subdomains */
6: SNES *subsnes; /* nonlinear solvers for each subdomain */
7: Vec *x; /* solution vectors */
8: Vec *xl; /* solution local vectors */
9: Vec *y; /* step vectors */
10: Vec *b; /* rhs vectors */
11: Vec weight; /* weighting for adding updates on overlaps, in global space */
12: VecScatter *oscatter; /* scatter from global space to the subdomain global space */
13: VecScatter *oscatter_copy; /* copy of the above */
14: VecScatter *iscatter; /* scatter from global space to the nonoverlapping subdomain space */
15: VecScatter *gscatter; /* scatter from global space to the subdomain local space */
16: PCASMType type; /* ASM type */
17: PetscBool usesdm; /* use the DM for setting up the subproblems */
18: PetscBool finaljacobian; /* compute the jacobian of the converged solution */
19: PetscReal damping; /* damping parameter for updates from the blocks */
20: PetscBool same_local_solves; /* flag to determine if the solvers have been individually modified */
21: PetscBool weight_set; /* use a weight in the overlap updates */
23: /* logging events */
24: PetscLogEvent eventrestrictinterp;
25: PetscLogEvent eventsubsolve;
27: PetscInt fjtype; /* type of computed jacobian */
28: Vec xinit; /* initial solution in case the final jacobian type is computed as first */
29: } SNES_NASM;
31: const char *const SNESNASMTypes[] = {"NONE","RESTRICT","INTERPOLATE","BASIC","PCASMType","PC_ASM_",0};
32: const char *const SNESNASMFJTypes[] = {"FINALOUTER","FINALINNER","INITIAL"};
34: static PetscErrorCode SNESReset_NASM(SNES snes)
35: {
36: SNES_NASM *nasm = (SNES_NASM*)snes->data;
38: PetscInt i;
41: for (i=0; i<nasm->n; i++) {
42: if (nasm->xl) { VecDestroy(&nasm->xl[i]); }
43: if (nasm->x) { VecDestroy(&nasm->x[i]); }
44: if (nasm->y) { VecDestroy(&nasm->y[i]); }
45: if (nasm->b) { VecDestroy(&nasm->b[i]); }
47: if (nasm->subsnes) { SNESDestroy(&nasm->subsnes[i]); }
48: if (nasm->oscatter) { VecScatterDestroy(&nasm->oscatter[i]); }
49: if (nasm->oscatter_copy) { VecScatterDestroy(&nasm->oscatter_copy[i]); }
50: if (nasm->iscatter) { VecScatterDestroy(&nasm->iscatter[i]); }
51: if (nasm->gscatter) { VecScatterDestroy(&nasm->gscatter[i]); }
52: }
54: PetscFree(nasm->x);
55: PetscFree(nasm->xl);
56: PetscFree(nasm->y);
57: PetscFree(nasm->b);
59: if (nasm->xinit) {VecDestroy(&nasm->xinit);}
61: PetscFree(nasm->subsnes);
62: PetscFree(nasm->oscatter);
63: PetscFree(nasm->oscatter_copy);
64: PetscFree(nasm->iscatter);
65: PetscFree(nasm->gscatter);
67: if (nasm->weight_set) {
68: VecDestroy(&nasm->weight);
69: }
71: nasm->eventrestrictinterp = 0;
72: nasm->eventsubsolve = 0;
73: return(0);
74: }
76: static PetscErrorCode SNESDestroy_NASM(SNES snes)
77: {
81: SNESReset_NASM(snes);
82: PetscFree(snes->data);
83: return(0);
84: }
86: static PetscErrorCode DMGlobalToLocalSubDomainDirichletHook_Private(DM dm,Vec g,InsertMode mode,Vec l,void *ctx)
87: {
89: Vec bcs = (Vec)ctx;
92: VecCopy(bcs,l);
93: return(0);
94: }
96: static PetscErrorCode SNESSetUp_NASM(SNES snes)
97: {
98: SNES_NASM *nasm = (SNES_NASM*)snes->data;
100: DM dm,subdm;
101: DM *subdms;
102: PetscInt i;
103: const char *optionsprefix;
104: Vec F;
105: PetscMPIInt size;
106: KSP ksp;
107: PC pc;
110: if (!nasm->subsnes) {
111: SNESGetDM(snes,&dm);
112: if (dm) {
113: nasm->usesdm = PETSC_TRUE;
114: DMCreateDomainDecomposition(dm,&nasm->n,NULL,NULL,NULL,&subdms);
115: if (!subdms) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM has no default decomposition defined. Set subsolves manually with SNESNASMSetSubdomains().");
116: DMCreateDomainDecompositionScatters(dm,nasm->n,subdms,&nasm->iscatter,&nasm->oscatter,&nasm->gscatter);
117: PetscMalloc1(nasm->n, &nasm->oscatter_copy);
118: for (i=0; i<nasm->n; i++) {
119: VecScatterCopy(nasm->oscatter[i], &nasm->oscatter_copy[i]);
120: }
122: SNESGetOptionsPrefix(snes, &optionsprefix);
123: PetscMalloc1(nasm->n,&nasm->subsnes);
124: for (i=0; i<nasm->n; i++) {
125: SNESCreate(PETSC_COMM_SELF,&nasm->subsnes[i]);
126: PetscObjectIncrementTabLevel((PetscObject)nasm->subsnes[i], (PetscObject)snes, 1);
127: SNESAppendOptionsPrefix(nasm->subsnes[i],optionsprefix);
128: SNESAppendOptionsPrefix(nasm->subsnes[i],"sub_");
129: SNESSetDM(nasm->subsnes[i],subdms[i]);
130: MPI_Comm_size(PetscObjectComm((PetscObject)nasm->subsnes[i]),&size);
131: if (size == 1) {
132: SNESGetKSP(nasm->subsnes[i],&ksp);
133: KSPGetPC(ksp,&pc);
134: KSPSetType(ksp,KSPPREONLY);
135: PCSetType(pc,PCLU);
136: }
137: SNESSetFromOptions(nasm->subsnes[i]);
138: DMDestroy(&subdms[i]);
139: }
140: PetscFree(subdms);
141: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Cannot construct local problems automatically without a DM!");
142: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set subproblems manually if there is no DM!");
143: /* allocate the global vectors */
144: if (!nasm->x) {
145: PetscCalloc1(nasm->n,&nasm->x);
146: }
147: if (!nasm->xl) {
148: PetscCalloc1(nasm->n,&nasm->xl);
149: }
150: if (!nasm->y) {
151: PetscCalloc1(nasm->n,&nasm->y);
152: }
153: if (!nasm->b) {
154: PetscCalloc1(nasm->n,&nasm->b);
155: }
157: for (i=0; i<nasm->n; i++) {
158: SNESGetFunction(nasm->subsnes[i],&F,NULL,NULL);
159: if (!nasm->x[i]) {VecDuplicate(F,&nasm->x[i]);}
160: if (!nasm->y[i]) {VecDuplicate(F,&nasm->y[i]);}
161: if (!nasm->b[i]) {VecDuplicate(F,&nasm->b[i]);}
162: if (!nasm->xl[i]) {
163: SNESGetDM(nasm->subsnes[i],&subdm);
164: DMCreateLocalVector(subdm,&nasm->xl[i]);
165: DMGlobalToLocalHookAdd(subdm,DMGlobalToLocalSubDomainDirichletHook_Private,NULL,nasm->xl[i]);
166: }
167: }
168: if (nasm->finaljacobian) {
169: SNESSetUpMatrices(snes);
170: if (nasm->fjtype == 2) {
171: VecDuplicate(snes->vec_sol,&nasm->xinit);
172: }
173: for (i=0; i<nasm->n;i++) {
174: SNESSetUpMatrices(nasm->subsnes[i]);
175: }
176: }
177: return(0);
178: }
180: static PetscErrorCode SNESSetFromOptions_NASM(PetscOptionItems *PetscOptionsObject,SNES snes)
181: {
182: PetscErrorCode ierr;
183: PCASMType asmtype;
184: PetscBool flg,monflg,subviewflg;
185: SNES_NASM *nasm = (SNES_NASM*)snes->data;
188: PetscOptionsHead(PetscOptionsObject,"Nonlinear Additive Schwartz options");
189: PetscOptionsEnum("-snes_nasm_type","Type of restriction/extension","",SNESNASMTypes,(PetscEnum)nasm->type,(PetscEnum*)&asmtype,&flg);
190: if (flg) {SNESNASMSetType(snes,asmtype);}
191: flg = PETSC_FALSE;
192: monflg = PETSC_TRUE;
193: PetscOptionsReal("-snes_nasm_damping","The new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)","SNESNASMSetDamping",nasm->damping,&nasm->damping,&flg);
194: if (flg) {SNESNASMSetDamping(snes,nasm->damping);}
195: subviewflg = PETSC_FALSE;
196: PetscOptionsBool("-snes_nasm_sub_view","Print detailed information for every processor when using -snes_view","",subviewflg,&subviewflg,&flg);
197: if (flg) {
198: nasm->same_local_solves = PETSC_FALSE;
199: if (!subviewflg) {
200: nasm->same_local_solves = PETSC_TRUE;
201: }
202: }
203: PetscOptionsBool("-snes_nasm_finaljacobian","Compute the global jacobian of the final iterate (for ASPIN)","",nasm->finaljacobian,&nasm->finaljacobian,NULL);
204: PetscOptionsEList("-snes_nasm_finaljacobian_type","The type of the final jacobian computed.","",SNESNASMFJTypes,3,SNESNASMFJTypes[0],&nasm->fjtype,NULL);
205: PetscOptionsBool("-snes_nasm_log","Log times for subSNES solves and restriction","",monflg,&monflg,&flg);
206: if (flg) {
207: PetscLogEventRegister("SNESNASMSubSolve",((PetscObject)snes)->classid,&nasm->eventsubsolve);
208: PetscLogEventRegister("SNESNASMRestrict",((PetscObject)snes)->classid,&nasm->eventrestrictinterp);
209: }
210: PetscOptionsTail();
211: return(0);
212: }
214: static PetscErrorCode SNESView_NASM(SNES snes, PetscViewer viewer)
215: {
216: SNES_NASM *nasm = (SNES_NASM*)snes->data;
218: PetscMPIInt rank,size;
219: PetscInt i,N,bsz;
220: PetscBool iascii,isstring;
221: PetscViewer sviewer;
222: MPI_Comm comm;
225: PetscObjectGetComm((PetscObject)snes,&comm);
226: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
227: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
228: MPI_Comm_rank(comm,&rank);
229: MPI_Comm_size(comm,&size);
230: MPIU_Allreduce(&nasm->n,&N,1,MPIU_INT,MPI_SUM,comm);
231: if (iascii) {
232: PetscViewerASCIIPrintf(viewer, " total subdomain blocks = %D\n",N);
233: if (nasm->same_local_solves) {
234: if (nasm->subsnes) {
235: PetscViewerASCIIPrintf(viewer," Local solve is the same for all blocks:\n");
236: PetscViewerASCIIPushTab(viewer);
237: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
238: if (!rank) {
239: PetscViewerASCIIPushTab(viewer);
240: SNESView(nasm->subsnes[0],sviewer);
241: PetscViewerASCIIPopTab(viewer);
242: }
243: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
244: PetscViewerASCIIPopTab(viewer);
245: }
246: } else {
247: /* print the solver on each block */
248: PetscViewerASCIIPushSynchronized(viewer);
249: PetscViewerASCIISynchronizedPrintf(viewer," [%d] number of local blocks = %D\n",(int)rank,nasm->n);
250: PetscViewerFlush(viewer);
251: PetscViewerASCIIPopSynchronized(viewer);
252: PetscViewerASCIIPrintf(viewer," Local solve info for each block is in the following SNES objects:\n");
253: PetscViewerASCIIPushTab(viewer);
254: PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
255: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
256: for (i=0; i<nasm->n; i++) {
257: VecGetLocalSize(nasm->x[i],&bsz);
258: PetscViewerASCIIPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);
259: SNESView(nasm->subsnes[i],sviewer);
260: PetscViewerASCIIPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
261: }
262: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
263: PetscViewerFlush(viewer);
264: PetscViewerASCIIPopTab(viewer);
265: }
266: } else if (isstring) {
267: PetscViewerStringSPrintf(viewer," blocks=%D,type=%s",N,SNESNASMTypes[nasm->type]);
268: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
269: if (nasm->subsnes && !rank) {SNESView(nasm->subsnes[0],sviewer);}
270: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
271: }
272: return(0);
273: }
275: /*@
276: SNESNASMSetType - Set the type of subdomain update used
278: Logically Collective on SNES
280: Input Parameters:
281: + SNES - the SNES context
282: - type - the type of update, PC_ASM_BASIC or PC_ASM_RESTRICT
284: Level: intermediate
286: .keywords: SNES, NASM
288: .seealso: SNESNASM, SNESNASMGetType(), PCASMSetType()
289: @*/
290: PetscErrorCode SNESNASMSetType(SNES snes,PCASMType type)
291: {
293: PetscErrorCode (*f)(SNES,PCASMType);
296: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetType_C",&f);
297: if (f) {(f)(snes,type);}
298: return(0);
299: }
301: static PetscErrorCode SNESNASMSetType_NASM(SNES snes,PCASMType type)
302: {
303: SNES_NASM *nasm = (SNES_NASM*)snes->data;
306: if (type != PC_ASM_BASIC && type != PC_ASM_RESTRICT) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"SNESNASM only supports basic and restrict types");
307: nasm->type = type;
308: return(0);
309: }
311: /*@
312: SNESNASMGetType - Get the type of subdomain update used
314: Logically Collective on SNES
316: Input Parameters:
317: . SNES - the SNES context
319: Output Parameters:
320: . type - the type of update
322: Level: intermediate
324: .keywords: SNES, NASM
326: .seealso: SNESNASM, SNESNASMSetType(), PCASMGetType()
327: @*/
328: PetscErrorCode SNESNASMGetType(SNES snes,PCASMType *type)
329: {
333: PetscUseMethod(snes,"SNESNASMGetType_C",(SNES,PCASMType*),(snes,type));
334: return(0);
335: }
337: static PetscErrorCode SNESNASMGetType_NASM(SNES snes,PCASMType *type)
338: {
339: SNES_NASM *nasm = (SNES_NASM*)snes->data;
342: *type = nasm->type;
343: return(0);
344: }
346: /*@
347: SNESNASMSetSubdomains - Manually Set the context required to restrict and solve subdomain problems.
349: Not Collective
351: Input Parameters:
352: + SNES - the SNES context
353: . n - the number of local subdomains
354: . subsnes - solvers defined on the local subdomains
355: . iscatter - scatters into the nonoverlapping portions of the local subdomains
356: . oscatter - scatters into the overlapping portions of the local subdomains
357: - gscatter - scatters into the (ghosted) local vector of the local subdomain
359: Level: intermediate
361: .keywords: SNES, NASM
363: .seealso: SNESNASM, SNESNASMGetSubdomains()
364: @*/
365: PetscErrorCode SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
366: {
368: PetscErrorCode (*f)(SNES,PetscInt,SNES*,VecScatter*,VecScatter*,VecScatter*);
371: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",&f);
372: if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
373: return(0);
374: }
376: static PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
377: {
378: PetscInt i;
380: SNES_NASM *nasm = (SNES_NASM*)snes->data;
383: if (snes->setupcalled) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"SNESNASMSetSubdomains() should be called before calling SNESSetUp().");
385: /* tear down the previously set things */
386: SNESReset(snes);
388: nasm->n = n;
389: if (oscatter) {
390: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)oscatter[i]);}
391: }
392: if (iscatter) {
393: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)iscatter[i]);}
394: }
395: if (gscatter) {
396: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)gscatter[i]);}
397: }
398: if (oscatter) {
399: PetscMalloc1(n,&nasm->oscatter);
400: PetscMalloc1(n,&nasm->oscatter_copy);
401: for (i=0; i<n; i++) {
402: nasm->oscatter[i] = oscatter[i];
403: VecScatterCopy(oscatter[i], &nasm->oscatter_copy[i]);
404: }
405: }
406: if (iscatter) {
407: PetscMalloc1(n,&nasm->iscatter);
408: for (i=0; i<n; i++) {
409: nasm->iscatter[i] = iscatter[i];
410: }
411: }
412: if (gscatter) {
413: PetscMalloc1(n,&nasm->gscatter);
414: for (i=0; i<n; i++) {
415: nasm->gscatter[i] = gscatter[i];
416: }
417: }
419: if (subsnes) {
420: PetscMalloc1(n,&nasm->subsnes);
421: for (i=0; i<n; i++) {
422: nasm->subsnes[i] = subsnes[i];
423: }
424: nasm->same_local_solves = PETSC_FALSE;
425: }
426: return(0);
427: }
429: /*@
430: SNESNASMGetSubdomains - Get the local subdomain context.
432: Not Collective
434: Input Parameters:
435: . SNES - the SNES context
437: Output Parameters:
438: + n - the number of local subdomains
439: . subsnes - solvers defined on the local subdomains
440: . iscatter - scatters into the nonoverlapping portions of the local subdomains
441: . oscatter - scatters into the overlapping portions of the local subdomains
442: - gscatter - scatters into the (ghosted) local vector of the local subdomain
444: Level: intermediate
446: .keywords: SNES, NASM
448: .seealso: SNESNASM, SNESNASMSetSubdomains()
449: @*/
450: PetscErrorCode SNESNASMGetSubdomains(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
451: {
453: PetscErrorCode (*f)(SNES,PetscInt*,SNES**,VecScatter**,VecScatter**,VecScatter**);
456: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",&f);
457: if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
458: return(0);
459: }
461: static PetscErrorCode SNESNASMGetSubdomains_NASM(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
462: {
463: SNES_NASM *nasm = (SNES_NASM*)snes->data;
466: if (n) *n = nasm->n;
467: if (oscatter) *oscatter = nasm->oscatter;
468: if (iscatter) *iscatter = nasm->iscatter;
469: if (gscatter) *gscatter = nasm->gscatter;
470: if (subsnes) {
471: *subsnes = nasm->subsnes;
472: nasm->same_local_solves = PETSC_FALSE;
473: }
474: return(0);
475: }
477: /*@
478: SNESNASMGetSubdomainVecs - Get the processor-local subdomain vectors
480: Not Collective
482: Input Parameters:
483: . SNES - the SNES context
485: Output Parameters:
486: + n - the number of local subdomains
487: . x - The subdomain solution vector
488: . y - The subdomain step vector
489: . b - The subdomain RHS vector
490: - xl - The subdomain local vectors (ghosted)
492: Level: developer
494: .keywords: SNES, NASM
496: .seealso: SNESNASM, SNESNASMGetSubdomains()
497: @*/
498: PetscErrorCode SNESNASMGetSubdomainVecs(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b, Vec **xl)
499: {
501: PetscErrorCode (*f)(SNES,PetscInt*,Vec**,Vec**,Vec**,Vec**);
504: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",&f);
505: if (f) {(f)(snes,n,x,y,b,xl);}
506: return(0);
507: }
509: static PetscErrorCode SNESNASMGetSubdomainVecs_NASM(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b,Vec **xl)
510: {
511: SNES_NASM *nasm = (SNES_NASM*)snes->data;
514: if (n) *n = nasm->n;
515: if (x) *x = nasm->x;
516: if (y) *y = nasm->y;
517: if (b) *b = nasm->b;
518: if (xl) *xl = nasm->xl;
519: return(0);
520: }
522: /*@
523: SNESNASMSetComputeFinalJacobian - Schedules the computation of the global and subdomain jacobians upon convergence
525: Collective on SNES
527: Input Parameters:
528: + SNES - the SNES context
529: - flg - indication of whether to compute the jacobians or not
531: Level: developer
533: Notes: This is used almost exclusively in the implementation of ASPIN, where the converged subdomain and global jacobian
534: is needed at each linear iteration.
536: .keywords: SNES, NASM, ASPIN
538: .seealso: SNESNASM, SNESNASMGetSubdomains()
539: @*/
540: PetscErrorCode SNESNASMSetComputeFinalJacobian(SNES snes,PetscBool flg)
541: {
542: PetscErrorCode (*f)(SNES,PetscBool);
546: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",&f);
547: if (f) {(f)(snes,flg);}
548: return(0);
549: }
551: static PetscErrorCode SNESNASMSetComputeFinalJacobian_NASM(SNES snes,PetscBool flg)
552: {
553: SNES_NASM *nasm = (SNES_NASM*)snes->data;
556: nasm->finaljacobian = flg;
557: if (flg) snes->usesksp = PETSC_TRUE;
558: return(0);
559: }
561: /*@
562: SNESNASMSetDamping - Sets the update damping for NASM
564: Logically collective on SNES
566: Input Parameters:
567: + SNES - the SNES context
568: - dmp - damping
570: Level: intermediate
572: Notes: The new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
574: .keywords: SNES, NASM, damping
576: .seealso: SNESNASM, SNESNASMGetDamping()
577: @*/
578: PetscErrorCode SNESNASMSetDamping(SNES snes,PetscReal dmp)
579: {
580: PetscErrorCode (*f)(SNES,PetscReal);
584: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetDamping_C",(void (**)(void))&f);
585: if (f) {(f)(snes,dmp);}
586: return(0);
587: }
589: static PetscErrorCode SNESNASMSetDamping_NASM(SNES snes,PetscReal dmp)
590: {
591: SNES_NASM *nasm = (SNES_NASM*)snes->data;
594: nasm->damping = dmp;
595: return(0);
596: }
598: /*@
599: SNESNASMGetDamping - Gets the update damping for NASM
601: Not Collective
603: Input Parameters:
604: + SNES - the SNES context
605: - dmp - damping
607: Level: intermediate
609: .keywords: SNES, NASM, damping
611: .seealso: SNESNASM, SNESNASMSetDamping()
612: @*/
613: PetscErrorCode SNESNASMGetDamping(SNES snes,PetscReal *dmp)
614: {
618: PetscUseMethod(snes,"SNESNASMGetDamping_C",(SNES,PetscReal*),(snes,dmp));
619: return(0);
620: }
622: static PetscErrorCode SNESNASMGetDamping_NASM(SNES snes,PetscReal *dmp)
623: {
624: SNES_NASM *nasm = (SNES_NASM*)snes->data;
627: *dmp = nasm->damping;
628: return(0);
629: }
632: /*
633: Input Parameters:
634: + snes - The solver
635: . B - The RHS vector
636: - X - The initial guess
638: Output Parameters:
639: . Y - The solution update
641: TODO: All scatters should be packed into one
642: */
643: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X)
644: {
645: SNES_NASM *nasm = (SNES_NASM*)snes->data;
646: SNES subsnes;
647: PetscInt i;
648: PetscReal dmp;
650: Vec Xl,Bl,Yl,Xlloc;
651: VecScatter iscat,oscat,gscat,oscat_copy;
652: DM dm,subdm;
653: PCASMType type;
656: SNESNASMGetType(snes,&type);
657: SNESGetDM(snes,&dm);
658: VecSet(Y,0);
659: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
660: for (i=0; i<nasm->n; i++) {
661: /* scatter the solution to the global solution and the local solution */
662: Xl = nasm->x[i];
663: Xlloc = nasm->xl[i];
664: oscat = nasm->oscatter[i];
665: oscat_copy = nasm->oscatter_copy[i];
666: gscat = nasm->gscatter[i];
667: VecScatterBegin(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
668: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
669: if (B) {
670: /* scatter the RHS to the local RHS */
671: Bl = nasm->b[i];
672: VecScatterBegin(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
673: }
674: }
675: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
678: if (nasm->eventsubsolve) {PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);}
679: for (i=0; i<nasm->n; i++) {
680: Xl = nasm->x[i];
681: Xlloc = nasm->xl[i];
682: Yl = nasm->y[i];
683: subsnes = nasm->subsnes[i];
684: SNESGetDM(subsnes,&subdm);
685: iscat = nasm->iscatter[i];
686: oscat = nasm->oscatter[i];
687: oscat_copy = nasm->oscatter_copy[i];
688: gscat = nasm->gscatter[i];
689: VecScatterEnd(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
690: VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
691: if (B) {
692: Bl = nasm->b[i];
693: VecScatterEnd(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
694: } else Bl = NULL;
696: DMSubDomainRestrict(dm,oscat,gscat,subdm);
697: VecCopy(Xl,Yl);
698: SNESSolve(subsnes,Bl,Xl);
699: VecAYPX(Yl,-1.0,Xl);
700: VecScale(Yl, nasm->damping);
701: if (type == PC_ASM_BASIC) {
702: VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
703: } else if (type == PC_ASM_RESTRICT) {
704: VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
705: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
706: }
707: if (nasm->eventsubsolve) {PetscLogEventEnd(nasm->eventsubsolve,snes,0,0,0);}
708: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
709: for (i=0; i<nasm->n; i++) {
710: Yl = nasm->y[i];
711: iscat = nasm->iscatter[i];
712: oscat = nasm->oscatter[i];
713: if (type == PC_ASM_BASIC) {
714: VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
715: } else if (type == PC_ASM_RESTRICT) {
716: VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
717: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
718: }
719: if (nasm->weight_set) {
720: VecPointwiseMult(Y,Y,nasm->weight);
721: }
722: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
723: SNESNASMGetDamping(snes,&dmp);
724: VecAXPY(X,dmp,Y);
725: return(0);
726: }
728: static PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
729: {
730: Vec X = Xfinal;
731: SNES_NASM *nasm = (SNES_NASM*)snes->data;
732: SNES subsnes;
733: PetscInt i,lag = 1;
735: Vec Xlloc,Xl,Fl,F;
736: VecScatter oscat,gscat;
737: DM dm,subdm;
740: if (nasm->fjtype == 2) X = nasm->xinit;
741: F = snes->vec_func;
742: if (snes->normschedule == SNES_NORM_NONE) {SNESComputeFunction(snes,X,F);}
743: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
744: SNESGetDM(snes,&dm);
745: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
746: if (nasm->fjtype != 1) {
747: for (i=0; i<nasm->n; i++) {
748: Xlloc = nasm->xl[i];
749: gscat = nasm->gscatter[i];
750: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
751: }
752: }
753: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
754: for (i=0; i<nasm->n; i++) {
755: Fl = nasm->subsnes[i]->vec_func;
756: Xl = nasm->x[i];
757: Xlloc = nasm->xl[i];
758: subsnes = nasm->subsnes[i];
759: oscat = nasm->oscatter[i];
760: gscat = nasm->gscatter[i];
761: if (nasm->fjtype != 1) {VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);}
762: SNESGetDM(subsnes,&subdm);
763: DMSubDomainRestrict(dm,oscat,gscat,subdm);
764: if (nasm->fjtype != 1) {
765: DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
766: DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
767: }
768: if (subsnes->lagjacobian == -1) subsnes->lagjacobian = -2;
769: else if (subsnes->lagjacobian > 1) lag = subsnes->lagjacobian;
770: SNESComputeFunction(subsnes,Xl,Fl);
771: SNESComputeJacobian(subsnes,Xl,subsnes->jacobian,subsnes->jacobian_pre);
772: if (lag > 1) subsnes->lagjacobian = lag;
773: }
774: return(0);
775: }
777: static PetscErrorCode SNESSolve_NASM(SNES snes)
778: {
779: Vec F;
780: Vec X;
781: Vec B;
782: Vec Y;
783: PetscInt i;
784: PetscReal fnorm = 0.0;
785: PetscErrorCode ierr;
786: SNESNormSchedule normschedule;
787: SNES_NASM *nasm = (SNES_NASM*)snes->data;
791: if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
793: PetscCitationsRegister(SNESCitation,&SNEScite);
794: X = snes->vec_sol;
795: Y = snes->vec_sol_update;
796: F = snes->vec_func;
797: B = snes->vec_rhs;
799: PetscObjectSAWsTakeAccess((PetscObject)snes);
800: snes->iter = 0;
801: snes->norm = 0.;
802: PetscObjectSAWsGrantAccess((PetscObject)snes);
803: snes->reason = SNES_CONVERGED_ITERATING;
804: SNESGetNormSchedule(snes, &normschedule);
805: if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
806: /* compute the initial function and preconditioned update delX */
807: if (!snes->vec_func_init_set) {
808: SNESComputeFunction(snes,X,F);
809: } else snes->vec_func_init_set = PETSC_FALSE;
811: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
812: SNESCheckFunctionNorm(snes,fnorm);
813: PetscObjectSAWsTakeAccess((PetscObject)snes);
814: snes->iter = 0;
815: snes->norm = fnorm;
816: PetscObjectSAWsGrantAccess((PetscObject)snes);
817: SNESLogConvergenceHistory(snes,snes->norm,0);
818: SNESMonitor(snes,0,snes->norm);
820: /* test convergence */
821: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
822: if (snes->reason) return(0);
823: } else {
824: PetscObjectSAWsGrantAccess((PetscObject)snes);
825: SNESLogConvergenceHistory(snes,snes->norm,0);
826: SNESMonitor(snes,0,snes->norm);
827: }
829: /* Call general purpose update function */
830: if (snes->ops->update) {
831: (*snes->ops->update)(snes, snes->iter);
832: }
833: /* copy the initial solution over for later */
834: if (nasm->fjtype == 2) {VecCopy(X,nasm->xinit);}
836: for (i=0; i < snes->max_its; i++) {
837: SNESNASMSolveLocal_Private(snes,B,Y,X);
838: if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
839: SNESComputeFunction(snes,X,F);
840: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
841: SNESCheckFunctionNorm(snes,fnorm);
842: }
843: /* Monitor convergence */
844: PetscObjectSAWsTakeAccess((PetscObject)snes);
845: snes->iter = i+1;
846: snes->norm = fnorm;
847: PetscObjectSAWsGrantAccess((PetscObject)snes);
848: SNESLogConvergenceHistory(snes,snes->norm,0);
849: SNESMonitor(snes,snes->iter,snes->norm);
850: /* Test for convergence */
851: if (normschedule == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);}
852: if (snes->reason) break;
853: /* Call general purpose update function */
854: if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
855: }
856: if (nasm->finaljacobian) {SNESNASMComputeFinalJacobian_Private(snes,X);}
857: if (normschedule == SNES_NORM_ALWAYS) {
858: if (i == snes->max_its) {
859: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
860: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
861: }
862: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */
863: return(0);
864: }
866: /*MC
867: SNESNASM - Nonlinear Additive Schwartz
869: Options Database:
870: + -snes_nasm_log - enable logging events for the communication and solve stages
871: . -snes_nasm_type <basic,restrict> - type of subdomain update used
872: . -snes_asm_damping <dmp> - the new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
873: . -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
874: . -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> - pick state the jacobian is calculated at
875: . -sub_snes_ - options prefix of the subdomain nonlinear solves
876: . -sub_ksp_ - options prefix of the subdomain Krylov solver
877: - -sub_pc_ - options prefix of the subdomain preconditioner
879: Level: advanced
881: References:
882: . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
883: SIAM Review, 57(4), 2015
885: .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types), SNESNASMSetType(), SNESNASMGetType(), SNESNASMSetSubdomains(), SNESNASMGetSubdomains(), SNESNASMGetSubdomainVecs(), SNESNASMSetComputeFinalJacobian(), SNESNASMSetDamping(), SNESNASMGetDamping()
886: M*/
888: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
889: {
890: SNES_NASM *nasm;
894: PetscNewLog(snes,&nasm);
895: snes->data = (void*)nasm;
897: nasm->n = PETSC_DECIDE;
898: nasm->subsnes = 0;
899: nasm->x = 0;
900: nasm->xl = 0;
901: nasm->y = 0;
902: nasm->b = 0;
903: nasm->oscatter = 0;
904: nasm->oscatter_copy = 0;
905: nasm->iscatter = 0;
906: nasm->gscatter = 0;
907: nasm->damping = 1.;
909: nasm->type = PC_ASM_BASIC;
910: nasm->finaljacobian = PETSC_FALSE;
911: nasm->same_local_solves = PETSC_TRUE;
912: nasm->weight_set = PETSC_FALSE;
914: snes->ops->destroy = SNESDestroy_NASM;
915: snes->ops->setup = SNESSetUp_NASM;
916: snes->ops->setfromoptions = SNESSetFromOptions_NASM;
917: snes->ops->view = SNESView_NASM;
918: snes->ops->solve = SNESSolve_NASM;
919: snes->ops->reset = SNESReset_NASM;
921: snes->usesksp = PETSC_FALSE;
922: snes->usesnpc = PETSC_FALSE;
924: snes->alwayscomputesfinalresidual = PETSC_FALSE;
926: nasm->fjtype = 0;
927: nasm->xinit = NULL;
928: nasm->eventrestrictinterp = 0;
929: nasm->eventsubsolve = 0;
931: if (!snes->tolerancesset) {
932: snes->max_its = 10000;
933: snes->max_funcs = 10000;
934: }
936: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetType_C",SNESNASMSetType_NASM);
937: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetType_C",SNESNASMGetType_NASM);
938: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
939: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
940: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
941: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
942: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
943: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
944: return(0);
945: }
947: /*@
948: SNESNASMGetSNES - Gets a subsolver
950: Not collective
952: Input Parameters:
953: + snes - the SNES context
954: - i - the number of the subsnes to get
956: Output Parameters:
957: . subsnes - the subsolver context
959: Level: intermediate
961: .keywords: SNES, NASM
963: .seealso: SNESNASM, SNESNASMGetNumber()
964: @*/
965: PetscErrorCode SNESNASMGetSNES(SNES snes,PetscInt i,SNES *subsnes)
966: {
967: SNES_NASM *nasm = (SNES_NASM*)snes->data;
970: if (i < 0 || i >= nasm->n) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"No such subsolver");
971: *subsnes = nasm->subsnes[i];
972: return(0);
973: }
975: /*@
976: SNESNASMGetNumber - Gets number of subsolvers
978: Not collective
980: Input Parameters:
981: . snes - the SNES context
983: Output Parameters:
984: . n - the number of subsolvers
986: Level: intermediate
988: .keywords: SNES, NASM
990: .seealso: SNESNASM, SNESNASMGetSNES()
991: @*/
992: PetscErrorCode SNESNASMGetNumber(SNES snes,PetscInt *n)
993: {
994: SNES_NASM *nasm = (SNES_NASM*)snes->data;
997: *n = nasm->n;
998: return(0);
999: }
1001: /*@
1002: SNESNASMSetWeight - Sets weight to use when adding overlapping updates
1004: Collective
1006: Input Parameters:
1007: + snes - the SNES context
1008: - weight - the weights to use (typically 1/N for each dof, where N is the number of patches it appears in)
1010: Level: intermediate
1012: .keywords: SNES, NASM
1014: .seealso: SNESNASM
1015: @*/
1016: PetscErrorCode SNESNASMSetWeight(SNES snes,Vec weight)
1017: {
1018: SNES_NASM *nasm = (SNES_NASM*)snes->data;
1023: VecDestroy(&nasm->weight);
1024: nasm->weight_set = PETSC_TRUE;
1025: nasm->weight = weight;
1026: PetscObjectReference((PetscObject)nasm->weight);
1028: return(0);
1029: }