Actual source code: nasm.c
petsc-3.10.5 2019-03-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: if (flg) snes->usesksp = PETSC_TRUE;
559: return(0);
560: }
562: /*@
563: SNESNASMSetDamping - Sets the update damping for NASM
565: Logically collective on SNES
567: Input Parameters:
568: + SNES - the SNES context
569: - dmp - damping
571: Level: intermediate
573: Notes:
574: The new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
576: .keywords: SNES, NASM, damping
578: .seealso: SNESNASM, SNESNASMGetDamping()
579: @*/
580: PetscErrorCode SNESNASMSetDamping(SNES snes,PetscReal dmp)
581: {
582: PetscErrorCode (*f)(SNES,PetscReal);
586: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetDamping_C",(void (**)(void))&f);
587: if (f) {(f)(snes,dmp);}
588: return(0);
589: }
591: static PetscErrorCode SNESNASMSetDamping_NASM(SNES snes,PetscReal dmp)
592: {
593: SNES_NASM *nasm = (SNES_NASM*)snes->data;
596: nasm->damping = dmp;
597: return(0);
598: }
600: /*@
601: SNESNASMGetDamping - Gets the update damping for NASM
603: Not Collective
605: Input Parameters:
606: + SNES - the SNES context
607: - dmp - damping
609: Level: intermediate
611: .keywords: SNES, NASM, damping
613: .seealso: SNESNASM, SNESNASMSetDamping()
614: @*/
615: PetscErrorCode SNESNASMGetDamping(SNES snes,PetscReal *dmp)
616: {
620: PetscUseMethod(snes,"SNESNASMGetDamping_C",(SNES,PetscReal*),(snes,dmp));
621: return(0);
622: }
624: static PetscErrorCode SNESNASMGetDamping_NASM(SNES snes,PetscReal *dmp)
625: {
626: SNES_NASM *nasm = (SNES_NASM*)snes->data;
629: *dmp = nasm->damping;
630: return(0);
631: }
634: /*
635: Input Parameters:
636: + snes - The solver
637: . B - The RHS vector
638: - X - The initial guess
640: Output Parameters:
641: . Y - The solution update
643: TODO: All scatters should be packed into one
644: */
645: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X)
646: {
647: SNES_NASM *nasm = (SNES_NASM*)snes->data;
648: SNES subsnes;
649: PetscInt i;
650: PetscReal dmp;
652: Vec Xl,Bl,Yl,Xlloc;
653: VecScatter iscat,oscat,gscat,oscat_copy;
654: DM dm,subdm;
655: PCASMType type;
658: SNESNASMGetType(snes,&type);
659: SNESGetDM(snes,&dm);
660: VecSet(Y,0);
661: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
662: for (i=0; i<nasm->n; i++) {
663: /* scatter the solution to the global solution and the local solution */
664: Xl = nasm->x[i];
665: Xlloc = nasm->xl[i];
666: oscat = nasm->oscatter[i];
667: oscat_copy = nasm->oscatter_copy[i];
668: gscat = nasm->gscatter[i];
669: VecScatterBegin(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
670: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
671: if (B) {
672: /* scatter the RHS to the local RHS */
673: Bl = nasm->b[i];
674: VecScatterBegin(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
675: }
676: }
677: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
680: if (nasm->eventsubsolve) {PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);}
681: for (i=0; i<nasm->n; i++) {
682: Xl = nasm->x[i];
683: Xlloc = nasm->xl[i];
684: Yl = nasm->y[i];
685: subsnes = nasm->subsnes[i];
686: SNESGetDM(subsnes,&subdm);
687: iscat = nasm->iscatter[i];
688: oscat = nasm->oscatter[i];
689: oscat_copy = nasm->oscatter_copy[i];
690: gscat = nasm->gscatter[i];
691: VecScatterEnd(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
692: VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
693: if (B) {
694: Bl = nasm->b[i];
695: VecScatterEnd(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
696: } else Bl = NULL;
698: DMSubDomainRestrict(dm,oscat,gscat,subdm);
699: VecCopy(Xl,Yl);
700: SNESSolve(subsnes,Bl,Xl);
701: VecAYPX(Yl,-1.0,Xl);
702: VecScale(Yl, nasm->damping);
703: if (type == PC_ASM_BASIC) {
704: VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
705: } else if (type == PC_ASM_RESTRICT) {
706: VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
707: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
708: }
709: if (nasm->eventsubsolve) {PetscLogEventEnd(nasm->eventsubsolve,snes,0,0,0);}
710: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
711: for (i=0; i<nasm->n; i++) {
712: Yl = nasm->y[i];
713: iscat = nasm->iscatter[i];
714: oscat = nasm->oscatter[i];
715: if (type == PC_ASM_BASIC) {
716: VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
717: } else if (type == PC_ASM_RESTRICT) {
718: VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
719: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
720: }
721: if (nasm->weight_set) {
722: VecPointwiseMult(Y,Y,nasm->weight);
723: }
724: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
725: SNESNASMGetDamping(snes,&dmp);
726: VecAXPY(X,dmp,Y);
727: return(0);
728: }
730: static PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
731: {
732: Vec X = Xfinal;
733: SNES_NASM *nasm = (SNES_NASM*)snes->data;
734: SNES subsnes;
735: PetscInt i,lag = 1;
737: Vec Xlloc,Xl,Fl,F;
738: VecScatter oscat,gscat;
739: DM dm,subdm;
742: if (nasm->fjtype == 2) X = nasm->xinit;
743: F = snes->vec_func;
744: if (snes->normschedule == SNES_NORM_NONE) {SNESComputeFunction(snes,X,F);}
745: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
746: SNESGetDM(snes,&dm);
747: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
748: if (nasm->fjtype != 1) {
749: for (i=0; i<nasm->n; i++) {
750: Xlloc = nasm->xl[i];
751: gscat = nasm->gscatter[i];
752: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
753: }
754: }
755: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
756: for (i=0; i<nasm->n; i++) {
757: Fl = nasm->subsnes[i]->vec_func;
758: Xl = nasm->x[i];
759: Xlloc = nasm->xl[i];
760: subsnes = nasm->subsnes[i];
761: oscat = nasm->oscatter[i];
762: gscat = nasm->gscatter[i];
763: if (nasm->fjtype != 1) {VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);}
764: SNESGetDM(subsnes,&subdm);
765: DMSubDomainRestrict(dm,oscat,gscat,subdm);
766: if (nasm->fjtype != 1) {
767: DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
768: DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
769: }
770: if (subsnes->lagjacobian == -1) subsnes->lagjacobian = -2;
771: else if (subsnes->lagjacobian > 1) lag = subsnes->lagjacobian;
772: SNESComputeFunction(subsnes,Xl,Fl);
773: SNESComputeJacobian(subsnes,Xl,subsnes->jacobian,subsnes->jacobian_pre);
774: if (lag > 1) subsnes->lagjacobian = lag;
775: }
776: return(0);
777: }
779: static PetscErrorCode SNESSolve_NASM(SNES snes)
780: {
781: Vec F;
782: Vec X;
783: Vec B;
784: Vec Y;
785: PetscInt i;
786: PetscReal fnorm = 0.0;
787: PetscErrorCode ierr;
788: SNESNormSchedule normschedule;
789: SNES_NASM *nasm = (SNES_NASM*)snes->data;
793: 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);
795: PetscCitationsRegister(SNESCitation,&SNEScite);
796: X = snes->vec_sol;
797: Y = snes->vec_sol_update;
798: F = snes->vec_func;
799: B = snes->vec_rhs;
801: PetscObjectSAWsTakeAccess((PetscObject)snes);
802: snes->iter = 0;
803: snes->norm = 0.;
804: PetscObjectSAWsGrantAccess((PetscObject)snes);
805: snes->reason = SNES_CONVERGED_ITERATING;
806: SNESGetNormSchedule(snes, &normschedule);
807: if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
808: /* compute the initial function and preconditioned update delX */
809: if (!snes->vec_func_init_set) {
810: SNESComputeFunction(snes,X,F);
811: } else snes->vec_func_init_set = PETSC_FALSE;
813: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
814: SNESCheckFunctionNorm(snes,fnorm);
815: PetscObjectSAWsTakeAccess((PetscObject)snes);
816: snes->iter = 0;
817: snes->norm = fnorm;
818: PetscObjectSAWsGrantAccess((PetscObject)snes);
819: SNESLogConvergenceHistory(snes,snes->norm,0);
820: SNESMonitor(snes,0,snes->norm);
822: /* test convergence */
823: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
824: if (snes->reason) return(0);
825: } else {
826: PetscObjectSAWsGrantAccess((PetscObject)snes);
827: SNESLogConvergenceHistory(snes,snes->norm,0);
828: SNESMonitor(snes,0,snes->norm);
829: }
831: /* Call general purpose update function */
832: if (snes->ops->update) {
833: (*snes->ops->update)(snes, snes->iter);
834: }
835: /* copy the initial solution over for later */
836: if (nasm->fjtype == 2) {VecCopy(X,nasm->xinit);}
838: for (i=0; i < snes->max_its; i++) {
839: SNESNASMSolveLocal_Private(snes,B,Y,X);
840: if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
841: SNESComputeFunction(snes,X,F);
842: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
843: SNESCheckFunctionNorm(snes,fnorm);
844: }
845: /* Monitor convergence */
846: PetscObjectSAWsTakeAccess((PetscObject)snes);
847: snes->iter = i+1;
848: snes->norm = fnorm;
849: PetscObjectSAWsGrantAccess((PetscObject)snes);
850: SNESLogConvergenceHistory(snes,snes->norm,0);
851: SNESMonitor(snes,snes->iter,snes->norm);
852: /* Test for convergence */
853: if (normschedule == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);}
854: if (snes->reason) break;
855: /* Call general purpose update function */
856: if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
857: }
858: if (nasm->finaljacobian) {SNESNASMComputeFinalJacobian_Private(snes,X);}
859: if (normschedule == SNES_NORM_ALWAYS) {
860: if (i == snes->max_its) {
861: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
862: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
863: }
864: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */
865: return(0);
866: }
868: /*MC
869: SNESNASM - Nonlinear Additive Schwartz
871: Options Database:
872: + -snes_nasm_log - enable logging events for the communication and solve stages
873: . -snes_nasm_type <basic,restrict> - type of subdomain update used
874: . -snes_asm_damping <dmp> - the new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
875: . -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
876: . -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> - pick state the jacobian is calculated at
877: . -sub_snes_ - options prefix of the subdomain nonlinear solves
878: . -sub_ksp_ - options prefix of the subdomain Krylov solver
879: - -sub_pc_ - options prefix of the subdomain preconditioner
881: Level: advanced
883: References:
884: . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
885: SIAM Review, 57(4), 2015
887: .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types), SNESNASMSetType(), SNESNASMGetType(), SNESNASMSetSubdomains(), SNESNASMGetSubdomains(), SNESNASMGetSubdomainVecs(), SNESNASMSetComputeFinalJacobian(), SNESNASMSetDamping(), SNESNASMGetDamping()
888: M*/
890: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
891: {
892: SNES_NASM *nasm;
896: PetscNewLog(snes,&nasm);
897: snes->data = (void*)nasm;
899: nasm->n = PETSC_DECIDE;
900: nasm->subsnes = 0;
901: nasm->x = 0;
902: nasm->xl = 0;
903: nasm->y = 0;
904: nasm->b = 0;
905: nasm->oscatter = 0;
906: nasm->oscatter_copy = 0;
907: nasm->iscatter = 0;
908: nasm->gscatter = 0;
909: nasm->damping = 1.;
911: nasm->type = PC_ASM_BASIC;
912: nasm->finaljacobian = PETSC_FALSE;
913: nasm->same_local_solves = PETSC_TRUE;
914: nasm->weight_set = PETSC_FALSE;
916: snes->ops->destroy = SNESDestroy_NASM;
917: snes->ops->setup = SNESSetUp_NASM;
918: snes->ops->setfromoptions = SNESSetFromOptions_NASM;
919: snes->ops->view = SNESView_NASM;
920: snes->ops->solve = SNESSolve_NASM;
921: snes->ops->reset = SNESReset_NASM;
923: snes->usesksp = PETSC_FALSE;
924: snes->usesnpc = PETSC_FALSE;
926: snes->alwayscomputesfinalresidual = PETSC_FALSE;
928: nasm->fjtype = 0;
929: nasm->xinit = NULL;
930: nasm->eventrestrictinterp = 0;
931: nasm->eventsubsolve = 0;
933: if (!snes->tolerancesset) {
934: snes->max_its = 10000;
935: snes->max_funcs = 10000;
936: }
938: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetType_C",SNESNASMSetType_NASM);
939: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetType_C",SNESNASMGetType_NASM);
940: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
941: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
942: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
943: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
944: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
945: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
946: return(0);
947: }
949: /*@
950: SNESNASMGetSNES - Gets a subsolver
952: Not collective
954: Input Parameters:
955: + snes - the SNES context
956: - i - the number of the subsnes to get
958: Output Parameters:
959: . subsnes - the subsolver context
961: Level: intermediate
963: .keywords: SNES, NASM
965: .seealso: SNESNASM, SNESNASMGetNumber()
966: @*/
967: PetscErrorCode SNESNASMGetSNES(SNES snes,PetscInt i,SNES *subsnes)
968: {
969: SNES_NASM *nasm = (SNES_NASM*)snes->data;
972: if (i < 0 || i >= nasm->n) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"No such subsolver");
973: *subsnes = nasm->subsnes[i];
974: return(0);
975: }
977: /*@
978: SNESNASMGetNumber - Gets number of subsolvers
980: Not collective
982: Input Parameters:
983: . snes - the SNES context
985: Output Parameters:
986: . n - the number of subsolvers
988: Level: intermediate
990: .keywords: SNES, NASM
992: .seealso: SNESNASM, SNESNASMGetSNES()
993: @*/
994: PetscErrorCode SNESNASMGetNumber(SNES snes,PetscInt *n)
995: {
996: SNES_NASM *nasm = (SNES_NASM*)snes->data;
999: *n = nasm->n;
1000: return(0);
1001: }
1003: /*@
1004: SNESNASMSetWeight - Sets weight to use when adding overlapping updates
1006: Collective
1008: Input Parameters:
1009: + snes - the SNES context
1010: - weight - the weights to use (typically 1/N for each dof, where N is the number of patches it appears in)
1012: Level: intermediate
1014: .keywords: SNES, NASM
1016: .seealso: SNESNASM
1017: @*/
1018: PetscErrorCode SNESNASMSetWeight(SNES snes,Vec weight)
1019: {
1020: SNES_NASM *nasm = (SNES_NASM*)snes->data;
1025: VecDestroy(&nasm->weight);
1026: nasm->weight_set = PETSC_TRUE;
1027: nasm->weight = weight;
1028: PetscObjectReference((PetscObject)nasm->weight);
1030: return(0);
1031: }