Actual source code: nasm.c
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 weight_set; /* use a weight in the overlap updates */
22: /* logging events */
23: PetscLogEvent eventrestrictinterp;
24: PetscLogEvent eventsubsolve;
26: PetscInt fjtype; /* type of computed jacobian */
27: Vec xinit; /* initial solution in case the final jacobian type is computed as first */
28: } SNES_NASM;
30: const char *const SNESNASMTypes[] = {"NONE","RESTRICT","INTERPOLATE","BASIC","PCASMType","PC_ASM_",NULL};
31: const char *const SNESNASMFJTypes[] = {"FINALOUTER","FINALINNER","INITIAL"};
33: static PetscErrorCode SNESReset_NASM(SNES snes)
34: {
35: SNES_NASM *nasm = (SNES_NASM*)snes->data;
36: PetscInt i;
38: for (i=0; i<nasm->n; i++) {
39: if (nasm->xl) VecDestroy(&nasm->xl[i]);
40: if (nasm->x) VecDestroy(&nasm->x[i]);
41: if (nasm->y) VecDestroy(&nasm->y[i]);
42: if (nasm->b) VecDestroy(&nasm->b[i]);
44: if (nasm->subsnes) SNESDestroy(&nasm->subsnes[i]);
45: if (nasm->oscatter) VecScatterDestroy(&nasm->oscatter[i]);
46: if (nasm->oscatter_copy) VecScatterDestroy(&nasm->oscatter_copy[i]);
47: if (nasm->iscatter) VecScatterDestroy(&nasm->iscatter[i]);
48: if (nasm->gscatter) VecScatterDestroy(&nasm->gscatter[i]);
49: }
51: PetscFree(nasm->x);
52: PetscFree(nasm->xl);
53: PetscFree(nasm->y);
54: PetscFree(nasm->b);
56: if (nasm->xinit) VecDestroy(&nasm->xinit);
58: PetscFree(nasm->subsnes);
59: PetscFree(nasm->oscatter);
60: PetscFree(nasm->oscatter_copy);
61: PetscFree(nasm->iscatter);
62: PetscFree(nasm->gscatter);
64: if (nasm->weight_set) {
65: VecDestroy(&nasm->weight);
66: }
68: nasm->eventrestrictinterp = 0;
69: nasm->eventsubsolve = 0;
70: return 0;
71: }
73: static PetscErrorCode SNESDestroy_NASM(SNES snes)
74: {
75: SNESReset_NASM(snes);
76: PetscFree(snes->data);
77: return 0;
78: }
80: static PetscErrorCode DMGlobalToLocalSubDomainDirichletHook_Private(DM dm,Vec g,InsertMode mode,Vec l,void *ctx)
81: {
82: Vec bcs = (Vec)ctx;
84: VecCopy(bcs,l);
85: return 0;
86: }
88: static PetscErrorCode SNESSetUp_NASM(SNES snes)
89: {
90: SNES_NASM *nasm = (SNES_NASM*)snes->data;
91: DM dm,subdm;
92: DM *subdms;
93: PetscInt i;
94: const char *optionsprefix;
95: Vec F;
96: PetscMPIInt size;
97: KSP ksp;
98: PC pc;
100: if (!nasm->subsnes) {
101: SNESGetDM(snes,&dm);
102: if (dm) {
103: nasm->usesdm = PETSC_TRUE;
104: DMCreateDomainDecomposition(dm,&nasm->n,NULL,NULL,NULL,&subdms);
106: DMCreateDomainDecompositionScatters(dm,nasm->n,subdms,&nasm->iscatter,&nasm->oscatter,&nasm->gscatter);
107: PetscMalloc1(nasm->n, &nasm->oscatter_copy);
108: for (i=0; i<nasm->n; i++) {
109: VecScatterCopy(nasm->oscatter[i], &nasm->oscatter_copy[i]);
110: }
112: SNESGetOptionsPrefix(snes, &optionsprefix);
113: PetscMalloc1(nasm->n,&nasm->subsnes);
114: for (i=0; i<nasm->n; i++) {
115: SNESCreate(PETSC_COMM_SELF,&nasm->subsnes[i]);
116: PetscObjectIncrementTabLevel((PetscObject)nasm->subsnes[i], (PetscObject)snes, 1);
117: SNESAppendOptionsPrefix(nasm->subsnes[i],optionsprefix);
118: SNESAppendOptionsPrefix(nasm->subsnes[i],"sub_");
119: SNESSetDM(nasm->subsnes[i],subdms[i]);
120: MPI_Comm_size(PetscObjectComm((PetscObject)nasm->subsnes[i]),&size);
121: if (size == 1) {
122: SNESGetKSP(nasm->subsnes[i],&ksp);
123: KSPGetPC(ksp,&pc);
124: KSPSetType(ksp,KSPPREONLY);
125: PCSetType(pc,PCLU);
126: }
127: SNESSetFromOptions(nasm->subsnes[i]);
128: DMDestroy(&subdms[i]);
129: }
130: PetscFree(subdms);
131: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Cannot construct local problems automatically without a DM!");
132: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set subproblems manually if there is no DM!");
133: /* allocate the global vectors */
134: if (!nasm->x) {
135: PetscCalloc1(nasm->n,&nasm->x);
136: }
137: if (!nasm->xl) {
138: PetscCalloc1(nasm->n,&nasm->xl);
139: }
140: if (!nasm->y) {
141: PetscCalloc1(nasm->n,&nasm->y);
142: }
143: if (!nasm->b) {
144: PetscCalloc1(nasm->n,&nasm->b);
145: }
147: for (i=0; i<nasm->n; i++) {
148: SNESGetFunction(nasm->subsnes[i],&F,NULL,NULL);
149: if (!nasm->x[i]) VecDuplicate(F,&nasm->x[i]);
150: if (!nasm->y[i]) VecDuplicate(F,&nasm->y[i]);
151: if (!nasm->b[i]) VecDuplicate(F,&nasm->b[i]);
152: if (!nasm->xl[i]) {
153: SNESGetDM(nasm->subsnes[i],&subdm);
154: DMCreateLocalVector(subdm,&nasm->xl[i]);
155: DMGlobalToLocalHookAdd(subdm,DMGlobalToLocalSubDomainDirichletHook_Private,NULL,nasm->xl[i]);
156: }
157: }
158: if (nasm->finaljacobian) {
159: SNESSetUpMatrices(snes);
160: if (nasm->fjtype == 2) {
161: VecDuplicate(snes->vec_sol,&nasm->xinit);
162: }
163: for (i=0; i<nasm->n;i++) {
164: SNESSetUpMatrices(nasm->subsnes[i]);
165: }
166: }
167: return 0;
168: }
170: static PetscErrorCode SNESSetFromOptions_NASM(PetscOptionItems *PetscOptionsObject,SNES snes)
171: {
172: PCASMType asmtype;
173: PetscBool flg,monflg;
174: SNES_NASM *nasm = (SNES_NASM*)snes->data;
176: PetscOptionsHead(PetscOptionsObject,"Nonlinear Additive Schwarz options");
177: PetscOptionsEnum("-snes_nasm_type","Type of restriction/extension","",SNESNASMTypes,(PetscEnum)nasm->type,(PetscEnum*)&asmtype,&flg);
178: if (flg) SNESNASMSetType(snes,asmtype);
179: flg = PETSC_FALSE;
180: monflg = PETSC_TRUE;
181: 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);
182: if (flg) SNESNASMSetDamping(snes,nasm->damping);
183: PetscOptionsDeprecated("-snes_nasm_sub_view",NULL,"3.15","Use -snes_view ::ascii_info_detail");
184: PetscOptionsBool("-snes_nasm_finaljacobian","Compute the global jacobian of the final iterate (for ASPIN)","",nasm->finaljacobian,&nasm->finaljacobian,NULL);
185: PetscOptionsEList("-snes_nasm_finaljacobian_type","The type of the final jacobian computed.","",SNESNASMFJTypes,3,SNESNASMFJTypes[0],&nasm->fjtype,NULL);
186: PetscOptionsBool("-snes_nasm_log","Log times for subSNES solves and restriction","",monflg,&monflg,&flg);
187: if (flg) {
188: PetscLogEventRegister("SNESNASMSubSolve",((PetscObject)snes)->classid,&nasm->eventsubsolve);
189: PetscLogEventRegister("SNESNASMRestrict",((PetscObject)snes)->classid,&nasm->eventrestrictinterp);
190: }
191: PetscOptionsTail();
192: return 0;
193: }
195: static PetscErrorCode SNESView_NASM(SNES snes, PetscViewer viewer)
196: {
197: SNES_NASM *nasm = (SNES_NASM*)snes->data;
198: PetscMPIInt rank,size;
199: PetscInt i,N,bsz;
200: PetscBool iascii,isstring;
201: PetscViewer sviewer;
202: MPI_Comm comm;
203: PetscViewerFormat format;
204: const char *prefix;
206: PetscObjectGetComm((PetscObject)snes,&comm);
207: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
208: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
209: MPI_Comm_rank(comm,&rank);
210: MPI_Comm_size(comm,&size);
211: MPIU_Allreduce(&nasm->n,&N,1,MPIU_INT,MPI_SUM,comm);
212: if (iascii) {
213: PetscViewerASCIIPrintf(viewer, " total subdomain blocks = %D\n",N);
214: PetscViewerGetFormat(viewer,&format);
215: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
216: if (nasm->subsnes) {
217: PetscViewerASCIIPrintf(viewer," Local solver information for first block on rank 0:\n");
218: SNESGetOptionsPrefix(snes,&prefix);
219: PetscViewerASCIIPrintf(viewer," Use -%ssnes_view ::ascii_info_detail to display information for all blocks\n",prefix?prefix:"");
220: PetscViewerASCIIPushTab(viewer);
221: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
222: if (rank == 0) {
223: PetscViewerASCIIPushTab(viewer);
224: SNESView(nasm->subsnes[0],sviewer);
225: PetscViewerASCIIPopTab(viewer);
226: }
227: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
228: PetscViewerASCIIPopTab(viewer);
229: }
230: } else {
231: /* print the solver on each block */
232: PetscViewerASCIIPushSynchronized(viewer);
233: PetscViewerASCIISynchronizedPrintf(viewer," [%d] number of local blocks = %D\n",(int)rank,nasm->n);
234: PetscViewerFlush(viewer);
235: PetscViewerASCIIPopSynchronized(viewer);
236: PetscViewerASCIIPrintf(viewer," Local solver information for each block is in the following SNES objects:\n");
237: PetscViewerASCIIPushTab(viewer);
238: PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
239: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
240: for (i=0; i<nasm->n; i++) {
241: VecGetLocalSize(nasm->x[i],&bsz);
242: PetscViewerASCIIPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);
243: SNESView(nasm->subsnes[i],sviewer);
244: PetscViewerASCIIPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
245: }
246: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
247: PetscViewerFlush(viewer);
248: PetscViewerASCIIPopTab(viewer);
249: }
250: } else if (isstring) {
251: PetscViewerStringSPrintf(viewer," blocks=%D,type=%s",N,SNESNASMTypes[nasm->type]);
252: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
253: if (nasm->subsnes && rank == 0) SNESView(nasm->subsnes[0],sviewer);
254: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
255: }
256: return 0;
257: }
259: /*@
260: SNESNASMSetType - Set the type of subdomain update used
262: Logically Collective on SNES
264: Input Parameters:
265: + SNES - the SNES context
266: - type - the type of update, PC_ASM_BASIC or PC_ASM_RESTRICT
268: Level: intermediate
270: .seealso: SNESNASM, SNESNASMGetType(), PCASMSetType()
271: @*/
272: PetscErrorCode SNESNASMSetType(SNES snes,PCASMType type)
273: {
274: PetscErrorCode (*f)(SNES,PCASMType);
276: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetType_C",&f);
277: if (f) (f)(snes,type);
278: return 0;
279: }
281: static PetscErrorCode SNESNASMSetType_NASM(SNES snes,PCASMType type)
282: {
283: SNES_NASM *nasm = (SNES_NASM*)snes->data;
286: nasm->type = type;
287: return 0;
288: }
290: /*@
291: SNESNASMGetType - Get the type of subdomain update used
293: Logically Collective on SNES
295: Input Parameters:
296: . SNES - the SNES context
298: Output Parameters:
299: . type - the type of update
301: Level: intermediate
303: .seealso: SNESNASM, SNESNASMSetType(), PCASMGetType()
304: @*/
305: PetscErrorCode SNESNASMGetType(SNES snes,PCASMType *type)
306: {
307: PetscUseMethod(snes,"SNESNASMGetType_C",(SNES,PCASMType*),(snes,type));
308: return 0;
309: }
311: static PetscErrorCode SNESNASMGetType_NASM(SNES snes,PCASMType *type)
312: {
313: SNES_NASM *nasm = (SNES_NASM*)snes->data;
315: *type = nasm->type;
316: return 0;
317: }
319: /*@
320: SNESNASMSetSubdomains - Manually Set the context required to restrict and solve subdomain problems.
322: Not Collective
324: Input Parameters:
325: + SNES - the SNES context
326: . n - the number of local subdomains
327: . subsnes - solvers defined on the local subdomains
328: . iscatter - scatters into the nonoverlapping portions of the local subdomains
329: . oscatter - scatters into the overlapping portions of the local subdomains
330: - gscatter - scatters into the (ghosted) local vector of the local subdomain
332: Level: intermediate
334: .seealso: SNESNASM, SNESNASMGetSubdomains()
335: @*/
336: PetscErrorCode SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
337: {
338: PetscErrorCode (*f)(SNES,PetscInt,SNES*,VecScatter*,VecScatter*,VecScatter*);
340: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",&f);
341: if (f) (f)(snes,n,subsnes,iscatter,oscatter,gscatter);
342: return 0;
343: }
345: static PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
346: {
347: PetscInt i;
348: SNES_NASM *nasm = (SNES_NASM*)snes->data;
352: /* tear down the previously set things */
353: SNESReset(snes);
355: nasm->n = n;
356: if (oscatter) {
357: for (i=0; i<n; i++) PetscObjectReference((PetscObject)oscatter[i]);
358: }
359: if (iscatter) {
360: for (i=0; i<n; i++) PetscObjectReference((PetscObject)iscatter[i]);
361: }
362: if (gscatter) {
363: for (i=0; i<n; i++) PetscObjectReference((PetscObject)gscatter[i]);
364: }
365: if (oscatter) {
366: PetscMalloc1(n,&nasm->oscatter);
367: PetscMalloc1(n,&nasm->oscatter_copy);
368: for (i=0; i<n; i++) {
369: nasm->oscatter[i] = oscatter[i];
370: VecScatterCopy(oscatter[i], &nasm->oscatter_copy[i]);
371: }
372: }
373: if (iscatter) {
374: PetscMalloc1(n,&nasm->iscatter);
375: for (i=0; i<n; i++) {
376: nasm->iscatter[i] = iscatter[i];
377: }
378: }
379: if (gscatter) {
380: PetscMalloc1(n,&nasm->gscatter);
381: for (i=0; i<n; i++) {
382: nasm->gscatter[i] = gscatter[i];
383: }
384: }
386: if (subsnes) {
387: PetscMalloc1(n,&nasm->subsnes);
388: for (i=0; i<n; i++) {
389: nasm->subsnes[i] = subsnes[i];
390: }
391: }
392: return 0;
393: }
395: /*@
396: SNESNASMGetSubdomains - Get the local subdomain context.
398: Not Collective
400: Input Parameter:
401: . SNES - the SNES context
403: Output Parameters:
404: + n - the number of local subdomains
405: . subsnes - solvers defined on the local subdomains
406: . iscatter - scatters into the nonoverlapping portions of the local subdomains
407: . oscatter - scatters into the overlapping portions of the local subdomains
408: - gscatter - scatters into the (ghosted) local vector of the local subdomain
410: Level: intermediate
412: .seealso: SNESNASM, SNESNASMSetSubdomains()
413: @*/
414: PetscErrorCode SNESNASMGetSubdomains(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
415: {
416: PetscErrorCode (*f)(SNES,PetscInt*,SNES**,VecScatter**,VecScatter**,VecScatter**);
418: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",&f);
419: if (f) (f)(snes,n,subsnes,iscatter,oscatter,gscatter);
420: return 0;
421: }
423: static PetscErrorCode SNESNASMGetSubdomains_NASM(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
424: {
425: SNES_NASM *nasm = (SNES_NASM*)snes->data;
427: if (n) *n = nasm->n;
428: if (oscatter) *oscatter = nasm->oscatter;
429: if (iscatter) *iscatter = nasm->iscatter;
430: if (gscatter) *gscatter = nasm->gscatter;
431: if (subsnes) *subsnes = nasm->subsnes;
432: return 0;
433: }
435: /*@
436: SNESNASMGetSubdomainVecs - Get the processor-local subdomain vectors
438: Not Collective
440: Input Parameter:
441: . SNES - the SNES context
443: Output Parameters:
444: + n - the number of local subdomains
445: . x - The subdomain solution vector
446: . y - The subdomain step vector
447: . b - The subdomain RHS vector
448: - xl - The subdomain local vectors (ghosted)
450: Level: developer
452: .seealso: SNESNASM, SNESNASMGetSubdomains()
453: @*/
454: PetscErrorCode SNESNASMGetSubdomainVecs(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b, Vec **xl)
455: {
456: PetscErrorCode (*f)(SNES,PetscInt*,Vec**,Vec**,Vec**,Vec**);
458: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",&f);
459: if (f) (f)(snes,n,x,y,b,xl);
460: return 0;
461: }
463: static PetscErrorCode SNESNASMGetSubdomainVecs_NASM(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b,Vec **xl)
464: {
465: SNES_NASM *nasm = (SNES_NASM*)snes->data;
467: if (n) *n = nasm->n;
468: if (x) *x = nasm->x;
469: if (y) *y = nasm->y;
470: if (b) *b = nasm->b;
471: if (xl) *xl = nasm->xl;
472: return 0;
473: }
475: /*@
476: SNESNASMSetComputeFinalJacobian - Schedules the computation of the global and subdomain Jacobians upon convergence
478: Collective on SNES
480: Input Parameters:
481: + SNES - the SNES context
482: - flg - indication of whether to compute the Jacobians or not
484: Level: developer
486: Notes:
487: This is used almost exclusively in the implementation of ASPIN, where the converged subdomain and global Jacobian
488: is needed at each linear iteration.
490: .seealso: SNESNASM, SNESNASMGetSubdomains()
491: @*/
492: PetscErrorCode SNESNASMSetComputeFinalJacobian(SNES snes,PetscBool flg)
493: {
494: PetscErrorCode (*f)(SNES,PetscBool);
496: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",&f);
497: if (f) (f)(snes,flg);
498: return 0;
499: }
501: static PetscErrorCode SNESNASMSetComputeFinalJacobian_NASM(SNES snes,PetscBool flg)
502: {
503: SNES_NASM *nasm = (SNES_NASM*)snes->data;
505: nasm->finaljacobian = flg;
506: return 0;
507: }
509: /*@
510: SNESNASMSetDamping - Sets the update damping for NASM
512: Logically collective on SNES
514: Input Parameters:
515: + SNES - the SNES context
516: - dmp - damping
518: Level: intermediate
520: Notes:
521: The new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
523: .seealso: SNESNASM, SNESNASMGetDamping()
524: @*/
525: PetscErrorCode SNESNASMSetDamping(SNES snes,PetscReal dmp)
526: {
527: PetscErrorCode (*f)(SNES,PetscReal);
529: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetDamping_C",(void (**)(void))&f);
530: if (f) (f)(snes,dmp);
531: return 0;
532: }
534: static PetscErrorCode SNESNASMSetDamping_NASM(SNES snes,PetscReal dmp)
535: {
536: SNES_NASM *nasm = (SNES_NASM*)snes->data;
538: nasm->damping = dmp;
539: return 0;
540: }
542: /*@
543: SNESNASMGetDamping - Gets the update damping for NASM
545: Not Collective
547: Input Parameters:
548: + SNES - the SNES context
549: - dmp - damping
551: Level: intermediate
553: .seealso: SNESNASM, SNESNASMSetDamping()
554: @*/
555: PetscErrorCode SNESNASMGetDamping(SNES snes,PetscReal *dmp)
556: {
557: PetscUseMethod(snes,"SNESNASMGetDamping_C",(SNES,PetscReal*),(snes,dmp));
558: return 0;
559: }
561: static PetscErrorCode SNESNASMGetDamping_NASM(SNES snes,PetscReal *dmp)
562: {
563: SNES_NASM *nasm = (SNES_NASM*)snes->data;
565: *dmp = nasm->damping;
566: return 0;
567: }
569: /*
570: Input Parameters:
571: + snes - The solver
572: . B - The RHS vector
573: - X - The initial guess
575: Output Parameters:
576: . Y - The solution update
578: TODO: All scatters should be packed into one
579: */
580: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X)
581: {
582: SNES_NASM *nasm = (SNES_NASM*)snes->data;
583: SNES subsnes;
584: PetscInt i;
585: PetscReal dmp;
586: Vec Xl,Bl,Yl,Xlloc;
587: VecScatter iscat,oscat,gscat,oscat_copy;
588: DM dm,subdm;
589: PCASMType type;
591: SNESNASMGetType(snes,&type);
592: SNESGetDM(snes,&dm);
593: VecSet(Y,0);
594: if (nasm->eventrestrictinterp) PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);
595: for (i=0; i<nasm->n; i++) {
596: /* scatter the solution to the global solution and the local solution */
597: Xl = nasm->x[i];
598: Xlloc = nasm->xl[i];
599: oscat = nasm->oscatter[i];
600: oscat_copy = nasm->oscatter_copy[i];
601: gscat = nasm->gscatter[i];
602: VecScatterBegin(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
603: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
604: if (B) {
605: /* scatter the RHS to the local RHS */
606: Bl = nasm->b[i];
607: VecScatterBegin(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
608: }
609: }
610: if (nasm->eventrestrictinterp) PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);
612: if (nasm->eventsubsolve) PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);
613: for (i=0; i<nasm->n; i++) {
614: Xl = nasm->x[i];
615: Xlloc = nasm->xl[i];
616: Yl = nasm->y[i];
617: subsnes = nasm->subsnes[i];
618: SNESGetDM(subsnes,&subdm);
619: iscat = nasm->iscatter[i];
620: oscat = nasm->oscatter[i];
621: oscat_copy = nasm->oscatter_copy[i];
622: gscat = nasm->gscatter[i];
623: VecScatterEnd(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
624: VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
625: if (B) {
626: Bl = nasm->b[i];
627: VecScatterEnd(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
628: } else Bl = NULL;
630: DMSubDomainRestrict(dm,oscat,gscat,subdm);
631: VecCopy(Xl,Yl);
632: SNESSolve(subsnes,Bl,Xl);
633: VecAYPX(Yl,-1.0,Xl);
634: VecScale(Yl, nasm->damping);
635: if (type == PC_ASM_BASIC) {
636: VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
637: VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
638: } else if (type == PC_ASM_RESTRICT) {
639: VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
640: VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
641: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
642: }
643: if (nasm->eventsubsolve) PetscLogEventEnd(nasm->eventsubsolve,snes,0,0,0);
644: if (nasm->eventrestrictinterp) PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);
645: if (nasm->weight_set) {
646: VecPointwiseMult(Y,Y,nasm->weight);
647: }
648: if (nasm->eventrestrictinterp) PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);
649: SNESNASMGetDamping(snes,&dmp);
650: VecAXPY(X,dmp,Y);
651: return 0;
652: }
654: static PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
655: {
656: Vec X = Xfinal;
657: SNES_NASM *nasm = (SNES_NASM*)snes->data;
658: SNES subsnes;
659: PetscInt i,lag = 1;
660: Vec Xlloc,Xl,Fl,F;
661: VecScatter oscat,gscat;
662: DM dm,subdm;
664: if (nasm->fjtype == 2) X = nasm->xinit;
665: F = snes->vec_func;
666: if (snes->normschedule == SNES_NORM_NONE) SNESComputeFunction(snes,X,F);
667: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
668: SNESGetDM(snes,&dm);
669: if (nasm->eventrestrictinterp) PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);
670: if (nasm->fjtype != 1) {
671: for (i=0; i<nasm->n; i++) {
672: Xlloc = nasm->xl[i];
673: gscat = nasm->gscatter[i];
674: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
675: }
676: }
677: if (nasm->eventrestrictinterp) PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);
678: for (i=0; i<nasm->n; i++) {
679: Fl = nasm->subsnes[i]->vec_func;
680: Xl = nasm->x[i];
681: Xlloc = nasm->xl[i];
682: subsnes = nasm->subsnes[i];
683: oscat = nasm->oscatter[i];
684: gscat = nasm->gscatter[i];
685: if (nasm->fjtype != 1) VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
686: SNESGetDM(subsnes,&subdm);
687: DMSubDomainRestrict(dm,oscat,gscat,subdm);
688: if (nasm->fjtype != 1) {
689: DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
690: DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
691: }
692: if (subsnes->lagjacobian == -1) subsnes->lagjacobian = -2;
693: else if (subsnes->lagjacobian > 1) lag = subsnes->lagjacobian;
694: SNESComputeFunction(subsnes,Xl,Fl);
695: SNESComputeJacobian(subsnes,Xl,subsnes->jacobian,subsnes->jacobian_pre);
696: if (lag > 1) subsnes->lagjacobian = lag;
697: }
698: return 0;
699: }
701: static PetscErrorCode SNESSolve_NASM(SNES snes)
702: {
703: Vec F;
704: Vec X;
705: Vec B;
706: Vec Y;
707: PetscInt i;
708: PetscReal fnorm = 0.0;
709: SNESNormSchedule normschedule;
710: SNES_NASM *nasm = (SNES_NASM*)snes->data;
715: PetscCitationsRegister(SNESCitation,&SNEScite);
716: X = snes->vec_sol;
717: Y = snes->vec_sol_update;
718: F = snes->vec_func;
719: B = snes->vec_rhs;
721: PetscObjectSAWsTakeAccess((PetscObject)snes);
722: snes->iter = 0;
723: snes->norm = 0.;
724: PetscObjectSAWsGrantAccess((PetscObject)snes);
725: snes->reason = SNES_CONVERGED_ITERATING;
726: SNESGetNormSchedule(snes, &normschedule);
727: if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
728: /* compute the initial function and preconditioned update delX */
729: if (!snes->vec_func_init_set) {
730: SNESComputeFunction(snes,X,F);
731: } else snes->vec_func_init_set = PETSC_FALSE;
733: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
734: SNESCheckFunctionNorm(snes,fnorm);
735: PetscObjectSAWsTakeAccess((PetscObject)snes);
736: snes->iter = 0;
737: snes->norm = fnorm;
738: PetscObjectSAWsGrantAccess((PetscObject)snes);
739: SNESLogConvergenceHistory(snes,snes->norm,0);
740: SNESMonitor(snes,0,snes->norm);
742: /* test convergence */
743: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
744: if (snes->reason) return 0;
745: } else {
746: PetscObjectSAWsGrantAccess((PetscObject)snes);
747: SNESLogConvergenceHistory(snes,snes->norm,0);
748: SNESMonitor(snes,0,snes->norm);
749: }
751: /* Call general purpose update function */
752: if (snes->ops->update) {
753: (*snes->ops->update)(snes, snes->iter);
754: }
755: /* copy the initial solution over for later */
756: if (nasm->fjtype == 2) VecCopy(X,nasm->xinit);
758: for (i=0; i < snes->max_its; i++) {
759: SNESNASMSolveLocal_Private(snes,B,Y,X);
760: if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
761: SNESComputeFunction(snes,X,F);
762: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
763: SNESCheckFunctionNorm(snes,fnorm);
764: }
765: /* Monitor convergence */
766: PetscObjectSAWsTakeAccess((PetscObject)snes);
767: snes->iter = i+1;
768: snes->norm = fnorm;
769: PetscObjectSAWsGrantAccess((PetscObject)snes);
770: SNESLogConvergenceHistory(snes,snes->norm,0);
771: SNESMonitor(snes,snes->iter,snes->norm);
772: /* Test for convergence */
773: if (normschedule == SNES_NORM_ALWAYS) (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
774: if (snes->reason) break;
775: /* Call general purpose update function */
776: if (snes->ops->update) (*snes->ops->update)(snes, snes->iter);
777: }
778: if (nasm->finaljacobian) {
779: SNESNASMComputeFinalJacobian_Private(snes,X);
780: SNESCheckJacobianDomainerror(snes);
781: }
782: if (normschedule == SNES_NORM_ALWAYS) {
783: if (i == snes->max_its) {
784: PetscInfo(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
785: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
786: }
787: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */
788: return 0;
789: }
791: /*MC
792: SNESNASM - Nonlinear Additive Schwarz
794: Options Database:
795: + -snes_nasm_log - enable logging events for the communication and solve stages
796: . -snes_nasm_type <basic,restrict> - type of subdomain update used
797: . -snes_asm_damping <dmp> - the new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
798: . -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
799: . -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> - pick state the jacobian is calculated at
800: . -sub_snes_ - options prefix of the subdomain nonlinear solves
801: . -sub_ksp_ - options prefix of the subdomain Krylov solver
802: - -sub_pc_ - options prefix of the subdomain preconditioner
804: Level: advanced
806: 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
807: false and SNESView() and -snes_view do not display a KSP object. However, if the flag nasm->finaljacobian is set (for example, if
808: NASM is used as a nonlinear preconditioner for KSPASPIN) then SNESSetUpMatrices() is called to generate the Jacobian (needed by KSPASPIN)
809: 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
810: used by SNESASPIN they share the same Jacobian matrices because SNESSetUp() (called on the outer SNES KSPASPIN) causes the inner SNES
811: object (in this case SNESNASM) to inherit the outer Jacobian matrices.
813: References:
814: . * - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
815: SIAM Review, 57(4), 2015
817: .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types), SNESNASMSetType(), SNESNASMGetType(), SNESNASMSetSubdomains(), SNESNASMGetSubdomains(), SNESNASMGetSubdomainVecs(), SNESNASMSetComputeFinalJacobian(), SNESNASMSetDamping(), SNESNASMGetDamping()
818: M*/
820: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
821: {
822: SNES_NASM *nasm;
824: PetscNewLog(snes,&nasm);
825: snes->data = (void*)nasm;
827: nasm->n = PETSC_DECIDE;
828: nasm->subsnes = NULL;
829: nasm->x = NULL;
830: nasm->xl = NULL;
831: nasm->y = NULL;
832: nasm->b = NULL;
833: nasm->oscatter = NULL;
834: nasm->oscatter_copy = NULL;
835: nasm->iscatter = NULL;
836: nasm->gscatter = NULL;
837: nasm->damping = 1.;
839: nasm->type = PC_ASM_BASIC;
840: nasm->finaljacobian = PETSC_FALSE;
841: nasm->weight_set = PETSC_FALSE;
843: snes->ops->destroy = SNESDestroy_NASM;
844: snes->ops->setup = SNESSetUp_NASM;
845: snes->ops->setfromoptions = SNESSetFromOptions_NASM;
846: snes->ops->view = SNESView_NASM;
847: snes->ops->solve = SNESSolve_NASM;
848: snes->ops->reset = SNESReset_NASM;
850: snes->usesksp = PETSC_FALSE;
851: snes->usesnpc = PETSC_FALSE;
853: snes->alwayscomputesfinalresidual = PETSC_FALSE;
855: nasm->fjtype = 0;
856: nasm->xinit = NULL;
857: nasm->eventrestrictinterp = 0;
858: nasm->eventsubsolve = 0;
860: if (!snes->tolerancesset) {
861: snes->max_its = 10000;
862: snes->max_funcs = 10000;
863: }
865: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetType_C",SNESNASMSetType_NASM);
866: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetType_C",SNESNASMGetType_NASM);
867: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
868: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
869: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
870: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
871: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
872: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
873: return 0;
874: }
876: /*@
877: SNESNASMGetSNES - Gets a subsolver
879: Not collective
881: Input Parameters:
882: + snes - the SNES context
883: - i - the number of the subsnes to get
885: Output Parameters:
886: . subsnes - the subsolver context
888: Level: intermediate
890: .seealso: SNESNASM, SNESNASMGetNumber()
891: @*/
892: PetscErrorCode SNESNASMGetSNES(SNES snes,PetscInt i,SNES *subsnes)
893: {
894: SNES_NASM *nasm = (SNES_NASM*)snes->data;
897: *subsnes = nasm->subsnes[i];
898: return 0;
899: }
901: /*@
902: SNESNASMGetNumber - Gets number of subsolvers
904: Not collective
906: Input Parameters:
907: . snes - the SNES context
909: Output Parameters:
910: . n - the number of subsolvers
912: Level: intermediate
914: .seealso: SNESNASM, SNESNASMGetSNES()
915: @*/
916: PetscErrorCode SNESNASMGetNumber(SNES snes,PetscInt *n)
917: {
918: SNES_NASM *nasm = (SNES_NASM*)snes->data;
920: *n = nasm->n;
921: return 0;
922: }
924: /*@
925: SNESNASMSetWeight - Sets weight to use when adding overlapping updates
927: Collective
929: Input Parameters:
930: + snes - the SNES context
931: - weight - the weights to use (typically 1/N for each dof, where N is the number of patches it appears in)
933: Level: intermediate
935: .seealso: SNESNASM
936: @*/
937: PetscErrorCode SNESNASMSetWeight(SNES snes,Vec weight)
938: {
939: SNES_NASM *nasm = (SNES_NASM*)snes->data;
942: VecDestroy(&nasm->weight);
943: nasm->weight_set = PETSC_TRUE;
944: nasm->weight = weight;
945: PetscObjectReference((PetscObject)nasm->weight);
947: return 0;
948: }