Actual source code: nasm.c
petsc-3.7.5 2017-01-01
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: PetscMalloc1(nasm->n,&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: PetscCalloc1(nasm->n,&nasm->x);
140: }
141: if (!nasm->xl) {
142: PetscCalloc1(nasm->n,&nasm->xl);
143: }
144: if (!nasm->y) {
145: PetscCalloc1(nasm->n,&nasm->y);
146: }
147: if (!nasm->b) {
148: PetscCalloc1(nasm->n,&nasm->b);
149: }
151: for (i=0; i<nasm->n; i++) {
152: SNESGetFunction(nasm->subsnes[i],&F,NULL,NULL);
153: if (!nasm->x[i]) {VecDuplicate(F,&nasm->x[i]);}
154: if (!nasm->y[i]) {VecDuplicate(F,&nasm->y[i]);}
155: if (!nasm->b[i]) {VecDuplicate(F,&nasm->b[i]);}
156: if (!nasm->xl[i]) {
157: SNESGetDM(nasm->subsnes[i],&subdm);
158: DMCreateLocalVector(subdm,&nasm->xl[i]);
159: DMGlobalToLocalHookAdd(subdm,DMGlobalToLocalSubDomainDirichletHook_Private,NULL,nasm->xl[i]);
160: }
161: }
162: if (nasm->finaljacobian) {
163: SNESSetUpMatrices(snes);
164: if (nasm->fjtype == 2) {
165: VecDuplicate(snes->vec_sol,&nasm->xinit);
166: }
167: for (i=0; i<nasm->n;i++) {
168: SNESSetUpMatrices(nasm->subsnes[i]);
169: }
170: }
171: return(0);
172: }
176: PetscErrorCode SNESSetFromOptions_NASM(PetscOptionItems *PetscOptionsObject,SNES snes)
177: {
178: PetscErrorCode ierr;
179: PCASMType asmtype;
180: PetscBool flg,monflg,subviewflg;
181: SNES_NASM *nasm = (SNES_NASM*)snes->data;
184: PetscOptionsHead(PetscOptionsObject,"Nonlinear Additive Schwartz options");
185: PetscOptionsEnum("-snes_nasm_type","Type of restriction/extension","",SNESNASMTypes,(PetscEnum)nasm->type,(PetscEnum*)&asmtype,&flg);
186: if (flg) {SNESNASMSetType(snes,asmtype);}
187: flg = PETSC_FALSE;
188: monflg = PETSC_TRUE;
189: 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);
190: if (flg) {SNESNASMSetDamping(snes,nasm->damping);}
191: subviewflg = PETSC_FALSE;
192: PetscOptionsBool("-snes_nasm_sub_view","Print detailed information for every processor when using -snes_view","",subviewflg,&subviewflg,&flg);
193: if (flg) {
194: nasm->same_local_solves = PETSC_FALSE;
195: if (!subviewflg) {
196: nasm->same_local_solves = PETSC_TRUE;
197: }
198: }
199: PetscOptionsBool("-snes_nasm_finaljacobian","Compute the global jacobian of the final iterate (for ASPIN)","",nasm->finaljacobian,&nasm->finaljacobian,NULL);
200: PetscOptionsEList("-snes_nasm_finaljacobian_type","The type of the final jacobian computed.","",SNESNASMFJTypes,3,SNESNASMFJTypes[0],&nasm->fjtype,NULL);
201: PetscOptionsBool("-snes_nasm_log","Log times for subSNES solves and restriction","",monflg,&monflg,&flg);
202: if (flg) {
203: PetscLogEventRegister("SNESNASMSubSolve",((PetscObject)snes)->classid,&nasm->eventsubsolve);
204: PetscLogEventRegister("SNESNASMRestrict",((PetscObject)snes)->classid,&nasm->eventrestrictinterp);
205: }
206: PetscOptionsTail();
207: return(0);
208: }
212: PetscErrorCode SNESView_NASM(SNES snes, PetscViewer viewer)
213: {
214: SNES_NASM *nasm = (SNES_NASM*)snes->data;
216: PetscMPIInt rank,size;
217: PetscInt i,N,bsz;
218: PetscBool iascii,isstring;
219: PetscViewer sviewer;
220: MPI_Comm comm;
223: PetscObjectGetComm((PetscObject)snes,&comm);
224: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
225: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
226: MPI_Comm_rank(comm,&rank);
227: MPI_Comm_size(comm,&size);
228: MPIU_Allreduce(&nasm->n,&N,1,MPIU_INT,MPI_SUM,comm);
229: if (iascii) {
230: PetscViewerASCIIPrintf(viewer, " Nonlinear Additive Schwarz: total subdomain blocks = %D\n",N);
231: if (nasm->same_local_solves) {
232: if (nasm->subsnes) {
233: PetscViewerASCIIPrintf(viewer," Local solve is the same for all blocks:\n");
234: PetscViewerASCIIPushTab(viewer);
235: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
236: if (!rank) {
237: PetscViewerASCIIPushTab(viewer);
238: SNESView(nasm->subsnes[0],sviewer);
239: PetscViewerASCIIPopTab(viewer);
240: }
241: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
242: PetscViewerASCIIPopTab(viewer);
243: }
244: } else {
245: /* print the solver on each block */
246: PetscViewerASCIIPushSynchronized(viewer);
247: PetscViewerASCIISynchronizedPrintf(viewer," [%d] number of local blocks = %D\n",(int)rank,nasm->n);
248: PetscViewerFlush(viewer);
249: PetscViewerASCIIPopSynchronized(viewer);
250: PetscViewerASCIIPrintf(viewer," Local solve info for each block is in the following SNES objects:\n");
251: PetscViewerASCIIPushTab(viewer);
252: PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
253: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
254: for (i=0; i<nasm->n; i++) {
255: VecGetLocalSize(nasm->x[i],&bsz);
256: PetscViewerASCIIPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);
257: SNESView(nasm->subsnes[i],sviewer);
258: PetscViewerASCIIPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
259: }
260: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
261: PetscViewerFlush(viewer);
262: PetscViewerASCIIPopTab(viewer);
263: }
264: } else if (isstring) {
265: PetscViewerStringSPrintf(viewer," blocks=%D,type=%s",N,SNESNASMTypes[nasm->type]);
266: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
267: if (nasm->subsnes && !rank) {SNESView(nasm->subsnes[0],sviewer);}
268: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
269: }
270: return(0);
271: }
275: /*@
276: SNESNASMSetType - Set the type of subdomain update used
278: Logically Collective on SNES
280: Input Parameters:
281: + SNES - the SNES context
282: - type - the type of update, PC_ASM_BASIC or PC_ASM_RESTRICT
284: Level: intermediate
286: .keywords: SNES, NASM
288: .seealso: SNESNASM, SNESNASMGetType(), PCASMSetType()
289: @*/
290: PetscErrorCode SNESNASMSetType(SNES snes,PCASMType type)
291: {
293: PetscErrorCode (*f)(SNES,PCASMType);
296: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetType_C",&f);
297: if (f) {(f)(snes,type);}
298: return(0);
299: }
303: PetscErrorCode SNESNASMSetType_NASM(SNES snes,PCASMType type)
304: {
305: SNES_NASM *nasm = (SNES_NASM*)snes->data;
308: if (type != PC_ASM_BASIC && type != PC_ASM_RESTRICT) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"SNESNASM only supports basic and restrict types");
309: nasm->type = type;
310: return(0);
311: }
315: /*@
316: SNESNASMGetType - Get the type of subdomain update used
318: Logically Collective on SNES
320: Input Parameters:
321: . SNES - the SNES context
323: Output Parameters:
324: . type - the type of update
326: Level: intermediate
328: .keywords: SNES, NASM
330: .seealso: SNESNASM, SNESNASMSetType(), PCASMGetType()
331: @*/
332: PetscErrorCode SNESNASMGetType(SNES snes,PCASMType *type)
333: {
335: PetscErrorCode (*f)(SNES,PCASMType*);
338: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetType_C",&f);
339: if (f) {(f)(snes,type);}
340: return(0);
341: }
345: PetscErrorCode SNESNASMGetType_NASM(SNES snes,PCASMType *type)
346: {
347: SNES_NASM *nasm = (SNES_NASM*)snes->data;
350: *type = nasm->type;
351: return(0);
352: }
356: /*@
357: SNESNASMSetSubdomains - Manually Set the context required to restrict and solve subdomain problems.
359: Not Collective
361: Input Parameters:
362: + SNES - the SNES context
363: . n - the number of local subdomains
364: . subsnes - solvers defined on the local subdomains
365: . iscatter - scatters into the nonoverlapping portions of the local subdomains
366: . oscatter - scatters into the overlapping portions of the local subdomains
367: - gscatter - scatters into the (ghosted) local vector of the local subdomain
369: Level: intermediate
371: .keywords: SNES, NASM
373: .seealso: SNESNASM, SNESNASMGetSubdomains()
374: @*/
375: PetscErrorCode SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
376: {
378: PetscErrorCode (*f)(SNES,PetscInt,SNES*,VecScatter*,VecScatter*,VecScatter*);
381: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",&f);
382: if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
383: return(0);
384: }
388: PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
389: {
390: PetscInt i;
392: SNES_NASM *nasm = (SNES_NASM*)snes->data;
395: if (snes->setupcalled) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"SNESNASMSetSubdomains() should be called before calling SNESSetUp().");
397: /* tear down the previously set things */
398: SNESReset(snes);
400: nasm->n = n;
401: if (oscatter) {
402: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)oscatter[i]);}
403: }
404: if (iscatter) {
405: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)iscatter[i]);}
406: }
407: if (gscatter) {
408: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)gscatter[i]);}
409: }
410: if (oscatter) {
411: PetscMalloc1(n,&nasm->oscatter);
412: for (i=0; i<n; i++) {
413: nasm->oscatter[i] = oscatter[i];
414: }
415: }
416: if (iscatter) {
417: PetscMalloc1(n,&nasm->iscatter);
418: for (i=0; i<n; i++) {
419: nasm->iscatter[i] = iscatter[i];
420: }
421: }
422: if (gscatter) {
423: PetscMalloc1(n,&nasm->gscatter);
424: for (i=0; i<n; i++) {
425: nasm->gscatter[i] = gscatter[i];
426: }
427: }
429: if (subsnes) {
430: PetscMalloc1(n,&nasm->subsnes);
431: for (i=0; i<n; i++) {
432: nasm->subsnes[i] = subsnes[i];
433: }
434: nasm->same_local_solves = PETSC_FALSE;
435: }
436: return(0);
437: }
441: /*@
442: SNESNASMGetSubdomains - Get the local subdomain context.
444: Not Collective
446: Input Parameters:
447: . SNES - the SNES context
449: Output Parameters:
450: + n - the number of local subdomains
451: . subsnes - solvers defined on the local subdomains
452: . iscatter - scatters into the nonoverlapping portions of the local subdomains
453: . oscatter - scatters into the overlapping portions of the local subdomains
454: - gscatter - scatters into the (ghosted) local vector of the local subdomain
456: Level: intermediate
458: .keywords: SNES, NASM
460: .seealso: SNESNASM, SNESNASMSetSubdomains()
461: @*/
462: PetscErrorCode SNESNASMGetSubdomains(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
463: {
465: PetscErrorCode (*f)(SNES,PetscInt*,SNES**,VecScatter**,VecScatter**,VecScatter**);
468: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",&f);
469: if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
470: return(0);
471: }
475: PetscErrorCode SNESNASMGetSubdomains_NASM(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
476: {
477: SNES_NASM *nasm = (SNES_NASM*)snes->data;
480: if (n) *n = nasm->n;
481: if (oscatter) *oscatter = nasm->oscatter;
482: if (iscatter) *iscatter = nasm->iscatter;
483: if (gscatter) *gscatter = nasm->gscatter;
484: if (subsnes) {
485: *subsnes = nasm->subsnes;
486: nasm->same_local_solves = PETSC_FALSE;
487: }
488: return(0);
489: }
493: /*@
494: SNESNASMGetSubdomainVecs - Get the processor-local subdomain vectors
496: Not Collective
498: Input Parameters:
499: . SNES - the SNES context
501: Output Parameters:
502: + n - the number of local subdomains
503: . x - The subdomain solution vector
504: . y - The subdomain step vector
505: . b - The subdomain RHS vector
506: - xl - The subdomain local vectors (ghosted)
508: Level: developer
510: .keywords: SNES, NASM
512: .seealso: SNESNASM, SNESNASMGetSubdomains()
513: @*/
514: PetscErrorCode SNESNASMGetSubdomainVecs(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b, Vec **xl)
515: {
517: PetscErrorCode (*f)(SNES,PetscInt*,Vec**,Vec**,Vec**,Vec**);
520: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",&f);
521: if (f) {(f)(snes,n,x,y,b,xl);}
522: return(0);
523: }
527: PetscErrorCode SNESNASMGetSubdomainVecs_NASM(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b,Vec **xl)
528: {
529: SNES_NASM *nasm = (SNES_NASM*)snes->data;
532: if (n) *n = nasm->n;
533: if (x) *x = nasm->x;
534: if (y) *y = nasm->y;
535: if (b) *b = nasm->b;
536: if (xl) *xl = nasm->xl;
537: return(0);
538: }
542: /*@
543: SNESNASMSetComputeFinalJacobian - Schedules the computation of the global and subdomain jacobians upon convergence
545: Collective on SNES
547: Input Parameters:
548: + SNES - the SNES context
549: - flg - indication of whether to compute the jacobians or not
551: Level: developer
553: Notes: This is used almost exclusively in the implementation of ASPIN, where the converged subdomain and global jacobian
554: is needed at each linear iteration.
556: .keywords: SNES, NASM, ASPIN
558: .seealso: SNESNASM, SNESNASMGetSubdomains()
559: @*/
560: PetscErrorCode SNESNASMSetComputeFinalJacobian(SNES snes,PetscBool flg)
561: {
562: PetscErrorCode (*f)(SNES,PetscBool);
566: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",&f);
567: if (f) {(f)(snes,flg);}
568: return(0);
569: }
573: PetscErrorCode SNESNASMSetComputeFinalJacobian_NASM(SNES snes,PetscBool flg)
574: {
575: SNES_NASM *nasm = (SNES_NASM*)snes->data;
578: nasm->finaljacobian = flg;
579: if (flg) snes->usesksp = PETSC_TRUE;
580: return(0);
581: }
585: /*@
586: SNESNASMSetDamping - Sets the update damping for NASM
588: Logically collective on SNES
590: Input Parameters:
591: + SNES - the SNES context
592: - dmp - damping
594: Level: intermediate
596: Notes: The new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
598: .keywords: SNES, NASM, damping
600: .seealso: SNESNASM, SNESNASMGetDamping()
601: @*/
602: PetscErrorCode SNESNASMSetDamping(SNES snes,PetscReal dmp)
603: {
604: PetscErrorCode (*f)(SNES,PetscReal);
608: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetDamping_C",(void (**)(void))&f);
609: if (f) {(f)(snes,dmp);}
610: return(0);
611: }
615: PetscErrorCode SNESNASMSetDamping_NASM(SNES snes,PetscReal dmp)
616: {
617: SNES_NASM *nasm = (SNES_NASM*)snes->data;
620: nasm->damping = dmp;
621: return(0);
622: }
626: /*@
627: SNESNASMGetDamping - Gets the update damping for NASM
629: Not Collective
631: Input Parameters:
632: + SNES - the SNES context
633: - dmp - damping
635: Level: intermediate
637: .keywords: SNES, NASM, damping
639: .seealso: SNESNASM, SNESNASMSetDamping()
640: @*/
641: PetscErrorCode SNESNASMGetDamping(SNES snes,PetscReal *dmp)
642: {
643: PetscErrorCode (*f)(SNES,PetscReal*);
647: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetDamping_C",(void (**)(void))&f);
648: if (f) {(f)(snes,dmp);}
649: return(0);
650: }
654: PetscErrorCode SNESNASMGetDamping_NASM(SNES snes,PetscReal *dmp)
655: {
656: SNES_NASM *nasm = (SNES_NASM*)snes->data;
659: *dmp = nasm->damping;
660: return(0);
661: }
666: /*
667: Input Parameters:
668: + snes - The solver
669: . B - The RHS vector
670: - X - The initial guess
672: Output Parameters:
673: . Y - The solution update
675: TODO: All scatters should be packed into one
676: */
677: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X)
678: {
679: SNES_NASM *nasm = (SNES_NASM*)snes->data;
680: SNES subsnes;
681: PetscInt i;
682: PetscReal dmp;
684: Vec Xlloc,Xl,Bl,Yl;
685: VecScatter iscat,oscat,gscat;
686: DM dm,subdm;
687: PCASMType type;
690: SNESNASMGetType(snes,&type);
691: SNESGetDM(snes,&dm);
692: VecSet(Y,0);
693: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
694: for (i=0; i<nasm->n; i++) {
695: /* scatter the solution to the local solution */
696: Xlloc = nasm->xl[i];
697: gscat = nasm->gscatter[i];
698: oscat = nasm->oscatter[i];
699: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
700: if (B) {
701: /* scatter the RHS to the local RHS */
702: Bl = nasm->b[i];
703: VecScatterBegin(oscat,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
704: }
705: }
706: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
709: if (nasm->eventsubsolve) {PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);}
710: for (i=0; i<nasm->n; i++) {
711: Xl = nasm->x[i];
712: Xlloc = nasm->xl[i];
713: Yl = nasm->y[i];
714: subsnes = nasm->subsnes[i];
715: SNESGetDM(subsnes,&subdm);
716: iscat = nasm->iscatter[i];
717: oscat = nasm->oscatter[i];
718: gscat = nasm->gscatter[i];
719: VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
720: if (B) {
721: Bl = nasm->b[i];
722: VecScatterEnd(oscat,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
723: } else Bl = NULL;
724: DMSubDomainRestrict(dm,oscat,gscat,subdm);
725: /* Could scatter directly from X */
726: DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
727: DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
728: VecCopy(Xl,Yl);
729: SNESSolve(subsnes,Bl,Xl);
730: VecAYPX(Yl,-1.0,Xl);
731: if (type == PC_ASM_BASIC) {
732: VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
733: } else if (type == PC_ASM_RESTRICT) {
734: VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
735: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
736: }
737: if (nasm->eventsubsolve) {PetscLogEventEnd(nasm->eventsubsolve,snes,0,0,0);}
738: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
739: for (i=0; i<nasm->n; i++) {
740: Yl = nasm->y[i];
741: iscat = nasm->iscatter[i];
742: oscat = nasm->oscatter[i];
743: if (type == PC_ASM_BASIC) {
744: VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
745: } else if (type == PC_ASM_RESTRICT) {
746: VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
747: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
748: }
749: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
750: SNESNASMGetDamping(snes,&dmp);
751: VecAXPY(X,dmp,Y);
752: return(0);
753: }
757: PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
758: {
759: Vec X = Xfinal;
760: SNES_NASM *nasm = (SNES_NASM*)snes->data;
761: SNES subsnes;
762: PetscInt i,lag = 1;
764: Vec Xlloc,Xl,Fl,F;
765: VecScatter oscat,gscat;
766: DM dm,subdm;
769: if (nasm->fjtype == 2) X = nasm->xinit;
770: F = snes->vec_func;
771: if (snes->normschedule == SNES_NORM_NONE) {SNESComputeFunction(snes,X,F);}
772: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
773: SNESGetDM(snes,&dm);
774: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
775: if (nasm->fjtype != 1) {
776: for (i=0; i<nasm->n; i++) {
777: Xlloc = nasm->xl[i];
778: gscat = nasm->gscatter[i];
779: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
780: }
781: }
782: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
783: for (i=0; i<nasm->n; i++) {
784: Fl = nasm->subsnes[i]->vec_func;
785: Xl = nasm->x[i];
786: Xlloc = nasm->xl[i];
787: subsnes = nasm->subsnes[i];
788: oscat = nasm->oscatter[i];
789: gscat = nasm->gscatter[i];
790: if (nasm->fjtype != 1) {VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);}
791: SNESGetDM(subsnes,&subdm);
792: DMSubDomainRestrict(dm,oscat,gscat,subdm);
793: if (nasm->fjtype != 1) {
794: DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
795: DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
796: }
797: if (subsnes->lagjacobian == -1) subsnes->lagjacobian = -2;
798: else if (subsnes->lagjacobian > 1) lag = subsnes->lagjacobian;
799: SNESComputeFunction(subsnes,Xl,Fl);
800: SNESComputeJacobian(subsnes,Xl,subsnes->jacobian,subsnes->jacobian_pre);
801: if (lag > 1) subsnes->lagjacobian = lag;
802: }
803: return(0);
804: }
808: PetscErrorCode SNESSolve_NASM(SNES snes)
809: {
810: Vec F;
811: Vec X;
812: Vec B;
813: Vec Y;
814: PetscInt i;
815: PetscReal fnorm = 0.0;
816: PetscErrorCode ierr;
817: SNESNormSchedule normschedule;
818: SNES_NASM *nasm = (SNES_NASM*)snes->data;
822: 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);
824: PetscCitationsRegister(SNESCitation,&SNEScite);
825: X = snes->vec_sol;
826: Y = snes->vec_sol_update;
827: F = snes->vec_func;
828: B = snes->vec_rhs;
830: PetscObjectSAWsTakeAccess((PetscObject)snes);
831: snes->iter = 0;
832: snes->norm = 0.;
833: PetscObjectSAWsGrantAccess((PetscObject)snes);
834: snes->reason = SNES_CONVERGED_ITERATING;
835: SNESGetNormSchedule(snes, &normschedule);
836: if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
837: /* compute the initial function and preconditioned update delX */
838: if (!snes->vec_func_init_set) {
839: SNESComputeFunction(snes,X,F);
840: } else snes->vec_func_init_set = PETSC_FALSE;
842: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
843: SNESCheckFunctionNorm(snes,fnorm);
844: PetscObjectSAWsTakeAccess((PetscObject)snes);
845: snes->iter = 0;
846: snes->norm = fnorm;
847: PetscObjectSAWsGrantAccess((PetscObject)snes);
848: SNESLogConvergenceHistory(snes,snes->norm,0);
849: SNESMonitor(snes,0,snes->norm);
851: /* test convergence */
852: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
853: if (snes->reason) return(0);
854: } else {
855: PetscObjectSAWsGrantAccess((PetscObject)snes);
856: SNESLogConvergenceHistory(snes,snes->norm,0);
857: SNESMonitor(snes,0,snes->norm);
858: }
860: /* Call general purpose update function */
861: if (snes->ops->update) {
862: (*snes->ops->update)(snes, snes->iter);
863: }
864: /* copy the initial solution over for later */
865: if (nasm->fjtype == 2) {VecCopy(X,nasm->xinit);}
867: for (i = 0; i < snes->max_its; i++) {
868: SNESNASMSolveLocal_Private(snes,B,Y,X);
869: if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
870: SNESComputeFunction(snes,X,F);
871: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
872: SNESCheckFunctionNorm(snes,fnorm);
873: }
874: /* Monitor convergence */
875: PetscObjectSAWsTakeAccess((PetscObject)snes);
876: snes->iter = i+1;
877: snes->norm = fnorm;
878: PetscObjectSAWsGrantAccess((PetscObject)snes);
879: SNESLogConvergenceHistory(snes,snes->norm,0);
880: SNESMonitor(snes,snes->iter,snes->norm);
881: /* Test for convergence */
882: if (normschedule == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);}
883: if (snes->reason) break;
884: /* Call general purpose update function */
885: if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
886: }
887: if (nasm->finaljacobian) {SNESNASMComputeFinalJacobian_Private(snes,X);}
888: if (normschedule == SNES_NORM_ALWAYS) {
889: if (i == snes->max_its) {
890: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
891: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
892: }
893: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */
894: return(0);
895: }
897: /*MC
898: SNESNASM - Nonlinear Additive Schwartz
900: Options Database:
901: + -snes_nasm_log - enable logging events for the communication and solve stages
902: . -snes_nasm_type <basic,restrict> - type of subdomain update used
903: . -snes_asm_damping <dmp> - the new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
904: . -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
905: . -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> - pick state the jacobian is calculated at
906: . -sub_snes_ - options prefix of the subdomain nonlinear solves
907: . -sub_ksp_ - options prefix of the subdomain Krylov solver
908: - -sub_pc_ - options prefix of the subdomain preconditioner
910: Level: advanced
912: References:
913: . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
914: SIAM Review, 57(4), 2015
916: .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types), SNESNASMSetType(), SNESNASMGetType(), SNESNASMSetSubdomains(), SNESNASMGetSubdomains(), SNESNASMGetSubdomainVecs(), SNESNASMSetComputeFinalJacobian(), SNESNASMSetDamping(), SNESNASMGetDamping()
917: M*/
921: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
922: {
923: SNES_NASM *nasm;
927: PetscNewLog(snes,&nasm);
928: snes->data = (void*)nasm;
930: nasm->n = PETSC_DECIDE;
931: nasm->subsnes = 0;
932: nasm->x = 0;
933: nasm->xl = 0;
934: nasm->y = 0;
935: nasm->b = 0;
936: nasm->oscatter = 0;
937: nasm->iscatter = 0;
938: nasm->gscatter = 0;
939: nasm->damping = 1.;
941: nasm->type = PC_ASM_BASIC;
942: nasm->finaljacobian = PETSC_FALSE;
943: nasm->same_local_solves = PETSC_TRUE;
945: snes->ops->destroy = SNESDestroy_NASM;
946: snes->ops->setup = SNESSetUp_NASM;
947: snes->ops->setfromoptions = SNESSetFromOptions_NASM;
948: snes->ops->view = SNESView_NASM;
949: snes->ops->solve = SNESSolve_NASM;
950: snes->ops->reset = SNESReset_NASM;
952: snes->usesksp = PETSC_FALSE;
953: snes->usespc = PETSC_FALSE;
955: nasm->fjtype = 0;
956: nasm->xinit = NULL;
957: nasm->eventrestrictinterp = 0;
958: nasm->eventsubsolve = 0;
960: if (!snes->tolerancesset) {
961: snes->max_its = 10000;
962: snes->max_funcs = 10000;
963: }
965: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetType_C",SNESNASMSetType_NASM);
966: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetType_C",SNESNASMGetType_NASM);
967: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
968: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
969: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
970: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
971: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
972: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
973: return(0);
974: }