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;
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->oscatter_copy) { VecScatterDestroy(&nasm->oscatter_copy[i]); }
49: if (nasm->iscatter) { VecScatterDestroy(&nasm->iscatter[i]); }
50: if (nasm->gscatter) { VecScatterDestroy(&nasm->gscatter[i]); }
51: }
53: PetscFree(nasm->x);
54: PetscFree(nasm->xl);
55: PetscFree(nasm->y);
56: PetscFree(nasm->b);
58: if (nasm->xinit) {VecDestroy(&nasm->xinit);}
60: PetscFree(nasm->subsnes);
61: PetscFree(nasm->oscatter);
62: PetscFree(nasm->oscatter_copy);
63: PetscFree(nasm->iscatter);
64: PetscFree(nasm->gscatter);
66: if (nasm->weight_set) {
67: VecDestroy(&nasm->weight);
68: }
70: nasm->eventrestrictinterp = 0;
71: nasm->eventsubsolve = 0;
72: return(0);
73: }
75: static PetscErrorCode SNESDestroy_NASM(SNES snes)
76: {
80: SNESReset_NASM(snes);
81: PetscFree(snes->data);
82: return(0);
83: }
85: static PetscErrorCode DMGlobalToLocalSubDomainDirichletHook_Private(DM dm,Vec g,InsertMode mode,Vec l,void *ctx)
86: {
88: Vec bcs = (Vec)ctx;
91: VecCopy(bcs,l);
92: return(0);
93: }
95: static 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);
116: PetscMalloc1(nasm->n, &nasm->oscatter_copy);
117: for (i=0; i<nasm->n; i++) {
118: VecScatterCopy(nasm->oscatter[i], &nasm->oscatter_copy[i]);
119: }
121: SNESGetOptionsPrefix(snes, &optionsprefix);
122: PetscMalloc1(nasm->n,&nasm->subsnes);
123: for (i=0; i<nasm->n; i++) {
124: SNESCreate(PETSC_COMM_SELF,&nasm->subsnes[i]);
125: PetscObjectIncrementTabLevel((PetscObject)nasm->subsnes[i], (PetscObject)snes, 1);
126: SNESAppendOptionsPrefix(nasm->subsnes[i],optionsprefix);
127: SNESAppendOptionsPrefix(nasm->subsnes[i],"sub_");
128: SNESSetDM(nasm->subsnes[i],subdms[i]);
129: MPI_Comm_size(PetscObjectComm((PetscObject)nasm->subsnes[i]),&size);
130: if (size == 1) {
131: SNESGetKSP(nasm->subsnes[i],&ksp);
132: KSPGetPC(ksp,&pc);
133: KSPSetType(ksp,KSPPREONLY);
134: PCSetType(pc,PCLU);
135: }
136: SNESSetFromOptions(nasm->subsnes[i]);
137: DMDestroy(&subdms[i]);
138: }
139: PetscFree(subdms);
140: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Cannot construct local problems automatically without a DM!");
141: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set subproblems manually if there is no DM!");
142: /* allocate the global vectors */
143: if (!nasm->x) {
144: PetscCalloc1(nasm->n,&nasm->x);
145: }
146: if (!nasm->xl) {
147: PetscCalloc1(nasm->n,&nasm->xl);
148: }
149: if (!nasm->y) {
150: PetscCalloc1(nasm->n,&nasm->y);
151: }
152: if (!nasm->b) {
153: PetscCalloc1(nasm->n,&nasm->b);
154: }
156: for (i=0; i<nasm->n; i++) {
157: SNESGetFunction(nasm->subsnes[i],&F,NULL,NULL);
158: if (!nasm->x[i]) {VecDuplicate(F,&nasm->x[i]);}
159: if (!nasm->y[i]) {VecDuplicate(F,&nasm->y[i]);}
160: if (!nasm->b[i]) {VecDuplicate(F,&nasm->b[i]);}
161: if (!nasm->xl[i]) {
162: SNESGetDM(nasm->subsnes[i],&subdm);
163: DMCreateLocalVector(subdm,&nasm->xl[i]);
164: DMGlobalToLocalHookAdd(subdm,DMGlobalToLocalSubDomainDirichletHook_Private,NULL,nasm->xl[i]);
165: }
166: }
167: if (nasm->finaljacobian) {
168: SNESSetUpMatrices(snes);
169: if (nasm->fjtype == 2) {
170: VecDuplicate(snes->vec_sol,&nasm->xinit);
171: }
172: for (i=0; i<nasm->n;i++) {
173: SNESSetUpMatrices(nasm->subsnes[i]);
174: }
175: }
176: return(0);
177: }
179: static PetscErrorCode SNESSetFromOptions_NASM(PetscOptionItems *PetscOptionsObject,SNES snes)
180: {
181: PetscErrorCode ierr;
182: PCASMType asmtype;
183: PetscBool flg,monflg;
184: SNES_NASM *nasm = (SNES_NASM*)snes->data;
187: PetscOptionsHead(PetscOptionsObject,"Nonlinear Additive Schwarz options");
188: PetscOptionsEnum("-snes_nasm_type","Type of restriction/extension","",SNESNASMTypes,(PetscEnum)nasm->type,(PetscEnum*)&asmtype,&flg);
189: if (flg) {SNESNASMSetType(snes,asmtype);}
190: flg = PETSC_FALSE;
191: monflg = PETSC_TRUE;
192: 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);
193: if (flg) {SNESNASMSetDamping(snes,nasm->damping);}
194: PetscOptionsDeprecated("-snes_nasm_sub_view",NULL,"3.15","Use -snes_view ::ascii_info_detail");
195: PetscOptionsBool("-snes_nasm_finaljacobian","Compute the global jacobian of the final iterate (for ASPIN)","",nasm->finaljacobian,&nasm->finaljacobian,NULL);
196: PetscOptionsEList("-snes_nasm_finaljacobian_type","The type of the final jacobian computed.","",SNESNASMFJTypes,3,SNESNASMFJTypes[0],&nasm->fjtype,NULL);
197: PetscOptionsBool("-snes_nasm_log","Log times for subSNES solves and restriction","",monflg,&monflg,&flg);
198: if (flg) {
199: PetscLogEventRegister("SNESNASMSubSolve",((PetscObject)snes)->classid,&nasm->eventsubsolve);
200: PetscLogEventRegister("SNESNASMRestrict",((PetscObject)snes)->classid,&nasm->eventrestrictinterp);
201: }
202: PetscOptionsTail();
203: return(0);
204: }
206: static PetscErrorCode SNESView_NASM(SNES snes, PetscViewer viewer)
207: {
208: SNES_NASM *nasm = (SNES_NASM*)snes->data;
209: PetscErrorCode ierr;
210: PetscMPIInt rank,size;
211: PetscInt i,N,bsz;
212: PetscBool iascii,isstring;
213: PetscViewer sviewer;
214: MPI_Comm comm;
215: PetscViewerFormat format;
216: const char *prefix;
219: PetscObjectGetComm((PetscObject)snes,&comm);
220: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
221: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
222: MPI_Comm_rank(comm,&rank);
223: MPI_Comm_size(comm,&size);
224: MPIU_Allreduce(&nasm->n,&N,1,MPIU_INT,MPI_SUM,comm);
225: if (iascii) {
226: PetscViewerASCIIPrintf(viewer, " total subdomain blocks = %D\n",N);
227: PetscViewerGetFormat(viewer,&format);
228: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
229: if (nasm->subsnes) {
230: PetscViewerASCIIPrintf(viewer," Local solver information for first block on rank 0:\n");
231: SNESGetOptionsPrefix(snes,&prefix);
232: PetscViewerASCIIPrintf(viewer," Use -%ssnes_view ::ascii_info_detail to display information for all blocks\n",prefix?prefix:"");
233: PetscViewerASCIIPushTab(viewer);
234: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
235: if (!rank) {
236: PetscViewerASCIIPushTab(viewer);
237: SNESView(nasm->subsnes[0],sviewer);
238: PetscViewerASCIIPopTab(viewer);
239: }
240: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
241: PetscViewerASCIIPopTab(viewer);
242: }
243: } else {
244: /* print the solver on each block */
245: PetscViewerASCIIPushSynchronized(viewer);
246: PetscViewerASCIISynchronizedPrintf(viewer," [%d] number of local blocks = %D\n",(int)rank,nasm->n);
247: PetscViewerFlush(viewer);
248: PetscViewerASCIIPopSynchronized(viewer);
249: PetscViewerASCIIPrintf(viewer," Local solver information for each block is in the following SNES objects:\n");
250: PetscViewerASCIIPushTab(viewer);
251: PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
252: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
253: for (i=0; i<nasm->n; i++) {
254: VecGetLocalSize(nasm->x[i],&bsz);
255: PetscViewerASCIIPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);
256: SNESView(nasm->subsnes[i],sviewer);
257: PetscViewerASCIIPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
258: }
259: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
260: PetscViewerFlush(viewer);
261: PetscViewerASCIIPopTab(viewer);
262: }
263: } else if (isstring) {
264: PetscViewerStringSPrintf(viewer," blocks=%D,type=%s",N,SNESNASMTypes[nasm->type]);
265: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
266: if (nasm->subsnes && !rank) {SNESView(nasm->subsnes[0],sviewer);}
267: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
268: }
269: return(0);
270: }
272: /*@
273: SNESNASMSetType - Set the type of subdomain update used
275: Logically Collective on SNES
277: Input Parameters:
278: + SNES - the SNES context
279: - type - the type of update, PC_ASM_BASIC or PC_ASM_RESTRICT
281: Level: intermediate
283: .seealso: SNESNASM, SNESNASMGetType(), PCASMSetType()
284: @*/
285: PetscErrorCode SNESNASMSetType(SNES snes,PCASMType type)
286: {
288: PetscErrorCode (*f)(SNES,PCASMType);
291: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetType_C",&f);
292: if (f) {(f)(snes,type);}
293: return(0);
294: }
296: static PetscErrorCode SNESNASMSetType_NASM(SNES snes,PCASMType type)
297: {
298: SNES_NASM *nasm = (SNES_NASM*)snes->data;
301: if (type != PC_ASM_BASIC && type != PC_ASM_RESTRICT) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"SNESNASM only supports basic and restrict types");
302: nasm->type = type;
303: return(0);
304: }
306: /*@
307: SNESNASMGetType - Get the type of subdomain update used
309: Logically Collective on SNES
311: Input Parameters:
312: . SNES - the SNES context
314: Output Parameters:
315: . type - the type of update
317: Level: intermediate
319: .seealso: SNESNASM, SNESNASMSetType(), PCASMGetType()
320: @*/
321: PetscErrorCode SNESNASMGetType(SNES snes,PCASMType *type)
322: {
326: PetscUseMethod(snes,"SNESNASMGetType_C",(SNES,PCASMType*),(snes,type));
327: return(0);
328: }
330: static PetscErrorCode SNESNASMGetType_NASM(SNES snes,PCASMType *type)
331: {
332: SNES_NASM *nasm = (SNES_NASM*)snes->data;
335: *type = nasm->type;
336: return(0);
337: }
339: /*@
340: SNESNASMSetSubdomains - Manually Set the context required to restrict and solve subdomain problems.
342: Not Collective
344: Input Parameters:
345: + SNES - the SNES context
346: . n - the number of local subdomains
347: . subsnes - solvers defined on the local subdomains
348: . iscatter - scatters into the nonoverlapping portions of the local subdomains
349: . oscatter - scatters into the overlapping portions of the local subdomains
350: - gscatter - scatters into the (ghosted) local vector of the local subdomain
352: Level: intermediate
354: .seealso: SNESNASM, SNESNASMGetSubdomains()
355: @*/
356: PetscErrorCode SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
357: {
359: PetscErrorCode (*f)(SNES,PetscInt,SNES*,VecScatter*,VecScatter*,VecScatter*);
362: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",&f);
363: if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
364: return(0);
365: }
367: static PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
368: {
369: PetscInt i;
371: SNES_NASM *nasm = (SNES_NASM*)snes->data;
374: if (snes->setupcalled) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"SNESNASMSetSubdomains() should be called before calling SNESSetUp().");
376: /* tear down the previously set things */
377: SNESReset(snes);
379: nasm->n = n;
380: if (oscatter) {
381: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)oscatter[i]);}
382: }
383: if (iscatter) {
384: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)iscatter[i]);}
385: }
386: if (gscatter) {
387: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)gscatter[i]);}
388: }
389: if (oscatter) {
390: PetscMalloc1(n,&nasm->oscatter);
391: PetscMalloc1(n,&nasm->oscatter_copy);
392: for (i=0; i<n; i++) {
393: nasm->oscatter[i] = oscatter[i];
394: VecScatterCopy(oscatter[i], &nasm->oscatter_copy[i]);
395: }
396: }
397: if (iscatter) {
398: PetscMalloc1(n,&nasm->iscatter);
399: for (i=0; i<n; i++) {
400: nasm->iscatter[i] = iscatter[i];
401: }
402: }
403: if (gscatter) {
404: PetscMalloc1(n,&nasm->gscatter);
405: for (i=0; i<n; i++) {
406: nasm->gscatter[i] = gscatter[i];
407: }
408: }
410: if (subsnes) {
411: PetscMalloc1(n,&nasm->subsnes);
412: for (i=0; i<n; i++) {
413: nasm->subsnes[i] = subsnes[i];
414: }
415: }
416: return(0);
417: }
419: /*@
420: SNESNASMGetSubdomains - Get the local subdomain context.
422: Not Collective
424: Input Parameters:
425: . SNES - the SNES context
427: Output Parameters:
428: + n - the number of local subdomains
429: . subsnes - solvers defined on the local subdomains
430: . iscatter - scatters into the nonoverlapping portions of the local subdomains
431: . oscatter - scatters into the overlapping portions of the local subdomains
432: - gscatter - scatters into the (ghosted) local vector of the local subdomain
434: Level: intermediate
436: .seealso: SNESNASM, SNESNASMSetSubdomains()
437: @*/
438: PetscErrorCode SNESNASMGetSubdomains(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
439: {
441: PetscErrorCode (*f)(SNES,PetscInt*,SNES**,VecScatter**,VecScatter**,VecScatter**);
444: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",&f);
445: if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
446: return(0);
447: }
449: static PetscErrorCode SNESNASMGetSubdomains_NASM(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
450: {
451: SNES_NASM *nasm = (SNES_NASM*)snes->data;
454: if (n) *n = nasm->n;
455: if (oscatter) *oscatter = nasm->oscatter;
456: if (iscatter) *iscatter = nasm->iscatter;
457: if (gscatter) *gscatter = nasm->gscatter;
458: if (subsnes) *subsnes = nasm->subsnes;
459: return(0);
460: }
462: /*@
463: SNESNASMGetSubdomainVecs - Get the processor-local subdomain vectors
465: Not Collective
467: Input Parameters:
468: . SNES - the SNES context
470: Output Parameters:
471: + n - the number of local subdomains
472: . x - The subdomain solution vector
473: . y - The subdomain step vector
474: . b - The subdomain RHS vector
475: - xl - The subdomain local vectors (ghosted)
477: Level: developer
479: .seealso: SNESNASM, SNESNASMGetSubdomains()
480: @*/
481: PetscErrorCode SNESNASMGetSubdomainVecs(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b, Vec **xl)
482: {
484: PetscErrorCode (*f)(SNES,PetscInt*,Vec**,Vec**,Vec**,Vec**);
487: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",&f);
488: if (f) {(f)(snes,n,x,y,b,xl);}
489: return(0);
490: }
492: static PetscErrorCode SNESNASMGetSubdomainVecs_NASM(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b,Vec **xl)
493: {
494: SNES_NASM *nasm = (SNES_NASM*)snes->data;
497: if (n) *n = nasm->n;
498: if (x) *x = nasm->x;
499: if (y) *y = nasm->y;
500: if (b) *b = nasm->b;
501: if (xl) *xl = nasm->xl;
502: return(0);
503: }
505: /*@
506: SNESNASMSetComputeFinalJacobian - Schedules the computation of the global and subdomain Jacobians upon convergence
508: Collective on SNES
510: Input Parameters:
511: + SNES - the SNES context
512: - flg - indication of whether to compute the Jacobians or not
514: Level: developer
516: Notes:
517: This is used almost exclusively in the implementation of ASPIN, where the converged subdomain and global Jacobian
518: is needed at each linear iteration.
520: .seealso: SNESNASM, SNESNASMGetSubdomains()
521: @*/
522: PetscErrorCode SNESNASMSetComputeFinalJacobian(SNES snes,PetscBool flg)
523: {
524: PetscErrorCode (*f)(SNES,PetscBool);
528: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",&f);
529: if (f) {(f)(snes,flg);}
530: return(0);
531: }
533: static PetscErrorCode SNESNASMSetComputeFinalJacobian_NASM(SNES snes,PetscBool flg)
534: {
535: SNES_NASM *nasm = (SNES_NASM*)snes->data;
538: nasm->finaljacobian = flg;
539: return(0);
540: }
542: /*@
543: SNESNASMSetDamping - Sets the update damping for NASM
545: Logically collective on SNES
547: Input Parameters:
548: + SNES - the SNES context
549: - dmp - damping
551: Level: intermediate
553: Notes:
554: The new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
556: .seealso: SNESNASM, SNESNASMGetDamping()
557: @*/
558: PetscErrorCode SNESNASMSetDamping(SNES snes,PetscReal dmp)
559: {
560: PetscErrorCode (*f)(SNES,PetscReal);
564: PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetDamping_C",(void (**)(void))&f);
565: if (f) {(f)(snes,dmp);}
566: return(0);
567: }
569: static PetscErrorCode SNESNASMSetDamping_NASM(SNES snes,PetscReal dmp)
570: {
571: SNES_NASM *nasm = (SNES_NASM*)snes->data;
574: nasm->damping = dmp;
575: return(0);
576: }
578: /*@
579: SNESNASMGetDamping - Gets the update damping for NASM
581: Not Collective
583: Input Parameters:
584: + SNES - the SNES context
585: - dmp - damping
587: Level: intermediate
589: .seealso: SNESNASM, SNESNASMSetDamping()
590: @*/
591: PetscErrorCode SNESNASMGetDamping(SNES snes,PetscReal *dmp)
592: {
596: PetscUseMethod(snes,"SNESNASMGetDamping_C",(SNES,PetscReal*),(snes,dmp));
597: return(0);
598: }
600: static PetscErrorCode SNESNASMGetDamping_NASM(SNES snes,PetscReal *dmp)
601: {
602: SNES_NASM *nasm = (SNES_NASM*)snes->data;
605: *dmp = nasm->damping;
606: return(0);
607: }
610: /*
611: Input Parameters:
612: + snes - The solver
613: . B - The RHS vector
614: - X - The initial guess
616: Output Parameters:
617: . Y - The solution update
619: TODO: All scatters should be packed into one
620: */
621: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X)
622: {
623: SNES_NASM *nasm = (SNES_NASM*)snes->data;
624: SNES subsnes;
625: PetscInt i;
626: PetscReal dmp;
628: Vec Xl,Bl,Yl,Xlloc;
629: VecScatter iscat,oscat,gscat,oscat_copy;
630: DM dm,subdm;
631: PCASMType type;
634: SNESNASMGetType(snes,&type);
635: SNESGetDM(snes,&dm);
636: VecSet(Y,0);
637: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
638: for (i=0; i<nasm->n; i++) {
639: /* scatter the solution to the global solution and the local solution */
640: Xl = nasm->x[i];
641: Xlloc = nasm->xl[i];
642: oscat = nasm->oscatter[i];
643: oscat_copy = nasm->oscatter_copy[i];
644: gscat = nasm->gscatter[i];
645: VecScatterBegin(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
646: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
647: if (B) {
648: /* scatter the RHS to the local RHS */
649: Bl = nasm->b[i];
650: VecScatterBegin(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
651: }
652: }
653: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
656: if (nasm->eventsubsolve) {PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);}
657: for (i=0; i<nasm->n; i++) {
658: Xl = nasm->x[i];
659: Xlloc = nasm->xl[i];
660: Yl = nasm->y[i];
661: subsnes = nasm->subsnes[i];
662: SNESGetDM(subsnes,&subdm);
663: iscat = nasm->iscatter[i];
664: oscat = nasm->oscatter[i];
665: oscat_copy = nasm->oscatter_copy[i];
666: gscat = nasm->gscatter[i];
667: VecScatterEnd(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
668: VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
669: if (B) {
670: Bl = nasm->b[i];
671: VecScatterEnd(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
672: } else Bl = NULL;
674: DMSubDomainRestrict(dm,oscat,gscat,subdm);
675: VecCopy(Xl,Yl);
676: SNESSolve(subsnes,Bl,Xl);
677: VecAYPX(Yl,-1.0,Xl);
678: VecScale(Yl, nasm->damping);
679: if (type == PC_ASM_BASIC) {
680: VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
681: VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
682: } else if (type == PC_ASM_RESTRICT) {
683: VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
684: VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
685: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
686: }
687: if (nasm->eventsubsolve) {PetscLogEventEnd(nasm->eventsubsolve,snes,0,0,0);}
688: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
689: if (nasm->weight_set) {
690: VecPointwiseMult(Y,Y,nasm->weight);
691: }
692: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
693: SNESNASMGetDamping(snes,&dmp);
694: VecAXPY(X,dmp,Y);
695: return(0);
696: }
698: static PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
699: {
700: Vec X = Xfinal;
701: SNES_NASM *nasm = (SNES_NASM*)snes->data;
702: SNES subsnes;
703: PetscInt i,lag = 1;
705: Vec Xlloc,Xl,Fl,F;
706: VecScatter oscat,gscat;
707: DM dm,subdm;
710: if (nasm->fjtype == 2) X = nasm->xinit;
711: F = snes->vec_func;
712: if (snes->normschedule == SNES_NORM_NONE) {SNESComputeFunction(snes,X,F);}
713: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
714: SNESGetDM(snes,&dm);
715: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
716: if (nasm->fjtype != 1) {
717: for (i=0; i<nasm->n; i++) {
718: Xlloc = nasm->xl[i];
719: gscat = nasm->gscatter[i];
720: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
721: }
722: }
723: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
724: for (i=0; i<nasm->n; i++) {
725: Fl = nasm->subsnes[i]->vec_func;
726: Xl = nasm->x[i];
727: Xlloc = nasm->xl[i];
728: subsnes = nasm->subsnes[i];
729: oscat = nasm->oscatter[i];
730: gscat = nasm->gscatter[i];
731: if (nasm->fjtype != 1) {VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);}
732: SNESGetDM(subsnes,&subdm);
733: DMSubDomainRestrict(dm,oscat,gscat,subdm);
734: if (nasm->fjtype != 1) {
735: DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
736: DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
737: }
738: if (subsnes->lagjacobian == -1) subsnes->lagjacobian = -2;
739: else if (subsnes->lagjacobian > 1) lag = subsnes->lagjacobian;
740: SNESComputeFunction(subsnes,Xl,Fl);
741: SNESComputeJacobian(subsnes,Xl,subsnes->jacobian,subsnes->jacobian_pre);
742: if (lag > 1) subsnes->lagjacobian = lag;
743: }
744: return(0);
745: }
747: static PetscErrorCode SNESSolve_NASM(SNES snes)
748: {
749: Vec F;
750: Vec X;
751: Vec B;
752: Vec Y;
753: PetscInt i;
754: PetscReal fnorm = 0.0;
755: PetscErrorCode ierr;
756: SNESNormSchedule normschedule;
757: SNES_NASM *nasm = (SNES_NASM*)snes->data;
761: 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);
763: PetscCitationsRegister(SNESCitation,&SNEScite);
764: X = snes->vec_sol;
765: Y = snes->vec_sol_update;
766: F = snes->vec_func;
767: B = snes->vec_rhs;
769: PetscObjectSAWsTakeAccess((PetscObject)snes);
770: snes->iter = 0;
771: snes->norm = 0.;
772: PetscObjectSAWsGrantAccess((PetscObject)snes);
773: snes->reason = SNES_CONVERGED_ITERATING;
774: SNESGetNormSchedule(snes, &normschedule);
775: if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
776: /* compute the initial function and preconditioned update delX */
777: if (!snes->vec_func_init_set) {
778: SNESComputeFunction(snes,X,F);
779: } else snes->vec_func_init_set = PETSC_FALSE;
781: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
782: SNESCheckFunctionNorm(snes,fnorm);
783: PetscObjectSAWsTakeAccess((PetscObject)snes);
784: snes->iter = 0;
785: snes->norm = fnorm;
786: PetscObjectSAWsGrantAccess((PetscObject)snes);
787: SNESLogConvergenceHistory(snes,snes->norm,0);
788: SNESMonitor(snes,0,snes->norm);
790: /* test convergence */
791: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
792: if (snes->reason) return(0);
793: } else {
794: PetscObjectSAWsGrantAccess((PetscObject)snes);
795: SNESLogConvergenceHistory(snes,snes->norm,0);
796: SNESMonitor(snes,0,snes->norm);
797: }
799: /* Call general purpose update function */
800: if (snes->ops->update) {
801: (*snes->ops->update)(snes, snes->iter);
802: }
803: /* copy the initial solution over for later */
804: if (nasm->fjtype == 2) {VecCopy(X,nasm->xinit);}
806: for (i=0; i < snes->max_its; i++) {
807: SNESNASMSolveLocal_Private(snes,B,Y,X);
808: if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
809: SNESComputeFunction(snes,X,F);
810: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
811: SNESCheckFunctionNorm(snes,fnorm);
812: }
813: /* Monitor convergence */
814: PetscObjectSAWsTakeAccess((PetscObject)snes);
815: snes->iter = i+1;
816: snes->norm = fnorm;
817: PetscObjectSAWsGrantAccess((PetscObject)snes);
818: SNESLogConvergenceHistory(snes,snes->norm,0);
819: SNESMonitor(snes,snes->iter,snes->norm);
820: /* Test for convergence */
821: if (normschedule == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);}
822: if (snes->reason) break;
823: /* Call general purpose update function */
824: if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
825: }
826: if (nasm->finaljacobian) {
827: SNESNASMComputeFinalJacobian_Private(snes,X);
828: SNESCheckJacobianDomainerror(snes);
829: }
830: if (normschedule == SNES_NORM_ALWAYS) {
831: if (i == snes->max_its) {
832: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
833: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
834: }
835: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */
836: return(0);
837: }
839: /*MC
840: SNESNASM - Nonlinear Additive Schwarz
842: Options Database:
843: + -snes_nasm_log - enable logging events for the communication and solve stages
844: . -snes_nasm_type <basic,restrict> - type of subdomain update used
845: . -snes_asm_damping <dmp> - the new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
846: . -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
847: . -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> - pick state the jacobian is calculated at
848: . -sub_snes_ - options prefix of the subdomain nonlinear solves
849: . -sub_ksp_ - options prefix of the subdomain Krylov solver
850: - -sub_pc_ - options prefix of the subdomain preconditioner
852: Level: advanced
854: 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
855: false and SNESView() and -snes_view do not display a KSP object. However, if the flag nasm->finaljacobian is set (for example, if
856: NASM is used as a nonlinear preconditioner for KSPASPIN) then SNESSetUpMatrices() is called to generate the Jacobian (needed by KSPASPIN)
857: 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
858: used by SNESASPIN they share the same Jacobian matrices because SNESSetUp() (called on the outer SNES KSPASPIN) causes the inner SNES
859: object (in this case SNESNASM) to inherit the outer Jacobian matrices.
861: References:
862: . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
863: SIAM Review, 57(4), 2015
865: .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types), SNESNASMSetType(), SNESNASMGetType(), SNESNASMSetSubdomains(), SNESNASMGetSubdomains(), SNESNASMGetSubdomainVecs(), SNESNASMSetComputeFinalJacobian(), SNESNASMSetDamping(), SNESNASMGetDamping()
866: M*/
868: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
869: {
870: SNES_NASM *nasm;
874: PetscNewLog(snes,&nasm);
875: snes->data = (void*)nasm;
877: nasm->n = PETSC_DECIDE;
878: nasm->subsnes = NULL;
879: nasm->x = NULL;
880: nasm->xl = NULL;
881: nasm->y = NULL;
882: nasm->b = NULL;
883: nasm->oscatter = NULL;
884: nasm->oscatter_copy = NULL;
885: nasm->iscatter = NULL;
886: nasm->gscatter = NULL;
887: nasm->damping = 1.;
889: nasm->type = PC_ASM_BASIC;
890: nasm->finaljacobian = PETSC_FALSE;
891: nasm->weight_set = PETSC_FALSE;
893: snes->ops->destroy = SNESDestroy_NASM;
894: snes->ops->setup = SNESSetUp_NASM;
895: snes->ops->setfromoptions = SNESSetFromOptions_NASM;
896: snes->ops->view = SNESView_NASM;
897: snes->ops->solve = SNESSolve_NASM;
898: snes->ops->reset = SNESReset_NASM;
900: snes->usesksp = PETSC_FALSE;
901: snes->usesnpc = PETSC_FALSE;
903: snes->alwayscomputesfinalresidual = PETSC_FALSE;
905: nasm->fjtype = 0;
906: nasm->xinit = NULL;
907: nasm->eventrestrictinterp = 0;
908: nasm->eventsubsolve = 0;
910: if (!snes->tolerancesset) {
911: snes->max_its = 10000;
912: snes->max_funcs = 10000;
913: }
915: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetType_C",SNESNASMSetType_NASM);
916: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetType_C",SNESNASMGetType_NASM);
917: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
918: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
919: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
920: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
921: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
922: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
923: return(0);
924: }
926: /*@
927: SNESNASMGetSNES - Gets a subsolver
929: Not collective
931: Input Parameters:
932: + snes - the SNES context
933: - i - the number of the subsnes to get
935: Output Parameters:
936: . subsnes - the subsolver context
938: Level: intermediate
940: .seealso: SNESNASM, SNESNASMGetNumber()
941: @*/
942: PetscErrorCode SNESNASMGetSNES(SNES snes,PetscInt i,SNES *subsnes)
943: {
944: SNES_NASM *nasm = (SNES_NASM*)snes->data;
947: if (i < 0 || i >= nasm->n) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"No such subsolver");
948: *subsnes = nasm->subsnes[i];
949: return(0);
950: }
952: /*@
953: SNESNASMGetNumber - Gets number of subsolvers
955: Not collective
957: Input Parameters:
958: . snes - the SNES context
960: Output Parameters:
961: . n - the number of subsolvers
963: Level: intermediate
965: .seealso: SNESNASM, SNESNASMGetSNES()
966: @*/
967: PetscErrorCode SNESNASMGetNumber(SNES snes,PetscInt *n)
968: {
969: SNES_NASM *nasm = (SNES_NASM*)snes->data;
972: *n = nasm->n;
973: return(0);
974: }
976: /*@
977: SNESNASMSetWeight - Sets weight to use when adding overlapping updates
979: Collective
981: Input Parameters:
982: + snes - the SNES context
983: - weight - the weights to use (typically 1/N for each dof, where N is the number of patches it appears in)
985: Level: intermediate
987: .seealso: SNESNASM
988: @*/
989: PetscErrorCode SNESNASMSetWeight(SNES snes,Vec weight)
990: {
991: SNES_NASM *nasm = (SNES_NASM*)snes->data;
996: VecDestroy(&nasm->weight);
997: nasm->weight_set = PETSC_TRUE;
998: nasm->weight = weight;
999: PetscObjectReference((PetscObject)nasm->weight);
1001: return(0);
1002: }