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 == 0) {
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 == 0) {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 Parameter:
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 Parameter:
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: }
609: /*
610: Input Parameters:
611: + snes - The solver
612: . B - The RHS vector
613: - X - The initial guess
615: Output Parameters:
616: . Y - The solution update
618: TODO: All scatters should be packed into one
619: */
620: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X)
621: {
622: SNES_NASM *nasm = (SNES_NASM*)snes->data;
623: SNES subsnes;
624: PetscInt i;
625: PetscReal dmp;
627: Vec Xl,Bl,Yl,Xlloc;
628: VecScatter iscat,oscat,gscat,oscat_copy;
629: DM dm,subdm;
630: PCASMType type;
633: SNESNASMGetType(snes,&type);
634: SNESGetDM(snes,&dm);
635: VecSet(Y,0);
636: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
637: for (i=0; i<nasm->n; i++) {
638: /* scatter the solution to the global solution and the local solution */
639: Xl = nasm->x[i];
640: Xlloc = nasm->xl[i];
641: oscat = nasm->oscatter[i];
642: oscat_copy = nasm->oscatter_copy[i];
643: gscat = nasm->gscatter[i];
644: VecScatterBegin(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
645: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
646: if (B) {
647: /* scatter the RHS to the local RHS */
648: Bl = nasm->b[i];
649: VecScatterBegin(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
650: }
651: }
652: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
654: if (nasm->eventsubsolve) {PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);}
655: for (i=0; i<nasm->n; i++) {
656: Xl = nasm->x[i];
657: Xlloc = nasm->xl[i];
658: Yl = nasm->y[i];
659: subsnes = nasm->subsnes[i];
660: SNESGetDM(subsnes,&subdm);
661: iscat = nasm->iscatter[i];
662: oscat = nasm->oscatter[i];
663: oscat_copy = nasm->oscatter_copy[i];
664: gscat = nasm->gscatter[i];
665: VecScatterEnd(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
666: VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
667: if (B) {
668: Bl = nasm->b[i];
669: VecScatterEnd(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
670: } else Bl = NULL;
672: DMSubDomainRestrict(dm,oscat,gscat,subdm);
673: VecCopy(Xl,Yl);
674: SNESSolve(subsnes,Bl,Xl);
675: VecAYPX(Yl,-1.0,Xl);
676: VecScale(Yl, nasm->damping);
677: if (type == PC_ASM_BASIC) {
678: VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
679: VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
680: } else if (type == PC_ASM_RESTRICT) {
681: VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
682: VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
683: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
684: }
685: if (nasm->eventsubsolve) {PetscLogEventEnd(nasm->eventsubsolve,snes,0,0,0);}
686: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
687: if (nasm->weight_set) {
688: VecPointwiseMult(Y,Y,nasm->weight);
689: }
690: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
691: SNESNASMGetDamping(snes,&dmp);
692: VecAXPY(X,dmp,Y);
693: return(0);
694: }
696: static PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
697: {
698: Vec X = Xfinal;
699: SNES_NASM *nasm = (SNES_NASM*)snes->data;
700: SNES subsnes;
701: PetscInt i,lag = 1;
703: Vec Xlloc,Xl,Fl,F;
704: VecScatter oscat,gscat;
705: DM dm,subdm;
708: if (nasm->fjtype == 2) X = nasm->xinit;
709: F = snes->vec_func;
710: if (snes->normschedule == SNES_NORM_NONE) {SNESComputeFunction(snes,X,F);}
711: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
712: SNESGetDM(snes,&dm);
713: if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
714: if (nasm->fjtype != 1) {
715: for (i=0; i<nasm->n; i++) {
716: Xlloc = nasm->xl[i];
717: gscat = nasm->gscatter[i];
718: VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
719: }
720: }
721: if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
722: for (i=0; i<nasm->n; i++) {
723: Fl = nasm->subsnes[i]->vec_func;
724: Xl = nasm->x[i];
725: Xlloc = nasm->xl[i];
726: subsnes = nasm->subsnes[i];
727: oscat = nasm->oscatter[i];
728: gscat = nasm->gscatter[i];
729: if (nasm->fjtype != 1) {VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);}
730: SNESGetDM(subsnes,&subdm);
731: DMSubDomainRestrict(dm,oscat,gscat,subdm);
732: if (nasm->fjtype != 1) {
733: DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
734: DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
735: }
736: if (subsnes->lagjacobian == -1) subsnes->lagjacobian = -2;
737: else if (subsnes->lagjacobian > 1) lag = subsnes->lagjacobian;
738: SNESComputeFunction(subsnes,Xl,Fl);
739: SNESComputeJacobian(subsnes,Xl,subsnes->jacobian,subsnes->jacobian_pre);
740: if (lag > 1) subsnes->lagjacobian = lag;
741: }
742: return(0);
743: }
745: static PetscErrorCode SNESSolve_NASM(SNES snes)
746: {
747: Vec F;
748: Vec X;
749: Vec B;
750: Vec Y;
751: PetscInt i;
752: PetscReal fnorm = 0.0;
753: PetscErrorCode ierr;
754: SNESNormSchedule normschedule;
755: SNES_NASM *nasm = (SNES_NASM*)snes->data;
759: 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);
761: PetscCitationsRegister(SNESCitation,&SNEScite);
762: X = snes->vec_sol;
763: Y = snes->vec_sol_update;
764: F = snes->vec_func;
765: B = snes->vec_rhs;
767: PetscObjectSAWsTakeAccess((PetscObject)snes);
768: snes->iter = 0;
769: snes->norm = 0.;
770: PetscObjectSAWsGrantAccess((PetscObject)snes);
771: snes->reason = SNES_CONVERGED_ITERATING;
772: SNESGetNormSchedule(snes, &normschedule);
773: if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
774: /* compute the initial function and preconditioned update delX */
775: if (!snes->vec_func_init_set) {
776: SNESComputeFunction(snes,X,F);
777: } else snes->vec_func_init_set = PETSC_FALSE;
779: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
780: SNESCheckFunctionNorm(snes,fnorm);
781: PetscObjectSAWsTakeAccess((PetscObject)snes);
782: snes->iter = 0;
783: snes->norm = fnorm;
784: PetscObjectSAWsGrantAccess((PetscObject)snes);
785: SNESLogConvergenceHistory(snes,snes->norm,0);
786: SNESMonitor(snes,0,snes->norm);
788: /* test convergence */
789: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
790: if (snes->reason) return(0);
791: } else {
792: PetscObjectSAWsGrantAccess((PetscObject)snes);
793: SNESLogConvergenceHistory(snes,snes->norm,0);
794: SNESMonitor(snes,0,snes->norm);
795: }
797: /* Call general purpose update function */
798: if (snes->ops->update) {
799: (*snes->ops->update)(snes, snes->iter);
800: }
801: /* copy the initial solution over for later */
802: if (nasm->fjtype == 2) {VecCopy(X,nasm->xinit);}
804: for (i=0; i < snes->max_its; i++) {
805: SNESNASMSolveLocal_Private(snes,B,Y,X);
806: if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
807: SNESComputeFunction(snes,X,F);
808: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
809: SNESCheckFunctionNorm(snes,fnorm);
810: }
811: /* Monitor convergence */
812: PetscObjectSAWsTakeAccess((PetscObject)snes);
813: snes->iter = i+1;
814: snes->norm = fnorm;
815: PetscObjectSAWsGrantAccess((PetscObject)snes);
816: SNESLogConvergenceHistory(snes,snes->norm,0);
817: SNESMonitor(snes,snes->iter,snes->norm);
818: /* Test for convergence */
819: if (normschedule == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);}
820: if (snes->reason) break;
821: /* Call general purpose update function */
822: if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
823: }
824: if (nasm->finaljacobian) {
825: SNESNASMComputeFinalJacobian_Private(snes,X);
826: SNESCheckJacobianDomainerror(snes);
827: }
828: if (normschedule == SNES_NORM_ALWAYS) {
829: if (i == snes->max_its) {
830: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
831: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
832: }
833: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */
834: return(0);
835: }
837: /*MC
838: SNESNASM - Nonlinear Additive Schwarz
840: Options Database:
841: + -snes_nasm_log - enable logging events for the communication and solve stages
842: . -snes_nasm_type <basic,restrict> - type of subdomain update used
843: . -snes_asm_damping <dmp> - the new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
844: . -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
845: . -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> - pick state the jacobian is calculated at
846: . -sub_snes_ - options prefix of the subdomain nonlinear solves
847: . -sub_ksp_ - options prefix of the subdomain Krylov solver
848: - -sub_pc_ - options prefix of the subdomain preconditioner
850: Level: advanced
852: 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
853: false and SNESView() and -snes_view do not display a KSP object. However, if the flag nasm->finaljacobian is set (for example, if
854: NASM is used as a nonlinear preconditioner for KSPASPIN) then SNESSetUpMatrices() is called to generate the Jacobian (needed by KSPASPIN)
855: 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
856: used by SNESASPIN they share the same Jacobian matrices because SNESSetUp() (called on the outer SNES KSPASPIN) causes the inner SNES
857: object (in this case SNESNASM) to inherit the outer Jacobian matrices.
859: References:
860: . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
861: SIAM Review, 57(4), 2015
863: .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types), SNESNASMSetType(), SNESNASMGetType(), SNESNASMSetSubdomains(), SNESNASMGetSubdomains(), SNESNASMGetSubdomainVecs(), SNESNASMSetComputeFinalJacobian(), SNESNASMSetDamping(), SNESNASMGetDamping()
864: M*/
866: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
867: {
868: SNES_NASM *nasm;
872: PetscNewLog(snes,&nasm);
873: snes->data = (void*)nasm;
875: nasm->n = PETSC_DECIDE;
876: nasm->subsnes = NULL;
877: nasm->x = NULL;
878: nasm->xl = NULL;
879: nasm->y = NULL;
880: nasm->b = NULL;
881: nasm->oscatter = NULL;
882: nasm->oscatter_copy = NULL;
883: nasm->iscatter = NULL;
884: nasm->gscatter = NULL;
885: nasm->damping = 1.;
887: nasm->type = PC_ASM_BASIC;
888: nasm->finaljacobian = PETSC_FALSE;
889: nasm->weight_set = PETSC_FALSE;
891: snes->ops->destroy = SNESDestroy_NASM;
892: snes->ops->setup = SNESSetUp_NASM;
893: snes->ops->setfromoptions = SNESSetFromOptions_NASM;
894: snes->ops->view = SNESView_NASM;
895: snes->ops->solve = SNESSolve_NASM;
896: snes->ops->reset = SNESReset_NASM;
898: snes->usesksp = PETSC_FALSE;
899: snes->usesnpc = PETSC_FALSE;
901: snes->alwayscomputesfinalresidual = PETSC_FALSE;
903: nasm->fjtype = 0;
904: nasm->xinit = NULL;
905: nasm->eventrestrictinterp = 0;
906: nasm->eventsubsolve = 0;
908: if (!snes->tolerancesset) {
909: snes->max_its = 10000;
910: snes->max_funcs = 10000;
911: }
913: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetType_C",SNESNASMSetType_NASM);
914: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetType_C",SNESNASMGetType_NASM);
915: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
916: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
917: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
918: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
919: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
920: PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
921: return(0);
922: }
924: /*@
925: SNESNASMGetSNES - Gets a subsolver
927: Not collective
929: Input Parameters:
930: + snes - the SNES context
931: - i - the number of the subsnes to get
933: Output Parameters:
934: . subsnes - the subsolver context
936: Level: intermediate
938: .seealso: SNESNASM, SNESNASMGetNumber()
939: @*/
940: PetscErrorCode SNESNASMGetSNES(SNES snes,PetscInt i,SNES *subsnes)
941: {
942: SNES_NASM *nasm = (SNES_NASM*)snes->data;
945: if (i < 0 || i >= nasm->n) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"No such subsolver");
946: *subsnes = nasm->subsnes[i];
947: return(0);
948: }
950: /*@
951: SNESNASMGetNumber - Gets number of subsolvers
953: Not collective
955: Input Parameters:
956: . snes - the SNES context
958: Output Parameters:
959: . n - the number of subsolvers
961: Level: intermediate
963: .seealso: SNESNASM, SNESNASMGetSNES()
964: @*/
965: PetscErrorCode SNESNASMGetNumber(SNES snes,PetscInt *n)
966: {
967: SNES_NASM *nasm = (SNES_NASM*)snes->data;
970: *n = nasm->n;
971: return(0);
972: }
974: /*@
975: SNESNASMSetWeight - Sets weight to use when adding overlapping updates
977: Collective
979: Input Parameters:
980: + snes - the SNES context
981: - weight - the weights to use (typically 1/N for each dof, where N is the number of patches it appears in)
983: Level: intermediate
985: .seealso: SNESNASM
986: @*/
987: PetscErrorCode SNESNASMSetWeight(SNES snes,Vec weight)
988: {
989: SNES_NASM *nasm = (SNES_NASM*)snes->data;
994: VecDestroy(&nasm->weight);
995: nasm->weight_set = PETSC_TRUE;
996: nasm->weight = weight;
997: PetscObjectReference((PetscObject)nasm->weight);
999: return(0);
1000: }