Actual source code: nasm.c
1: #include <petsc/private/snesimpl.h>
2: #include <petscdm.h>
4: typedef struct {
5: PetscInt n; /* local subdomains */
6: SNES *subsnes; /* nonlinear solvers for each subdomain */
7: Vec *x; /* solution vectors */
8: Vec *xl; /* solution local vectors */
9: Vec *y; /* step vectors */
10: Vec *b; /* rhs vectors */
11: Vec weight; /* weighting for adding updates on overlaps, in global space */
12: VecScatter *oscatter; /* scatter from global space to the subdomain global space */
13: VecScatter *oscatter_copy; /* copy of the above */
14: VecScatter *iscatter; /* scatter from global space to the nonoverlapping subdomain space */
15: VecScatter *gscatter; /* scatter from global space to the subdomain local space */
16: PCASMType type; /* ASM type */
17: PetscBool usesdm; /* use the DM for setting up the subproblems */
18: PetscBool finaljacobian; /* compute the jacobian of the converged solution */
19: PetscReal damping; /* damping parameter for updates from the blocks */
20: PetscBool weight_set; /* use a weight in the overlap updates */
22: /* logging events */
23: PetscLogEvent eventrestrictinterp;
24: PetscLogEvent eventsubsolve;
26: PetscInt fjtype; /* type of computed jacobian */
27: Vec xinit; /* initial solution in case the final jacobian type is computed as first */
28: } SNES_NASM;
30: const char *const SNESNASMTypes[] = {"NONE", "RESTRICT", "INTERPOLATE", "BASIC", "PCASMType", "PC_ASM_", NULL};
31: const char *const SNESNASMFJTypes[] = {"FINALOUTER", "FINALINNER", "INITIAL"};
33: static PetscErrorCode SNESReset_NASM(SNES snes)
34: {
35: SNES_NASM *nasm = (SNES_NASM *)snes->data;
36: PetscInt i;
38: for (i = 0; i < nasm->n; i++) {
39: if (nasm->xl) VecDestroy(&nasm->xl[i]);
40: if (nasm->x) VecDestroy(&nasm->x[i]);
41: if (nasm->y) VecDestroy(&nasm->y[i]);
42: if (nasm->b) VecDestroy(&nasm->b[i]);
44: if (nasm->subsnes) SNESDestroy(&nasm->subsnes[i]);
45: if (nasm->oscatter) VecScatterDestroy(&nasm->oscatter[i]);
46: if (nasm->oscatter_copy) VecScatterDestroy(&nasm->oscatter_copy[i]);
47: if (nasm->iscatter) VecScatterDestroy(&nasm->iscatter[i]);
48: if (nasm->gscatter) VecScatterDestroy(&nasm->gscatter[i]);
49: }
51: PetscFree(nasm->x);
52: PetscFree(nasm->xl);
53: PetscFree(nasm->y);
54: PetscFree(nasm->b);
56: if (nasm->xinit) VecDestroy(&nasm->xinit);
58: PetscFree(nasm->subsnes);
59: PetscFree(nasm->oscatter);
60: PetscFree(nasm->oscatter_copy);
61: PetscFree(nasm->iscatter);
62: PetscFree(nasm->gscatter);
64: if (nasm->weight_set) VecDestroy(&nasm->weight);
66: nasm->eventrestrictinterp = 0;
67: nasm->eventsubsolve = 0;
68: return 0;
69: }
71: static PetscErrorCode SNESDestroy_NASM(SNES snes)
72: {
73: SNESReset_NASM(snes);
74: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetType_C", NULL);
75: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetType_C", NULL);
76: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetSubdomains_C", NULL);
77: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetSubdomains_C", NULL);
78: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetDamping_C", NULL);
79: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetDamping_C", NULL);
80: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetSubdomainVecs_C", NULL);
81: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetComputeFinalJacobian_C", NULL);
82: PetscFree(snes->data);
83: return 0;
84: }
86: static PetscErrorCode DMGlobalToLocalSubDomainDirichletHook_Private(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
87: {
88: Vec bcs = (Vec)ctx;
90: VecCopy(bcs, l);
91: return 0;
92: }
94: static PetscErrorCode SNESSetUp_NASM(SNES snes)
95: {
96: SNES_NASM *nasm = (SNES_NASM *)snes->data;
97: DM dm, subdm;
98: DM *subdms;
99: PetscInt i;
100: const char *optionsprefix;
101: Vec F;
102: PetscMPIInt size;
103: KSP ksp;
104: PC pc;
106: if (!nasm->subsnes) {
107: SNESGetDM(snes, &dm);
108: if (dm) {
109: nasm->usesdm = PETSC_TRUE;
110: DMCreateDomainDecomposition(dm, &nasm->n, NULL, NULL, NULL, &subdms);
112: DMCreateDomainDecompositionScatters(dm, nasm->n, subdms, &nasm->iscatter, &nasm->oscatter, &nasm->gscatter);
113: PetscMalloc1(nasm->n, &nasm->oscatter_copy);
114: for (i = 0; i < nasm->n; i++) VecScatterCopy(nasm->oscatter[i], &nasm->oscatter_copy[i]);
116: SNESGetOptionsPrefix(snes, &optionsprefix);
117: PetscMalloc1(nasm->n, &nasm->subsnes);
118: for (i = 0; i < nasm->n; i++) {
119: SNESCreate(PETSC_COMM_SELF, &nasm->subsnes[i]);
120: PetscObjectIncrementTabLevel((PetscObject)nasm->subsnes[i], (PetscObject)snes, 1);
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) PetscCalloc1(nasm->n, &nasm->x);
139: if (!nasm->xl) PetscCalloc1(nasm->n, &nasm->xl);
140: if (!nasm->y) PetscCalloc1(nasm->n, &nasm->y);
141: if (!nasm->b) PetscCalloc1(nasm->n, &nasm->b);
143: for (i = 0; i < nasm->n; i++) {
144: SNESGetFunction(nasm->subsnes[i], &F, NULL, NULL);
145: if (!nasm->x[i]) VecDuplicate(F, &nasm->x[i]);
146: if (!nasm->y[i]) VecDuplicate(F, &nasm->y[i]);
147: if (!nasm->b[i]) VecDuplicate(F, &nasm->b[i]);
148: if (!nasm->xl[i]) {
149: SNESGetDM(nasm->subsnes[i], &subdm);
150: DMCreateLocalVector(subdm, &nasm->xl[i]);
151: DMGlobalToLocalHookAdd(subdm, DMGlobalToLocalSubDomainDirichletHook_Private, NULL, nasm->xl[i]);
152: }
153: }
154: if (nasm->finaljacobian) {
155: SNESSetUpMatrices(snes);
156: if (nasm->fjtype == 2) VecDuplicate(snes->vec_sol, &nasm->xinit);
157: for (i = 0; i < nasm->n; i++) SNESSetUpMatrices(nasm->subsnes[i]);
158: }
159: return 0;
160: }
162: static PetscErrorCode SNESSetFromOptions_NASM(SNES snes, PetscOptionItems *PetscOptionsObject)
163: {
164: PCASMType asmtype;
165: PetscBool flg, monflg;
166: SNES_NASM *nasm = (SNES_NASM *)snes->data;
168: PetscOptionsHeadBegin(PetscOptionsObject, "Nonlinear Additive Schwarz options");
169: PetscOptionsEnum("-snes_nasm_type", "Type of restriction/extension", "", SNESNASMTypes, (PetscEnum)nasm->type, (PetscEnum *)&asmtype, &flg);
170: if (flg) SNESNASMSetType(snes, asmtype);
171: flg = PETSC_FALSE;
172: monflg = PETSC_TRUE;
173: 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);
174: if (flg) SNESNASMSetDamping(snes, nasm->damping);
175: PetscOptionsDeprecated("-snes_nasm_sub_view", NULL, "3.15", "Use -snes_view ::ascii_info_detail");
176: PetscOptionsBool("-snes_nasm_finaljacobian", "Compute the global jacobian of the final iterate (for ASPIN)", "", nasm->finaljacobian, &nasm->finaljacobian, NULL);
177: PetscOptionsEList("-snes_nasm_finaljacobian_type", "The type of the final jacobian computed.", "", SNESNASMFJTypes, 3, SNESNASMFJTypes[0], &nasm->fjtype, NULL);
178: PetscOptionsBool("-snes_nasm_log", "Log times for subSNES solves and restriction", "", monflg, &monflg, &flg);
179: if (flg) {
180: PetscLogEventRegister("SNESNASMSubSolve", ((PetscObject)snes)->classid, &nasm->eventsubsolve);
181: PetscLogEventRegister("SNESNASMRestrict", ((PetscObject)snes)->classid, &nasm->eventrestrictinterp);
182: }
183: PetscOptionsHeadEnd();
184: return 0;
185: }
187: static PetscErrorCode SNESView_NASM(SNES snes, PetscViewer viewer)
188: {
189: SNES_NASM *nasm = (SNES_NASM *)snes->data;
190: PetscMPIInt rank, size;
191: PetscInt i, N, bsz;
192: PetscBool iascii, isstring;
193: PetscViewer sviewer;
194: MPI_Comm comm;
195: PetscViewerFormat format;
196: const char *prefix;
198: PetscObjectGetComm((PetscObject)snes, &comm);
199: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
200: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);
201: MPI_Comm_rank(comm, &rank);
202: MPI_Comm_size(comm, &size);
203: MPIU_Allreduce(&nasm->n, &N, 1, MPIU_INT, MPI_SUM, comm);
204: if (iascii) {
205: PetscViewerASCIIPrintf(viewer, " total subdomain blocks = %" PetscInt_FMT "\n", N);
206: PetscViewerGetFormat(viewer, &format);
207: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
208: if (nasm->subsnes) {
209: PetscViewerASCIIPrintf(viewer, " Local solver information for first block on rank 0:\n");
210: SNESGetOptionsPrefix(snes, &prefix);
211: PetscViewerASCIIPrintf(viewer, " Use -%ssnes_view ::ascii_info_detail to display information for all blocks\n", prefix ? prefix : "");
212: PetscViewerASCIIPushTab(viewer);
213: PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
214: if (rank == 0) {
215: PetscViewerASCIIPushTab(viewer);
216: SNESView(nasm->subsnes[0], sviewer);
217: PetscViewerASCIIPopTab(viewer);
218: }
219: PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
220: PetscViewerASCIIPopTab(viewer);
221: }
222: } else {
223: /* print the solver on each block */
224: PetscViewerASCIIPushSynchronized(viewer);
225: PetscViewerASCIISynchronizedPrintf(viewer, " [%d] number of local blocks = %" PetscInt_FMT "\n", (int)rank, nasm->n);
226: PetscViewerFlush(viewer);
227: PetscViewerASCIIPopSynchronized(viewer);
228: PetscViewerASCIIPrintf(viewer, " Local solver information for each block is in the following SNES objects:\n");
229: PetscViewerASCIIPushTab(viewer);
230: PetscViewerASCIIPrintf(viewer, "- - - - - - - - - - - - - - - - - -\n");
231: PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
232: for (i = 0; i < nasm->n; i++) {
233: VecGetLocalSize(nasm->x[i], &bsz);
234: PetscViewerASCIIPrintf(sviewer, "[%d] local block number %" PetscInt_FMT ", size = %" PetscInt_FMT "\n", (int)rank, i, bsz);
235: SNESView(nasm->subsnes[i], sviewer);
236: PetscViewerASCIIPrintf(sviewer, "- - - - - - - - - - - - - - - - - -\n");
237: }
238: PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
239: PetscViewerFlush(viewer);
240: PetscViewerASCIIPopTab(viewer);
241: }
242: } else if (isstring) {
243: PetscViewerStringSPrintf(viewer, " blocks=%" PetscInt_FMT ",type=%s", N, SNESNASMTypes[nasm->type]);
244: PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
245: if (nasm->subsnes && rank == 0) SNESView(nasm->subsnes[0], sviewer);
246: PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
247: }
248: return 0;
249: }
251: /*@
252: SNESNASMSetType - Set the type of subdomain update used for the nonlinear additive Schwarz solver
254: Logically Collective
256: Input Parameters:
257: + snes - the `SNES` context
258: - type - the type of update, `PC_ASM_BASIC` or `PC_ASM_RESTRICT`
260: Level: intermediate
262: .seealso: `SNESNASM`, `SNESNASMGetType()`, `PCASMSetType()`, `PC_ASM_BASIC`, `PC_ASM_RESTRICT`, `PCASMType`
263: @*/
264: PetscErrorCode SNESNASMSetType(SNES snes, PCASMType type)
265: {
266: PetscErrorCode (*f)(SNES, PCASMType);
268: PetscObjectQueryFunction((PetscObject)snes, "SNESNASMSetType_C", &f);
269: if (f) (f)(snes, type);
270: return 0;
271: }
273: static PetscErrorCode SNESNASMSetType_NASM(SNES snes, PCASMType type)
274: {
275: SNES_NASM *nasm = (SNES_NASM *)snes->data;
278: nasm->type = type;
279: return 0;
280: }
282: /*@
283: SNESNASMGetType - Get the type of subdomain update used for the nonlinear additive Schwarz solver
285: Logically Collective
287: Input Parameters:
288: . snes - the `SNES` context
290: Output Parameters:
291: . type - the type of update
293: Level: intermediate
295: .seealso: `SNESNASM`, `SNESNASMSetType()`, `PCASMGetType()`, `PC_ASM_BASIC`, `PC_ASM_RESTRICT`, `PCASMType`
296: @*/
297: PetscErrorCode SNESNASMGetType(SNES snes, PCASMType *type)
298: {
299: PetscUseMethod(snes, "SNESNASMGetType_C", (SNES, PCASMType *), (snes, type));
300: return 0;
301: }
303: static PetscErrorCode SNESNASMGetType_NASM(SNES snes, PCASMType *type)
304: {
305: SNES_NASM *nasm = (SNES_NASM *)snes->data;
307: *type = nasm->type;
308: return 0;
309: }
311: /*@
312: SNESNASMSetSubdomains - Manually Set the context required to restrict and solve subdomain problems in the nonlinear additive Schwarz solver
314: Logically Collective
316: Input Parameters:
317: + snes - the SNES context
318: . n - the number of local subdomains
319: . subsnes - solvers defined on the local subdomains
320: . iscatter - scatters into the nonoverlapping portions of the local subdomains
321: . oscatter - scatters into the overlapping portions of the local subdomains
322: - gscatter - scatters into the (ghosted) local vector of the local subdomain
324: Level: intermediate
326: .seealso: `SNESNASM`, `SNESNASMGetSubdomains()`
327: @*/
328: PetscErrorCode SNESNASMSetSubdomains(SNES snes, PetscInt n, SNES subsnes[], VecScatter iscatter[], VecScatter oscatter[], VecScatter gscatter[])
329: {
330: PetscErrorCode (*f)(SNES, PetscInt, SNES *, VecScatter *, VecScatter *, VecScatter *);
332: PetscObjectQueryFunction((PetscObject)snes, "SNESNASMSetSubdomains_C", &f);
333: if (f) (f)(snes, n, subsnes, iscatter, oscatter, gscatter);
334: return 0;
335: }
337: static PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes, PetscInt n, SNES subsnes[], VecScatter iscatter[], VecScatter oscatter[], VecScatter gscatter[])
338: {
339: PetscInt i;
340: SNES_NASM *nasm = (SNES_NASM *)snes->data;
344: /* tear down the previously set things */
345: SNESReset(snes);
347: nasm->n = n;
348: if (oscatter) {
349: for (i = 0; i < n; i++) PetscObjectReference((PetscObject)oscatter[i]);
350: }
351: if (iscatter) {
352: for (i = 0; i < n; i++) PetscObjectReference((PetscObject)iscatter[i]);
353: }
354: if (gscatter) {
355: for (i = 0; i < n; i++) PetscObjectReference((PetscObject)gscatter[i]);
356: }
357: if (oscatter) {
358: PetscMalloc1(n, &nasm->oscatter);
359: PetscMalloc1(n, &nasm->oscatter_copy);
360: for (i = 0; i < n; i++) {
361: nasm->oscatter[i] = oscatter[i];
362: VecScatterCopy(oscatter[i], &nasm->oscatter_copy[i]);
363: }
364: }
365: if (iscatter) {
366: PetscMalloc1(n, &nasm->iscatter);
367: for (i = 0; i < n; i++) nasm->iscatter[i] = iscatter[i];
368: }
369: if (gscatter) {
370: PetscMalloc1(n, &nasm->gscatter);
371: for (i = 0; i < n; i++) nasm->gscatter[i] = gscatter[i];
372: }
374: if (subsnes) {
375: PetscMalloc1(n, &nasm->subsnes);
376: for (i = 0; i < n; i++) nasm->subsnes[i] = subsnes[i];
377: }
378: return 0;
379: }
381: /*@
382: SNESNASMGetSubdomains - Get the local subdomain contexts for the nonlinear additive Schwarz solver
384: Not Collective but some of the objects returned will be parallel
386: Input Parameter:
387: . snes - the `SNES` context
389: Output Parameters:
390: + n - the number of local subdomains
391: . subsnes - solvers defined on the local subdomains
392: . iscatter - scatters into the nonoverlapping portions of the local subdomains
393: . oscatter - scatters into the overlapping portions of the local subdomains
394: - gscatter - scatters into the (ghosted) local vector of the local subdomain
396: Level: intermediate
398: .seealso: `SNESNASM`, `SNESNASMSetSubdomains()`
399: @*/
400: PetscErrorCode SNESNASMGetSubdomains(SNES snes, PetscInt *n, SNES *subsnes[], VecScatter *iscatter[], VecScatter *oscatter[], VecScatter *gscatter[])
401: {
402: PetscErrorCode (*f)(SNES, PetscInt *, SNES **, VecScatter **, VecScatter **, VecScatter **);
404: PetscObjectQueryFunction((PetscObject)snes, "SNESNASMGetSubdomains_C", &f);
405: if (f) (f)(snes, n, subsnes, iscatter, oscatter, gscatter);
406: return 0;
407: }
409: static PetscErrorCode SNESNASMGetSubdomains_NASM(SNES snes, PetscInt *n, SNES *subsnes[], VecScatter *iscatter[], VecScatter *oscatter[], VecScatter *gscatter[])
410: {
411: SNES_NASM *nasm = (SNES_NASM *)snes->data;
413: if (n) *n = nasm->n;
414: if (oscatter) *oscatter = nasm->oscatter;
415: if (iscatter) *iscatter = nasm->iscatter;
416: if (gscatter) *gscatter = nasm->gscatter;
417: if (subsnes) *subsnes = nasm->subsnes;
418: return 0;
419: }
421: /*@
422: SNESNASMGetSubdomainVecs - Get the processor-local subdomain vectors for the nonlinear additive Schwarz solver
424: Not Collective
426: Input Parameter:
427: . snes - the `SNES` context
429: Output Parameters:
430: + n - the number of local subdomains
431: . x - The subdomain solution vector
432: . y - The subdomain step vector
433: . b - The subdomain RHS vector
434: - xl - The subdomain local vectors (ghosted)
436: Level: developer
438: .seealso: `SNESNASM`, `SNESNASMGetSubdomains()`
439: @*/
440: PetscErrorCode SNESNASMGetSubdomainVecs(SNES snes, PetscInt *n, Vec **x, Vec **y, Vec **b, Vec **xl)
441: {
442: PetscErrorCode (*f)(SNES, PetscInt *, Vec **, Vec **, Vec **, Vec **);
444: PetscObjectQueryFunction((PetscObject)snes, "SNESNASMGetSubdomainVecs_C", &f);
445: if (f) (f)(snes, n, x, y, b, xl);
446: return 0;
447: }
449: static PetscErrorCode SNESNASMGetSubdomainVecs_NASM(SNES snes, PetscInt *n, Vec **x, Vec **y, Vec **b, Vec **xl)
450: {
451: SNES_NASM *nasm = (SNES_NASM *)snes->data;
453: if (n) *n = nasm->n;
454: if (x) *x = nasm->x;
455: if (y) *y = nasm->y;
456: if (b) *b = nasm->b;
457: if (xl) *xl = nasm->xl;
458: return 0;
459: }
461: /*@
462: SNESNASMSetComputeFinalJacobian - Schedules the computation of the global and subdomain Jacobians upon convergence for the
463: nonlinear additive Schwarz solver
465: Collective
467: Input Parameters:
468: + snes - the SNES context
469: - flg - `PETSC_TRUE` to compute the Jacobians
471: Level: developer
473: Notes:
474: This is used almost exclusively in the implementation of `SNESASPIN`, where the converged subdomain and global Jacobian
475: is needed at each linear iteration.
477: .seealso: `SNESNASM`, `SNESNASMGetSubdomains()`
478: @*/
479: PetscErrorCode SNESNASMSetComputeFinalJacobian(SNES snes, PetscBool flg)
480: {
481: PetscErrorCode (*f)(SNES, PetscBool);
483: PetscObjectQueryFunction((PetscObject)snes, "SNESNASMSetComputeFinalJacobian_C", &f);
484: if (f) (f)(snes, flg);
485: return 0;
486: }
488: static PetscErrorCode SNESNASMSetComputeFinalJacobian_NASM(SNES snes, PetscBool flg)
489: {
490: SNES_NASM *nasm = (SNES_NASM *)snes->data;
492: nasm->finaljacobian = flg;
493: return 0;
494: }
496: /*@
497: SNESNASMSetDamping - Sets the update damping for `SNESNASM` the nonlinear additive Schwarz solver
499: Logically collective
501: Input Parameters:
502: + snes - the `SNES` context
503: - dmp - damping
505: Level: intermediate
507: Notes:
508: The new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
510: .seealso: `SNESNASM`, `SNESNASMGetDamping()`
511: @*/
512: PetscErrorCode SNESNASMSetDamping(SNES snes, PetscReal dmp)
513: {
514: PetscErrorCode (*f)(SNES, PetscReal);
516: PetscObjectQueryFunction((PetscObject)snes, "SNESNASMSetDamping_C", (void (**)(void)) & f);
517: if (f) (f)(snes, dmp);
518: return 0;
519: }
521: static PetscErrorCode SNESNASMSetDamping_NASM(SNES snes, PetscReal dmp)
522: {
523: SNES_NASM *nasm = (SNES_NASM *)snes->data;
525: nasm->damping = dmp;
526: return 0;
527: }
529: /*@
530: SNESNASMGetDamping - Gets the update damping for `SNESNASM` the nonlinear additive Schwarz solver
532: Not Collective
534: Input Parameter:
535: . snes - the SNES context
537: Output Parameter:
538: . dmp - damping
540: Level: intermediate
542: .seealso: `SNESNASM`, `SNESNASMSetDamping()`
543: @*/
544: PetscErrorCode SNESNASMGetDamping(SNES snes, PetscReal *dmp)
545: {
546: PetscUseMethod(snes, "SNESNASMGetDamping_C", (SNES, PetscReal *), (snes, dmp));
547: return 0;
548: }
550: static PetscErrorCode SNESNASMGetDamping_NASM(SNES snes, PetscReal *dmp)
551: {
552: SNES_NASM *nasm = (SNES_NASM *)snes->data;
554: *dmp = nasm->damping;
555: return 0;
556: }
558: /*
559: Input Parameters:
560: + snes - The solver
561: . B - The RHS vector
562: - X - The initial guess
564: Output Parameters:
565: . Y - The solution update
567: TODO: All scatters should be packed into one
568: */
569: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes, Vec B, Vec Y, Vec X)
570: {
571: SNES_NASM *nasm = (SNES_NASM *)snes->data;
572: SNES subsnes;
573: PetscInt i;
574: PetscReal dmp;
575: Vec Xl, Bl, Yl, Xlloc;
576: VecScatter iscat, oscat, gscat, oscat_copy;
577: DM dm, subdm;
578: PCASMType type;
580: SNESNASMGetType(snes, &type);
581: SNESGetDM(snes, &dm);
582: VecSet(Y, 0);
583: if (nasm->eventrestrictinterp) PetscLogEventBegin(nasm->eventrestrictinterp, snes, 0, 0, 0);
584: for (i = 0; i < nasm->n; i++) {
585: /* scatter the solution to the global solution and the local solution */
586: Xl = nasm->x[i];
587: Xlloc = nasm->xl[i];
588: oscat = nasm->oscatter[i];
589: oscat_copy = nasm->oscatter_copy[i];
590: gscat = nasm->gscatter[i];
591: VecScatterBegin(oscat, X, Xl, INSERT_VALUES, SCATTER_FORWARD);
592: VecScatterBegin(gscat, X, Xlloc, INSERT_VALUES, SCATTER_FORWARD);
593: if (B) {
594: /* scatter the RHS to the local RHS */
595: Bl = nasm->b[i];
596: VecScatterBegin(oscat_copy, B, Bl, INSERT_VALUES, SCATTER_FORWARD);
597: }
598: }
599: if (nasm->eventrestrictinterp) PetscLogEventEnd(nasm->eventrestrictinterp, snes, 0, 0, 0);
601: if (nasm->eventsubsolve) PetscLogEventBegin(nasm->eventsubsolve, snes, 0, 0, 0);
602: for (i = 0; i < nasm->n; i++) {
603: Xl = nasm->x[i];
604: Xlloc = nasm->xl[i];
605: Yl = nasm->y[i];
606: subsnes = nasm->subsnes[i];
607: SNESGetDM(subsnes, &subdm);
608: iscat = nasm->iscatter[i];
609: oscat = nasm->oscatter[i];
610: oscat_copy = nasm->oscatter_copy[i];
611: gscat = nasm->gscatter[i];
612: VecScatterEnd(oscat, X, Xl, INSERT_VALUES, SCATTER_FORWARD);
613: VecScatterEnd(gscat, X, Xlloc, INSERT_VALUES, SCATTER_FORWARD);
614: if (B) {
615: Bl = nasm->b[i];
616: VecScatterEnd(oscat_copy, B, Bl, INSERT_VALUES, SCATTER_FORWARD);
617: } else Bl = NULL;
619: DMSubDomainRestrict(dm, oscat, gscat, subdm);
620: VecCopy(Xl, Yl);
621: SNESSolve(subsnes, Bl, Xl);
622: VecAYPX(Yl, -1.0, Xl);
623: VecScale(Yl, nasm->damping);
624: if (type == PC_ASM_BASIC) {
625: VecScatterBegin(oscat, Yl, Y, ADD_VALUES, SCATTER_REVERSE);
626: VecScatterEnd(oscat, Yl, Y, ADD_VALUES, SCATTER_REVERSE);
627: } else if (type == PC_ASM_RESTRICT) {
628: VecScatterBegin(iscat, Yl, Y, ADD_VALUES, SCATTER_REVERSE);
629: VecScatterEnd(iscat, Yl, Y, ADD_VALUES, SCATTER_REVERSE);
630: } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Only basic and restrict types are supported for SNESNASM");
631: }
632: if (nasm->eventsubsolve) PetscLogEventEnd(nasm->eventsubsolve, snes, 0, 0, 0);
633: if (nasm->eventrestrictinterp) PetscLogEventBegin(nasm->eventrestrictinterp, snes, 0, 0, 0);
634: if (nasm->weight_set) VecPointwiseMult(Y, Y, nasm->weight);
635: if (nasm->eventrestrictinterp) PetscLogEventEnd(nasm->eventrestrictinterp, snes, 0, 0, 0);
636: SNESNASMGetDamping(snes, &dmp);
637: VecAXPY(X, dmp, Y);
638: return 0;
639: }
641: static PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
642: {
643: Vec X = Xfinal;
644: SNES_NASM *nasm = (SNES_NASM *)snes->data;
645: SNES subsnes;
646: PetscInt i, lag = 1;
647: Vec Xlloc, Xl, Fl, F;
648: VecScatter oscat, gscat;
649: DM dm, subdm;
651: if (nasm->fjtype == 2) X = nasm->xinit;
652: F = snes->vec_func;
653: if (snes->normschedule == SNES_NORM_NONE) SNESComputeFunction(snes, X, F);
654: SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre);
655: SNESGetDM(snes, &dm);
656: if (nasm->eventrestrictinterp) PetscLogEventBegin(nasm->eventrestrictinterp, snes, 0, 0, 0);
657: if (nasm->fjtype != 1) {
658: for (i = 0; i < nasm->n; i++) {
659: Xlloc = nasm->xl[i];
660: gscat = nasm->gscatter[i];
661: VecScatterBegin(gscat, X, Xlloc, INSERT_VALUES, SCATTER_FORWARD);
662: }
663: }
664: if (nasm->eventrestrictinterp) PetscLogEventEnd(nasm->eventrestrictinterp, snes, 0, 0, 0);
665: for (i = 0; i < nasm->n; i++) {
666: Fl = nasm->subsnes[i]->vec_func;
667: Xl = nasm->x[i];
668: Xlloc = nasm->xl[i];
669: subsnes = nasm->subsnes[i];
670: oscat = nasm->oscatter[i];
671: gscat = nasm->gscatter[i];
672: if (nasm->fjtype != 1) VecScatterEnd(gscat, X, Xlloc, INSERT_VALUES, SCATTER_FORWARD);
673: SNESGetDM(subsnes, &subdm);
674: DMSubDomainRestrict(dm, oscat, gscat, subdm);
675: if (nasm->fjtype != 1) {
676: DMLocalToGlobalBegin(subdm, Xlloc, INSERT_VALUES, Xl);
677: DMLocalToGlobalEnd(subdm, Xlloc, INSERT_VALUES, Xl);
678: }
679: if (subsnes->lagjacobian == -1) subsnes->lagjacobian = -2;
680: else if (subsnes->lagjacobian > 1) lag = subsnes->lagjacobian;
681: SNESComputeFunction(subsnes, Xl, Fl);
682: SNESComputeJacobian(subsnes, Xl, subsnes->jacobian, subsnes->jacobian_pre);
683: if (lag > 1) subsnes->lagjacobian = lag;
684: }
685: return 0;
686: }
688: static PetscErrorCode SNESSolve_NASM(SNES snes)
689: {
690: Vec F;
691: Vec X;
692: Vec B;
693: Vec Y;
694: PetscInt i;
695: PetscReal fnorm = 0.0;
696: SNESNormSchedule normschedule;
697: SNES_NASM *nasm = (SNES_NASM *)snes->data;
702: PetscCitationsRegister(SNESCitation, &SNEScite);
703: X = snes->vec_sol;
704: Y = snes->vec_sol_update;
705: F = snes->vec_func;
706: B = snes->vec_rhs;
708: PetscObjectSAWsTakeAccess((PetscObject)snes);
709: snes->iter = 0;
710: snes->norm = 0.;
711: PetscObjectSAWsGrantAccess((PetscObject)snes);
712: snes->reason = SNES_CONVERGED_ITERATING;
713: SNESGetNormSchedule(snes, &normschedule);
714: if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
715: /* compute the initial function and preconditioned update delX */
716: if (!snes->vec_func_init_set) {
717: SNESComputeFunction(snes, X, F);
718: } else snes->vec_func_init_set = PETSC_FALSE;
720: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
721: SNESCheckFunctionNorm(snes, fnorm);
722: PetscObjectSAWsTakeAccess((PetscObject)snes);
723: snes->iter = 0;
724: snes->norm = fnorm;
725: PetscObjectSAWsGrantAccess((PetscObject)snes);
726: SNESLogConvergenceHistory(snes, snes->norm, 0);
727: SNESMonitor(snes, 0, snes->norm);
729: /* test convergence */
730: PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
731: if (snes->reason) return 0;
732: } else {
733: PetscObjectSAWsGrantAccess((PetscObject)snes);
734: SNESLogConvergenceHistory(snes, snes->norm, 0);
735: SNESMonitor(snes, 0, snes->norm);
736: }
738: /* Call general purpose update function */
739: PetscTryTypeMethod(snes, update, snes->iter);
740: /* copy the initial solution over for later */
741: if (nasm->fjtype == 2) VecCopy(X, nasm->xinit);
743: for (i = 0; i < snes->max_its; i++) {
744: SNESNASMSolveLocal_Private(snes, B, Y, X);
745: if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
746: SNESComputeFunction(snes, X, F);
747: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
748: SNESCheckFunctionNorm(snes, fnorm);
749: }
750: /* Monitor convergence */
751: PetscObjectSAWsTakeAccess((PetscObject)snes);
752: snes->iter = i + 1;
753: snes->norm = fnorm;
754: PetscObjectSAWsGrantAccess((PetscObject)snes);
755: SNESLogConvergenceHistory(snes, snes->norm, 0);
756: SNESMonitor(snes, snes->iter, snes->norm);
757: /* Test for convergence */
758: if (normschedule == SNES_NORM_ALWAYS) PetscUseTypeMethod(snes, converged, snes->iter, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
759: if (snes->reason) break;
760: /* Call general purpose update function */
761: PetscTryTypeMethod(snes, update, snes->iter);
762: }
763: if (nasm->finaljacobian) {
764: SNESNASMComputeFinalJacobian_Private(snes, X);
765: SNESCheckJacobianDomainerror(snes);
766: }
767: if (normschedule == SNES_NORM_ALWAYS) {
768: if (i == snes->max_its) {
769: PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", snes->max_its);
770: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
771: }
772: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a nonlinear preconditioner */
773: return 0;
774: }
776: /*MC
777: SNESNASM - Nonlinear Additive Schwarz solver
779: Options Database Keys:
780: + -snes_nasm_log - enable logging events for the communication and solve stages
781: . -snes_nasm_type <basic,restrict> - type of subdomain update used
782: . -snes_asm_damping <dmp> - the new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
783: . -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
784: . -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> - pick state the jacobian is calculated at
785: . -sub_snes_ - options prefix of the subdomain nonlinear solves
786: . -sub_ksp_ - options prefix of the subdomain Krylov solver
787: - -sub_pc_ - options prefix of the subdomain preconditioner
789: Level: advanced
791: Note:
792: This is not often used directly as a solver, it converges too slowly. However it works well as a nonlinear preconditioner for
793: the `SNESASPIN` solver
795: Developer Note:
796: This is a non-Newton based nonlinear solver that does not directly require a Jacobian; hence the flag snes->usesksp is set to
797: false and `SNESView()` and -snes_view do not display a `KSP` object. However, if the flag nasm->finaljacobian is set (for example, if
798: `SNESNASM` is used as a nonlinear preconditioner for `SNESASPIN`) then `SNESSetUpMatrices()` is called to generate the
799: Jacobian (needed by `SNESASPIN`)
800: and this utilizes the inner `KSP` object for storing the matrices, but the `KSP` is never used for solving a linear system. When `SNESNASM` is
801: used by `SNESASPIN` they share the same Jacobian matrices because `SNESSetUp()` (called on the outer `SNESASPIN`) causes the inner `SNES`
802: object (in this case `SNESNASM`) to inherit the outer Jacobian matrices.
804: References:
805: . * - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
806: SIAM Review, 57(4), 2015
808: .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESType`, `SNESNASMSetType()`, `SNESNASMGetType()`, `SNESNASMSetSubdomains()`, `SNESNASMGetSubdomains()`, `SNESNASMGetSubdomainVecs()`, `SNESNASMSetComputeFinalJacobian()`, `SNESNASMSetDamping()`, `SNESNASMGetDamping()`
809: M*/
811: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
812: {
813: SNES_NASM *nasm;
815: PetscNew(&nasm);
816: snes->data = (void *)nasm;
818: nasm->n = PETSC_DECIDE;
819: nasm->subsnes = NULL;
820: nasm->x = NULL;
821: nasm->xl = NULL;
822: nasm->y = NULL;
823: nasm->b = NULL;
824: nasm->oscatter = NULL;
825: nasm->oscatter_copy = NULL;
826: nasm->iscatter = NULL;
827: nasm->gscatter = NULL;
828: nasm->damping = 1.;
830: nasm->type = PC_ASM_BASIC;
831: nasm->finaljacobian = PETSC_FALSE;
832: nasm->weight_set = PETSC_FALSE;
834: snes->ops->destroy = SNESDestroy_NASM;
835: snes->ops->setup = SNESSetUp_NASM;
836: snes->ops->setfromoptions = SNESSetFromOptions_NASM;
837: snes->ops->view = SNESView_NASM;
838: snes->ops->solve = SNESSolve_NASM;
839: snes->ops->reset = SNESReset_NASM;
841: snes->usesksp = PETSC_FALSE;
842: snes->usesnpc = PETSC_FALSE;
844: snes->alwayscomputesfinalresidual = PETSC_FALSE;
846: nasm->fjtype = 0;
847: nasm->xinit = NULL;
848: nasm->eventrestrictinterp = 0;
849: nasm->eventsubsolve = 0;
851: if (!snes->tolerancesset) {
852: snes->max_its = 10000;
853: snes->max_funcs = 10000;
854: }
856: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetType_C", SNESNASMSetType_NASM);
857: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetType_C", SNESNASMGetType_NASM);
858: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetSubdomains_C", SNESNASMSetSubdomains_NASM);
859: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetSubdomains_C", SNESNASMGetSubdomains_NASM);
860: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetDamping_C", SNESNASMSetDamping_NASM);
861: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetDamping_C", SNESNASMGetDamping_NASM);
862: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetSubdomainVecs_C", SNESNASMGetSubdomainVecs_NASM);
863: PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetComputeFinalJacobian_C", SNESNASMSetComputeFinalJacobian_NASM);
864: return 0;
865: }
867: /*@
868: SNESNASMGetSNES - Gets a subsolver
870: Not collective
872: Input Parameters:
873: + snes - the SNES context
874: - i - the number of the subsnes to get
876: Output Parameters:
877: . subsnes - the subsolver context
879: Level: intermediate
881: .seealso: `SNESNASM`, `SNESNASMGetNumber()`
882: @*/
883: PetscErrorCode SNESNASMGetSNES(SNES snes, PetscInt i, SNES *subsnes)
884: {
885: SNES_NASM *nasm = (SNES_NASM *)snes->data;
888: *subsnes = nasm->subsnes[i];
889: return 0;
890: }
892: /*@
893: SNESNASMGetNumber - Gets number of subsolvers
895: Not collective
897: Input Parameters:
898: . snes - the SNES context
900: Output Parameters:
901: . n - the number of subsolvers
903: Level: intermediate
905: .seealso: `SNESNASM`, `SNESNASMGetSNES()`
906: @*/
907: PetscErrorCode SNESNASMGetNumber(SNES snes, PetscInt *n)
908: {
909: SNES_NASM *nasm = (SNES_NASM *)snes->data;
911: *n = nasm->n;
912: return 0;
913: }
915: /*@
916: SNESNASMSetWeight - Sets weight to use when adding overlapping updates
918: Collective
920: Input Parameters:
921: + snes - the SNES context
922: - weight - the weights to use (typically 1/N for each dof, where N is the number of patches it appears in)
924: Level: intermediate
926: .seealso: `SNESNASM`
927: @*/
928: PetscErrorCode SNESNASMSetWeight(SNES snes, Vec weight)
929: {
930: SNES_NASM *nasm = (SNES_NASM *)snes->data;
933: VecDestroy(&nasm->weight);
934: nasm->weight_set = PETSC_TRUE;
935: nasm->weight = weight;
936: PetscObjectReference((PetscObject)nasm->weight);
938: return 0;
939: }