Actual source code: snescomposite.c
2: /*
3: Defines a SNES that can consist of a collection of SNESes
4: */
5: #include <petsc/private/snesimpl.h>
6: #include <petscblaslapack.h>
8: const char *const SNESCompositeTypes[] = {"ADDITIVE", "MULTIPLICATIVE", "ADDITIVEOPTIMAL", "SNESCompositeType", "SNES_COMPOSITE", NULL};
10: typedef struct _SNES_CompositeLink *SNES_CompositeLink;
11: struct _SNES_CompositeLink {
12: SNES snes;
13: PetscReal dmp;
14: Vec X;
15: SNES_CompositeLink next;
16: SNES_CompositeLink previous;
17: };
19: typedef struct {
20: SNES_CompositeLink head;
21: PetscInt nsnes;
22: SNESCompositeType type;
23: Vec Xorig;
24: PetscInt innerFailures; /* the number of inner failures we've seen */
26: /* context for ADDITIVEOPTIMAL */
27: Vec *Xes, *Fes; /* solution and residual vectors for the subsolvers */
28: PetscReal *fnorms; /* norms of the solutions */
29: PetscScalar *h; /* the matrix formed as q_ij = (rdot_i, rdot_j) */
30: PetscScalar *g; /* the dotproducts of the previous function with the candidate functions */
31: PetscBLASInt n; /* matrix dimension -- nsnes */
32: PetscBLASInt nrhs; /* the number of right hand sides */
33: PetscBLASInt lda; /* the padded matrix dimension */
34: PetscBLASInt ldb; /* the padded vector dimension */
35: PetscReal *s; /* the singular values */
36: PetscScalar *beta; /* the RHS and combination */
37: PetscReal rcond; /* the exit condition */
38: PetscBLASInt rank; /* the effective rank */
39: PetscScalar *work; /* the work vector */
40: PetscReal *rwork; /* the real work vector used for complex */
41: PetscBLASInt lwork; /* the size of the work vector */
42: PetscBLASInt info; /* the output condition */
44: PetscReal rtol; /* restart tolerance for accepting the combination */
45: PetscReal stol; /* restart tolerance for the combination */
46: } SNES_Composite;
48: static PetscErrorCode SNESCompositeApply_Multiplicative(SNES snes, Vec X, Vec B, Vec F, PetscReal *fnorm)
49: {
50: SNES_Composite *jac = (SNES_Composite *)snes->data;
51: SNES_CompositeLink next = jac->head;
52: Vec FSub;
53: SNESConvergedReason reason;
55: PetscFunctionBegin;
56: PetscCheck(next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
57: if (snes->normschedule == SNES_NORM_ALWAYS) PetscCall(SNESSetInitialFunction(next->snes, F));
58: PetscCall(SNESSolve(next->snes, B, X));
59: PetscCall(SNESGetConvergedReason(next->snes, &reason));
60: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
61: jac->innerFailures++;
62: if (jac->innerFailures >= snes->maxFailures) {
63: snes->reason = SNES_DIVERGED_INNER;
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
66: }
68: while (next->next) {
69: /* only copy the function over in the case where the functions correspond */
70: if (next->snes->npcside == PC_RIGHT && next->snes->normschedule != SNES_NORM_NONE) {
71: PetscCall(SNESGetFunction(next->snes, &FSub, NULL, NULL));
72: next = next->next;
73: PetscCall(SNESSetInitialFunction(next->snes, FSub));
74: } else {
75: next = next->next;
76: }
77: PetscCall(SNESSolve(next->snes, B, X));
78: PetscCall(SNESGetConvergedReason(next->snes, &reason));
79: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
80: jac->innerFailures++;
81: if (jac->innerFailures >= snes->maxFailures) {
82: snes->reason = SNES_DIVERGED_INNER;
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
85: }
86: }
87: if (next->snes->npcside == PC_RIGHT) {
88: PetscCall(SNESGetFunction(next->snes, &FSub, NULL, NULL));
89: PetscCall(VecCopy(FSub, F));
90: if (fnorm) {
91: if (snes->xl && snes->xu) {
92: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
93: } else {
94: PetscCall(VecNorm(F, NORM_2, fnorm));
95: }
96: SNESCheckFunctionNorm(snes, *fnorm);
97: }
98: } else if (snes->normschedule == SNES_NORM_ALWAYS) {
99: PetscCall(SNESComputeFunction(snes, X, F));
100: if (fnorm) {
101: if (snes->xl && snes->xu) {
102: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
103: } else {
104: PetscCall(VecNorm(F, NORM_2, fnorm));
105: }
106: SNESCheckFunctionNorm(snes, *fnorm);
107: }
108: }
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
112: static PetscErrorCode SNESCompositeApply_Additive(SNES snes, Vec X, Vec B, Vec F, PetscReal *fnorm)
113: {
114: SNES_Composite *jac = (SNES_Composite *)snes->data;
115: SNES_CompositeLink next = jac->head;
116: Vec Y, Xorig;
117: SNESConvergedReason reason;
119: PetscFunctionBegin;
120: Y = snes->vec_sol_update;
121: if (!jac->Xorig) PetscCall(VecDuplicate(X, &jac->Xorig));
122: Xorig = jac->Xorig;
123: PetscCall(VecCopy(X, Xorig));
124: PetscCheck(next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
125: if (snes->normschedule == SNES_NORM_ALWAYS) {
126: PetscCall(SNESSetInitialFunction(next->snes, F));
127: while (next->next) {
128: next = next->next;
129: PetscCall(SNESSetInitialFunction(next->snes, F));
130: }
131: }
132: next = jac->head;
133: PetscCall(VecCopy(Xorig, Y));
134: PetscCall(SNESSolve(next->snes, B, Y));
135: PetscCall(SNESGetConvergedReason(next->snes, &reason));
136: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
137: jac->innerFailures++;
138: if (jac->innerFailures >= snes->maxFailures) {
139: snes->reason = SNES_DIVERGED_INNER;
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
142: }
143: PetscCall(VecAXPY(Y, -1.0, Xorig));
144: PetscCall(VecAXPY(X, next->dmp, Y));
145: while (next->next) {
146: next = next->next;
147: PetscCall(VecCopy(Xorig, Y));
148: PetscCall(SNESSolve(next->snes, B, Y));
149: PetscCall(SNESGetConvergedReason(next->snes, &reason));
150: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
151: jac->innerFailures++;
152: if (jac->innerFailures >= snes->maxFailures) {
153: snes->reason = SNES_DIVERGED_INNER;
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
156: }
157: PetscCall(VecAXPY(Y, -1.0, Xorig));
158: PetscCall(VecAXPY(X, next->dmp, Y));
159: }
160: if (snes->normschedule == SNES_NORM_ALWAYS) {
161: PetscCall(SNESComputeFunction(snes, X, F));
162: if (fnorm) {
163: if (snes->xl && snes->xu) {
164: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
165: } else {
166: PetscCall(VecNorm(F, NORM_2, fnorm));
167: }
168: SNESCheckFunctionNorm(snes, *fnorm);
169: }
170: }
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: /* approximately solve the overdetermined system:
176: 2*F(x_i)\cdot F(\x_j)\alpha_i = 0
177: \alpha_i = 1
179: Which minimizes the L2 norm of the linearization of:
180: ||F(\sum_i \alpha_i*x_i)||^2
182: With the constraint that \sum_i\alpha_i = 1
183: Where x_i is the solution from the ith subsolver.
184: */
185: static PetscErrorCode SNESCompositeApply_AdditiveOptimal(SNES snes, Vec X, Vec B, Vec F, PetscReal *fnorm)
186: {
187: SNES_Composite *jac = (SNES_Composite *)snes->data;
188: SNES_CompositeLink next = jac->head;
189: Vec *Xes = jac->Xes, *Fes = jac->Fes;
190: PetscInt i, j;
191: PetscScalar tot, total, ftf;
192: PetscReal min_fnorm;
193: PetscInt min_i;
194: SNESConvergedReason reason;
196: PetscFunctionBegin;
197: PetscCheck(next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
199: if (snes->normschedule == SNES_NORM_ALWAYS) {
200: next = jac->head;
201: PetscCall(SNESSetInitialFunction(next->snes, F));
202: while (next->next) {
203: next = next->next;
204: PetscCall(SNESSetInitialFunction(next->snes, F));
205: }
206: }
208: next = jac->head;
209: i = 0;
210: PetscCall(VecCopy(X, Xes[i]));
211: PetscCall(SNESSolve(next->snes, B, Xes[i]));
212: PetscCall(SNESGetConvergedReason(next->snes, &reason));
213: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
214: jac->innerFailures++;
215: if (jac->innerFailures >= snes->maxFailures) {
216: snes->reason = SNES_DIVERGED_INNER;
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
219: }
220: while (next->next) {
221: i++;
222: next = next->next;
223: PetscCall(VecCopy(X, Xes[i]));
224: PetscCall(SNESSolve(next->snes, B, Xes[i]));
225: PetscCall(SNESGetConvergedReason(next->snes, &reason));
226: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
227: jac->innerFailures++;
228: if (jac->innerFailures >= snes->maxFailures) {
229: snes->reason = SNES_DIVERGED_INNER;
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
232: }
233: }
235: /* all the solutions are collected; combine optimally */
236: for (i = 0; i < jac->n; i++) {
237: for (j = 0; j < i + 1; j++) PetscCall(VecDotBegin(Fes[i], Fes[j], &jac->h[i + j * jac->n]));
238: PetscCall(VecDotBegin(Fes[i], F, &jac->g[i]));
239: }
241: for (i = 0; i < jac->n; i++) {
242: for (j = 0; j < i + 1; j++) {
243: PetscCall(VecDotEnd(Fes[i], Fes[j], &jac->h[i + j * jac->n]));
244: if (i == j) jac->fnorms[i] = PetscSqrtReal(PetscRealPart(jac->h[i + j * jac->n]));
245: }
246: PetscCall(VecDotEnd(Fes[i], F, &jac->g[i]));
247: }
249: ftf = (*fnorm) * (*fnorm);
251: for (i = 0; i < jac->n; i++) {
252: for (j = i + 1; j < jac->n; j++) jac->h[i + j * jac->n] = jac->h[j + i * jac->n];
253: }
255: for (i = 0; i < jac->n; i++) {
256: for (j = 0; j < jac->n; j++) jac->h[i + j * jac->n] = jac->h[i + j * jac->n] - jac->g[j] - jac->g[i] + ftf;
257: jac->beta[i] = ftf - jac->g[i];
258: }
260: jac->info = 0;
261: jac->rcond = -1.;
262: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
263: #if defined(PETSC_USE_COMPLEX)
264: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&jac->n, &jac->n, &jac->nrhs, jac->h, &jac->lda, jac->beta, &jac->lda, jac->s, &jac->rcond, &jac->rank, jac->work, &jac->lwork, jac->rwork, &jac->info));
265: #else
266: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&jac->n, &jac->n, &jac->nrhs, jac->h, &jac->lda, jac->beta, &jac->lda, jac->s, &jac->rcond, &jac->rank, jac->work, &jac->lwork, &jac->info));
267: #endif
268: PetscCall(PetscFPTrapPop());
269: PetscCheck(jac->info >= 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_LIB, "Bad argument to GELSS");
270: PetscCheck(jac->info <= 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_LIB, "SVD failed to converge");
271: tot = 0.;
272: total = 0.;
273: for (i = 0; i < jac->n; i++) {
274: PetscCheck(!snes->errorifnotconverged || !PetscIsInfOrNanScalar(jac->beta[i]), PetscObjectComm((PetscObject)snes), PETSC_ERR_LIB, "SVD generated inconsistent output");
275: PetscCall(PetscInfo(snes, "%" PetscInt_FMT ": %g\n", i, (double)PetscRealPart(jac->beta[i])));
276: tot += jac->beta[i];
277: total += PetscAbsScalar(jac->beta[i]);
278: }
279: PetscCall(VecScale(X, (1. - tot)));
280: PetscCall(VecMAXPY(X, jac->n, jac->beta, Xes));
281: PetscCall(SNESComputeFunction(snes, X, F));
283: if (snes->xl && snes->xu) {
284: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
285: } else {
286: PetscCall(VecNorm(F, NORM_2, fnorm));
287: }
289: /* take the minimum-normed candidate if it beats the combination by a factor of rtol or the combination has stagnated */
290: min_fnorm = jac->fnorms[0];
291: min_i = 0;
292: for (i = 0; i < jac->n; i++) {
293: if (jac->fnorms[i] < min_fnorm) {
294: min_fnorm = jac->fnorms[i];
295: min_i = i;
296: }
297: }
299: /* stagnation or divergence restart to the solution of the solver that failed the least */
300: if (PetscRealPart(total) < jac->stol || min_fnorm * jac->rtol < *fnorm) {
301: PetscCall(VecCopy(jac->Xes[min_i], X));
302: PetscCall(VecCopy(jac->Fes[min_i], F));
303: *fnorm = min_fnorm;
304: }
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }
308: static PetscErrorCode SNESSetUp_Composite(SNES snes)
309: {
310: DM dm;
311: SNES_Composite *jac = (SNES_Composite *)snes->data;
312: SNES_CompositeLink next = jac->head;
313: PetscInt n = 0, i;
314: Vec F;
316: PetscFunctionBegin;
317: PetscCall(SNESGetDM(snes, &dm));
319: if (snes->ops->computevariablebounds) {
320: /* SNESVI only ever calls computevariablebounds once, so calling it once here is justified */
321: if (!snes->xl) PetscCall(VecDuplicate(snes->vec_sol, &snes->xl));
322: if (!snes->xu) PetscCall(VecDuplicate(snes->vec_sol, &snes->xu));
323: PetscUseTypeMethod(snes, computevariablebounds, snes->xl, snes->xu);
324: }
326: while (next) {
327: n++;
328: PetscCall(SNESSetDM(next->snes, dm));
329: PetscCall(SNESSetApplicationContext(next->snes, snes->user));
330: if (snes->xl && snes->xu) {
331: if (snes->ops->computevariablebounds) {
332: PetscCall(SNESVISetComputeVariableBounds(next->snes, snes->ops->computevariablebounds));
333: } else {
334: PetscCall(SNESVISetVariableBounds(next->snes, snes->xl, snes->xu));
335: }
336: }
338: next = next->next;
339: }
340: jac->nsnes = n;
341: PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
342: if (jac->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
343: PetscCall(VecDuplicateVecs(F, jac->nsnes, &jac->Xes));
344: PetscCall(PetscMalloc1(n, &jac->Fes));
345: PetscCall(PetscMalloc1(n, &jac->fnorms));
346: next = jac->head;
347: i = 0;
348: while (next) {
349: PetscCall(SNESGetFunction(next->snes, &F, NULL, NULL));
350: jac->Fes[i] = F;
351: PetscCall(PetscObjectReference((PetscObject)F));
352: next = next->next;
353: i++;
354: }
355: /* allocate the subspace direct solve area */
356: jac->nrhs = 1;
357: jac->lda = jac->nsnes;
358: jac->ldb = jac->nsnes;
359: jac->n = jac->nsnes;
361: PetscCall(PetscMalloc1(jac->n * jac->n, &jac->h));
362: PetscCall(PetscMalloc1(jac->n, &jac->beta));
363: PetscCall(PetscMalloc1(jac->n, &jac->s));
364: PetscCall(PetscMalloc1(jac->n, &jac->g));
365: jac->lwork = 12 * jac->n;
366: #if defined(PETSC_USE_COMPLEX)
367: PetscCall(PetscMalloc1(jac->lwork, &jac->rwork));
368: #endif
369: PetscCall(PetscMalloc1(jac->lwork, &jac->work));
370: }
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: static PetscErrorCode SNESReset_Composite(SNES snes)
376: {
377: SNES_Composite *jac = (SNES_Composite *)snes->data;
378: SNES_CompositeLink next = jac->head;
380: PetscFunctionBegin;
381: while (next) {
382: PetscCall(SNESReset(next->snes));
383: next = next->next;
384: }
385: PetscCall(VecDestroy(&jac->Xorig));
386: if (jac->Xes) PetscCall(VecDestroyVecs(jac->nsnes, &jac->Xes));
387: if (jac->Fes) PetscCall(VecDestroyVecs(jac->nsnes, &jac->Fes));
388: PetscCall(PetscFree(jac->fnorms));
389: PetscCall(PetscFree(jac->h));
390: PetscCall(PetscFree(jac->s));
391: PetscCall(PetscFree(jac->g));
392: PetscCall(PetscFree(jac->beta));
393: PetscCall(PetscFree(jac->work));
394: PetscCall(PetscFree(jac->rwork));
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
398: static PetscErrorCode SNESDestroy_Composite(SNES snes)
399: {
400: SNES_Composite *jac = (SNES_Composite *)snes->data;
401: SNES_CompositeLink next = jac->head, next_tmp;
403: PetscFunctionBegin;
404: PetscCall(SNESReset_Composite(snes));
405: while (next) {
406: PetscCall(SNESDestroy(&next->snes));
407: next_tmp = next;
408: next = next->next;
409: PetscCall(PetscFree(next_tmp));
410: }
411: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetType_C", NULL));
412: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeAddSNES_C", NULL));
413: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeGetSNES_C", NULL));
414: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetDamping_C", NULL));
415: PetscCall(PetscFree(snes->data));
416: PetscFunctionReturn(PETSC_SUCCESS);
417: }
419: static PetscErrorCode SNESSetFromOptions_Composite(SNES snes, PetscOptionItems *PetscOptionsObject)
420: {
421: SNES_Composite *jac = (SNES_Composite *)snes->data;
422: PetscInt nmax = 8, i;
423: SNES_CompositeLink next;
424: char *sneses[8];
425: PetscReal dmps[8];
426: PetscBool flg;
428: PetscFunctionBegin;
429: PetscOptionsHeadBegin(PetscOptionsObject, "Composite preconditioner options");
430: PetscCall(PetscOptionsEnum("-snes_composite_type", "Type of composition", "SNESCompositeSetType", SNESCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&jac->type, &flg));
431: if (flg) PetscCall(SNESCompositeSetType(snes, jac->type));
432: PetscCall(PetscOptionsStringArray("-snes_composite_sneses", "List of composite solvers", "SNESCompositeAddSNES", sneses, &nmax, &flg));
433: if (flg) {
434: for (i = 0; i < nmax; i++) {
435: PetscCall(SNESCompositeAddSNES(snes, sneses[i]));
436: PetscCall(PetscFree(sneses[i])); /* deallocate string sneses[i], which is allocated in PetscOptionsStringArray() */
437: }
438: }
439: PetscCall(PetscOptionsRealArray("-snes_composite_damping", "Damping of the additive composite solvers", "SNESCompositeSetDamping", dmps, &nmax, &flg));
440: if (flg) {
441: for (i = 0; i < nmax; i++) PetscCall(SNESCompositeSetDamping(snes, i, dmps[i]));
442: }
443: PetscCall(PetscOptionsReal("-snes_composite_stol", "Step tolerance for restart on the additive composite solvers", "", jac->stol, &jac->stol, NULL));
444: PetscCall(PetscOptionsReal("-snes_composite_rtol", "Residual tolerance for the additive composite solvers", "", jac->rtol, &jac->rtol, NULL));
445: PetscOptionsHeadEnd();
447: next = jac->head;
448: while (next) {
449: PetscCall(SNESSetFromOptions(next->snes));
450: next = next->next;
451: }
452: PetscFunctionReturn(PETSC_SUCCESS);
453: }
455: static PetscErrorCode SNESView_Composite(SNES snes, PetscViewer viewer)
456: {
457: SNES_Composite *jac = (SNES_Composite *)snes->data;
458: SNES_CompositeLink next = jac->head;
459: PetscBool iascii;
461: PetscFunctionBegin;
462: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
463: if (iascii) {
464: PetscCall(PetscViewerASCIIPrintf(viewer, " type - %s\n", SNESCompositeTypes[jac->type]));
465: PetscCall(PetscViewerASCIIPrintf(viewer, " SNESes on composite preconditioner follow\n"));
466: PetscCall(PetscViewerASCIIPrintf(viewer, " ---------------------------------\n"));
467: }
468: if (iascii) PetscCall(PetscViewerASCIIPushTab(viewer));
469: while (next) {
470: PetscCall(SNESView(next->snes, viewer));
471: next = next->next;
472: }
473: if (iascii) {
474: PetscCall(PetscViewerASCIIPopTab(viewer));
475: PetscCall(PetscViewerASCIIPrintf(viewer, " ---------------------------------\n"));
476: }
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }
480: static PetscErrorCode SNESCompositeSetType_Composite(SNES snes, SNESCompositeType type)
481: {
482: SNES_Composite *jac = (SNES_Composite *)snes->data;
484: PetscFunctionBegin;
485: jac->type = type;
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
489: static PetscErrorCode SNESCompositeAddSNES_Composite(SNES snes, SNESType type)
490: {
491: SNES_Composite *jac;
492: SNES_CompositeLink next, ilink;
493: PetscInt cnt = 0;
494: const char *prefix;
495: char newprefix[20];
496: DM dm;
498: PetscFunctionBegin;
499: PetscCall(PetscNew(&ilink));
500: ilink->next = NULL;
501: PetscCall(SNESCreate(PetscObjectComm((PetscObject)snes), &ilink->snes));
502: PetscCall(SNESGetDM(snes, &dm));
503: PetscCall(SNESSetDM(ilink->snes, dm));
504: PetscCall(SNESSetTolerances(ilink->snes, snes->abstol, snes->rtol, snes->stol, 1, snes->max_funcs));
505: PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)ilink->snes));
506: jac = (SNES_Composite *)snes->data;
507: next = jac->head;
508: if (!next) {
509: jac->head = ilink;
510: ilink->previous = NULL;
511: } else {
512: cnt++;
513: while (next->next) {
514: next = next->next;
515: cnt++;
516: }
517: next->next = ilink;
518: ilink->previous = next;
519: }
520: PetscCall(SNESGetOptionsPrefix(snes, &prefix));
521: PetscCall(SNESSetOptionsPrefix(ilink->snes, prefix));
522: PetscCall(PetscSNPrintf(newprefix, sizeof(newprefix), "sub_%d_", (int)cnt));
523: PetscCall(SNESAppendOptionsPrefix(ilink->snes, newprefix));
524: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->snes, (PetscObject)snes, 1));
525: PetscCall(SNESSetType(ilink->snes, type));
526: PetscCall(SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY));
528: ilink->dmp = 1.0;
529: jac->nsnes++;
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: static PetscErrorCode SNESCompositeGetSNES_Composite(SNES snes, PetscInt n, SNES *subsnes)
534: {
535: SNES_Composite *jac;
536: SNES_CompositeLink next;
537: PetscInt i;
539: PetscFunctionBegin;
540: jac = (SNES_Composite *)snes->data;
541: next = jac->head;
542: for (i = 0; i < n; i++) {
543: PetscCheck(next->next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Not enough SNESes in composite preconditioner");
544: next = next->next;
545: }
546: *subsnes = next->snes;
547: PetscFunctionReturn(PETSC_SUCCESS);
548: }
550: /*@C
551: SNESCompositeSetType - Sets the type of composite preconditioner.
553: Logically Collective
555: Input Parameters:
556: + snes - the preconditioner context
557: - type - `SNES_COMPOSITE_ADDITIVE` (default), `SNES_COMPOSITE_MULTIPLICATIVE`
559: Options Database Key:
560: . -snes_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
562: Level: Developer
564: .seealso: `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCOMPOSITE`
565: @*/
566: PetscErrorCode SNESCompositeSetType(SNES snes, SNESCompositeType type)
567: {
568: PetscFunctionBegin;
571: PetscTryMethod(snes, "SNESCompositeSetType_C", (SNES, SNESCompositeType), (snes, type));
572: PetscFunctionReturn(PETSC_SUCCESS);
573: }
575: /*@C
576: SNESCompositeAddSNES - Adds another `SNES` to the `SNESCOMPOSITE`
578: Collective
580: Input Parameters:
581: + snes - the snes context of type `SNESCOMPOSITE`
582: - type - the type of the new solver
584: Level: Developer
586: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeGetSNES()`
587: @*/
588: PetscErrorCode SNESCompositeAddSNES(SNES snes, SNESType type)
589: {
590: PetscFunctionBegin;
592: PetscTryMethod(snes, "SNESCompositeAddSNES_C", (SNES, SNESType), (snes, type));
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
596: /*@
597: SNESCompositeGetSNES - Gets one of the `SNES` objects in the composite `SNESCOMPOSITE`
599: Not Collective
601: Input Parameters:
602: + snes - the snes context
603: - n - the number of the composed snes requested
605: Output Parameter:
606: . subsnes - the `SNES` requested
608: Level: Developer
610: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`
611: @*/
612: PetscErrorCode SNESCompositeGetSNES(SNES snes, PetscInt n, SNES *subsnes)
613: {
614: PetscFunctionBegin;
617: PetscUseMethod(snes, "SNESCompositeGetSNES_C", (SNES, PetscInt, SNES *), (snes, n, subsnes));
618: PetscFunctionReturn(PETSC_SUCCESS);
619: }
621: /*@
622: SNESCompositeGetNumber - Get the number of subsolvers in the `SNESCOMPOSITE`
624: Logically Collective
626: Input Parameter:
627: . snes - the snes context
629: Output Parameter:
630: . n - the number of subsolvers
632: Level: Developer
634: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`
635: @*/
636: PetscErrorCode SNESCompositeGetNumber(SNES snes, PetscInt *n)
637: {
638: SNES_Composite *jac;
639: SNES_CompositeLink next;
641: PetscFunctionBegin;
642: jac = (SNES_Composite *)snes->data;
643: next = jac->head;
645: *n = 0;
646: while (next) {
647: *n = *n + 1;
648: next = next->next;
649: }
650: PetscFunctionReturn(PETSC_SUCCESS);
651: }
653: static PetscErrorCode SNESCompositeSetDamping_Composite(SNES snes, PetscInt n, PetscReal dmp)
654: {
655: SNES_Composite *jac;
656: SNES_CompositeLink next;
657: PetscInt i;
659: PetscFunctionBegin;
660: jac = (SNES_Composite *)snes->data;
661: next = jac->head;
662: for (i = 0; i < n; i++) {
663: PetscCheck(next->next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Not enough SNESes in composite preconditioner");
664: next = next->next;
665: }
666: next->dmp = dmp;
667: PetscFunctionReturn(PETSC_SUCCESS);
668: }
670: /*@
671: SNESCompositeSetDamping - Sets the damping of a subsolver when using `SNES_COMPOSITE_ADDITIVE` `SNESCOMPOSITE`
673: Not Collective
675: Input Parameters:
676: + snes - the snes context
677: . n - the number of the sub-`SNES` object requested
678: - dmp - the damping
680: Level: intermediate
682: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`,
683: `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCompositeSetType()`
684: @*/
685: PetscErrorCode SNESCompositeSetDamping(SNES snes, PetscInt n, PetscReal dmp)
686: {
687: PetscFunctionBegin;
689: PetscUseMethod(snes, "SNESCompositeSetDamping_C", (SNES, PetscInt, PetscReal), (snes, n, dmp));
690: PetscFunctionReturn(PETSC_SUCCESS);
691: }
693: static PetscErrorCode SNESSolve_Composite(SNES snes)
694: {
695: Vec F, X, B, Y;
696: PetscInt i;
697: PetscReal fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
698: SNESNormSchedule normtype;
699: SNES_Composite *comp = (SNES_Composite *)snes->data;
701: PetscFunctionBegin;
702: X = snes->vec_sol;
703: F = snes->vec_func;
704: B = snes->vec_rhs;
705: Y = snes->vec_sol_update;
707: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
708: snes->iter = 0;
709: snes->norm = 0.;
710: comp->innerFailures = 0;
711: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
712: snes->reason = SNES_CONVERGED_ITERATING;
713: PetscCall(SNESGetNormSchedule(snes, &normtype));
714: if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
715: if (!snes->vec_func_init_set) {
716: PetscCall(SNESComputeFunction(snes, X, F));
717: } else snes->vec_func_init_set = PETSC_FALSE;
719: if (snes->xl && snes->xu) {
720: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm));
721: } else {
722: PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F|| */
723: }
724: SNESCheckFunctionNorm(snes, fnorm);
725: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
726: snes->iter = 0;
727: snes->norm = fnorm;
728: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
729: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
730: PetscCall(SNESMonitor(snes, 0, snes->norm));
732: /* test convergence */
733: PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
734: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
735: } else {
736: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
737: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
738: PetscCall(SNESMonitor(snes, 0, snes->norm));
739: }
741: for (i = 0; i < snes->max_its; i++) {
742: /* Call general purpose update function */
743: PetscTryTypeMethod(snes, update, snes->iter);
745: /* Copy the state before modification by application of the composite solver;
746: we will subtract the new state after application */
747: PetscCall(VecCopy(X, Y));
749: if (comp->type == SNES_COMPOSITE_ADDITIVE) {
750: PetscCall(SNESCompositeApply_Additive(snes, X, B, F, &fnorm));
751: } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
752: PetscCall(SNESCompositeApply_Multiplicative(snes, X, B, F, &fnorm));
753: } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
754: PetscCall(SNESCompositeApply_AdditiveOptimal(snes, X, B, F, &fnorm));
755: } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Unsupported SNESComposite type");
756: if (snes->reason < 0) break;
758: /* Compute the solution update for convergence testing */
759: PetscCall(VecAYPX(Y, -1.0, X));
761: if ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY)) {
762: PetscCall(SNESComputeFunction(snes, X, F));
764: if (snes->xl && snes->xu) {
765: PetscCall(VecNormBegin(X, NORM_2, &xnorm));
766: PetscCall(VecNormBegin(Y, NORM_2, &snorm));
767: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm));
768: PetscCall(VecNormEnd(X, NORM_2, &xnorm));
769: PetscCall(VecNormEnd(Y, NORM_2, &snorm));
770: } else {
771: PetscCall(VecNormBegin(F, NORM_2, &fnorm));
772: PetscCall(VecNormBegin(X, NORM_2, &xnorm));
773: PetscCall(VecNormBegin(Y, NORM_2, &snorm));
775: PetscCall(VecNormEnd(F, NORM_2, &fnorm));
776: PetscCall(VecNormEnd(X, NORM_2, &xnorm));
777: PetscCall(VecNormEnd(Y, NORM_2, &snorm));
778: }
779: SNESCheckFunctionNorm(snes, fnorm);
780: } else if (normtype == SNES_NORM_ALWAYS) {
781: PetscCall(VecNormBegin(X, NORM_2, &xnorm));
782: PetscCall(VecNormBegin(Y, NORM_2, &snorm));
783: PetscCall(VecNormEnd(X, NORM_2, &xnorm));
784: PetscCall(VecNormEnd(Y, NORM_2, &snorm));
785: }
786: /* Monitor convergence */
787: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
788: snes->iter = i + 1;
789: snes->norm = fnorm;
790: snes->xnorm = xnorm;
791: snes->ynorm = snorm;
792: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
793: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
794: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
795: /* Test for convergence */
796: if (normtype == SNES_NORM_ALWAYS) PetscUseTypeMethod(snes, converged, snes->iter, xnorm, snorm, fnorm, &snes->reason, snes->cnvP);
797: if (snes->reason) break;
798: }
799: if (normtype == SNES_NORM_ALWAYS) {
800: if (i == snes->max_its) {
801: PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", snes->max_its));
802: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
803: }
804: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
805: PetscFunctionReturn(PETSC_SUCCESS);
806: }
808: /*MC
809: SNESCOMPOSITE - Build a preconditioner by composing together several nonlinear solvers
811: Options Database Keys:
812: + -snes_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
813: - -snes_composite_sneses - <snes0,snes1,...> list of `SNES` to compose
815: Level: intermediate
817: References:
818: . * - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
819: SIAM Review, 57(4), 2015
821: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`,
822: `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCompositeSetType()`,
823: `SNESCompositeSetDamping()`
824: M*/
826: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
827: {
828: SNES_Composite *jac;
830: PetscFunctionBegin;
831: PetscCall(PetscNew(&jac));
833: snes->ops->solve = SNESSolve_Composite;
834: snes->ops->setup = SNESSetUp_Composite;
835: snes->ops->reset = SNESReset_Composite;
836: snes->ops->destroy = SNESDestroy_Composite;
837: snes->ops->setfromoptions = SNESSetFromOptions_Composite;
838: snes->ops->view = SNESView_Composite;
840: snes->usesksp = PETSC_FALSE;
842: snes->alwayscomputesfinalresidual = PETSC_FALSE;
844: snes->data = (void *)jac;
845: jac->type = SNES_COMPOSITE_ADDITIVEOPTIMAL;
846: jac->Fes = NULL;
847: jac->Xes = NULL;
848: jac->fnorms = NULL;
849: jac->nsnes = 0;
850: jac->head = NULL;
851: jac->stol = 0.1;
852: jac->rtol = 1.1;
854: jac->h = NULL;
855: jac->s = NULL;
856: jac->beta = NULL;
857: jac->work = NULL;
858: jac->rwork = NULL;
860: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetType_C", SNESCompositeSetType_Composite));
861: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeAddSNES_C", SNESCompositeAddSNES_Composite));
862: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeGetSNES_C", SNESCompositeGetSNES_Composite));
863: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetDamping_C", SNESCompositeSetDamping_Composite));
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }