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;
56: if (snes->normschedule == SNES_NORM_ALWAYS) SNESSetInitialFunction(next->snes, F);
57: SNESSolve(next->snes, B, X);
58: 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: return 0;
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: SNESGetFunction(next->snes, &FSub, NULL, NULL);
71: next = next->next;
72: SNESSetInitialFunction(next->snes, FSub);
73: } else {
74: next = next->next;
75: }
76: SNESSolve(next->snes, B, X);
77: 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: return 0;
83: }
84: }
85: }
86: if (next->snes->npcside == PC_RIGHT) {
87: SNESGetFunction(next->snes, &FSub, NULL, NULL);
88: VecCopy(FSub, F);
89: if (fnorm) {
90: if (snes->xl && snes->xu) {
91: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
92: } else {
93: VecNorm(F, NORM_2, fnorm);
94: }
95: SNESCheckFunctionNorm(snes, *fnorm);
96: }
97: } else if (snes->normschedule == SNES_NORM_ALWAYS) {
98: SNESComputeFunction(snes, X, F);
99: if (fnorm) {
100: if (snes->xl && snes->xu) {
101: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
102: } else {
103: VecNorm(F, NORM_2, fnorm);
104: }
105: SNESCheckFunctionNorm(snes, *fnorm);
106: }
107: }
108: return 0;
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: Y = snes->vec_sol_update;
119: if (!jac->Xorig) VecDuplicate(X, &jac->Xorig);
120: Xorig = jac->Xorig;
121: VecCopy(X, Xorig);
123: if (snes->normschedule == SNES_NORM_ALWAYS) {
124: SNESSetInitialFunction(next->snes, F);
125: while (next->next) {
126: next = next->next;
127: SNESSetInitialFunction(next->snes, F);
128: }
129: }
130: next = jac->head;
131: VecCopy(Xorig, Y);
132: SNESSolve(next->snes, B, Y);
133: 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: return 0;
139: }
140: }
141: VecAXPY(Y, -1.0, Xorig);
142: VecAXPY(X, next->dmp, Y);
143: while (next->next) {
144: next = next->next;
145: VecCopy(Xorig, Y);
146: SNESSolve(next->snes, B, Y);
147: 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: return 0;
153: }
154: }
155: VecAXPY(Y, -1.0, Xorig);
156: VecAXPY(X, next->dmp, Y);
157: }
158: if (snes->normschedule == SNES_NORM_ALWAYS) {
159: SNESComputeFunction(snes, X, F);
160: if (fnorm) {
161: if (snes->xl && snes->xu) {
162: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
163: } else {
164: VecNorm(F, NORM_2, fnorm);
165: }
166: SNESCheckFunctionNorm(snes, *fnorm);
167: }
168: }
169: return 0;
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, j;
189: PetscScalar tot, total, ftf;
190: PetscReal min_fnorm;
191: PetscInt min_i;
192: SNESConvergedReason reason;
196: if (snes->normschedule == SNES_NORM_ALWAYS) {
197: next = jac->head;
198: SNESSetInitialFunction(next->snes, F);
199: while (next->next) {
200: next = next->next;
201: SNESSetInitialFunction(next->snes, F);
202: }
203: }
205: next = jac->head;
206: i = 0;
207: VecCopy(X, Xes[i]);
208: SNESSolve(next->snes, B, Xes[i]);
209: SNESGetConvergedReason(next->snes, &reason);
210: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
211: jac->innerFailures++;
212: if (jac->innerFailures >= snes->maxFailures) {
213: snes->reason = SNES_DIVERGED_INNER;
214: return 0;
215: }
216: }
217: while (next->next) {
218: i++;
219: next = next->next;
220: VecCopy(X, Xes[i]);
221: SNESSolve(next->snes, B, Xes[i]);
222: SNESGetConvergedReason(next->snes, &reason);
223: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
224: jac->innerFailures++;
225: if (jac->innerFailures >= snes->maxFailures) {
226: snes->reason = SNES_DIVERGED_INNER;
227: return 0;
228: }
229: }
230: }
232: /* all the solutions are collected; combine optimally */
233: for (i = 0; i < jac->n; i++) {
234: for (j = 0; j < i + 1; j++) VecDotBegin(Fes[i], Fes[j], &jac->h[i + j * jac->n]);
235: VecDotBegin(Fes[i], F, &jac->g[i]);
236: }
238: for (i = 0; i < jac->n; i++) {
239: for (j = 0; j < i + 1; j++) {
240: VecDotEnd(Fes[i], Fes[j], &jac->h[i + j * jac->n]);
241: if (i == j) jac->fnorms[i] = PetscSqrtReal(PetscRealPart(jac->h[i + j * jac->n]));
242: }
243: VecDotEnd(Fes[i], F, &jac->g[i]);
244: }
246: ftf = (*fnorm) * (*fnorm);
248: for (i = 0; i < jac->n; i++) {
249: for (j = i + 1; j < jac->n; j++) jac->h[i + j * jac->n] = jac->h[j + i * jac->n];
250: }
252: for (i = 0; i < jac->n; i++) {
253: 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;
254: jac->beta[i] = ftf - jac->g[i];
255: }
257: jac->info = 0;
258: jac->rcond = -1.;
259: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
260: #if defined(PETSC_USE_COMPLEX)
261: 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));
262: #else
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->info));
264: #endif
265: PetscFPTrapPop();
268: tot = 0.;
269: total = 0.;
270: for (i = 0; i < jac->n; i++) {
272: PetscInfo(snes, "%" PetscInt_FMT ": %g\n", i, (double)PetscRealPart(jac->beta[i]));
273: tot += jac->beta[i];
274: total += PetscAbsScalar(jac->beta[i]);
275: }
276: VecScale(X, (1. - tot));
277: VecMAXPY(X, jac->n, jac->beta, Xes);
278: SNESComputeFunction(snes, X, F);
280: if (snes->xl && snes->xu) {
281: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
282: } else {
283: VecNorm(F, NORM_2, fnorm);
284: }
286: /* take the minimum-normed candidate if it beats the combination by a factor of rtol or the combination has stagnated */
287: min_fnorm = jac->fnorms[0];
288: min_i = 0;
289: for (i = 0; i < jac->n; i++) {
290: if (jac->fnorms[i] < min_fnorm) {
291: min_fnorm = jac->fnorms[i];
292: min_i = i;
293: }
294: }
296: /* stagnation or divergence restart to the solution of the solver that failed the least */
297: if (PetscRealPart(total) < jac->stol || min_fnorm * jac->rtol < *fnorm) {
298: VecCopy(jac->Xes[min_i], X);
299: VecCopy(jac->Fes[min_i], F);
300: *fnorm = min_fnorm;
301: }
302: return 0;
303: }
305: static PetscErrorCode SNESSetUp_Composite(SNES snes)
306: {
307: DM dm;
308: SNES_Composite *jac = (SNES_Composite *)snes->data;
309: SNES_CompositeLink next = jac->head;
310: PetscInt n = 0, i;
311: Vec F;
313: 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) VecDuplicate(snes->vec_sol, &snes->xl);
318: if (!snes->xu) VecDuplicate(snes->vec_sol, &snes->xu);
319: PetscUseTypeMethod(snes, computevariablebounds, snes->xl, snes->xu);
320: }
322: while (next) {
323: n++;
324: SNESSetDM(next->snes, dm);
325: SNESSetApplicationContext(next->snes, snes->user);
326: if (snes->xl && snes->xu) {
327: if (snes->ops->computevariablebounds) {
328: SNESVISetComputeVariableBounds(next->snes, snes->ops->computevariablebounds);
329: } else {
330: SNESVISetVariableBounds(next->snes, snes->xl, snes->xu);
331: }
332: }
334: next = next->next;
335: }
336: jac->nsnes = n;
337: SNESGetFunction(snes, &F, NULL, NULL);
338: if (jac->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
339: VecDuplicateVecs(F, jac->nsnes, &jac->Xes);
340: PetscMalloc1(n, &jac->Fes);
341: PetscMalloc1(n, &jac->fnorms);
342: next = jac->head;
343: i = 0;
344: while (next) {
345: SNESGetFunction(next->snes, &F, NULL, NULL);
346: jac->Fes[i] = F;
347: PetscObjectReference((PetscObject)F);
348: next = next->next;
349: i++;
350: }
351: /* allocate the subspace direct solve area */
352: jac->nrhs = 1;
353: jac->lda = jac->nsnes;
354: jac->ldb = jac->nsnes;
355: jac->n = jac->nsnes;
357: PetscMalloc1(jac->n * jac->n, &jac->h);
358: PetscMalloc1(jac->n, &jac->beta);
359: PetscMalloc1(jac->n, &jac->s);
360: PetscMalloc1(jac->n, &jac->g);
361: jac->lwork = 12 * jac->n;
362: #if defined(PETSC_USE_COMPLEX)
363: PetscMalloc1(jac->lwork, &jac->rwork);
364: #endif
365: PetscMalloc1(jac->lwork, &jac->work);
366: }
368: return 0;
369: }
371: static PetscErrorCode SNESReset_Composite(SNES snes)
372: {
373: SNES_Composite *jac = (SNES_Composite *)snes->data;
374: SNES_CompositeLink next = jac->head;
376: while (next) {
377: SNESReset(next->snes);
378: next = next->next;
379: }
380: VecDestroy(&jac->Xorig);
381: if (jac->Xes) VecDestroyVecs(jac->nsnes, &jac->Xes);
382: if (jac->Fes) VecDestroyVecs(jac->nsnes, &jac->Fes);
383: PetscFree(jac->fnorms);
384: PetscFree(jac->h);
385: PetscFree(jac->s);
386: PetscFree(jac->g);
387: PetscFree(jac->beta);
388: PetscFree(jac->work);
389: PetscFree(jac->rwork);
390: return 0;
391: }
393: static PetscErrorCode SNESDestroy_Composite(SNES snes)
394: {
395: SNES_Composite *jac = (SNES_Composite *)snes->data;
396: SNES_CompositeLink next = jac->head, next_tmp;
398: SNESReset_Composite(snes);
399: while (next) {
400: SNESDestroy(&next->snes);
401: next_tmp = next;
402: next = next->next;
403: PetscFree(next_tmp);
404: }
405: PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetType_C", NULL);
406: PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeAddSNES_C", NULL);
407: PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeGetSNES_C", NULL);
408: PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetDamping_C", NULL);
409: PetscFree(snes->data);
410: return 0;
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: PetscOptionsHeadBegin(PetscOptionsObject, "Composite preconditioner options");
423: PetscOptionsEnum("-snes_composite_type", "Type of composition", "SNESCompositeSetType", SNESCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&jac->type, &flg);
424: if (flg) SNESCompositeSetType(snes, jac->type);
425: PetscOptionsStringArray("-snes_composite_sneses", "List of composite solvers", "SNESCompositeAddSNES", sneses, &nmax, &flg);
426: if (flg) {
427: for (i = 0; i < nmax; i++) {
428: SNESCompositeAddSNES(snes, sneses[i]);
429: PetscFree(sneses[i]); /* deallocate string sneses[i], which is allocated in PetscOptionsStringArray() */
430: }
431: }
432: PetscOptionsRealArray("-snes_composite_damping", "Damping of the additive composite solvers", "SNESCompositeSetDamping", dmps, &nmax, &flg);
433: if (flg) {
434: for (i = 0; i < nmax; i++) SNESCompositeSetDamping(snes, i, dmps[i]);
435: }
436: PetscOptionsReal("-snes_composite_stol", "Step tolerance for restart on the additive composite solvers", "", jac->stol, &jac->stol, NULL);
437: PetscOptionsReal("-snes_composite_rtol", "Residual tolerance for the additive composite solvers", "", jac->rtol, &jac->rtol, NULL);
438: PetscOptionsHeadEnd();
440: next = jac->head;
441: while (next) {
442: SNESSetFromOptions(next->snes);
443: next = next->next;
444: }
445: return 0;
446: }
448: static PetscErrorCode SNESView_Composite(SNES snes, PetscViewer viewer)
449: {
450: SNES_Composite *jac = (SNES_Composite *)snes->data;
451: SNES_CompositeLink next = jac->head;
452: PetscBool iascii;
454: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
455: if (iascii) {
456: PetscViewerASCIIPrintf(viewer, " type - %s\n", SNESCompositeTypes[jac->type]);
457: PetscViewerASCIIPrintf(viewer, " SNESes on composite preconditioner follow\n");
458: PetscViewerASCIIPrintf(viewer, " ---------------------------------\n");
459: }
460: if (iascii) PetscViewerASCIIPushTab(viewer);
461: while (next) {
462: SNESView(next->snes, viewer);
463: next = next->next;
464: }
465: if (iascii) {
466: PetscViewerASCIIPopTab(viewer);
467: PetscViewerASCIIPrintf(viewer, " ---------------------------------\n");
468: }
469: return 0;
470: }
472: static PetscErrorCode SNESCompositeSetType_Composite(SNES snes, SNESCompositeType type)
473: {
474: SNES_Composite *jac = (SNES_Composite *)snes->data;
476: jac->type = type;
477: return 0;
478: }
480: static PetscErrorCode SNESCompositeAddSNES_Composite(SNES snes, SNESType type)
481: {
482: SNES_Composite *jac;
483: SNES_CompositeLink next, ilink;
484: PetscInt cnt = 0;
485: const char *prefix;
486: char newprefix[20];
487: DM dm;
489: PetscNew(&ilink);
490: ilink->next = NULL;
491: SNESCreate(PetscObjectComm((PetscObject)snes), &ilink->snes);
492: SNESGetDM(snes, &dm);
493: SNESSetDM(ilink->snes, dm);
494: SNESSetTolerances(ilink->snes, snes->abstol, snes->rtol, snes->stol, 1, snes->max_funcs);
495: PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)ilink->snes);
496: jac = (SNES_Composite *)snes->data;
497: next = jac->head;
498: if (!next) {
499: jac->head = ilink;
500: ilink->previous = NULL;
501: } else {
502: cnt++;
503: while (next->next) {
504: next = next->next;
505: cnt++;
506: }
507: next->next = ilink;
508: ilink->previous = next;
509: }
510: SNESGetOptionsPrefix(snes, &prefix);
511: SNESSetOptionsPrefix(ilink->snes, prefix);
512: PetscSNPrintf(newprefix, sizeof(newprefix), "sub_%d_", (int)cnt);
513: SNESAppendOptionsPrefix(ilink->snes, newprefix);
514: PetscObjectIncrementTabLevel((PetscObject)ilink->snes, (PetscObject)snes, 1);
515: SNESSetType(ilink->snes, type);
516: SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY);
518: ilink->dmp = 1.0;
519: jac->nsnes++;
520: return 0;
521: }
523: static PetscErrorCode SNESCompositeGetSNES_Composite(SNES snes, PetscInt n, SNES *subsnes)
524: {
525: SNES_Composite *jac;
526: SNES_CompositeLink next;
527: PetscInt i;
529: jac = (SNES_Composite *)snes->data;
530: next = jac->head;
531: for (i = 0; i < n; i++) {
533: next = next->next;
534: }
535: *subsnes = next->snes;
536: return 0;
537: }
539: /*@C
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`
548: Options Database Key:
549: . -snes_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
551: Level: Developer
553: .seealso: `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCOMPOSITE`
554: @*/
555: PetscErrorCode SNESCompositeSetType(SNES snes, SNESCompositeType type)
556: {
559: PetscTryMethod(snes, "SNESCompositeSetType_C", (SNES, SNESCompositeType), (snes, type));
560: return 0;
561: }
563: /*@C
564: SNESCompositeAddSNES - Adds another `SNES` to the `SNESCOMPOSITE`
566: Collective
568: Input Parameters:
569: + snes - the snes context of type `SNESCOMPOSITE`
570: - type - the type of the new solver
572: Level: Developer
574: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeGetSNES()`
575: @*/
576: PetscErrorCode SNESCompositeAddSNES(SNES snes, SNESType type)
577: {
579: PetscTryMethod(snes, "SNESCompositeAddSNES_C", (SNES, SNESType), (snes, type));
580: return 0;
581: }
583: /*@
584: SNESCompositeGetSNES - Gets one of the `SNES` objects in the composite `SNESCOMPOSITE`
586: Not Collective
588: Input Parameters:
589: + snes - the snes context
590: - n - the number of the composed snes requested
592: Output Parameters:
593: . subsnes - the `SNES` requested
595: Level: Developer
597: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`
598: @*/
599: PetscErrorCode SNESCompositeGetSNES(SNES snes, PetscInt n, SNES *subsnes)
600: {
603: PetscUseMethod(snes, "SNESCompositeGetSNES_C", (SNES, PetscInt, SNES *), (snes, n, subsnes));
604: return 0;
605: }
607: /*@
608: SNESCompositeGetNumber - Get the number of subsolvers in the `SNESCOMPOSITE`
610: Logically Collective
612: Input Parameter:
613: snes - the snes context
615: Output Parameter:
616: n - the number of subsolvers
618: Level: Developer
620: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`
621: @*/
622: PetscErrorCode SNESCompositeGetNumber(SNES snes, PetscInt *n)
623: {
624: SNES_Composite *jac;
625: SNES_CompositeLink next;
627: jac = (SNES_Composite *)snes->data;
628: next = jac->head;
630: *n = 0;
631: while (next) {
632: *n = *n + 1;
633: next = next->next;
634: }
635: return 0;
636: }
638: static PetscErrorCode SNESCompositeSetDamping_Composite(SNES snes, PetscInt n, PetscReal dmp)
639: {
640: SNES_Composite *jac;
641: SNES_CompositeLink next;
642: PetscInt i;
644: jac = (SNES_Composite *)snes->data;
645: next = jac->head;
646: for (i = 0; i < n; i++) {
648: next = next->next;
649: }
650: next->dmp = dmp;
651: return 0;
652: }
654: /*@
655: SNESCompositeSetDamping - Sets the damping of a subsolver when using `SNES_COMPOSITE_ADDITIVE` `SNESCOMPOSITE`
657: Not Collective
659: Input Parameters:
660: + snes - the snes context
661: . n - the number of the snes requested
662: - dmp - the damping
664: Level: intermediate
666: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`,
667: `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCompositeSetType()`
668: @*/
669: PetscErrorCode SNESCompositeSetDamping(SNES snes, PetscInt n, PetscReal dmp)
670: {
672: PetscUseMethod(snes, "SNESCompositeSetDamping_C", (SNES, PetscInt, PetscReal), (snes, n, dmp));
673: return 0;
674: }
676: static PetscErrorCode SNESSolve_Composite(SNES snes)
677: {
678: Vec F, X, B, Y;
679: PetscInt i;
680: PetscReal fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
681: SNESNormSchedule normtype;
682: SNES_Composite *comp = (SNES_Composite *)snes->data;
684: X = snes->vec_sol;
685: F = snes->vec_func;
686: B = snes->vec_rhs;
687: Y = snes->vec_sol_update;
689: PetscObjectSAWsTakeAccess((PetscObject)snes);
690: snes->iter = 0;
691: snes->norm = 0.;
692: comp->innerFailures = 0;
693: PetscObjectSAWsGrantAccess((PetscObject)snes);
694: snes->reason = SNES_CONVERGED_ITERATING;
695: SNESGetNormSchedule(snes, &normtype);
696: if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
697: if (!snes->vec_func_init_set) {
698: SNESComputeFunction(snes, X, F);
699: } else snes->vec_func_init_set = PETSC_FALSE;
701: if (snes->xl && snes->xu) {
702: SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
703: } else {
704: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
705: }
706: SNESCheckFunctionNorm(snes, fnorm);
707: PetscObjectSAWsTakeAccess((PetscObject)snes);
708: snes->iter = 0;
709: snes->norm = fnorm;
710: PetscObjectSAWsGrantAccess((PetscObject)snes);
711: SNESLogConvergenceHistory(snes, snes->norm, 0);
712: SNESMonitor(snes, 0, snes->norm);
714: /* test convergence */
715: PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
716: if (snes->reason) return 0;
717: } else {
718: PetscObjectSAWsGrantAccess((PetscObject)snes);
719: SNESLogConvergenceHistory(snes, snes->norm, 0);
720: SNESMonitor(snes, 0, snes->norm);
721: }
723: for (i = 0; i < snes->max_its; i++) {
724: /* Call general purpose update function */
725: PetscTryTypeMethod(snes, update, snes->iter);
727: /* Copy the state before modification by application of the composite solver;
728: we will subtract the new state after application */
729: VecCopy(X, Y);
731: if (comp->type == SNES_COMPOSITE_ADDITIVE) {
732: SNESCompositeApply_Additive(snes, X, B, F, &fnorm);
733: } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
734: SNESCompositeApply_Multiplicative(snes, X, B, F, &fnorm);
735: } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
736: SNESCompositeApply_AdditiveOptimal(snes, X, B, F, &fnorm);
737: } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Unsupported SNESComposite type");
738: if (snes->reason < 0) break;
740: /* Compute the solution update for convergence testing */
741: VecAYPX(Y, -1.0, X);
743: if ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY)) {
744: SNESComputeFunction(snes, X, F);
746: if (snes->xl && snes->xu) {
747: VecNormBegin(X, NORM_2, &xnorm);
748: VecNormBegin(Y, NORM_2, &snorm);
749: SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
750: VecNormEnd(X, NORM_2, &xnorm);
751: VecNormEnd(Y, NORM_2, &snorm);
752: } else {
753: VecNormBegin(F, NORM_2, &fnorm);
754: VecNormBegin(X, NORM_2, &xnorm);
755: VecNormBegin(Y, NORM_2, &snorm);
757: VecNormEnd(F, NORM_2, &fnorm);
758: VecNormEnd(X, NORM_2, &xnorm);
759: VecNormEnd(Y, NORM_2, &snorm);
760: }
761: SNESCheckFunctionNorm(snes, fnorm);
762: } else if (normtype == SNES_NORM_ALWAYS) {
763: VecNormBegin(X, NORM_2, &xnorm);
764: VecNormBegin(Y, NORM_2, &snorm);
765: VecNormEnd(X, NORM_2, &xnorm);
766: VecNormEnd(Y, NORM_2, &snorm);
767: }
768: /* Monitor convergence */
769: PetscObjectSAWsTakeAccess((PetscObject)snes);
770: snes->iter = i + 1;
771: snes->norm = fnorm;
772: snes->xnorm = xnorm;
773: snes->ynorm = snorm;
774: PetscObjectSAWsGrantAccess((PetscObject)snes);
775: SNESLogConvergenceHistory(snes, snes->norm, 0);
776: SNESMonitor(snes, snes->iter, snes->norm);
777: /* Test for convergence */
778: if (normtype == SNES_NORM_ALWAYS) PetscUseTypeMethod(snes, converged, snes->iter, xnorm, snorm, fnorm, &snes->reason, snes->cnvP);
779: if (snes->reason) break;
780: }
781: if (normtype == SNES_NORM_ALWAYS) {
782: if (i == snes->max_its) {
783: PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", snes->max_its);
784: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
785: }
786: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
787: return 0;
788: }
790: /*MC
791: SNESCOMPOSITE - Build a preconditioner by composing together several nonlinear solvers
793: Options Database Keys:
794: + -snes_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
795: - -snes_composite_sneses - <snes0,snes1,...> list of SNESes to compose
797: Level: intermediate
799: References:
800: . * - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
801: SIAM Review, 57(4), 2015
803: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`,
804: `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCompositeSetType()`,
805: `SNESCompositeSetDamping()`
806: M*/
808: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
809: {
810: SNES_Composite *jac;
812: PetscNew(&jac);
814: snes->ops->solve = SNESSolve_Composite;
815: snes->ops->setup = SNESSetUp_Composite;
816: snes->ops->reset = SNESReset_Composite;
817: snes->ops->destroy = SNESDestroy_Composite;
818: snes->ops->setfromoptions = SNESSetFromOptions_Composite;
819: snes->ops->view = SNESView_Composite;
821: snes->usesksp = PETSC_FALSE;
823: snes->alwayscomputesfinalresidual = PETSC_FALSE;
825: snes->data = (void *)jac;
826: jac->type = SNES_COMPOSITE_ADDITIVEOPTIMAL;
827: jac->Fes = NULL;
828: jac->Xes = NULL;
829: jac->fnorms = NULL;
830: jac->nsnes = 0;
831: jac->head = NULL;
832: jac->stol = 0.1;
833: jac->rtol = 1.1;
835: jac->h = NULL;
836: jac->s = NULL;
837: jac->beta = NULL;
838: jac->work = NULL;
839: jac->rwork = NULL;
841: PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetType_C", SNESCompositeSetType_Composite);
842: PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeAddSNES_C", SNESCompositeAddSNES_Composite);
843: PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeGetSNES_C", SNESCompositeGetSNES_Composite);
844: PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetDamping_C", SNESCompositeSetDamping_Composite);
845: return 0;
846: }