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