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 */
42: PetscReal rtol; /* restart tolerance for accepting the combination */
43: PetscReal stol; /* restart tolerance for the combination */
44: } SNES_Composite;
46: static PetscErrorCode SNESCompositeApply_Multiplicative(SNES snes, Vec X, Vec B, Vec F, PetscReal *fnorm)
47: {
48: SNES_Composite *jac = (SNES_Composite *)snes->data;
49: SNES_CompositeLink next = jac->head;
50: Vec FSub;
51: SNESConvergedReason reason;
53: PetscFunctionBegin;
54: PetscCheck(next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
55: if (snes->normschedule == SNES_NORM_ALWAYS) PetscCall(SNESSetInitialFunction(next->snes, F));
56: PetscCall(SNESSolve(next->snes, B, X));
57: PetscCall(SNESGetConvergedReason(next->snes, &reason));
58: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
59: jac->innerFailures++;
60: if (jac->innerFailures >= snes->maxFailures) {
61: snes->reason = SNES_DIVERGED_INNER;
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
64: }
66: while (next->next) {
67: /* only copy the function over in the case where the functions correspond */
68: if (next->snes->npcside == PC_RIGHT && next->snes->normschedule != SNES_NORM_NONE) {
69: PetscCall(SNESGetFunction(next->snes, &FSub, NULL, NULL));
70: next = next->next;
71: PetscCall(SNESSetInitialFunction(next->snes, FSub));
72: } else {
73: next = next->next;
74: }
75: PetscCall(SNESSolve(next->snes, B, X));
76: PetscCall(SNESGetConvergedReason(next->snes, &reason));
77: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
78: jac->innerFailures++;
79: if (jac->innerFailures >= snes->maxFailures) {
80: snes->reason = SNES_DIVERGED_INNER;
81: PetscFunctionReturn(PETSC_SUCCESS);
82: }
83: }
84: }
85: if (next->snes->npcside == PC_RIGHT) {
86: PetscCall(SNESGetFunction(next->snes, &FSub, NULL, NULL));
87: PetscCall(VecCopy(FSub, F));
88: if (fnorm) {
89: if (snes->xl && snes->xu) {
90: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
91: } else {
92: PetscCall(VecNorm(F, NORM_2, fnorm));
93: }
94: SNESCheckFunctionDomainError(snes, *fnorm);
95: }
96: } else if (snes->normschedule == SNES_NORM_ALWAYS) {
97: PetscCall(SNESComputeFunction(snes, X, F));
98: if (fnorm) {
99: if (snes->xl && snes->xu) {
100: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
101: } else {
102: PetscCall(VecNorm(F, NORM_2, fnorm));
103: }
104: SNESCheckFunctionDomainError(snes, *fnorm);
105: }
106: }
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: static PetscErrorCode SNESCompositeApply_Additive(SNES snes, Vec X, Vec B, Vec F, PetscReal *fnorm)
111: {
112: SNES_Composite *jac = (SNES_Composite *)snes->data;
113: SNES_CompositeLink next = jac->head;
114: Vec Y, Xorig;
115: SNESConvergedReason reason;
117: PetscFunctionBegin;
118: Y = snes->vec_sol_update;
119: if (!jac->Xorig) PetscCall(VecDuplicate(X, &jac->Xorig));
120: Xorig = jac->Xorig;
121: PetscCall(VecCopy(X, Xorig));
122: PetscCheck(next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
123: if (snes->normschedule == SNES_NORM_ALWAYS) {
124: PetscCall(SNESSetInitialFunction(next->snes, F));
125: while (next->next) {
126: next = next->next;
127: PetscCall(SNESSetInitialFunction(next->snes, F));
128: }
129: }
130: next = jac->head;
131: PetscCall(VecCopy(Xorig, Y));
132: PetscCall(SNESSolve(next->snes, B, Y));
133: PetscCall(SNESGetConvergedReason(next->snes, &reason));
134: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
135: jac->innerFailures++;
136: if (jac->innerFailures >= snes->maxFailures) {
137: snes->reason = SNES_DIVERGED_INNER;
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
140: }
141: PetscCall(VecAXPY(Y, -1.0, Xorig));
142: PetscCall(VecAXPY(X, next->dmp, Y));
143: while (next->next) {
144: next = next->next;
145: PetscCall(VecCopy(Xorig, Y));
146: PetscCall(SNESSolve(next->snes, B, Y));
147: PetscCall(SNESGetConvergedReason(next->snes, &reason));
148: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
149: jac->innerFailures++;
150: if (jac->innerFailures >= snes->maxFailures) {
151: snes->reason = SNES_DIVERGED_INNER;
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
154: }
155: PetscCall(VecAXPY(Y, -1.0, Xorig));
156: PetscCall(VecAXPY(X, next->dmp, Y));
157: }
158: if (snes->normschedule == SNES_NORM_ALWAYS) {
159: PetscCall(SNESComputeFunction(snes, X, F));
160: if (fnorm) {
161: if (snes->xl && snes->xu) {
162: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
163: } else {
164: PetscCall(VecNorm(F, NORM_2, fnorm));
165: }
166: SNESCheckFunctionDomainError(snes, *fnorm);
167: }
168: }
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: /* approximately solve the overdetermined system:
174: 2*F(x_i)\cdot F(\x_j)\alpha_i = 0
175: \alpha_i = 1
177: Which minimizes the L2 norm of the linearization of:
178: ||F(\sum_i \alpha_i*x_i)||^2
180: With the constraint that \sum_i\alpha_i = 1
181: Where x_i is the solution from the ith subsolver.
182: */
183: static PetscErrorCode SNESCompositeApply_AdditiveOptimal(SNES snes, Vec X, Vec B, Vec F, PetscReal *fnorm)
184: {
185: SNES_Composite *jac = (SNES_Composite *)snes->data;
186: SNES_CompositeLink next = jac->head;
187: Vec *Xes = jac->Xes, *Fes = jac->Fes;
188: PetscInt i;
189: PetscScalar tot, total, ftf;
190: PetscReal min_fnorm;
191: PetscInt min_i;
192: SNESConvergedReason reason;
194: PetscFunctionBegin;
195: PetscCheck(next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
197: if (snes->normschedule == SNES_NORM_ALWAYS) {
198: next = jac->head;
199: PetscCall(SNESSetInitialFunction(next->snes, F));
200: while (next->next) {
201: next = next->next;
202: PetscCall(SNESSetInitialFunction(next->snes, F));
203: }
204: }
206: next = jac->head;
207: i = 0;
208: PetscCall(VecCopy(X, Xes[i]));
209: PetscCall(SNESSolve(next->snes, B, Xes[i]));
210: PetscCall(SNESGetConvergedReason(next->snes, &reason));
211: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
212: jac->innerFailures++;
213: if (jac->innerFailures >= snes->maxFailures) {
214: snes->reason = SNES_DIVERGED_INNER;
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
217: }
218: while (next->next) {
219: i++;
220: next = next->next;
221: PetscCall(VecCopy(X, Xes[i]));
222: PetscCall(SNESSolve(next->snes, B, Xes[i]));
223: PetscCall(SNESGetConvergedReason(next->snes, &reason));
224: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
225: jac->innerFailures++;
226: if (jac->innerFailures >= snes->maxFailures) {
227: snes->reason = SNES_DIVERGED_INNER;
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
230: }
231: }
233: /* all the solutions are collected; combine optimally */
234: for (i = 0; i < jac->n; i++) {
235: for (PetscInt j = 0; j < i + 1; j++) PetscCall(VecDotBegin(Fes[i], Fes[j], &jac->h[i + j * jac->n]));
236: PetscCall(VecDotBegin(Fes[i], F, &jac->g[i]));
237: }
239: for (i = 0; i < jac->n; i++) {
240: for (PetscInt j = 0; j < i + 1; j++) {
241: PetscCall(VecDotEnd(Fes[i], Fes[j], &jac->h[i + j * jac->n]));
242: if (i == j) jac->fnorms[i] = PetscSqrtReal(PetscRealPart(jac->h[i + j * jac->n]));
243: }
244: PetscCall(VecDotEnd(Fes[i], F, &jac->g[i]));
245: }
247: ftf = (*fnorm) * (*fnorm);
249: for (i = 0; i < jac->n; i++) {
250: for (PetscInt j = i + 1; j < jac->n; j++) jac->h[i + j * jac->n] = jac->h[j + i * jac->n];
251: }
253: for (i = 0; i < jac->n; i++) {
254: for (PetscInt 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;
255: jac->beta[i] = ftf - jac->g[i];
256: }
258: jac->rcond = -1.;
259: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
260: #if defined(PETSC_USE_COMPLEX)
261: PetscCallLAPACKInfo("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, &info));
262: #else
263: PetscCallLAPACKInfo("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, &info));
264: #endif
265: PetscCall(PetscFPTrapPop());
266: tot = 0.;
267: total = 0.;
268: for (i = 0; i < jac->n; i++) {
269: PetscCheck(!snes->errorifnotconverged || !PetscIsInfOrNanScalar(jac->beta[i]), PetscObjectComm((PetscObject)snes), PETSC_ERR_LIB, "SVD generated inconsistent output");
270: PetscCall(PetscInfo(snes, "%" PetscInt_FMT ": %g\n", i, (double)PetscRealPart(jac->beta[i])));
271: tot += jac->beta[i];
272: total += PetscAbsScalar(jac->beta[i]);
273: }
274: PetscCall(VecScale(X, 1. - tot));
275: PetscCall(VecMAXPY(X, jac->n, jac->beta, Xes));
276: PetscCall(SNESComputeFunction(snes, X, F));
278: if (snes->xl && snes->xu) {
279: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
280: } else {
281: PetscCall(VecNorm(F, NORM_2, fnorm));
282: }
283: SNESCheckFunctionDomainError(snes, *fnorm);
285: /* take the minimum-normed candidate if it beats the combination by a factor of rtol or the combination has stagnated */
286: min_fnorm = jac->fnorms[0];
287: min_i = 0;
288: for (i = 0; i < jac->n; i++) {
289: if (jac->fnorms[i] < min_fnorm) {
290: min_fnorm = jac->fnorms[i];
291: min_i = i;
292: }
293: }
295: /* stagnation or divergence restart to the solution of the solver that failed the least */
296: if (PetscRealPart(total) < jac->stol || min_fnorm * jac->rtol < *fnorm) {
297: PetscCall(VecCopy(jac->Xes[min_i], X));
298: PetscCall(VecCopy(jac->Fes[min_i], F));
299: *fnorm = min_fnorm;
300: }
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: static PetscErrorCode SNESSetUp_Composite(SNES snes)
305: {
306: DM dm;
307: SNES_Composite *jac = (SNES_Composite *)snes->data;
308: SNES_CompositeLink next = jac->head;
309: PetscInt n = 0, i;
310: Vec F;
312: PetscFunctionBegin;
313: PetscCall(SNESGetDM(snes, &dm));
315: if (snes->ops->computevariablebounds) {
316: /* SNESVI only ever calls computevariablebounds once, so calling it once here is justified */
317: if (!snes->xl) PetscCall(VecDuplicate(snes->vec_sol, &snes->xl));
318: if (!snes->xu) PetscCall(VecDuplicate(snes->vec_sol, &snes->xu));
319: PetscUseTypeMethod(snes, computevariablebounds, snes->xl, snes->xu);
320: }
322: while (next) {
323: n++;
324: PetscCall(SNESSetDM(next->snes, dm));
325: PetscCall(SNESSetJacobian(next->snes, snes->jacobian, snes->jacobian_pre, NULL, NULL));
326: PetscCall(SNESSetApplicationContext(next->snes, snes->ctx));
327: if (snes->xl && snes->xu) {
328: if (snes->ops->computevariablebounds) {
329: PetscCall(SNESVISetComputeVariableBounds(next->snes, snes->ops->computevariablebounds));
330: } else {
331: PetscCall(SNESVISetVariableBounds(next->snes, snes->xl, snes->xu));
332: }
333: }
335: next = next->next;
336: }
337: jac->nsnes = n;
338: PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
339: if (jac->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
340: PetscCall(VecDuplicateVecs(F, jac->nsnes, &jac->Xes));
341: PetscCall(PetscMalloc1(n, &jac->Fes));
342: PetscCall(PetscMalloc1(n, &jac->fnorms));
343: next = jac->head;
344: i = 0;
345: while (next) {
346: PetscCall(SNESGetFunction(next->snes, &F, NULL, NULL));
347: jac->Fes[i] = F;
348: PetscCall(PetscObjectReference((PetscObject)F));
349: next = next->next;
350: i++;
351: }
352: /* allocate the subspace direct solve area */
353: jac->nrhs = 1;
354: PetscCall(PetscBLASIntCast(jac->nsnes, &jac->lda));
355: PetscCall(PetscBLASIntCast(jac->nsnes, &jac->ldb));
356: PetscCall(PetscBLASIntCast(jac->nsnes, &jac->n));
358: PetscCall(PetscMalloc4(jac->n * jac->n, &jac->h, jac->n, &jac->beta, jac->n, &jac->s, jac->n, &jac->g));
359: jac->lwork = 12 * jac->n;
360: #if defined(PETSC_USE_COMPLEX)
361: PetscCall(PetscMalloc1(jac->lwork, &jac->rwork));
362: #endif
363: PetscCall(PetscMalloc1(jac->lwork, &jac->work));
364: }
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: static PetscErrorCode SNESReset_Composite(SNES snes)
369: {
370: SNES_Composite *jac = (SNES_Composite *)snes->data;
371: SNES_CompositeLink next = jac->head;
373: PetscFunctionBegin;
374: while (next) {
375: PetscCall(SNESReset(next->snes));
376: next = next->next;
377: }
378: PetscCall(VecDestroy(&jac->Xorig));
379: if (jac->Xes) PetscCall(VecDestroyVecs(jac->nsnes, &jac->Xes));
380: if (jac->Fes) PetscCall(VecDestroyVecs(jac->nsnes, &jac->Fes));
381: PetscCall(PetscFree(jac->fnorms));
382: PetscCall(PetscFree4(jac->h, jac->beta, jac->s, jac->g));
383: PetscCall(PetscFree(jac->work));
384: PetscCall(PetscFree(jac->rwork));
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: static PetscErrorCode SNESDestroy_Composite(SNES snes)
389: {
390: SNES_Composite *jac = (SNES_Composite *)snes->data;
391: SNES_CompositeLink next = jac->head, next_tmp;
393: PetscFunctionBegin;
394: PetscCall(SNESReset_Composite(snes));
395: while (next) {
396: PetscCall(SNESDestroy(&next->snes));
397: next_tmp = next;
398: next = next->next;
399: PetscCall(PetscFree(next_tmp));
400: }
401: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetType_C", NULL));
402: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeAddSNES_C", NULL));
403: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeGetSNES_C", NULL));
404: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetDamping_C", NULL));
405: PetscCall(PetscFree(snes->data));
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: static PetscErrorCode SNESSetFromOptions_Composite(SNES snes, PetscOptionItems PetscOptionsObject)
410: {
411: SNES_Composite *jac = (SNES_Composite *)snes->data;
412: PetscInt nmax = 8, i;
413: SNES_CompositeLink next;
414: char *sneses[8];
415: PetscReal dmps[8];
416: PetscBool flg;
418: PetscFunctionBegin;
419: PetscOptionsHeadBegin(PetscOptionsObject, "Composite preconditioner options");
420: PetscCall(PetscOptionsEnum("-snes_composite_type", "Type of composition", "SNESCompositeSetType", SNESCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&jac->type, &flg));
421: if (flg) PetscCall(SNESCompositeSetType(snes, jac->type));
422: PetscCall(PetscOptionsStringArray("-snes_composite_sneses", "List of composite solvers", "SNESCompositeAddSNES", sneses, &nmax, &flg));
423: if (flg) {
424: for (i = 0; i < nmax; i++) {
425: PetscCall(SNESCompositeAddSNES(snes, sneses[i]));
426: PetscCall(PetscFree(sneses[i])); /* deallocate string sneses[i], which is allocated in PetscOptionsStringArray() */
427: }
428: }
429: PetscCall(PetscOptionsRealArray("-snes_composite_damping", "Damping of the additive composite solvers", "SNESCompositeSetDamping", dmps, &nmax, &flg));
430: if (flg) {
431: for (i = 0; i < nmax; i++) PetscCall(SNESCompositeSetDamping(snes, i, dmps[i]));
432: }
433: PetscCall(PetscOptionsReal("-snes_composite_stol", "Step tolerance for restart on the additive composite solvers", "", jac->stol, &jac->stol, NULL));
434: PetscCall(PetscOptionsReal("-snes_composite_rtol", "Residual tolerance for the additive composite solvers", "", jac->rtol, &jac->rtol, NULL));
435: PetscOptionsHeadEnd();
437: next = jac->head;
438: while (next) {
439: PetscCall(SNESSetFromOptions(next->snes));
440: next = next->next;
441: }
442: PetscFunctionReturn(PETSC_SUCCESS);
443: }
445: static PetscErrorCode SNESView_Composite(SNES snes, PetscViewer viewer)
446: {
447: SNES_Composite *jac = (SNES_Composite *)snes->data;
448: SNES_CompositeLink next = jac->head;
449: PetscBool isascii;
451: PetscFunctionBegin;
452: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
453: if (isascii) {
454: PetscCall(PetscViewerASCIIPrintf(viewer, " type - %s\n", SNESCompositeTypes[jac->type]));
455: PetscCall(PetscViewerASCIIPrintf(viewer, " SNESes on composite preconditioner follow\n"));
456: PetscCall(PetscViewerASCIIPrintf(viewer, " ---------------------------------\n"));
457: }
458: if (isascii) PetscCall(PetscViewerASCIIPushTab(viewer));
459: while (next) {
460: PetscCall(SNESView(next->snes, viewer));
461: next = next->next;
462: }
463: if (isascii) {
464: PetscCall(PetscViewerASCIIPopTab(viewer));
465: PetscCall(PetscViewerASCIIPrintf(viewer, " ---------------------------------\n"));
466: }
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
470: static PetscErrorCode SNESCompositeSetType_Composite(SNES snes, SNESCompositeType type)
471: {
472: SNES_Composite *jac = (SNES_Composite *)snes->data;
474: PetscFunctionBegin;
475: jac->type = type;
476: PetscFunctionReturn(PETSC_SUCCESS);
477: }
479: static PetscErrorCode SNESCompositeAddSNES_Composite(SNES snes, SNESType type)
480: {
481: SNES_Composite *jac;
482: SNES_CompositeLink next, ilink;
483: PetscInt cnt = 0;
484: const char *prefix;
485: char newprefix[20];
486: DM dm;
488: PetscFunctionBegin;
489: PetscCall(PetscNew(&ilink));
490: ilink->next = NULL;
491: PetscCall(SNESCreate(PetscObjectComm((PetscObject)snes), &ilink->snes));
492: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->snes, (PetscObject)snes, 1));
493: PetscCall(SNESGetDM(snes, &dm));
494: PetscCall(SNESSetDM(ilink->snes, dm));
495: PetscCall(SNESSetTolerances(ilink->snes, snes->abstol, snes->rtol, snes->stol, 1, snes->max_funcs));
496: PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)ilink->snes));
497: jac = (SNES_Composite *)snes->data;
498: next = jac->head;
499: if (!next) {
500: jac->head = ilink;
501: ilink->previous = NULL;
502: } else {
503: cnt++;
504: while (next->next) {
505: next = next->next;
506: cnt++;
507: }
508: next->next = ilink;
509: ilink->previous = next;
510: }
511: PetscCall(SNESGetOptionsPrefix(snes, &prefix));
512: PetscCall(SNESSetOptionsPrefix(ilink->snes, prefix));
513: PetscCall(PetscSNPrintf(newprefix, sizeof(newprefix), "sub_%" PetscInt_FMT "_", cnt));
514: PetscCall(SNESAppendOptionsPrefix(ilink->snes, newprefix));
515: PetscCall(SNESSetType(ilink->snes, type));
516: PetscCall(SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY));
518: ilink->dmp = 1.0;
519: jac->nsnes++;
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: static PetscErrorCode SNESCompositeGetSNES_Composite(SNES snes, PetscInt n, SNES *subsnes)
524: {
525: SNES_Composite *jac;
526: SNES_CompositeLink next;
528: PetscFunctionBegin;
529: jac = (SNES_Composite *)snes->data;
530: next = jac->head;
531: for (PetscInt i = 0; i < n; i++) {
532: PetscCheck(next->next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Not enough SNESes in composite preconditioner");
533: next = next->next;
534: }
535: *subsnes = next->snes;
536: PetscFunctionReturn(PETSC_SUCCESS);
537: }
539: /*@
540: SNESCompositeSetType - Sets the type of composite preconditioner.
542: Logically Collective
544: Input Parameters:
545: + snes - the preconditioner context
546: - type - `SNES_COMPOSITE_ADDITIVE` (default), `SNES_COMPOSITE_MULTIPLICATIVE`, or `SNES_COMPOSITE_ADDITIVEOPTIMAL`
548: Options Database Key:
549: . -snes_composite_type (multiplicative|additive|additive_optimal) - Sets composite preconditioner type
551: Level: developer
553: .seealso: [](ch_snes), `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCOMPOSITE`, `SNES_COMPOSITE_ADDITIVEOPTIMAL`,
554: `PCCompositeType`
555: @*/
556: PetscErrorCode SNESCompositeSetType(SNES snes, SNESCompositeType type)
557: {
558: PetscFunctionBegin;
561: PetscTryMethod(snes, "SNESCompositeSetType_C", (SNES, SNESCompositeType), (snes, type));
562: PetscFunctionReturn(PETSC_SUCCESS);
563: }
565: /*@
566: SNESCompositeAddSNES - Adds another `SNES` to the `SNESCOMPOSITE`
568: Collective
570: Input Parameters:
571: + snes - the `SNES` context of type `SNESCOMPOSITE`
572: - type - the `SNESType` of the new solver
574: Level: developer
576: .seealso: [](ch_snes), `SNES`, `SNESCOMPOSITE`, `SNESCompositeGetSNES()`
577: @*/
578: PetscErrorCode SNESCompositeAddSNES(SNES snes, SNESType type)
579: {
580: PetscFunctionBegin;
582: PetscTryMethod(snes, "SNESCompositeAddSNES_C", (SNES, SNESType), (snes, type));
583: PetscFunctionReturn(PETSC_SUCCESS);
584: }
586: /*@
587: SNESCompositeGetSNES - Gets one of the `SNES` objects in the `SNES` of `SNESType` `SNESCOMPOSITE`
589: Not Collective
591: Input Parameters:
592: + snes - the `SNES` context
593: - n - the number of the composed `SNES` requested
595: Output Parameter:
596: . subsnes - the `SNES` requested
598: Level: developer
600: .seealso: [](ch_snes), `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetNumber()`
601: @*/
602: PetscErrorCode SNESCompositeGetSNES(SNES snes, PetscInt n, SNES *subsnes)
603: {
604: PetscFunctionBegin;
606: PetscAssertPointer(subsnes, 3);
607: PetscUseMethod(snes, "SNESCompositeGetSNES_C", (SNES, PetscInt, SNES *), (snes, n, subsnes));
608: PetscFunctionReturn(PETSC_SUCCESS);
609: }
611: /*@
612: SNESCompositeGetNumber - Get the number of subsolvers in the `SNESCOMPOSITE`
614: Logically Collective
616: Input Parameter:
617: . snes - the `SNES` context
619: Output Parameter:
620: . n - the number of subsolvers
622: Level: developer
624: .seealso: [](ch_snes), `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`
625: @*/
626: PetscErrorCode SNESCompositeGetNumber(SNES snes, PetscInt *n)
627: {
628: SNES_Composite *jac;
629: SNES_CompositeLink next;
631: PetscFunctionBegin;
632: jac = (SNES_Composite *)snes->data;
633: next = jac->head;
635: *n = 0;
636: while (next) {
637: *n = *n + 1;
638: next = next->next;
639: }
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: static PetscErrorCode SNESCompositeSetDamping_Composite(SNES snes, PetscInt n, PetscReal dmp)
644: {
645: SNES_Composite *jac;
646: SNES_CompositeLink next;
648: PetscFunctionBegin;
649: jac = (SNES_Composite *)snes->data;
650: next = jac->head;
651: for (PetscInt i = 0; i < n; i++) {
652: PetscCheck(next->next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Not enough SNESes in composite preconditioner");
653: next = next->next;
654: }
655: next->dmp = dmp;
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: /*@
660: SNESCompositeSetDamping - Sets the damping of a subsolver when using `SNES_COMPOSITE_ADDITIVE` with a `SNES` of `SNESType` `SNESCOMPOSITE`
662: Not Collective
664: Input Parameters:
665: + snes - the `SNES` context
666: . n - the number of the sub-`SNES` object requested
667: - dmp - the damping
669: Level: intermediate
671: .seealso: [](ch_snes), `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`,
672: `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCompositeSetType()`
673: @*/
674: PetscErrorCode SNESCompositeSetDamping(SNES snes, PetscInt n, PetscReal dmp)
675: {
676: PetscFunctionBegin;
678: PetscUseMethod(snes, "SNESCompositeSetDamping_C", (SNES, PetscInt, PetscReal), (snes, n, dmp));
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: static PetscErrorCode SNESSolve_Composite(SNES snes)
683: {
684: Vec F, X, B, Y;
685: PetscReal fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
686: SNESNormSchedule normtype;
687: SNES_Composite *comp = (SNES_Composite *)snes->data;
689: PetscFunctionBegin;
690: X = snes->vec_sol;
691: F = snes->vec_func;
692: B = snes->vec_rhs;
693: Y = snes->vec_sol_update;
695: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
696: snes->iter = 0;
697: snes->norm = 0.;
698: comp->innerFailures = 0;
699: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
700: snes->reason = SNES_CONVERGED_ITERATING;
701: PetscCall(SNESGetNormSchedule(snes, &normtype));
702: if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
703: if (!snes->vec_func_init_set) PetscCall(SNESComputeFunction(snes, X, F));
704: else snes->vec_func_init_set = PETSC_FALSE;
706: if (snes->xl && snes->xu) {
707: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm));
708: } else {
709: PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F|| */
710: }
711: SNESCheckFunctionDomainError(snes, fnorm);
712: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
713: snes->iter = 0;
714: snes->norm = fnorm;
715: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
716: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
718: /* test convergence */
719: PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
720: PetscCall(SNESMonitor(snes, 0, snes->norm));
721: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
722: } else {
723: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
724: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
725: PetscCall(SNESMonitor(snes, 0, snes->norm));
726: }
728: for (PetscInt i = 0; i < snes->max_its; i++) {
729: /* Call general purpose update function */
730: PetscTryTypeMethod(snes, update, snes->iter);
732: /* Copy the state before modification by application of the composite solver;
733: we will subtract the new state after application */
734: PetscCall(VecCopy(X, Y));
736: if (comp->type == SNES_COMPOSITE_ADDITIVE) {
737: PetscCall(SNESCompositeApply_Additive(snes, X, B, F, &fnorm));
738: } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
739: PetscCall(SNESCompositeApply_Multiplicative(snes, X, B, F, &fnorm));
740: } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
741: PetscCall(SNESCompositeApply_AdditiveOptimal(snes, X, B, F, &fnorm));
742: } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Unsupported SNESComposite type");
743: if (snes->reason < 0) break;
745: /* Compute the solution update for convergence testing */
746: PetscCall(VecAYPX(Y, -1.0, X));
748: if ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY)) {
749: PetscCall(SNESComputeFunction(snes, X, F));
751: if (snes->xl && snes->xu) {
752: PetscCall(VecNormBegin(X, NORM_2, &xnorm));
753: PetscCall(VecNormBegin(Y, NORM_2, &snorm));
754: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm));
755: PetscCall(VecNormEnd(X, NORM_2, &xnorm));
756: PetscCall(VecNormEnd(Y, NORM_2, &snorm));
757: } else {
758: PetscCall(VecNormBegin(F, NORM_2, &fnorm));
759: PetscCall(VecNormBegin(X, NORM_2, &xnorm));
760: PetscCall(VecNormBegin(Y, NORM_2, &snorm));
762: PetscCall(VecNormEnd(F, NORM_2, &fnorm));
763: PetscCall(VecNormEnd(X, NORM_2, &xnorm));
764: PetscCall(VecNormEnd(Y, NORM_2, &snorm));
765: }
766: SNESCheckFunctionDomainError(snes, fnorm);
767: } else if (normtype == SNES_NORM_ALWAYS) {
768: PetscCall(VecNormBegin(X, NORM_2, &xnorm));
769: PetscCall(VecNormBegin(Y, NORM_2, &snorm));
770: PetscCall(VecNormEnd(X, NORM_2, &xnorm));
771: PetscCall(VecNormEnd(Y, NORM_2, &snorm));
772: }
773: /* Monitor convergence */
774: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
775: snes->iter = i + 1;
776: snes->norm = fnorm;
777: snes->xnorm = xnorm;
778: snes->ynorm = snorm;
779: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
780: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
781: /* Test for convergence */
782: PetscCall(SNESConverged(snes, snes->iter, xnorm, snorm, fnorm));
783: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
784: if (snes->reason) break;
785: }
786: PetscFunctionReturn(PETSC_SUCCESS);
787: }
789: /*MC
790: SNESCOMPOSITE - Builds a nonlinear solver/preconditioner by composing together several `SNES` nonlinear solvers {cite}`bruneknepleysmithtu15`
792: Options Database Keys:
793: + -snes_composite_type (multiplicative|additive|additiveoptimal) - Sets composite preconditioner type
794: - -snes_composite_sneses snes0,snes1,... - list of `SNES` to compose
796: Level: intermediate
798: .seealso: [](ch_snes), `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`,
799: `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_ADDITIVEOPTIMAL`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`,
800: `SNESCompositeSetType()`, `SNESCompositeSetDamping()`, `SNESCompositeGetNumber()`, `PCCompositeType`
801: M*/
803: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
804: {
805: SNES_Composite *jac;
807: PetscFunctionBegin;
808: PetscCall(PetscNew(&jac));
810: snes->ops->solve = SNESSolve_Composite;
811: snes->ops->setup = SNESSetUp_Composite;
812: snes->ops->reset = SNESReset_Composite;
813: snes->ops->destroy = SNESDestroy_Composite;
814: snes->ops->setfromoptions = SNESSetFromOptions_Composite;
815: snes->ops->view = SNESView_Composite;
817: snes->usesksp = PETSC_FALSE;
819: snes->alwayscomputesfinalresidual = PETSC_FALSE;
821: PetscCall(SNESParametersInitialize(snes));
823: snes->data = (void *)jac;
824: jac->type = SNES_COMPOSITE_ADDITIVEOPTIMAL;
825: jac->Fes = NULL;
826: jac->Xes = NULL;
827: jac->fnorms = NULL;
828: jac->nsnes = 0;
829: jac->head = NULL;
830: jac->stol = 0.1;
831: jac->rtol = 1.1;
833: jac->h = NULL;
834: jac->s = NULL;
835: jac->beta = NULL;
836: jac->work = NULL;
837: jac->rwork = NULL;
839: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetType_C", SNESCompositeSetType_Composite));
840: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeAddSNES_C", SNESCompositeAddSNES_Composite));
841: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeGetSNES_C", SNESCompositeGetSNES_Composite));
842: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetDamping_C", SNESCompositeSetDamping_Composite));
843: PetscFunctionReturn(PETSC_SUCCESS);
844: }