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) {
57: SNESSetInitialFunction(next->snes,F);
58: }
59: SNESSolve(next->snes,B,X);
60: SNESGetConvergedReason(next->snes,&reason);
61: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
62: jac->innerFailures++;
63: if (jac->innerFailures >= snes->maxFailures) {
64: snes->reason = SNES_DIVERGED_INNER;
65: return 0;
66: }
67: }
69: while (next->next) {
70: /* only copy the function over in the case where the functions correspond */
71: if (next->snes->npcside== PC_RIGHT && next->snes->normschedule != SNES_NORM_NONE) {
72: SNESGetFunction(next->snes,&FSub,NULL,NULL);
73: next = next->next;
74: SNESSetInitialFunction(next->snes,FSub);
75: } else {
76: next = next->next;
77: }
78: SNESSolve(next->snes,B,X);
79: SNESGetConvergedReason(next->snes,&reason);
80: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
81: jac->innerFailures++;
82: if (jac->innerFailures >= snes->maxFailures) {
83: snes->reason = SNES_DIVERGED_INNER;
84: return 0;
85: }
86: }
87: }
88: if (next->snes->npcside== PC_RIGHT) {
89: SNESGetFunction(next->snes,&FSub,NULL,NULL);
90: VecCopy(FSub,F);
91: if (fnorm) {
92: if (snes->xl && snes->xu) {
93: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
94: } else {
95: VecNorm(F, NORM_2, fnorm);
96: }
97: SNESCheckFunctionNorm(snes,*fnorm);
98: }
99: } else if (snes->normschedule == SNES_NORM_ALWAYS) {
100: SNESComputeFunction(snes,X,F);
101: if (fnorm) {
102: if (snes->xl && snes->xu) {
103: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
104: } else {
105: VecNorm(F, NORM_2, fnorm);
106: }
107: SNESCheckFunctionNorm(snes,*fnorm);
108: }
109: }
110: return 0;
111: }
113: static PetscErrorCode SNESCompositeApply_Additive(SNES snes,Vec X,Vec B,Vec F,PetscReal *fnorm)
114: {
115: SNES_Composite *jac = (SNES_Composite*)snes->data;
116: SNES_CompositeLink next = jac->head;
117: Vec Y,Xorig;
118: SNESConvergedReason reason;
120: Y = snes->vec_sol_update;
121: if (!jac->Xorig) VecDuplicate(X,&jac->Xorig);
122: Xorig = jac->Xorig;
123: VecCopy(X,Xorig);
125: if (snes->normschedule == SNES_NORM_ALWAYS) {
126: SNESSetInitialFunction(next->snes,F);
127: while (next->next) {
128: next = next->next;
129: SNESSetInitialFunction(next->snes,F);
130: }
131: }
132: next = jac->head;
133: VecCopy(Xorig,Y);
134: SNESSolve(next->snes,B,Y);
135: 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: return 0;
141: }
142: }
143: VecAXPY(Y,-1.0,Xorig);
144: VecAXPY(X,next->dmp,Y);
145: while (next->next) {
146: next = next->next;
147: VecCopy(Xorig,Y);
148: SNESSolve(next->snes,B,Y);
149: 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: return 0;
155: }
156: }
157: VecAXPY(Y,-1.0,Xorig);
158: VecAXPY(X,next->dmp,Y);
159: }
160: if (snes->normschedule == SNES_NORM_ALWAYS) {
161: SNESComputeFunction(snes,X,F);
162: if (fnorm) {
163: if (snes->xl && snes->xu) {
164: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
165: } else {
166: VecNorm(F, NORM_2, fnorm);
167: }
168: SNESCheckFunctionNorm(snes,*fnorm);
169: }
170: }
171: return 0;
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;
198: if (snes->normschedule == SNES_NORM_ALWAYS) {
199: next = jac->head;
200: SNESSetInitialFunction(next->snes,F);
201: while (next->next) {
202: next = next->next;
203: SNESSetInitialFunction(next->snes,F);
204: }
205: }
207: next = jac->head;
208: i = 0;
209: VecCopy(X,Xes[i]);
210: SNESSolve(next->snes,B,Xes[i]);
211: 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: return 0;
217: }
218: }
219: while (next->next) {
220: i++;
221: next = next->next;
222: VecCopy(X,Xes[i]);
223: SNESSolve(next->snes,B,Xes[i]);
224: 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: return 0;
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++) {
237: VecDotBegin(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
238: }
239: VecDotBegin(Fes[i],F,&jac->g[i]);
240: }
242: for (i=0;i<jac->n;i++) {
243: for (j=0;j<i+1;j++) {
244: VecDotEnd(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
245: if (i == j) jac->fnorms[i] = PetscSqrtReal(PetscRealPart(jac->h[i + j*jac->n]));
246: }
247: VecDotEnd(Fes[i],F,&jac->g[i]);
248: }
250: ftf = (*fnorm)*(*fnorm);
252: for (i=0; i<jac->n; i++) {
253: for (j=i+1;j<jac->n;j++) {
254: jac->h[i + j*jac->n] = jac->h[j + i*jac->n];
255: }
256: }
258: for (i=0; i<jac->n; i++) {
259: for (j=0; j<jac->n; j++) {
260: jac->h[i + j*jac->n] = jac->h[i + j*jac->n] - jac->g[j] - jac->g[i] + ftf;
261: }
262: jac->beta[i] = ftf - jac->g[i];
263: }
265: jac->info = 0;
266: jac->rcond = -1.;
267: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
268: #if defined(PETSC_USE_COMPLEX)
269: PetscStackCallBLAS("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));
270: #else
271: PetscStackCallBLAS("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));
272: #endif
273: PetscFPTrapPop();
276: tot = 0.;
277: total = 0.;
278: for (i=0; i<jac->n; i++) {
280: PetscInfo(snes,"%D: %g\n",i,(double)PetscRealPart(jac->beta[i]));
281: tot += jac->beta[i];
282: total += PetscAbsScalar(jac->beta[i]);
283: }
284: VecScale(X,(1. - tot));
285: VecMAXPY(X,jac->n,jac->beta,Xes);
286: SNESComputeFunction(snes,X,F);
288: if (snes->xl && snes->xu) {
289: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
290: } else {
291: VecNorm(F, NORM_2, fnorm);
292: }
294: /* take the minimum-normed candidate if it beats the combination by a factor of rtol or the combination has stagnated */
295: min_fnorm = jac->fnorms[0];
296: min_i = 0;
297: for (i=0; i<jac->n; i++) {
298: if (jac->fnorms[i] < min_fnorm) {
299: min_fnorm = jac->fnorms[i];
300: min_i = i;
301: }
302: }
304: /* stagnation or divergence restart to the solution of the solver that failed the least */
305: if (PetscRealPart(total) < jac->stol || min_fnorm*jac->rtol < *fnorm) {
306: VecCopy(jac->Xes[min_i],X);
307: VecCopy(jac->Fes[min_i],F);
308: *fnorm = min_fnorm;
309: }
310: return 0;
311: }
313: static PetscErrorCode SNESSetUp_Composite(SNES snes)
314: {
315: DM dm;
316: SNES_Composite *jac = (SNES_Composite*)snes->data;
317: SNES_CompositeLink next = jac->head;
318: PetscInt n=0,i;
319: Vec F;
321: SNESGetDM(snes,&dm);
323: if (snes->ops->computevariablebounds) {
324: /* SNESVI only ever calls computevariablebounds once, so calling it once here is justified */
325: if (!snes->xl) VecDuplicate(snes->vec_sol,&snes->xl);
326: if (!snes->xu) VecDuplicate(snes->vec_sol,&snes->xu);
327: (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);
328: }
330: while (next) {
331: n++;
332: SNESSetDM(next->snes,dm);
333: SNESSetApplicationContext(next->snes, snes->user);
334: if (snes->xl && snes->xu) {
335: if (snes->ops->computevariablebounds) {
336: SNESVISetComputeVariableBounds(next->snes, snes->ops->computevariablebounds);
337: } else {
338: SNESVISetVariableBounds(next->snes,snes->xl,snes->xu);
339: }
340: }
342: next = next->next;
343: }
344: jac->nsnes = n;
345: SNESGetFunction(snes,&F,NULL,NULL);
346: if (jac->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
347: VecDuplicateVecs(F,jac->nsnes,&jac->Xes);
348: PetscMalloc1(n,&jac->Fes);
349: PetscMalloc1(n,&jac->fnorms);
350: next = jac->head;
351: i = 0;
352: while (next) {
353: SNESGetFunction(next->snes,&F,NULL,NULL);
354: jac->Fes[i] = F;
355: PetscObjectReference((PetscObject)F);
356: next = next->next;
357: i++;
358: }
359: /* allocate the subspace direct solve area */
360: jac->nrhs = 1;
361: jac->lda = jac->nsnes;
362: jac->ldb = jac->nsnes;
363: jac->n = jac->nsnes;
365: PetscMalloc1(jac->n*jac->n,&jac->h);
366: PetscMalloc1(jac->n,&jac->beta);
367: PetscMalloc1(jac->n,&jac->s);
368: PetscMalloc1(jac->n,&jac->g);
369: jac->lwork = 12*jac->n;
370: #if defined(PETSC_USE_COMPLEX)
371: PetscMalloc1(jac->lwork,&jac->rwork);
372: #endif
373: PetscMalloc1(jac->lwork,&jac->work);
374: }
376: return 0;
377: }
379: static PetscErrorCode SNESReset_Composite(SNES snes)
380: {
381: SNES_Composite *jac = (SNES_Composite*)snes->data;
382: SNES_CompositeLink next = jac->head;
384: while (next) {
385: SNESReset(next->snes);
386: next = next->next;
387: }
388: VecDestroy(&jac->Xorig);
389: if (jac->Xes) VecDestroyVecs(jac->nsnes,&jac->Xes);
390: if (jac->Fes) VecDestroyVecs(jac->nsnes,&jac->Fes);
391: PetscFree(jac->fnorms);
392: PetscFree(jac->h);
393: PetscFree(jac->s);
394: PetscFree(jac->g);
395: PetscFree(jac->beta);
396: PetscFree(jac->work);
397: PetscFree(jac->rwork);
398: return 0;
399: }
401: static PetscErrorCode SNESDestroy_Composite(SNES snes)
402: {
403: SNES_Composite *jac = (SNES_Composite*)snes->data;
404: SNES_CompositeLink next = jac->head,next_tmp;
406: SNESReset_Composite(snes);
407: while (next) {
408: SNESDestroy(&next->snes);
409: next_tmp = next;
410: next = next->next;
411: PetscFree(next_tmp);
412: }
413: PetscFree(snes->data);
414: return 0;
415: }
417: static PetscErrorCode SNESSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,SNES snes)
418: {
419: SNES_Composite *jac = (SNES_Composite*)snes->data;
420: PetscInt nmax = 8,i;
421: SNES_CompositeLink next;
422: char *sneses[8];
423: PetscReal dmps[8];
424: PetscBool flg;
426: PetscOptionsHead(PetscOptionsObject,"Composite preconditioner options");
427: PetscOptionsEnum("-snes_composite_type","Type of composition","SNESCompositeSetType",SNESCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
428: if (flg) {
429: SNESCompositeSetType(snes,jac->type);
430: }
431: PetscOptionsStringArray("-snes_composite_sneses","List of composite solvers","SNESCompositeAddSNES",sneses,&nmax,&flg);
432: if (flg) {
433: for (i=0; i<nmax; i++) {
434: SNESCompositeAddSNES(snes,sneses[i]);
435: PetscFree(sneses[i]); /* deallocate string sneses[i], which is allocated in PetscOptionsStringArray() */
436: }
437: }
438: PetscOptionsRealArray("-snes_composite_damping","Damping of the additive composite solvers","SNESCompositeSetDamping",dmps,&nmax,&flg);
439: if (flg) {
440: for (i=0; i<nmax; i++) {
441: SNESCompositeSetDamping(snes,i,dmps[i]);
442: }
443: }
444: PetscOptionsReal("-snes_composite_stol","Step tolerance for restart on the additive composite solvers","",jac->stol,&jac->stol,NULL);
445: PetscOptionsReal("-snes_composite_rtol","Residual tolerance for the additive composite solvers","",jac->rtol,&jac->rtol,NULL);
446: PetscOptionsTail();
448: next = jac->head;
449: while (next) {
450: SNESSetFromOptions(next->snes);
451: next = next->next;
452: }
453: return 0;
454: }
456: static PetscErrorCode SNESView_Composite(SNES snes,PetscViewer viewer)
457: {
458: SNES_Composite *jac = (SNES_Composite*)snes->data;
459: SNES_CompositeLink next = jac->head;
460: PetscBool iascii;
462: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
463: if (iascii) {
464: PetscViewerASCIIPrintf(viewer," type - %s\n",SNESCompositeTypes[jac->type]);
465: PetscViewerASCIIPrintf(viewer," SNESes on composite preconditioner follow\n");
466: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
467: }
468: if (iascii) {
469: PetscViewerASCIIPushTab(viewer);
470: }
471: while (next) {
472: SNESView(next->snes,viewer);
473: next = next->next;
474: }
475: if (iascii) {
476: PetscViewerASCIIPopTab(viewer);
477: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
478: }
479: return 0;
480: }
482: /* ------------------------------------------------------------------------------*/
484: static PetscErrorCode SNESCompositeSetType_Composite(SNES snes,SNESCompositeType type)
485: {
486: SNES_Composite *jac = (SNES_Composite*)snes->data;
488: jac->type = type;
489: return 0;
490: }
492: static PetscErrorCode SNESCompositeAddSNES_Composite(SNES snes,SNESType type)
493: {
494: SNES_Composite *jac;
495: SNES_CompositeLink next,ilink;
496: PetscInt cnt = 0;
497: const char *prefix;
498: char newprefix[20];
499: DM dm;
501: PetscNewLog(snes,&ilink);
502: ilink->next = NULL;
503: SNESCreate(PetscObjectComm((PetscObject)snes),&ilink->snes);
504: PetscLogObjectParent((PetscObject)snes,(PetscObject)ilink->snes);
505: SNESGetDM(snes,&dm);
506: SNESSetDM(ilink->snes,dm);
507: SNESSetTolerances(ilink->snes,snes->abstol,snes->rtol,snes->stol,1,snes->max_funcs);
508: PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)ilink->snes);
509: jac = (SNES_Composite*)snes->data;
510: next = jac->head;
511: if (!next) {
512: jac->head = ilink;
513: ilink->previous = NULL;
514: } else {
515: cnt++;
516: while (next->next) {
517: next = next->next;
518: cnt++;
519: }
520: next->next = ilink;
521: ilink->previous = next;
522: }
523: SNESGetOptionsPrefix(snes,&prefix);
524: SNESSetOptionsPrefix(ilink->snes,prefix);
525: PetscSNPrintf(newprefix,sizeof(newprefix),"sub_%d_",(int)cnt);
526: SNESAppendOptionsPrefix(ilink->snes,newprefix);
527: PetscObjectIncrementTabLevel((PetscObject)ilink->snes,(PetscObject)snes,1);
528: SNESSetType(ilink->snes,type);
529: SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY);
531: ilink->dmp = 1.0;
532: jac->nsnes++;
533: return 0;
534: }
536: static PetscErrorCode SNESCompositeGetSNES_Composite(SNES snes,PetscInt n,SNES *subsnes)
537: {
538: SNES_Composite *jac;
539: SNES_CompositeLink next;
540: PetscInt i;
542: jac = (SNES_Composite*)snes->data;
543: next = jac->head;
544: for (i=0; i<n; i++) {
546: next = next->next;
547: }
548: *subsnes = next->snes;
549: return 0;
550: }
552: /* -------------------------------------------------------------------------------- */
553: /*@C
554: SNESCompositeSetType - Sets the type of composite preconditioner.
556: Logically Collective on SNES
558: Input Parameters:
559: + snes - the preconditioner context
560: - type - SNES_COMPOSITE_ADDITIVE (default), SNES_COMPOSITE_MULTIPLICATIVE
562: Options Database Key:
563: . -snes_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
565: Level: Developer
567: @*/
568: PetscErrorCode SNESCompositeSetType(SNES snes,SNESCompositeType type)
569: {
572: PetscTryMethod(snes,"SNESCompositeSetType_C",(SNES,SNESCompositeType),(snes,type));
573: return 0;
574: }
576: /*@C
577: SNESCompositeAddSNES - Adds another SNES to the composite SNES.
579: Collective on SNES
581: Input Parameters:
582: + snes - the preconditioner context
583: - type - the type of the new preconditioner
585: Level: Developer
587: @*/
588: PetscErrorCode SNESCompositeAddSNES(SNES snes,SNESType type)
589: {
591: PetscTryMethod(snes,"SNESCompositeAddSNES_C",(SNES,SNESType),(snes,type));
592: return 0;
593: }
595: /*@
596: SNESCompositeGetSNES - Gets one of the SNES objects in the composite SNES.
598: Not Collective
600: Input Parameters:
601: + snes - the preconditioner context
602: - n - the number of the snes requested
604: Output Parameters:
605: . subsnes - the SNES requested
607: Level: Developer
609: .seealso: SNESCompositeAddSNES()
610: @*/
611: PetscErrorCode SNESCompositeGetSNES(SNES snes,PetscInt n,SNES *subsnes)
612: {
615: PetscUseMethod(snes,"SNESCompositeGetSNES_C",(SNES,PetscInt,SNES*),(snes,n,subsnes));
616: return 0;
617: }
619: /*@
620: SNESCompositeGetNumber - Get the number of subsolvers in the composite SNES.
622: Logically Collective on SNES
624: Input Parameter:
625: snes - the preconditioner context
627: Output Parameter:
628: n - the number of subsolvers
630: Level: Developer
632: @*/
633: PetscErrorCode SNESCompositeGetNumber(SNES snes,PetscInt *n)
634: {
635: SNES_Composite *jac;
636: SNES_CompositeLink next;
638: jac = (SNES_Composite*)snes->data;
639: next = jac->head;
641: *n = 0;
642: while (next) {
643: *n = *n + 1;
644: next = next->next;
645: }
646: return 0;
647: }
649: static PetscErrorCode SNESCompositeSetDamping_Composite(SNES snes,PetscInt n,PetscReal dmp)
650: {
651: SNES_Composite *jac;
652: SNES_CompositeLink next;
653: PetscInt i;
655: jac = (SNES_Composite*)snes->data;
656: next = jac->head;
657: for (i=0; i<n; i++) {
659: next = next->next;
660: }
661: next->dmp = dmp;
662: return 0;
663: }
665: /*@
666: SNESCompositeSetDamping - Sets the damping of a subsolver when using additive composite SNES.
668: Not Collective
670: Input Parameters:
671: + snes - the preconditioner context
672: . n - the number of the snes requested
673: - dmp - the damping
675: Level: Developer
677: .seealso: SNESCompositeAddSNES()
678: @*/
679: PetscErrorCode SNESCompositeSetDamping(SNES snes,PetscInt n,PetscReal dmp)
680: {
682: PetscUseMethod(snes,"SNESCompositeSetDamping_C",(SNES,PetscInt,PetscReal),(snes,n,dmp));
683: return 0;
684: }
686: static PetscErrorCode SNESSolve_Composite(SNES snes)
687: {
688: Vec F,X,B,Y;
689: PetscInt i;
690: PetscReal fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
691: SNESNormSchedule normtype;
692: SNES_Composite *comp = (SNES_Composite*)snes->data;
694: X = snes->vec_sol;
695: F = snes->vec_func;
696: B = snes->vec_rhs;
697: Y = snes->vec_sol_update;
699: PetscObjectSAWsTakeAccess((PetscObject)snes);
700: snes->iter = 0;
701: snes->norm = 0.;
702: comp->innerFailures = 0;
703: PetscObjectSAWsGrantAccess((PetscObject)snes);
704: snes->reason = SNES_CONVERGED_ITERATING;
705: SNESGetNormSchedule(snes, &normtype);
706: if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
707: if (!snes->vec_func_init_set) {
708: SNESComputeFunction(snes,X,F);
709: } else snes->vec_func_init_set = PETSC_FALSE;
711: if (snes->xl && snes->xu) {
712: SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
713: } else {
714: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
715: }
716: SNESCheckFunctionNorm(snes,fnorm);
717: PetscObjectSAWsTakeAccess((PetscObject)snes);
718: snes->iter = 0;
719: snes->norm = fnorm;
720: PetscObjectSAWsGrantAccess((PetscObject)snes);
721: SNESLogConvergenceHistory(snes,snes->norm,0);
722: SNESMonitor(snes,0,snes->norm);
724: /* test convergence */
725: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
726: if (snes->reason) return 0;
727: } else {
728: PetscObjectSAWsGrantAccess((PetscObject)snes);
729: SNESLogConvergenceHistory(snes,snes->norm,0);
730: SNESMonitor(snes,0,snes->norm);
731: }
733: for (i = 0; i < snes->max_its; i++) {
734: /* Call general purpose update function */
735: if (snes->ops->update) {
736: (*snes->ops->update)(snes, snes->iter);
737: }
739: /* Copy the state before modification by application of the composite solver;
740: we will subtract the new state after application */
741: VecCopy(X, Y);
743: if (comp->type == SNES_COMPOSITE_ADDITIVE) {
744: SNESCompositeApply_Additive(snes,X,B,F,&fnorm);
745: } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
746: SNESCompositeApply_Multiplicative(snes,X,B,F,&fnorm);
747: } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
748: 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: VecAYPX(Y, -1.0, X);
755: if ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY)) {
756: SNESComputeFunction(snes,X,F);
758: if (snes->xl && snes->xu) {
759: VecNormBegin(X, NORM_2, &xnorm);
760: VecNormBegin(Y, NORM_2, &snorm);
761: SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
762: VecNormEnd(X, NORM_2, &xnorm);
763: VecNormEnd(Y, NORM_2, &snorm);
764: } else {
765: VecNormBegin(F, NORM_2, &fnorm);
766: VecNormBegin(X, NORM_2, &xnorm);
767: VecNormBegin(Y, NORM_2, &snorm);
769: VecNormEnd(F, NORM_2, &fnorm);
770: VecNormEnd(X, NORM_2, &xnorm);
771: VecNormEnd(Y, NORM_2, &snorm);
772: }
773: SNESCheckFunctionNorm(snes,fnorm);
774: } else if (normtype == SNES_NORM_ALWAYS) {
775: VecNormBegin(X, NORM_2, &xnorm);
776: VecNormBegin(Y, NORM_2, &snorm);
777: VecNormEnd(X, NORM_2, &xnorm);
778: VecNormEnd(Y, NORM_2, &snorm);
779: }
780: /* Monitor convergence */
781: PetscObjectSAWsTakeAccess((PetscObject)snes);
782: snes->iter = i+1;
783: snes->norm = fnorm;
784: snes->xnorm = xnorm;
785: snes->ynorm = snorm;
786: PetscObjectSAWsGrantAccess((PetscObject)snes);
787: SNESLogConvergenceHistory(snes,snes->norm,0);
788: SNESMonitor(snes,snes->iter,snes->norm);
789: /* Test for convergence */
790: if (normtype == SNES_NORM_ALWAYS) (*snes->ops->converged)(snes,snes->iter,xnorm,snorm,fnorm,&snes->reason,snes->cnvP);
791: if (snes->reason) break;
792: }
793: if (normtype == SNES_NORM_ALWAYS) {
794: if (i == snes->max_its) {
795: PetscInfo(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
796: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
797: }
798: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
799: return 0;
800: }
802: /* -------------------------------------------------------------------------------------------*/
804: /*MC
805: SNESCOMPOSITE - Build a preconditioner by composing together several nonlinear solvers
807: Options Database Keys:
808: + -snes_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
809: - -snes_composite_sneses - <snes0,snes1,...> list of SNESes to compose
811: Level: intermediate
813: .seealso: SNESCreate(), SNESSetType(), SNESType (for list of available types), SNES,
814: SNESSHELL, SNESCompositeSetType(), SNESCompositeSpecialSetAlpha(), SNESCompositeAddSNES(),
815: SNESCompositeGetSNES()
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: M*/
823: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
824: {
825: SNES_Composite *jac;
827: PetscNewLog(snes,&jac);
829: snes->ops->solve = SNESSolve_Composite;
830: snes->ops->setup = SNESSetUp_Composite;
831: snes->ops->reset = SNESReset_Composite;
832: snes->ops->destroy = SNESDestroy_Composite;
833: snes->ops->setfromoptions = SNESSetFromOptions_Composite;
834: snes->ops->view = SNESView_Composite;
836: snes->usesksp = PETSC_FALSE;
838: snes->alwayscomputesfinalresidual = PETSC_FALSE;
840: snes->data = (void*)jac;
841: jac->type = SNES_COMPOSITE_ADDITIVEOPTIMAL;
842: jac->Fes = NULL;
843: jac->Xes = NULL;
844: jac->fnorms = NULL;
845: jac->nsnes = 0;
846: jac->head = NULL;
847: jac->stol = 0.1;
848: jac->rtol = 1.1;
850: jac->h = NULL;
851: jac->s = NULL;
852: jac->beta = NULL;
853: jac->work = NULL;
854: jac->rwork = NULL;
856: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetType_C",SNESCompositeSetType_Composite);
857: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeAddSNES_C",SNESCompositeAddSNES_Composite);
858: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeGetSNES_C",SNESCompositeGetSNES_Composite);
859: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetDamping_C",SNESCompositeSetDamping_Composite);
860: return 0;
861: }