Actual source code: nasm.c
petsc-3.13.6 2020-09-29
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 Schwarz 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: .seealso: SNESNASM, SNESNASMGetType(), PCASMSetType()
287: @*/
288: PetscErrorCode SNESNASMSetType(SNES snes,PCASMType type)
289: {
291: PetscErrorCode (*f)(SNES,PCASMType);
294: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetType_C",&f);
295: if (f) {(f)(snes,type);}
296: return(0);
297: }
299: static PetscErrorCode SNESNASMSetType_NASM(SNES snes,PCASMType type)
300: {
301: SNES_NASM *nasm = (SNES_NASM*)snes->data;
304: if (type != PC_ASM_BASIC && type != PC_ASM_RESTRICT) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"SNESNASM only supports basic and restrict types");
305: nasm->type = type;
306: return(0);
307: }
309: /*@
310: SNESNASMGetType - Get the type of subdomain update used
312: Logically Collective on SNES
314: Input Parameters:
315: . SNES - the SNES context
317: Output Parameters:
318: . type - the type of update
320: Level: intermediate
322: .seealso: SNESNASM, SNESNASMSetType(), PCASMGetType()
323: @*/
324: PetscErrorCode SNESNASMGetType(SNES snes,PCASMType *type)
325: {
329: PetscUseMethod(snes,"SNESNASMGetType_C",(SNES,PCASMType*),(snes,type));
330: return(0);
331: }
333: static PetscErrorCode SNESNASMGetType_NASM(SNES snes,PCASMType *type)
334: {
335: SNES_NASM *nasm = (SNES_NASM*)snes->data;
338: *type = nasm->type;
339: return(0);
340: }
342: /*@
343: SNESNASMSetSubdomains - Manually Set the context required to restrict and solve subdomain problems.
345: Not Collective
347: Input Parameters:
348: + SNES - the SNES context
349: . n - the number of local subdomains
350: . subsnes - solvers defined on the local subdomains
351: . iscatter - scatters into the nonoverlapping portions of the local subdomains
352: . oscatter - scatters into the overlapping portions of the local subdomains
353: - gscatter - scatters into the (ghosted) local vector of the local subdomain
355: Level: intermediate
357: .seealso: SNESNASM, SNESNASMGetSubdomains()
358: @*/
359: PetscErrorCode SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
360: {
362: PetscErrorCode (*f)(SNES,PetscInt,SNES*,VecScatter*,VecScatter*,VecScatter*);
365: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",&f);
366: if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
367: return(0);
368: }
370: static PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
371: {
372: PetscInt i;
374: SNES_NASM *nasm = (SNES_NASM*)snes->data;
377: if (snes->setupcalled) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"SNESNASMSetSubdomains() should be called before calling SNESSetUp().");
379: /* tear down the previously set things */
380: SNESReset(snes);
382: nasm->n = n;
383: if (oscatter) {
384: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)oscatter[i]);}
385: }
386: if (iscatter) {
387: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)iscatter[i]);}
388: }
389: if (gscatter) {
390: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)gscatter[i]);}
391: }
392: if (oscatter) {
393: PetscMalloc1(n,&nasm->oscatter);
394: PetscMalloc1(n,&nasm->oscatter_copy);
395: for (i=0; i<n; i++) {
396: nasm->oscatter[i] = oscatter[i];
397: VecScatterCopy(oscatter[i], &nasm->oscatter_copy[i]);
398: }
399: }
400: if (iscatter) {
401: PetscMalloc1(n,&nasm->iscatter);
402: for (i=0; i<n; i++) {
403: nasm->iscatter[i] = iscatter[i];
404: }
405: }
406: if (gscatter) {
407: PetscMalloc1(n,&nasm->gscatter);
408: for (i=0; i<n; i++) {
409: nasm->gscatter[i] = gscatter[i];
410: }
411: }
413: if (subsnes) {
414: PetscMalloc1(n,&nasm->subsnes);
415: for (i=0; i<n; i++) {
416: nasm->subsnes[i] = subsnes[i];
417: }
418: nasm->same_local_solves = PETSC_FALSE;
419: }
420: return(0);
421: }
423: /*@
424: SNESNASMGetSubdomains - Get the local subdomain context.
426: Not Collective
428: Input Parameters:
429: . SNES - the SNES context
431: Output Parameters:
432: + n - the number of local subdomains
433: . subsnes - solvers defined on the local subdomains
434: . iscatter - scatters into the nonoverlapping portions of the local subdomains
435: . oscatter - scatters into the overlapping portions of the local subdomains
436: - gscatter - scatters into the (ghosted) local vector of the local subdomain
438: Level: intermediate
440: .seealso: SNESNASM, SNESNASMSetSubdomains()
441: @*/
442: PetscErrorCode SNESNASMGetSubdomains(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
443: {
445: PetscErrorCode (*f)(SNES,PetscInt*,SNES**,VecScatter**,VecScatter**,VecScatter**);
448: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",&f);
449: if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
450: return(0);
451: }
453: static PetscErrorCode SNESNASMGetSubdomains_NASM(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
454: {
455: SNES_NASM *nasm = (SNES_NASM*)snes->data;
458: if (n) *n = nasm->n;
459: if (oscatter) *oscatter = nasm->oscatter;
460: if (iscatter) *iscatter = nasm->iscatter;
461: if (gscatter) *gscatter = nasm->gscatter;
462: if (subsnes) {
463: *subsnes = nasm->subsnes;
464: nasm->same_local_solves = PETSC_FALSE;
465: }
466: return(0);
467: }
469: /*@
470: SNESNASMGetSubdomainVecs - Get the processor-local subdomain vectors
472: Not Collective
474: Input Parameters:
475: . SNES - the SNES context
477: Output Parameters:
478: + n - the number of local subdomains
479: . x - The subdomain solution vector
480: . y - The subdomain step vector
481: . b - The subdomain RHS vector
482: - xl - The subdomain local vectors (ghosted)
484: Level: developer
486: .seealso: SNESNASM, SNESNASMGetSubdomains()
487: @*/
488: PetscErrorCode SNESNASMGetSubdomainVecs(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b, Vec **xl)
489: {
491: PetscErrorCode (*f)(SNES,PetscInt*,Vec**,Vec**,Vec**,Vec**);
494: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",&f);
495: if (f) {(f)(snes,n,x,y,b,xl);}
496: return(0);
497: }
499: static PetscErrorCode SNESNASMGetSubdomainVecs_NASM(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b,Vec **xl)
500: {
501: SNES_NASM *nasm = (SNES_NASM*)snes->data;
504: if (n) *n = nasm->n;
505: if (x) *x = nasm->x;
506: if (y) *y = nasm->y;
507: if (b) *b = nasm->b;
508: if (xl) *xl = nasm->xl;
509: return(0);
510: }
512: /*@
513: SNESNASMSetComputeFinalJacobian - Schedules the computation of the global and subdomain Jacobians upon convergence
515: Collective on SNES
517: Input Parameters:
518: + SNES - the SNES context
519: - flg - indication of whether to compute the Jacobians or not
521: Level: developer
523: Notes:
524: This is used almost exclusively in the implementation of ASPIN, where the converged subdomain and global Jacobian
525: is needed at each linear iteration.
527: .seealso: SNESNASM, SNESNASMGetSubdomains()
528: @*/
529: PetscErrorCode SNESNASMSetComputeFinalJacobian(SNES snes,PetscBool flg)
530: {
531: PetscErrorCode (*f)(SNES,PetscBool);
535: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",&f);
536: if (f) {(f)(snes,flg);}
537: return(0);
538: }
540: static PetscErrorCode SNESNASMSetComputeFinalJacobian_NASM(SNES snes,PetscBool flg)
541: {
542: SNES_NASM *nasm = (SNES_NASM*)snes->data;
545: nasm->finaljacobian = flg;
546: return(0);
547: }
549: /*@
550: SNESNASMSetDamping - Sets the update damping for NASM
552: Logically collective on SNES
554: Input Parameters:
555: + SNES - the SNES context
556: - dmp - damping
558: Level: intermediate
560: Notes:
561: The new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
563: .seealso: SNESNASM, SNESNASMGetDamping()
564: @*/
565: PetscErrorCode SNESNASMSetDamping(SNES snes,PetscReal dmp)
566: {
567: PetscErrorCode (*f)(SNES,PetscReal);
571: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetDamping_C",(void (**)(void))&f);
572: if (f) {(f)(snes,dmp);}
573: return(0);
574: }
576: static PetscErrorCode SNESNASMSetDamping_NASM(SNES snes,PetscReal dmp)
577: {
578: SNES_NASM *nasm = (SNES_NASM*)snes->data;
581: nasm->damping = dmp;
582: return(0);
583: }
585: /*@
586: SNESNASMGetDamping - Gets the update damping for NASM
588: Not Collective
590: Input Parameters:
591: + SNES - the SNES context
592: - dmp - damping
594: Level: intermediate
596: .seealso: SNESNASM, SNESNASMSetDamping()
597: @*/
598: PetscErrorCode SNESNASMGetDamping(SNES snes,PetscReal *dmp)
599: {
603: PetscUseMethod(snes,"SNESNASMGetDamping_C",(SNES,PetscReal*),(snes,dmp));
604: return(0);
605: }
607: static PetscErrorCode SNESNASMGetDamping_NASM(SNES snes,PetscReal *dmp)
608: {
609: SNES_NASM *nasm = (SNES_NASM*)snes->data;
612: *dmp = nasm->damping;
613: return(0);
614: }
617: /*
618: Input Parameters:
619: + snes - The solver
620: . B - The RHS vector
621: - X - The initial guess
623: Output Parameters:
624: . Y - The solution update
626: TODO: All scatters should be packed into one
627: */
628: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X)
629: {
630: SNES_NASM *nasm = (SNES_NASM*)snes->data;
631: SNES subsnes;
632: PetscInt i;
633: PetscReal dmp;
635: Vec Xl,Bl,Yl,Xlloc;
636: VecScatter iscat,oscat,gscat,oscat_copy;
637: DM dm,subdm;
638: PCASMType type;
641: SNESNASMGetType(snes,&type);
642: SNESGetDM(snes,&dm);
643: VecSet(Y,0);
644: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
645: for (i=0; i<nasm->n; i++) {
646: /* scatter the solution to the global solution and the local solution */
647: Xl = nasm->x[i];
648: Xlloc = nasm->xl[i];
649: oscat = nasm->oscatter[i];
650: oscat_copy = nasm->oscatter_copy[i];
651: gscat = nasm->gscatter[i];
652: VecScatterBegin(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
653: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
654: if (B) {
655: /* scatter the RHS to the local RHS */
656: Bl = nasm->b[i];
657: VecScatterBegin(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
658: }
659: }
660: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
663: if (nasm->eventsubsolve) {PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);}
664: for (i=0; i<nasm->n; i++) {
665: Xl = nasm->x[i];
666: Xlloc = nasm->xl[i];
667: Yl = nasm->y[i];
668: subsnes = nasm->subsnes[i];
669: SNESGetDM(subsnes,&subdm);
670: iscat = nasm->iscatter[i];
671: oscat = nasm->oscatter[i];
672: oscat_copy = nasm->oscatter_copy[i];
673: gscat = nasm->gscatter[i];
674: VecScatterEnd(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
675: VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
676: if (B) {
677: Bl = nasm->b[i];
678: VecScatterEnd(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
679: } else Bl = NULL;
681: DMSubDomainRestrict(dm,oscat,gscat,subdm);
682: VecCopy(Xl,Yl);
683: SNESSolve(subsnes,Bl,Xl);
684: VecAYPX(Yl,-1.0,Xl);
685: VecScale(Yl, nasm->damping);
686: if (type == PC_ASM_BASIC) {
687: VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
688: VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
689: } else if (type == PC_ASM_RESTRICT) {
690: VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
691: VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
692: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
693: }
694: if (nasm->eventsubsolve) {PetscLogEventEnd(nasm->eventsubsolve,snes,0,0,0);}
695: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
696: if (nasm->weight_set) {
697: VecPointwiseMult(Y,Y,nasm->weight);
698: }
699: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
700: SNESNASMGetDamping(snes,&dmp);
701: VecAXPY(X,dmp,Y);
702: return(0);
703: }
705: static PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
706: {
707: Vec X = Xfinal;
708: SNES_NASM *nasm = (SNES_NASM*)snes->data;
709: SNES subsnes;
710: PetscInt i,lag = 1;
712: Vec Xlloc,Xl,Fl,F;
713: VecScatter oscat,gscat;
714: DM dm,subdm;
717: if (nasm->fjtype == 2) X = nasm->xinit;
718: F = snes->vec_func;
719: if (snes->normschedule == SNES_NORM_NONE) {SNESComputeFunction(snes,X,F);}
720: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
721: SNESGetDM(snes,&dm);
722: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
723: if (nasm->fjtype != 1) {
724: for (i=0; i<nasm->n; i++) {
725: Xlloc = nasm->xl[i];
726: gscat = nasm->gscatter[i];
727: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
728: }
729: }
730: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
731: for (i=0; i<nasm->n; i++) {
732: Fl = nasm->subsnes[i]->vec_func;
733: Xl = nasm->x[i];
734: Xlloc = nasm->xl[i];
735: subsnes = nasm->subsnes[i];
736: oscat = nasm->oscatter[i];
737: gscat = nasm->gscatter[i];
738: if (nasm->fjtype != 1) {VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);}
739: SNESGetDM(subsnes,&subdm);
740: DMSubDomainRestrict(dm,oscat,gscat,subdm);
741: if (nasm->fjtype != 1) {
742: DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
743: DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
744: }
745: if (subsnes->lagjacobian == -1) subsnes->lagjacobian = -2;
746: else if (subsnes->lagjacobian > 1) lag = subsnes->lagjacobian;
747: SNESComputeFunction(subsnes,Xl,Fl);
748: SNESComputeJacobian(subsnes,Xl,subsnes->jacobian,subsnes->jacobian_pre);
749: if (lag > 1) subsnes->lagjacobian = lag;
750: }
751: return(0);
752: }
754: static PetscErrorCode SNESSolve_NASM(SNES snes)
755: {
756: Vec F;
757: Vec X;
758: Vec B;
759: Vec Y;
760: PetscInt i;
761: PetscReal fnorm = 0.0;
762: PetscErrorCode ierr;
763: SNESNormSchedule normschedule;
764: SNES_NASM *nasm = (SNES_NASM*)snes->data;
768: 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);
770: PetscCitationsRegister(SNESCitation,&SNEScite);
771: X = snes->vec_sol;
772: Y = snes->vec_sol_update;
773: F = snes->vec_func;
774: B = snes->vec_rhs;
776: PetscObjectSAWsTakeAccess((PetscObject)snes);
777: snes->iter = 0;
778: snes->norm = 0.;
779: PetscObjectSAWsGrantAccess((PetscObject)snes);
780: snes->reason = SNES_CONVERGED_ITERATING;
781: SNESGetNormSchedule(snes, &normschedule);
782: if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
783: /* compute the initial function and preconditioned update delX */
784: if (!snes->vec_func_init_set) {
785: SNESComputeFunction(snes,X,F);
786: } else snes->vec_func_init_set = PETSC_FALSE;
788: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
789: SNESCheckFunctionNorm(snes,fnorm);
790: PetscObjectSAWsTakeAccess((PetscObject)snes);
791: snes->iter = 0;
792: snes->norm = fnorm;
793: PetscObjectSAWsGrantAccess((PetscObject)snes);
794: SNESLogConvergenceHistory(snes,snes->norm,0);
795: SNESMonitor(snes,0,snes->norm);
797: /* test convergence */
798: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
799: if (snes->reason) return(0);
800: } else {
801: PetscObjectSAWsGrantAccess((PetscObject)snes);
802: SNESLogConvergenceHistory(snes,snes->norm,0);
803: SNESMonitor(snes,0,snes->norm);
804: }
806: /* Call general purpose update function */
807: if (snes->ops->update) {
808: (*snes->ops->update)(snes, snes->iter);
809: }
810: /* copy the initial solution over for later */
811: if (nasm->fjtype == 2) {VecCopy(X,nasm->xinit);}
813: for (i=0; i < snes->max_its; i++) {
814: SNESNASMSolveLocal_Private(snes,B,Y,X);
815: if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
816: SNESComputeFunction(snes,X,F);
817: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
818: SNESCheckFunctionNorm(snes,fnorm);
819: }
820: /* Monitor convergence */
821: PetscObjectSAWsTakeAccess((PetscObject)snes);
822: snes->iter = i+1;
823: snes->norm = fnorm;
824: PetscObjectSAWsGrantAccess((PetscObject)snes);
825: SNESLogConvergenceHistory(snes,snes->norm,0);
826: SNESMonitor(snes,snes->iter,snes->norm);
827: /* Test for convergence */
828: if (normschedule == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);}
829: if (snes->reason) break;
830: /* Call general purpose update function */
831: if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
832: }
833: if (nasm->finaljacobian) {
834: SNESNASMComputeFinalJacobian_Private(snes,X);
835: SNESCheckJacobianDomainerror(snes);
836: }
837: if (normschedule == SNES_NORM_ALWAYS) {
838: if (i == snes->max_its) {
839: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
840: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
841: }
842: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */
843: return(0);
844: }
846: /*MC
847: SNESNASM - Nonlinear Additive Schwarz
849: Options Database:
850: + -snes_nasm_log - enable logging events for the communication and solve stages
851: . -snes_nasm_type <basic,restrict> - type of subdomain update used
852: . -snes_asm_damping <dmp> - the new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
853: . -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
854: . -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> - pick state the jacobian is calculated at
855: . -sub_snes_ - options prefix of the subdomain nonlinear solves
856: . -sub_ksp_ - options prefix of the subdomain Krylov solver
857: - -sub_pc_ - options prefix of the subdomain preconditioner
859: Level: advanced
861: 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
862: false and SNESView() and -snes_view do not display a KSP object. However, if the flag nasm->finaljacobian is set (for example, if
863: NASM is used as a nonlinear preconditioner for KSPASPIN) then SNESSetUpMatrices() is called to generate the Jacobian (needed by KSPASPIN)
864: 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
865: used by SNESASPIN they share the same Jacobian matrices because SNESSetUp() (called on the outer SNES KSPASPIN) causes the inner SNES
866: object (in this case SNESNASM) to inherit the outer Jacobian matrices.
868: References:
869: . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
870: SIAM Review, 57(4), 2015
872: .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types), SNESNASMSetType(), SNESNASMGetType(), SNESNASMSetSubdomains(), SNESNASMGetSubdomains(), SNESNASMGetSubdomainVecs(), SNESNASMSetComputeFinalJacobian(), SNESNASMSetDamping(), SNESNASMGetDamping()
873: M*/
875: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
876: {
877: SNES_NASM *nasm;
881: PetscNewLog(snes,&nasm);
882: snes->data = (void*)nasm;
884: nasm->n = PETSC_DECIDE;
885: nasm->subsnes = 0;
886: nasm->x = 0;
887: nasm->xl = 0;
888: nasm->y = 0;
889: nasm->b = 0;
890: nasm->oscatter = 0;
891: nasm->oscatter_copy = 0;
892: nasm->iscatter = 0;
893: nasm->gscatter = 0;
894: nasm->damping = 1.;
896: nasm->type = PC_ASM_BASIC;
897: nasm->finaljacobian = PETSC_FALSE;
898: nasm->same_local_solves = PETSC_TRUE;
899: nasm->weight_set = PETSC_FALSE;
901: snes->ops->destroy = SNESDestroy_NASM;
902: snes->ops->setup = SNESSetUp_NASM;
903: snes->ops->setfromoptions = SNESSetFromOptions_NASM;
904: snes->ops->view = SNESView_NASM;
905: snes->ops->solve = SNESSolve_NASM;
906: snes->ops->reset = SNESReset_NASM;
908: snes->usesksp = PETSC_FALSE;
909: snes->usesnpc = PETSC_FALSE;
911: snes->alwayscomputesfinalresidual = PETSC_FALSE;
913: nasm->fjtype = 0;
914: nasm->xinit = NULL;
915: nasm->eventrestrictinterp = 0;
916: nasm->eventsubsolve = 0;
918: if (!snes->tolerancesset) {
919: snes->max_its = 10000;
920: snes->max_funcs = 10000;
921: }
923: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetType_C",SNESNASMSetType_NASM);
924: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetType_C",SNESNASMGetType_NASM);
925: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
926: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
927: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
928: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
929: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
930: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
931: return(0);
932: }
934: /*@
935: SNESNASMGetSNES - Gets a subsolver
937: Not collective
939: Input Parameters:
940: + snes - the SNES context
941: - i - the number of the subsnes to get
943: Output Parameters:
944: . subsnes - the subsolver context
946: Level: intermediate
948: .seealso: SNESNASM, SNESNASMGetNumber()
949: @*/
950: PetscErrorCode SNESNASMGetSNES(SNES snes,PetscInt i,SNES *subsnes)
951: {
952: SNES_NASM *nasm = (SNES_NASM*)snes->data;
955: if (i < 0 || i >= nasm->n) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"No such subsolver");
956: *subsnes = nasm->subsnes[i];
957: return(0);
958: }
960: /*@
961: SNESNASMGetNumber - Gets number of subsolvers
963: Not collective
965: Input Parameters:
966: . snes - the SNES context
968: Output Parameters:
969: . n - the number of subsolvers
971: Level: intermediate
973: .seealso: SNESNASM, SNESNASMGetSNES()
974: @*/
975: PetscErrorCode SNESNASMGetNumber(SNES snes,PetscInt *n)
976: {
977: SNES_NASM *nasm = (SNES_NASM*)snes->data;
980: *n = nasm->n;
981: return(0);
982: }
984: /*@
985: SNESNASMSetWeight - Sets weight to use when adding overlapping updates
987: Collective
989: Input Parameters:
990: + snes - the SNES context
991: - weight - the weights to use (typically 1/N for each dof, where N is the number of patches it appears in)
993: Level: intermediate
995: .seealso: SNESNASM
996: @*/
997: PetscErrorCode SNESNASMSetWeight(SNES snes,Vec weight)
998: {
999: SNES_NASM *nasm = (SNES_NASM*)snes->data;
1004: VecDestroy(&nasm->weight);
1005: nasm->weight_set = PETSC_TRUE;
1006: nasm->weight = weight;
1007: PetscObjectReference((PetscObject)nasm->weight);
1009: return(0);
1010: }