Actual source code: snescomposite.c
1: /*
2: Defines a SNES that can consist of a collection of SNESes
3: */
4: #include <petsc/private/snesimpl.h>
5: #include <petscblaslapack.h>
7: const char *const SNESCompositeTypes[] = {"ADDITIVE", "MULTIPLICATIVE", "ADDITIVEOPTIMAL", "SNESCompositeType", "SNES_COMPOSITE", NULL};
9: typedef struct _SNES_CompositeLink *SNES_CompositeLink;
10: struct _SNES_CompositeLink {
11: SNES snes;
12: PetscReal dmp;
13: Vec X;
14: SNES_CompositeLink next;
15: SNES_CompositeLink previous;
16: };
18: typedef struct {
19: SNES_CompositeLink head;
20: PetscInt nsnes;
21: SNESCompositeType type;
22: Vec Xorig;
23: PetscInt innerFailures; /* the number of inner failures we've seen */
25: /* context for ADDITIVEOPTIMAL */
26: Vec *Xes, *Fes; /* solution and residual vectors for the subsolvers */
27: PetscReal *fnorms; /* norms of the solutions */
28: PetscScalar *h; /* the matrix formed as q_ij = (rdot_i, rdot_j) */
29: PetscScalar *g; /* the dotproducts of the previous function with the candidate functions */
30: PetscBLASInt n; /* matrix dimension -- nsnes */
31: PetscBLASInt nrhs; /* the number of right-hand sides */
32: PetscBLASInt lda; /* the padded matrix dimension */
33: PetscBLASInt ldb; /* the padded vector dimension */
34: PetscReal *s; /* the singular values */
35: PetscScalar *beta; /* the RHS and combination */
36: PetscReal rcond; /* the exit condition */
37: PetscBLASInt rank; /* the effective rank */
38: PetscScalar *work; /* the work vector */
39: PetscReal *rwork; /* the real work vector used for complex */
40: PetscBLASInt lwork; /* the size of the work vector */
41: PetscBLASInt info; /* the output condition */
43: PetscReal rtol; /* restart tolerance for accepting the combination */
44: PetscReal stol; /* restart tolerance for the combination */
45: } SNES_Composite;
47: static PetscErrorCode SNESCompositeApply_Multiplicative(SNES snes, Vec X, Vec B, Vec F, PetscReal *fnorm)
48: {
49: SNES_Composite *jac = (SNES_Composite *)snes->data;
50: SNES_CompositeLink next = jac->head;
51: Vec FSub;
52: SNESConvergedReason reason;
54: PetscFunctionBegin;
55: PetscCheck(next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
56: if (snes->normschedule == SNES_NORM_ALWAYS) PetscCall(SNESSetInitialFunction(next->snes, F));
57: PetscCall(SNESSolve(next->snes, B, X));
58: PetscCall(SNESGetConvergedReason(next->snes, &reason));
59: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
60: jac->innerFailures++;
61: if (jac->innerFailures >= snes->maxFailures) {
62: snes->reason = SNES_DIVERGED_INNER;
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
65: }
67: while (next->next) {
68: /* only copy the function over in the case where the functions correspond */
69: if (next->snes->npcside == PC_RIGHT && next->snes->normschedule != SNES_NORM_NONE) {
70: PetscCall(SNESGetFunction(next->snes, &FSub, NULL, NULL));
71: next = next->next;
72: PetscCall(SNESSetInitialFunction(next->snes, FSub));
73: } else {
74: next = next->next;
75: }
76: PetscCall(SNESSolve(next->snes, B, X));
77: PetscCall(SNESGetConvergedReason(next->snes, &reason));
78: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
79: jac->innerFailures++;
80: if (jac->innerFailures >= snes->maxFailures) {
81: snes->reason = SNES_DIVERGED_INNER;
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
84: }
85: }
86: if (next->snes->npcside == PC_RIGHT) {
87: PetscCall(SNESGetFunction(next->snes, &FSub, NULL, NULL));
88: PetscCall(VecCopy(FSub, F));
89: if (fnorm) {
90: if (snes->xl && snes->xu) {
91: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
92: } else {
93: PetscCall(VecNorm(F, NORM_2, fnorm));
94: }
95: SNESCheckFunctionDomainError(snes, *fnorm);
96: }
97: } else if (snes->normschedule == SNES_NORM_ALWAYS) {
98: PetscCall(SNESComputeFunction(snes, X, F));
99: if (fnorm) {
100: if (snes->xl && snes->xu) {
101: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
102: } else {
103: PetscCall(VecNorm(F, NORM_2, fnorm));
104: }
105: SNESCheckFunctionDomainError(snes, *fnorm);
106: }
107: }
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: static PetscErrorCode SNESCompositeApply_Additive(SNES snes, Vec X, Vec B, Vec F, PetscReal *fnorm)
112: {
113: SNES_Composite *jac = (SNES_Composite *)snes->data;
114: SNES_CompositeLink next = jac->head;
115: Vec Y, Xorig;
116: SNESConvergedReason reason;
118: PetscFunctionBegin;
119: Y = snes->vec_sol_update;
120: if (!jac->Xorig) PetscCall(VecDuplicate(X, &jac->Xorig));
121: Xorig = jac->Xorig;
122: PetscCall(VecCopy(X, Xorig));
123: PetscCheck(next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
124: if (snes->normschedule == SNES_NORM_ALWAYS) {
125: PetscCall(SNESSetInitialFunction(next->snes, F));
126: while (next->next) {
127: next = next->next;
128: PetscCall(SNESSetInitialFunction(next->snes, F));
129: }
130: }
131: next = jac->head;
132: PetscCall(VecCopy(Xorig, Y));
133: PetscCall(SNESSolve(next->snes, B, Y));
134: PetscCall(SNESGetConvergedReason(next->snes, &reason));
135: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
136: jac->innerFailures++;
137: if (jac->innerFailures >= snes->maxFailures) {
138: snes->reason = SNES_DIVERGED_INNER;
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
141: }
142: PetscCall(VecAXPY(Y, -1.0, Xorig));
143: PetscCall(VecAXPY(X, next->dmp, Y));
144: while (next->next) {
145: next = next->next;
146: PetscCall(VecCopy(Xorig, Y));
147: PetscCall(SNESSolve(next->snes, B, Y));
148: PetscCall(SNESGetConvergedReason(next->snes, &reason));
149: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
150: jac->innerFailures++;
151: if (jac->innerFailures >= snes->maxFailures) {
152: snes->reason = SNES_DIVERGED_INNER;
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
155: }
156: PetscCall(VecAXPY(Y, -1.0, Xorig));
157: PetscCall(VecAXPY(X, next->dmp, Y));
158: }
159: if (snes->normschedule == SNES_NORM_ALWAYS) {
160: PetscCall(SNESComputeFunction(snes, X, F));
161: if (fnorm) {
162: if (snes->xl && snes->xu) {
163: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
164: } else {
165: PetscCall(VecNorm(F, NORM_2, fnorm));
166: }
167: SNESCheckFunctionDomainError(snes, *fnorm);
168: }
169: }
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: /* approximately solve the overdetermined system:
175: 2*F(x_i)\cdot F(\x_j)\alpha_i = 0
176: \alpha_i = 1
178: Which minimizes the L2 norm of the linearization of:
179: ||F(\sum_i \alpha_i*x_i)||^2
181: With the constraint that \sum_i\alpha_i = 1
182: Where x_i is the solution from the ith subsolver.
183: */
184: static PetscErrorCode SNESCompositeApply_AdditiveOptimal(SNES snes, Vec X, Vec B, Vec F, PetscReal *fnorm)
185: {
186: SNES_Composite *jac = (SNES_Composite *)snes->data;
187: SNES_CompositeLink next = jac->head;
188: Vec *Xes = jac->Xes, *Fes = jac->Fes;
189: PetscInt i, j;
190: PetscScalar tot, total, ftf;
191: PetscReal min_fnorm;
192: PetscInt min_i;
193: SNESConvergedReason reason;
195: PetscFunctionBegin;
196: PetscCheck(next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
198: if (snes->normschedule == SNES_NORM_ALWAYS) {
199: next = jac->head;
200: PetscCall(SNESSetInitialFunction(next->snes, F));
201: while (next->next) {
202: next = next->next;
203: PetscCall(SNESSetInitialFunction(next->snes, F));
204: }
205: }
207: next = jac->head;
208: i = 0;
209: PetscCall(VecCopy(X, Xes[i]));
210: PetscCall(SNESSolve(next->snes, B, Xes[i]));
211: PetscCall(SNESGetConvergedReason(next->snes, &reason));
212: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
213: jac->innerFailures++;
214: if (jac->innerFailures >= snes->maxFailures) {
215: snes->reason = SNES_DIVERGED_INNER;
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
218: }
219: while (next->next) {
220: i++;
221: next = next->next;
222: PetscCall(VecCopy(X, Xes[i]));
223: PetscCall(SNESSolve(next->snes, B, Xes[i]));
224: PetscCall(SNESGetConvergedReason(next->snes, &reason));
225: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
226: jac->innerFailures++;
227: if (jac->innerFailures >= snes->maxFailures) {
228: snes->reason = SNES_DIVERGED_INNER;
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
231: }
232: }
234: /* all the solutions are collected; combine optimally */
235: for (i = 0; i < jac->n; i++) {
236: for (j = 0; j < i + 1; j++) PetscCall(VecDotBegin(Fes[i], Fes[j], &jac->h[i + j * jac->n]));
237: PetscCall(VecDotBegin(Fes[i], F, &jac->g[i]));
238: }
240: for (i = 0; i < jac->n; i++) {
241: for (j = 0; j < i + 1; j++) {
242: PetscCall(VecDotEnd(Fes[i], Fes[j], &jac->h[i + j * jac->n]));
243: if (i == j) jac->fnorms[i] = PetscSqrtReal(PetscRealPart(jac->h[i + j * jac->n]));
244: }
245: PetscCall(VecDotEnd(Fes[i], F, &jac->g[i]));
246: }
248: ftf = (*fnorm) * (*fnorm);
250: for (i = 0; i < jac->n; i++) {
251: for (j = i + 1; j < jac->n; j++) jac->h[i + j * jac->n] = jac->h[j + i * jac->n];
252: }
254: for (i = 0; i < jac->n; i++) {
255: 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;
256: jac->beta[i] = ftf - jac->g[i];
257: }
259: jac->info = 0;
260: jac->rcond = -1.;
261: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
262: #if defined(PETSC_USE_COMPLEX)
263: 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));
264: #else
265: 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));
266: #endif
267: PetscCall(PetscFPTrapPop());
268: PetscCheck(jac->info >= 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_LIB, "Bad argument to GELSS");
269: PetscCheck(jac->info <= 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_LIB, "SVD failed to converge");
270: tot = 0.;
271: total = 0.;
272: for (i = 0; i < jac->n; i++) {
273: PetscCheck(!snes->errorifnotconverged || !PetscIsInfOrNanScalar(jac->beta[i]), PetscObjectComm((PetscObject)snes), PETSC_ERR_LIB, "SVD generated inconsistent output");
274: PetscCall(PetscInfo(snes, "%" PetscInt_FMT ": %g\n", i, (double)PetscRealPart(jac->beta[i])));
275: tot += jac->beta[i];
276: total += PetscAbsScalar(jac->beta[i]);
277: }
278: PetscCall(VecScale(X, 1. - tot));
279: PetscCall(VecMAXPY(X, jac->n, jac->beta, Xes));
280: PetscCall(SNESComputeFunction(snes, X, F));
282: if (snes->xl && snes->xu) {
283: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
284: } else {
285: PetscCall(VecNorm(F, NORM_2, fnorm));
286: }
287: SNESCheckFunctionDomainError(snes, *fnorm);
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(SNESSetJacobian(next->snes, snes->jacobian, snes->jacobian_pre, NULL, NULL));
330: PetscCall(SNESSetApplicationContext(next->snes, snes->ctx));
331: if (snes->xl && snes->xu) {
332: if (snes->ops->computevariablebounds) {
333: PetscCall(SNESVISetComputeVariableBounds(next->snes, snes->ops->computevariablebounds));
334: } else {
335: PetscCall(SNESVISetVariableBounds(next->snes, snes->xl, snes->xu));
336: }
337: }
339: next = next->next;
340: }
341: jac->nsnes = n;
342: PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
343: if (jac->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
344: PetscCall(VecDuplicateVecs(F, jac->nsnes, &jac->Xes));
345: PetscCall(PetscMalloc1(n, &jac->Fes));
346: PetscCall(PetscMalloc1(n, &jac->fnorms));
347: next = jac->head;
348: i = 0;
349: while (next) {
350: PetscCall(SNESGetFunction(next->snes, &F, NULL, NULL));
351: jac->Fes[i] = F;
352: PetscCall(PetscObjectReference((PetscObject)F));
353: next = next->next;
354: i++;
355: }
356: /* allocate the subspace direct solve area */
357: jac->nrhs = 1;
358: PetscCall(PetscBLASIntCast(jac->nsnes, &jac->lda));
359: PetscCall(PetscBLASIntCast(jac->nsnes, &jac->ldb));
360: PetscCall(PetscBLASIntCast(jac->nsnes, &jac->n));
362: PetscCall(PetscMalloc4(jac->n * jac->n, &jac->h, jac->n, &jac->beta, jac->n, &jac->s, jac->n, &jac->g));
363: jac->lwork = 12 * jac->n;
364: #if defined(PETSC_USE_COMPLEX)
365: PetscCall(PetscMalloc1(jac->lwork, &jac->rwork));
366: #endif
367: PetscCall(PetscMalloc1(jac->lwork, &jac->work));
368: }
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: static PetscErrorCode SNESReset_Composite(SNES snes)
373: {
374: SNES_Composite *jac = (SNES_Composite *)snes->data;
375: SNES_CompositeLink next = jac->head;
377: PetscFunctionBegin;
378: while (next) {
379: PetscCall(SNESReset(next->snes));
380: next = next->next;
381: }
382: PetscCall(VecDestroy(&jac->Xorig));
383: if (jac->Xes) PetscCall(VecDestroyVecs(jac->nsnes, &jac->Xes));
384: if (jac->Fes) PetscCall(VecDestroyVecs(jac->nsnes, &jac->Fes));
385: PetscCall(PetscFree(jac->fnorms));
386: PetscCall(PetscFree4(jac->h, jac->beta, jac->s, jac->g));
387: PetscCall(PetscFree(jac->work));
388: PetscCall(PetscFree(jac->rwork));
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: static PetscErrorCode SNESDestroy_Composite(SNES snes)
393: {
394: SNES_Composite *jac = (SNES_Composite *)snes->data;
395: SNES_CompositeLink next = jac->head, next_tmp;
397: PetscFunctionBegin;
398: PetscCall(SNESReset_Composite(snes));
399: while (next) {
400: PetscCall(SNESDestroy(&next->snes));
401: next_tmp = next;
402: next = next->next;
403: PetscCall(PetscFree(next_tmp));
404: }
405: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetType_C", NULL));
406: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeAddSNES_C", NULL));
407: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeGetSNES_C", NULL));
408: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetDamping_C", NULL));
409: PetscCall(PetscFree(snes->data));
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: static PetscErrorCode SNESSetFromOptions_Composite(SNES snes, PetscOptionItems PetscOptionsObject)
414: {
415: SNES_Composite *jac = (SNES_Composite *)snes->data;
416: PetscInt nmax = 8, i;
417: SNES_CompositeLink next;
418: char *sneses[8];
419: PetscReal dmps[8];
420: PetscBool flg;
422: PetscFunctionBegin;
423: PetscOptionsHeadBegin(PetscOptionsObject, "Composite preconditioner options");
424: PetscCall(PetscOptionsEnum("-snes_composite_type", "Type of composition", "SNESCompositeSetType", SNESCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&jac->type, &flg));
425: if (flg) PetscCall(SNESCompositeSetType(snes, jac->type));
426: PetscCall(PetscOptionsStringArray("-snes_composite_sneses", "List of composite solvers", "SNESCompositeAddSNES", sneses, &nmax, &flg));
427: if (flg) {
428: for (i = 0; i < nmax; i++) {
429: PetscCall(SNESCompositeAddSNES(snes, sneses[i]));
430: PetscCall(PetscFree(sneses[i])); /* deallocate string sneses[i], which is allocated in PetscOptionsStringArray() */
431: }
432: }
433: PetscCall(PetscOptionsRealArray("-snes_composite_damping", "Damping of the additive composite solvers", "SNESCompositeSetDamping", dmps, &nmax, &flg));
434: if (flg) {
435: for (i = 0; i < nmax; i++) PetscCall(SNESCompositeSetDamping(snes, i, dmps[i]));
436: }
437: PetscCall(PetscOptionsReal("-snes_composite_stol", "Step tolerance for restart on the additive composite solvers", "", jac->stol, &jac->stol, NULL));
438: PetscCall(PetscOptionsReal("-snes_composite_rtol", "Residual tolerance for the additive composite solvers", "", jac->rtol, &jac->rtol, NULL));
439: PetscOptionsHeadEnd();
441: next = jac->head;
442: while (next) {
443: PetscCall(SNESSetFromOptions(next->snes));
444: next = next->next;
445: }
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: static PetscErrorCode SNESView_Composite(SNES snes, PetscViewer viewer)
450: {
451: SNES_Composite *jac = (SNES_Composite *)snes->data;
452: SNES_CompositeLink next = jac->head;
453: PetscBool isascii;
455: PetscFunctionBegin;
456: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
457: if (isascii) {
458: PetscCall(PetscViewerASCIIPrintf(viewer, " type - %s\n", SNESCompositeTypes[jac->type]));
459: PetscCall(PetscViewerASCIIPrintf(viewer, " SNESes on composite preconditioner follow\n"));
460: PetscCall(PetscViewerASCIIPrintf(viewer, " ---------------------------------\n"));
461: }
462: if (isascii) PetscCall(PetscViewerASCIIPushTab(viewer));
463: while (next) {
464: PetscCall(SNESView(next->snes, viewer));
465: next = next->next;
466: }
467: if (isascii) {
468: PetscCall(PetscViewerASCIIPopTab(viewer));
469: PetscCall(PetscViewerASCIIPrintf(viewer, " ---------------------------------\n"));
470: }
471: PetscFunctionReturn(PETSC_SUCCESS);
472: }
474: static PetscErrorCode SNESCompositeSetType_Composite(SNES snes, SNESCompositeType type)
475: {
476: SNES_Composite *jac = (SNES_Composite *)snes->data;
478: PetscFunctionBegin;
479: jac->type = type;
480: PetscFunctionReturn(PETSC_SUCCESS);
481: }
483: static PetscErrorCode SNESCompositeAddSNES_Composite(SNES snes, SNESType type)
484: {
485: SNES_Composite *jac;
486: SNES_CompositeLink next, ilink;
487: PetscInt cnt = 0;
488: const char *prefix;
489: char newprefix[20];
490: DM dm;
492: PetscFunctionBegin;
493: PetscCall(PetscNew(&ilink));
494: ilink->next = NULL;
495: PetscCall(SNESCreate(PetscObjectComm((PetscObject)snes), &ilink->snes));
496: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->snes, (PetscObject)snes, 1));
497: PetscCall(SNESGetDM(snes, &dm));
498: PetscCall(SNESSetDM(ilink->snes, dm));
499: PetscCall(SNESSetTolerances(ilink->snes, snes->abstol, snes->rtol, snes->stol, 1, snes->max_funcs));
500: PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)ilink->snes));
501: jac = (SNES_Composite *)snes->data;
502: next = jac->head;
503: if (!next) {
504: jac->head = ilink;
505: ilink->previous = NULL;
506: } else {
507: cnt++;
508: while (next->next) {
509: next = next->next;
510: cnt++;
511: }
512: next->next = ilink;
513: ilink->previous = next;
514: }
515: PetscCall(SNESGetOptionsPrefix(snes, &prefix));
516: PetscCall(SNESSetOptionsPrefix(ilink->snes, prefix));
517: PetscCall(PetscSNPrintf(newprefix, sizeof(newprefix), "sub_%" PetscInt_FMT "_", cnt));
518: PetscCall(SNESAppendOptionsPrefix(ilink->snes, newprefix));
519: PetscCall(SNESSetType(ilink->snes, type));
520: PetscCall(SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY));
522: ilink->dmp = 1.0;
523: jac->nsnes++;
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: static PetscErrorCode SNESCompositeGetSNES_Composite(SNES snes, PetscInt n, SNES *subsnes)
528: {
529: SNES_Composite *jac;
530: SNES_CompositeLink next;
531: PetscInt i;
533: PetscFunctionBegin;
534: jac = (SNES_Composite *)snes->data;
535: next = jac->head;
536: for (i = 0; i < n; i++) {
537: PetscCheck(next->next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Not enough SNESes in composite preconditioner");
538: next = next->next;
539: }
540: *subsnes = next->snes;
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: /*@
545: SNESCompositeSetType - Sets the type of composite preconditioner.
547: Logically Collective
549: Input Parameters:
550: + snes - the preconditioner context
551: - type - `SNES_COMPOSITE_ADDITIVE` (default), `SNES_COMPOSITE_MULTIPLICATIVE`, or `SNES_COMPOSITE_ADDITIVEOPTIMAL`
553: Options Database Key:
554: . -snes_composite_type <type: one of multiplicative, additive, additiveoptimal> - Sets composite preconditioner type
556: Level: developer
558: .seealso: [](ch_snes), `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCOMPOSITE`, `SNES_COMPOSITE_ADDITIVEOPTIMAL`,
559: `PCCompositeType`
560: @*/
561: PetscErrorCode SNESCompositeSetType(SNES snes, SNESCompositeType type)
562: {
563: PetscFunctionBegin;
566: PetscTryMethod(snes, "SNESCompositeSetType_C", (SNES, SNESCompositeType), (snes, type));
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: /*@
571: SNESCompositeAddSNES - Adds another `SNES` to the `SNESCOMPOSITE`
573: Collective
575: Input Parameters:
576: + snes - the `SNES` context of type `SNESCOMPOSITE`
577: - type - the `SNESType` of the new solver
579: Level: developer
581: .seealso: [](ch_snes), `SNES`, `SNESCOMPOSITE`, `SNESCompositeGetSNES()`
582: @*/
583: PetscErrorCode SNESCompositeAddSNES(SNES snes, SNESType type)
584: {
585: PetscFunctionBegin;
587: PetscTryMethod(snes, "SNESCompositeAddSNES_C", (SNES, SNESType), (snes, type));
588: PetscFunctionReturn(PETSC_SUCCESS);
589: }
591: /*@
592: SNESCompositeGetSNES - Gets one of the `SNES` objects in the `SNES` of `SNESType` `SNESCOMPOSITE`
594: Not Collective
596: Input Parameters:
597: + snes - the `SNES` context
598: - n - the number of the composed `SNES` requested
600: Output Parameter:
601: . subsnes - the `SNES` requested
603: Level: developer
605: .seealso: [](ch_snes), `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetNumber()`
606: @*/
607: PetscErrorCode SNESCompositeGetSNES(SNES snes, PetscInt n, SNES *subsnes)
608: {
609: PetscFunctionBegin;
611: PetscAssertPointer(subsnes, 3);
612: PetscUseMethod(snes, "SNESCompositeGetSNES_C", (SNES, PetscInt, SNES *), (snes, n, subsnes));
613: PetscFunctionReturn(PETSC_SUCCESS);
614: }
616: /*@
617: SNESCompositeGetNumber - Get the number of subsolvers in the `SNESCOMPOSITE`
619: Logically Collective
621: Input Parameter:
622: . snes - the `SNES` context
624: Output Parameter:
625: . n - the number of subsolvers
627: Level: developer
629: .seealso: [](ch_snes), `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`
630: @*/
631: PetscErrorCode SNESCompositeGetNumber(SNES snes, PetscInt *n)
632: {
633: SNES_Composite *jac;
634: SNES_CompositeLink next;
636: PetscFunctionBegin;
637: jac = (SNES_Composite *)snes->data;
638: next = jac->head;
640: *n = 0;
641: while (next) {
642: *n = *n + 1;
643: next = next->next;
644: }
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
648: static PetscErrorCode SNESCompositeSetDamping_Composite(SNES snes, PetscInt n, PetscReal dmp)
649: {
650: SNES_Composite *jac;
651: SNES_CompositeLink next;
652: PetscInt i;
654: PetscFunctionBegin;
655: jac = (SNES_Composite *)snes->data;
656: next = jac->head;
657: for (i = 0; i < n; i++) {
658: PetscCheck(next->next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Not enough SNESes in composite preconditioner");
659: next = next->next;
660: }
661: next->dmp = dmp;
662: PetscFunctionReturn(PETSC_SUCCESS);
663: }
665: /*@
666: SNESCompositeSetDamping - Sets the damping of a subsolver when using `SNES_COMPOSITE_ADDITIVE` with a `SNES` of `SNESType` `SNESCOMPOSITE`
668: Not Collective
670: Input Parameters:
671: + snes - the `SNES` context
672: . n - the number of the sub-`SNES` object requested
673: - dmp - the damping
675: Level: intermediate
677: .seealso: [](ch_snes), `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`,
678: `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCompositeSetType()`
679: @*/
680: PetscErrorCode SNESCompositeSetDamping(SNES snes, PetscInt n, PetscReal dmp)
681: {
682: PetscFunctionBegin;
684: PetscUseMethod(snes, "SNESCompositeSetDamping_C", (SNES, PetscInt, PetscReal), (snes, n, dmp));
685: PetscFunctionReturn(PETSC_SUCCESS);
686: }
688: static PetscErrorCode SNESSolve_Composite(SNES snes)
689: {
690: Vec F, X, B, Y;
691: PetscInt i;
692: PetscReal fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
693: SNESNormSchedule normtype;
694: SNES_Composite *comp = (SNES_Composite *)snes->data;
696: PetscFunctionBegin;
697: X = snes->vec_sol;
698: F = snes->vec_func;
699: B = snes->vec_rhs;
700: Y = snes->vec_sol_update;
702: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
703: snes->iter = 0;
704: snes->norm = 0.;
705: comp->innerFailures = 0;
706: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
707: snes->reason = SNES_CONVERGED_ITERATING;
708: PetscCall(SNESGetNormSchedule(snes, &normtype));
709: if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
710: if (!snes->vec_func_init_set) PetscCall(SNESComputeFunction(snes, X, F));
711: else snes->vec_func_init_set = PETSC_FALSE;
713: if (snes->xl && snes->xu) {
714: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm));
715: } else {
716: PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F|| */
717: }
718: SNESCheckFunctionDomainError(snes, fnorm);
719: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
720: snes->iter = 0;
721: snes->norm = fnorm;
722: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
723: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
725: /* test convergence */
726: PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
727: PetscCall(SNESMonitor(snes, 0, snes->norm));
728: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
729: } else {
730: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
731: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
732: PetscCall(SNESMonitor(snes, 0, snes->norm));
733: }
735: for (i = 0; i < snes->max_its; i++) {
736: /* Call general purpose update function */
737: PetscTryTypeMethod(snes, update, snes->iter);
739: /* Copy the state before modification by application of the composite solver;
740: we will subtract the new state after application */
741: PetscCall(VecCopy(X, Y));
743: if (comp->type == SNES_COMPOSITE_ADDITIVE) {
744: PetscCall(SNESCompositeApply_Additive(snes, X, B, F, &fnorm));
745: } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
746: PetscCall(SNESCompositeApply_Multiplicative(snes, X, B, F, &fnorm));
747: } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
748: PetscCall(SNESCompositeApply_AdditiveOptimal(snes, X, B, F, &fnorm));
749: } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Unsupported SNESComposite type");
750: if (snes->reason < 0) break;
752: /* Compute the solution update for convergence testing */
753: PetscCall(VecAYPX(Y, -1.0, X));
755: if ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY)) {
756: PetscCall(SNESComputeFunction(snes, X, F));
758: if (snes->xl && snes->xu) {
759: PetscCall(VecNormBegin(X, NORM_2, &xnorm));
760: PetscCall(VecNormBegin(Y, NORM_2, &snorm));
761: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm));
762: PetscCall(VecNormEnd(X, NORM_2, &xnorm));
763: PetscCall(VecNormEnd(Y, NORM_2, &snorm));
764: } else {
765: PetscCall(VecNormBegin(F, NORM_2, &fnorm));
766: PetscCall(VecNormBegin(X, NORM_2, &xnorm));
767: PetscCall(VecNormBegin(Y, NORM_2, &snorm));
769: PetscCall(VecNormEnd(F, NORM_2, &fnorm));
770: PetscCall(VecNormEnd(X, NORM_2, &xnorm));
771: PetscCall(VecNormEnd(Y, NORM_2, &snorm));
772: }
773: SNESCheckFunctionDomainError(snes, fnorm);
774: } else if (normtype == SNES_NORM_ALWAYS) {
775: PetscCall(VecNormBegin(X, NORM_2, &xnorm));
776: PetscCall(VecNormBegin(Y, NORM_2, &snorm));
777: PetscCall(VecNormEnd(X, NORM_2, &xnorm));
778: PetscCall(VecNormEnd(Y, NORM_2, &snorm));
779: }
780: /* Monitor convergence */
781: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
782: snes->iter = i + 1;
783: snes->norm = fnorm;
784: snes->xnorm = xnorm;
785: snes->ynorm = snorm;
786: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
787: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
788: /* Test for convergence */
789: PetscCall(SNESConverged(snes, snes->iter, xnorm, snorm, fnorm));
790: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
791: if (snes->reason) break;
792: }
793: PetscFunctionReturn(PETSC_SUCCESS);
794: }
796: /*MC
797: SNESCOMPOSITE - Builds a nonlinear solver/preconditioner by composing together several `SNES` nonlinear solvers {cite}`bruneknepleysmithtu15`
799: Options Database Keys:
800: + -snes_composite_type <type: one of multiplicative, additive, additiveoptimal> - Sets composite preconditioner type
801: - -snes_composite_sneses - <snes0,snes1,...> list of `SNES` to compose
803: Level: intermediate
805: .seealso: [](ch_snes), `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`,
806: `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_ADDITIVEOPTIMAL`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`,
807: `SNESCompositeSetType()`, `SNESCompositeSetDamping()`, `SNESCompositeGetNumber()`, `PCCompositeType`
808: M*/
810: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
811: {
812: SNES_Composite *jac;
814: PetscFunctionBegin;
815: PetscCall(PetscNew(&jac));
817: snes->ops->solve = SNESSolve_Composite;
818: snes->ops->setup = SNESSetUp_Composite;
819: snes->ops->reset = SNESReset_Composite;
820: snes->ops->destroy = SNESDestroy_Composite;
821: snes->ops->setfromoptions = SNESSetFromOptions_Composite;
822: snes->ops->view = SNESView_Composite;
824: snes->usesksp = PETSC_FALSE;
826: snes->alwayscomputesfinalresidual = PETSC_FALSE;
828: PetscCall(SNESParametersInitialize(snes));
830: snes->data = (void *)jac;
831: jac->type = SNES_COMPOSITE_ADDITIVEOPTIMAL;
832: jac->Fes = NULL;
833: jac->Xes = NULL;
834: jac->fnorms = NULL;
835: jac->nsnes = 0;
836: jac->head = NULL;
837: jac->stol = 0.1;
838: jac->rtol = 1.1;
840: jac->h = NULL;
841: jac->s = NULL;
842: jac->beta = NULL;
843: jac->work = NULL;
844: jac->rwork = NULL;
846: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetType_C", SNESCompositeSetType_Composite));
847: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeAddSNES_C", SNESCompositeAddSNES_Composite));
848: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeGetSNES_C", SNESCompositeGetSNES_Composite));
849: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetDamping_C", SNESCompositeSetDamping_Composite));
850: PetscFunctionReturn(PETSC_SUCCESS);
851: }