Actual source code: nasm.c
petsc-3.4.5 2014-06-29
1: #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/
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: VecScatter *oscatter; /* scatter from global space to the subdomain global space */
12: VecScatter *iscatter; /* scatter from global space to the nonoverlapping subdomain space */
13: VecScatter *gscatter; /* scatter from global space to the subdomain local space */
14: PCASMType type; /* ASM type */
15: PetscBool usesdm; /* use the DM for setting up the subproblems */
16: PetscBool finaljacobian; /* compute the jacobian of the converged solution */
17: PetscReal damping; /* damping parameter for updates from the blocks */
18: PetscBool same_local_solves; /* flag to determine if the solvers have been individually modified */
20: /* logging events */
21: PetscLogEvent eventrestrictinterp;
22: PetscLogEvent eventsubsolve;
24: PetscInt fjtype; /* type of computed jacobian */
25: Vec xinit; /* initial solution in case the final jacobian type is computed as first */
26: } SNES_NASM;
28: const char *const SNESNASMTypes[] = {"NONE","RESTRICT","INTERPOLATE","BASIC","PCASMType","PC_ASM_",0};
29: const char *const SNESNASMFJTypes[] = {"FINALOUTER","FINALINNER","INITIAL"};
33: PetscErrorCode SNESReset_NASM(SNES snes)
34: {
35: SNES_NASM *nasm = (SNES_NASM*)snes->data;
37: PetscInt i;
40: for (i=0; i<nasm->n; i++) {
41: if (nasm->xl) { VecDestroy(&nasm->xl[i]); }
42: if (nasm->x) { VecDestroy(&nasm->x[i]); }
43: if (nasm->y) { VecDestroy(&nasm->y[i]); }
44: if (nasm->b) { VecDestroy(&nasm->b[i]); }
46: if (nasm->subsnes) { SNESDestroy(&nasm->subsnes[i]); }
47: if (nasm->oscatter) { VecScatterDestroy(&nasm->oscatter[i]); }
48: if (nasm->iscatter) { VecScatterDestroy(&nasm->iscatter[i]); }
49: if (nasm->gscatter) { VecScatterDestroy(&nasm->gscatter[i]); }
50: }
52: if (nasm->x) {PetscFree(nasm->x);}
53: if (nasm->xl) {PetscFree(nasm->xl);}
54: if (nasm->y) {PetscFree(nasm->y);}
55: if (nasm->b) {PetscFree(nasm->b);}
57: if (nasm->xinit) {VecDestroy(&nasm->xinit);}
59: if (nasm->subsnes) {PetscFree(nasm->subsnes);}
60: if (nasm->oscatter) {PetscFree(nasm->oscatter);}
61: if (nasm->iscatter) {PetscFree(nasm->iscatter);}
62: if (nasm->gscatter) {PetscFree(nasm->gscatter);}
64: nasm->eventrestrictinterp = 0;
65: nasm->eventsubsolve = 0;
66: return(0);
67: }
71: PetscErrorCode SNESDestroy_NASM(SNES snes)
72: {
76: SNESReset_NASM(snes);
77: PetscFree(snes->data);
78: return(0);
79: }
83: PetscErrorCode DMGlobalToLocalSubDomainDirichletHook_Private(DM dm,Vec g,InsertMode mode,Vec l,void *ctx)
84: {
86: Vec bcs = (Vec)ctx;
89: VecCopy(bcs,l);
90: return(0);
91: }
95: PetscErrorCode SNESSetUp_NASM(SNES snes)
96: {
97: SNES_NASM *nasm = (SNES_NASM*)snes->data;
99: DM dm,subdm;
100: DM *subdms;
101: PetscInt i;
102: const char *optionsprefix;
103: Vec F;
104: PetscMPIInt size;
105: KSP ksp;
106: PC pc;
109: if (!nasm->subsnes) {
110: SNESGetDM(snes,&dm);
111: if (dm) {
112: nasm->usesdm = PETSC_TRUE;
113: DMCreateDomainDecomposition(dm,&nasm->n,NULL,NULL,NULL,&subdms);
114: if (!subdms) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM has no default decomposition defined. Set subsolves manually with SNESNASMSetSubdomains().");
115: DMCreateDomainDecompositionScatters(dm,nasm->n,subdms,&nasm->iscatter,&nasm->oscatter,&nasm->gscatter);
117: SNESGetOptionsPrefix(snes, &optionsprefix);
118: PetscMalloc(nasm->n*sizeof(SNES),&nasm->subsnes);
119: for (i=0; i<nasm->n; i++) {
120: SNESCreate(PETSC_COMM_SELF,&nasm->subsnes[i]);
121: SNESAppendOptionsPrefix(nasm->subsnes[i],optionsprefix);
122: SNESAppendOptionsPrefix(nasm->subsnes[i],"sub_");
123: SNESSetDM(nasm->subsnes[i],subdms[i]);
124: MPI_Comm_size(PetscObjectComm((PetscObject)nasm->subsnes[i]),&size);
125: if (size == 1) {
126: SNESGetKSP(nasm->subsnes[i],&ksp);
127: KSPGetPC(ksp,&pc);
128: KSPSetType(ksp,KSPPREONLY);
129: PCSetType(pc,PCLU);
130: }
131: SNESSetFromOptions(nasm->subsnes[i]);
132: DMDestroy(&subdms[i]);
133: }
134: PetscFree(subdms);
135: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Cannot construct local problems automatically without a DM!");
136: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set subproblems manually if there is no DM!");
137: /* allocate the global vectors */
138: if (!nasm->x) {
139: PetscMalloc(nasm->n*sizeof(Vec),&nasm->x);
140: PetscMemzero(nasm->x,nasm->n*sizeof(Vec));
141: }
142: if (!nasm->xl) {
143: PetscMalloc(nasm->n*sizeof(Vec),&nasm->xl);
144: PetscMemzero(nasm->xl,nasm->n*sizeof(Vec));
145: }
146: if (!nasm->y) {
147: PetscMalloc(nasm->n*sizeof(Vec),&nasm->y);
148: PetscMemzero(nasm->y,nasm->n*sizeof(Vec));
149: }
150: if (!nasm->b) {
151: PetscMalloc(nasm->n*sizeof(Vec),&nasm->b);
152: PetscMemzero(nasm->b,nasm->n*sizeof(Vec));
153: }
155: for (i=0; i<nasm->n; i++) {
156: SNESGetFunction(nasm->subsnes[i],&F,NULL,NULL);
157: if (!nasm->x[i]) {VecDuplicate(F,&nasm->x[i]);}
158: if (!nasm->y[i]) {VecDuplicate(F,&nasm->y[i]);}
159: if (!nasm->b[i]) {VecDuplicate(F,&nasm->b[i]);}
160: if (!nasm->xl[i]) {
161: SNESGetDM(nasm->subsnes[i],&subdm);
162: DMCreateLocalVector(subdm,&nasm->xl[i]);
163: }
164: DMGlobalToLocalHookAdd(subdm,DMGlobalToLocalSubDomainDirichletHook_Private,NULL,nasm->xl[i]);
165: }
166: if (nasm->finaljacobian) {
167: SNESSetUpMatrices(snes);
168: if (nasm->fjtype == 2) {
169: VecDuplicate(snes->vec_sol,&nasm->xinit);
170: }
171: for (i=0; i<nasm->n;i++) {
172: SNESSetUpMatrices(nasm->subsnes[i]);
173: }
174: }
175: return(0);
176: }
180: PetscErrorCode SNESSetFromOptions_NASM(SNES snes)
181: {
182: PetscErrorCode ierr;
183: PCASMType asmtype;
184: PetscBool flg,monflg,subviewflg;
185: SNES_NASM *nasm = (SNES_NASM*)snes->data;
188: PetscOptionsHead("Nonlinear Additive Schwartz options");
189: PetscOptionsEnum("-snes_nasm_type","Type of restriction/extension","",SNESNASMTypes,(PetscEnum)nasm->type,(PetscEnum*)&asmtype,&flg);
190: if (flg) nasm->type = asmtype;
191: flg = PETSC_FALSE;
192: monflg = PETSC_TRUE;
193: PetscOptionsReal("-snes_nasm_damping","Log times for subSNES solves and restriction","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: }
216: PetscErrorCode SNESView_NASM(SNES snes, PetscViewer viewer)
217: {
218: SNES_NASM *nasm = (SNES_NASM*)snes->data;
220: PetscMPIInt rank,size;
221: PetscInt i,j,N,bsz;
222: PetscBool iascii,isstring;
223: PetscViewer sviewer;
224: MPI_Comm comm;
227: PetscObjectGetComm((PetscObject)snes,&comm);
228: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
229: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
230: MPI_Comm_rank(comm,&rank);
231: MPI_Comm_size(comm,&size);
232: MPI_Allreduce(&nasm->n,&N,1,MPIU_INT,MPI_SUM,comm);
233: if (iascii) {
234: PetscViewerASCIIPrintf(viewer, " Nonlinear Additive Schwarz: total subdomain blocks = %D\n",N);
235: if (nasm->same_local_solves) {
236: if (nasm->subsnes) {
237: PetscViewerASCIIPrintf(viewer," Local solve is the same for all blocks:\n");
238: PetscViewerASCIIPushTab(viewer);
239: PetscViewerGetSingleton(viewer,&sviewer);
240: if (!rank) {
241: PetscViewerASCIIPushTab(viewer);
242: SNESView(nasm->subsnes[0],sviewer);
243: PetscViewerASCIIPopTab(viewer);
244: }
245: PetscViewerRestoreSingleton(viewer,&sviewer);
246: PetscViewerASCIIPopTab(viewer);
247: }
248: } else {
249: /* print the solver on each block */
250: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
251: PetscViewerASCIISynchronizedPrintf(viewer," [%d] number of local blocks = %D\n",(int)rank,nasm->n);
252: PetscViewerFlush(viewer);
253: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
254: PetscViewerASCIIPrintf(viewer," Local solve info for each block is in the following SNES objects:\n");
255: PetscViewerASCIIPushTab(viewer);
256: PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
257: for (j=0; j<size; j++) {
258: PetscViewerGetSingleton(viewer,&sviewer);
259: if (rank == j) {
260: for (i=0; i<nasm->n; i++) {
261: VecGetLocalSize(nasm->x[i],&bsz);
262: PetscViewerASCIIPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);
263: SNESView(nasm->subsnes[i],sviewer);
264: PetscViewerASCIIPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
265: }
266: }
267: PetscViewerRestoreSingleton(viewer,&sviewer);
268: PetscViewerFlush(viewer);
269: }
270: PetscViewerASCIIPopTab(viewer);
271: }
272: } else if (isstring) {
273: PetscViewerStringSPrintf(viewer," blocks=%D,type=%s",N,SNESNASMTypes[nasm->type]);
274: PetscViewerGetSingleton(viewer,&sviewer);
275: if (nasm->subsnes && !rank) {SNESView(nasm->subsnes[0],sviewer);}
276: PetscViewerRestoreSingleton(viewer,&sviewer);
277: }
278: return(0);
279: }
283: /*@
284: SNESNASMSetSubdomains - Manually Set the context required to restrict and solve subdomain problems.
286: Not Collective
288: Input Parameters:
289: + SNES - the SNES context
290: . n - the number of local subdomains
291: . subsnes - solvers defined on the local subdomains
292: . iscatter - scatters into the nonoverlapping portions of the local subdomains
293: . oscatter - scatters into the overlapping portions of the local subdomains
294: - gscatter - scatters into the (ghosted) local vector of the local subdomain
296: Level: intermediate
298: .keywords: SNES, NASM
300: .seealso: SNESNASM, SNESNASMGetSubdomains()
301: @*/
302: PetscErrorCode SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
303: {
305: PetscErrorCode (*f)(SNES,PetscInt,SNES*,VecScatter*,VecScatter*,VecScatter*);
308: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",&f);
309: if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
310: return(0);
311: }
315: PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
316: {
317: PetscInt i;
319: SNES_NASM *nasm = (SNES_NASM*)snes->data;
322: if (snes->setupcalled) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"SNESNASMSetSubdomains() should be called before calling SNESSetUp().");
324: /* tear down the previously set things */
325: SNESReset(snes);
327: nasm->n = n;
328: if (oscatter) {
329: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)oscatter[i]);}
330: }
331: if (iscatter) {
332: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)iscatter[i]);}
333: }
334: if (gscatter) {
335: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)gscatter[i]);}
336: }
337: if (oscatter) {
338: PetscMalloc(n*sizeof(IS),&nasm->oscatter);
339: for (i=0; i<n; i++) {
340: nasm->oscatter[i] = oscatter[i];
341: }
342: }
343: if (iscatter) {
344: PetscMalloc(n*sizeof(IS),&nasm->iscatter);
345: for (i=0; i<n; i++) {
346: nasm->iscatter[i] = iscatter[i];
347: }
348: }
349: if (gscatter) {
350: PetscMalloc(n*sizeof(IS),&nasm->gscatter);
351: for (i=0; i<n; i++) {
352: nasm->gscatter[i] = gscatter[i];
353: }
354: }
356: if (subsnes) {
357: PetscMalloc(n*sizeof(SNES),&nasm->subsnes);
358: for (i=0; i<n; i++) {
359: nasm->subsnes[i] = subsnes[i];
360: }
361: nasm->same_local_solves = PETSC_FALSE;
362: }
363: return(0);
364: }
368: /*@
369: SNESNASMGetSubdomains - Get the local subdomain context.
371: Not Collective
373: Input Parameters:
374: . SNES - the SNES context
376: Output Parameters:
377: + n - the number of local subdomains
378: . subsnes - solvers defined on the local subdomains
379: . iscatter - scatters into the nonoverlapping portions of the local subdomains
380: . oscatter - scatters into the overlapping portions of the local subdomains
381: - gscatter - scatters into the (ghosted) local vector of the local subdomain
383: Level: intermediate
385: .keywords: SNES, NASM
387: .seealso: SNESNASM, SNESNASMSetSubdomains()
388: @*/
389: PetscErrorCode SNESNASMGetSubdomains(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
390: {
392: PetscErrorCode (*f)(SNES,PetscInt*,SNES**,VecScatter**,VecScatter**,VecScatter**);
395: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",&f);
396: if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
397: return(0);
398: }
402: PetscErrorCode SNESNASMGetSubdomains_NASM(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
403: {
404: SNES_NASM *nasm = (SNES_NASM*)snes->data;
407: if (n) *n = nasm->n;
408: if (oscatter) *oscatter = nasm->oscatter;
409: if (iscatter) *iscatter = nasm->iscatter;
410: if (gscatter) *gscatter = nasm->gscatter;
411: if (subsnes) {
412: *subsnes = nasm->subsnes;
413: nasm->same_local_solves = PETSC_FALSE;
414: }
415: return(0);
416: }
420: /*@
421: SNESNASMGetSubdomainVecs - Get the processor-local subdomain vectors
423: Not Collective
425: Input Parameters:
426: . SNES - the SNES context
428: Output Parameters:
429: + n - the number of local subdomains
430: . x - The subdomain solution vector
431: . y - The subdomain step vector
432: . b - The subdomain RHS vector
433: - xl - The subdomain local vectors (ghosted)
435: Level: developer
437: .keywords: SNES, NASM
439: .seealso: SNESNASM, SNESNASMGetSubdomains()
440: @*/
441: PetscErrorCode SNESNASMGetSubdomainVecs(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b, Vec **xl)
442: {
444: PetscErrorCode (*f)(SNES,PetscInt*,Vec**,Vec**,Vec**,Vec**);
447: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",&f);
448: if (f) {(f)(snes,n,x,y,b,xl);}
449: return(0);
450: }
454: PetscErrorCode SNESNASMGetSubdomainVecs_NASM(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b,Vec **xl)
455: {
456: SNES_NASM *nasm = (SNES_NASM*)snes->data;
459: if (n) *n = nasm->n;
460: if (x) *x = nasm->x;
461: if (y) *y = nasm->y;
462: if (b) *b = nasm->b;
463: if (xl) *xl = nasm->xl;
464: return(0);
465: }
469: /*@
470: SNESNASMSetComputeFinalJacobian - Schedules the computation of the global and subdomain jacobians upon convergence
472: Collective on SNES
474: Input Parameters:
475: + SNES - the SNES context
476: - flg - indication of whether to compute the jacobians or not
478: Level: developer
480: Notes: This is used almost exclusively in the implementation of ASPIN, where the converged subdomain and global jacobian
481: is needed at each linear iteration.
483: .keywords: SNES, NASM, ASPIN
485: .seealso: SNESNASM, SNESNASMGetSubdomains()
486: @*/
487: PetscErrorCode SNESNASMSetComputeFinalJacobian(SNES snes,PetscBool flg)
488: {
489: PetscErrorCode (*f)(SNES,PetscBool);
493: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",&f);
494: if (f) {(f)(snes,flg);}
495: return(0);
496: }
500: PetscErrorCode SNESNASMSetComputeFinalJacobian_NASM(SNES snes,PetscBool flg)
501: {
502: SNES_NASM *nasm = (SNES_NASM*)snes->data;
505: nasm->finaljacobian = flg;
506: if (flg) snes->usesksp = PETSC_TRUE;
507: return(0);
508: }
512: /*@
513: SNESNASMSetDamping - Sets the update damping for NASM
515: Logically collective on SNES
517: Input Parameters:
518: + SNES - the SNES context
519: - dmp - damping
521: Level: intermediate
523: .keywords: SNES, NASM, damping
525: .seealso: SNESNASM, SNESNASMGetDamping()
526: @*/
527: PetscErrorCode SNESNASMSetDamping(SNES snes,PetscReal dmp)
528: {
529: PetscErrorCode (*f)(SNES,PetscReal);
533: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetDamping_C",(void (**)(void))&f);
534: if (f) {(f)(snes,dmp);}
535: return(0);
536: }
540: PetscErrorCode SNESNASMSetDamping_NASM(SNES snes,PetscReal dmp)
541: {
542: SNES_NASM *nasm = (SNES_NASM*)snes->data;
545: nasm->damping = dmp;
546: return(0);
547: }
551: /*@
552: SNESNASMGetDamping - Gets the update damping for NASM
554: Not Collective
556: Input Parameters:
557: + SNES - the SNES context
558: - dmp - damping
560: Level: intermediate
562: .keywords: SNES, NASM, damping
564: .seealso: SNESNASM, SNESNASMSetDamping()
565: @*/
566: PetscErrorCode SNESNASMGetDamping(SNES snes,PetscReal *dmp)
567: {
568: PetscErrorCode (*f)(SNES,PetscReal*);
572: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetDamping_C",(void (**)(void))&f);
573: if (f) {(f)(snes,dmp);}
574: return(0);
575: }
579: PetscErrorCode SNESNASMGetDamping_NASM(SNES snes,PetscReal *dmp)
580: {
581: SNES_NASM *nasm = (SNES_NASM*)snes->data;
584: *dmp = nasm->damping;
585: return(0);
586: }
591: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X)
592: {
593: SNES_NASM *nasm = (SNES_NASM*)snes->data;
594: SNES subsnes;
595: PetscInt i;
596: PetscReal dmp;
598: Vec Xlloc,Xl,Bl,Yl;
599: VecScatter iscat,oscat,gscat;
600: DM dm,subdm;
603: SNESGetDM(snes,&dm);
604: SNESNASMGetDamping(snes,&dmp);
605: VecSet(Y,0);
606: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
607: for (i=0; i<nasm->n; i++) {
608: /* scatter the solution to the local solution */
609: Xlloc = nasm->xl[i];
610: gscat = nasm->gscatter[i];
611: oscat = nasm->oscatter[i];
612: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
613: if (B) {
614: /* scatter the RHS to the local RHS */
615: Bl = nasm->b[i];
616: VecScatterBegin(oscat,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
617: }
618: }
619: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
622: if (nasm->eventsubsolve) {PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);}
623: for (i=0; i<nasm->n; i++) {
624: Xl = nasm->x[i];
625: Xlloc = nasm->xl[i];
626: Yl = nasm->y[i];
627: subsnes = nasm->subsnes[i];
628: SNESGetDM(subsnes,&subdm);
629: iscat = nasm->iscatter[i];
630: oscat = nasm->oscatter[i];
631: gscat = nasm->gscatter[i];
632: VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
633: if (B) {
634: Bl = nasm->b[i];
635: VecScatterEnd(oscat,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
636: } else Bl = NULL;
637: DMSubDomainRestrict(dm,oscat,gscat,subdm);
638: DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
639: DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
640: VecCopy(Xl,Yl);
641: SNESSolve(subsnes,Bl,Xl);
642: VecAYPX(Yl,-1.0,Xl);
643: if (nasm->type == PC_ASM_BASIC) {
644: VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
645: } else if (nasm->type == PC_ASM_RESTRICT) {
646: VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
647: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
648: }
649: if (nasm->eventsubsolve) {PetscLogEventEnd(nasm->eventsubsolve,snes,0,0,0);}
650: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
651: for (i=0; i<nasm->n; i++) {
652: Yl = nasm->y[i];
653: iscat = nasm->iscatter[i];
654: oscat = nasm->oscatter[i];
655: if (nasm->type == PC_ASM_BASIC) {
656: VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
657: } else if (nasm->type == PC_ASM_RESTRICT) {
658: VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
659: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
660: }
661: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
662: VecAXPY(X,dmp,Y);
663: return(0);
664: }
668: PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
669: {
670: Vec X = Xfinal;
671: SNES_NASM *nasm = (SNES_NASM*)snes->data;
672: SNES subsnes;
673: PetscInt i,lag = 1;
675: Vec Xlloc,Xl,Fl,F;
676: VecScatter oscat,gscat;
677: DM dm,subdm;
678: MatStructure flg = DIFFERENT_NONZERO_PATTERN;
680: if (nasm->fjtype == 2) X = nasm->xinit;
681: F = snes->vec_func;
682: if (snes->normtype == SNES_NORM_NONE) {SNESComputeFunction(snes,X,F);}
683: SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);
684: SNESGetDM(snes,&dm);
685: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
686: if (nasm->fjtype != 1) {
687: for (i=0; i<nasm->n; i++) {
688: Xlloc = nasm->xl[i];
689: gscat = nasm->gscatter[i];
690: oscat = nasm->oscatter[i];
691: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
692: }
693: }
694: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
695: for (i=0; i<nasm->n; i++) {
696: Fl = nasm->subsnes[i]->vec_func;
697: Xl = nasm->x[i];
698: Xlloc = nasm->xl[i];
699: subsnes = nasm->subsnes[i];
700: oscat = nasm->oscatter[i];
701: gscat = nasm->gscatter[i];
702: if (nasm->fjtype != 1) {VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);}
703: SNESGetDM(subsnes,&subdm);
704: DMSubDomainRestrict(dm,oscat,gscat,subdm);
705: if (nasm->fjtype != 1) {
706: DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
707: DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
708: }
709: if (subsnes->lagjacobian == -1) subsnes->lagjacobian = -2;
710: else if (subsnes->lagjacobian > 1) lag = subsnes->lagjacobian;
711: SNESComputeFunction(subsnes,Xl,Fl);
712: SNESComputeJacobian(subsnes,Xl,&subsnes->jacobian,&subsnes->jacobian_pre,&flg);
713: if (lag > 1) subsnes->lagjacobian = lag;
714: KSPSetOperators(subsnes->ksp,subsnes->jacobian,subsnes->jacobian_pre,flg);
715: }
716: return(0);
717: }
721: PetscErrorCode SNESSolve_NASM(SNES snes)
722: {
723: Vec F;
724: Vec X;
725: Vec B;
726: Vec Y;
727: PetscInt i;
728: PetscReal fnorm = 0.0;
730: SNESNormType normtype;
731: SNES_NASM *nasm = (SNES_NASM*)snes->data;
734: X = snes->vec_sol;
735: Y = snes->vec_sol_update;
736: F = snes->vec_func;
737: B = snes->vec_rhs;
739: PetscObjectAMSTakeAccess((PetscObject)snes);
740: snes->iter = 0;
741: snes->norm = 0.;
742: PetscObjectAMSGrantAccess((PetscObject)snes);
743: snes->reason = SNES_CONVERGED_ITERATING;
744: SNESGetNormType(snes, &normtype);
745: if (normtype == SNES_NORM_FUNCTION || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
746: /* compute the initial function and preconditioned update delX */
747: if (!snes->vec_func_init_set) {
748: SNESComputeFunction(snes,X,F);
749: if (snes->domainerror) {
750: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
751: return(0);
752: }
753: } else snes->vec_func_init_set = PETSC_FALSE;
755: /* convergence test */
756: if (!snes->norm_init_set) {
757: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
758: if (PetscIsInfOrNanReal(fnorm)) {
759: snes->reason = SNES_DIVERGED_FNORM_NAN;
760: return(0);
761: }
762: } else {
763: fnorm = snes->norm_init;
764: snes->norm_init_set = PETSC_FALSE;
765: }
766: PetscObjectAMSTakeAccess((PetscObject)snes);
767: snes->iter = 0;
768: snes->norm = fnorm;
769: PetscObjectAMSGrantAccess((PetscObject)snes);
770: SNESLogConvergenceHistory(snes,snes->norm,0);
771: SNESMonitor(snes,0,snes->norm);
773: /* set parameter for default relative tolerance convergence test */
774: snes->ttol = fnorm*snes->rtol;
776: /* test convergence */
777: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
778: if (snes->reason) return(0);
779: } else {
780: PetscObjectAMSGrantAccess((PetscObject)snes);
781: SNESLogConvergenceHistory(snes,snes->norm,0);
782: SNESMonitor(snes,0,snes->norm);
783: }
785: /* Call general purpose update function */
786: if (snes->ops->update) {
787: (*snes->ops->update)(snes, snes->iter);
788: }
789: /* copy the initial solution over for later */
790: if (nasm->fjtype == 2) {VecCopy(X,nasm->xinit);}
792: for (i = 0; i < snes->max_its; i++) {
793: SNESNASMSolveLocal_Private(snes,B,Y,X);
794: if (normtype == SNES_NORM_FUNCTION || ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY))) {
795: SNESComputeFunction(snes,X,F);
796: if (snes->domainerror) {
797: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
798: break;
799: }
800: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
801: if (PetscIsInfOrNanReal(fnorm)) {
802: snes->reason = SNES_DIVERGED_FNORM_NAN;
803: break;
804: }
805: }
806: /* Monitor convergence */
807: PetscObjectAMSTakeAccess((PetscObject)snes);
808: snes->iter = i+1;
809: snes->norm = fnorm;
810: PetscObjectAMSGrantAccess((PetscObject)snes);
811: SNESLogConvergenceHistory(snes,snes->norm,0);
812: SNESMonitor(snes,snes->iter,snes->norm);
813: /* Test for convergence */
814: if (normtype == SNES_NORM_FUNCTION) {(*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);}
815: if (snes->reason) break;
816: /* Call general purpose update function */
817: if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
818: }
819: if (nasm->finaljacobian) {SNESNASMComputeFinalJacobian_Private(snes,X);}
820: if (normtype == SNES_NORM_FUNCTION) {
821: if (i == snes->max_its) {
822: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
823: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
824: }
825: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */
826: return(0);
827: }
829: /*MC
830: SNESNASM - Nonlinear Additive Schwartz
832: Options Database:
833: + -snes_nasm_log - enable logging events for the communication and solve stages
834: . -snes_nasm_type <basic,restrict> - type of subdomain update used
835: . -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
836: . -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> pick state the jacobian is calculated at
837: . -sub_snes_ - options prefix of the subdomain nonlinear solves
838: . -sub_ksp_ - options prefix of the subdomain Krylov solver
839: - -sub_pc_ - options prefix of the subdomain preconditioner
841: Level: advanced
843: .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types)
844: M*/
848: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
849: {
850: SNES_NASM *nasm;
854: PetscNewLog(snes, SNES_NASM, &nasm);
855: snes->data = (void*)nasm;
857: nasm->n = PETSC_DECIDE;
858: nasm->subsnes = 0;
859: nasm->x = 0;
860: nasm->xl = 0;
861: nasm->y = 0;
862: nasm->b = 0;
863: nasm->oscatter = 0;
864: nasm->iscatter = 0;
865: nasm->gscatter = 0;
866: nasm->damping = 1.;
868: nasm->type = PC_ASM_BASIC;
869: nasm->finaljacobian = PETSC_FALSE;
870: nasm->same_local_solves = PETSC_TRUE;
872: snes->ops->destroy = SNESDestroy_NASM;
873: snes->ops->setup = SNESSetUp_NASM;
874: snes->ops->setfromoptions = SNESSetFromOptions_NASM;
875: snes->ops->view = SNESView_NASM;
876: snes->ops->solve = SNESSolve_NASM;
877: snes->ops->reset = SNESReset_NASM;
879: snes->usesksp = PETSC_FALSE;
880: snes->usespc = PETSC_FALSE;
882: nasm->fjtype = 0;
883: nasm->xinit = NULL;
884: nasm->eventrestrictinterp = 0;
885: nasm->eventsubsolve = 0;
887: if (!snes->tolerancesset) {
888: snes->max_its = 10000;
889: snes->max_funcs = 10000;
890: }
892: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
893: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
894: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
895: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
896: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
897: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
898: return(0);
899: }