Actual source code: snescomposite.c
petsc-3.5.4 2015-05-23
2: /*
3: Defines a SNES that can consist of a collection of SNESes
4: */
5: #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/
6: #include <petscblaslapack.h>
8: const char *const SNESCompositeTypes[] = {"ADDITIVE","MULTIPLICATIVE","ADDITIVEOPTIMAL","SNESCompositeType","SNES_COMPOSITE",0};
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;
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;
49: static PetscErrorCode SNESCompositeApply_Multiplicative(SNES snes,Vec X,Vec B,Vec F,PetscReal *fnorm)
50: {
51: PetscErrorCode ierr;
52: SNES_Composite *jac = (SNES_Composite*)snes->data;
53: SNES_CompositeLink next = jac->head;
54: Vec FSub;
55: SNESConvergedReason reason;
58: if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
59: if (snes->normschedule == SNES_NORM_ALWAYS) {
60: SNESSetInitialFunction(next->snes,F);
61: }
62: SNESSolve(next->snes,B,X);
63: SNESGetConvergedReason(next->snes,&reason);
64: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
65: snes->reason = SNES_DIVERGED_INNER;
66: return(0);
67: }
69: while (next->next) {
70: /* only copy the function over in the case where the functions correspond */
71: if (next->snes->pcside == 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: snes->reason = SNES_DIVERGED_INNER;
82: return(0);
83: }
84: }
85: if (next->snes->pcside == PC_RIGHT) {
86: SNESGetFunction(next->snes,&FSub,NULL,NULL);
87: VecCopy(FSub,F);
88: if (fnorm) {
89: VecNorm(F,NORM_2,fnorm);
90: if (PetscIsInfOrNanReal(*fnorm)) {
91: snes->reason = SNES_DIVERGED_FNORM_NAN;
92: return(0);
93: }
94: }
95: } else if (snes->normschedule == SNES_NORM_ALWAYS) {
96: SNESComputeFunction(snes,X,F);
97: if (snes->domainerror) {
98: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
99: return(0);
100: }
101: if (fnorm) {
102: VecNorm(F,NORM_2,fnorm);
103: if (PetscIsInfOrNanReal(*fnorm)) {
104: snes->reason = SNES_DIVERGED_FNORM_NAN;
105: return(0);
106: }
107: }
108: }
109: return(0);
110: }
114: static PetscErrorCode SNESCompositeApply_Additive(SNES snes,Vec X,Vec B,Vec F,PetscReal *fnorm)
115: {
116: PetscErrorCode ierr;
117: SNES_Composite *jac = (SNES_Composite*)snes->data;
118: SNES_CompositeLink next = jac->head;
119: Vec Y,Xorig;
120: SNESConvergedReason reason;
123: Y = snes->vec_sol_update;
124: if (!jac->Xorig) {VecDuplicate(X,&jac->Xorig);}
125: Xorig = jac->Xorig;
126: VecCopy(X,Xorig);
127: if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
128: if (snes->normschedule == SNES_NORM_ALWAYS) {
129: SNESSetInitialFunction(next->snes,F);
130: while (next->next) {
131: next = next->next;
132: SNESSetInitialFunction(next->snes,F);
133: }
134: }
135: next = jac->head;
136: VecCopy(Xorig,Y);
137: SNESSolve(next->snes,B,Y);
138: SNESGetConvergedReason(next->snes,&reason);
139: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
140: snes->reason = SNES_DIVERGED_INNER;
141: return(0);
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: snes->reason = SNES_DIVERGED_INNER;
152: return(0);
153: }
154: VecAXPY(Y,-1.0,Xorig);
155: VecAXPY(X,next->dmp,Y);
156: }
157: if (snes->normschedule == SNES_NORM_ALWAYS) {
158: SNESComputeFunction(snes,X,F);
159: if (fnorm) {VecNorm(F,NORM_2,fnorm);}
160: }
161: return(0);
162: }
166: /* approximately solve the overdetermined system:
168: 2*F(x_i)\cdot F(\x_j)\alpha_i = 0
169: \alpha_i = 1
171: Which minimizes the L2 norm of the linearization of:
172: ||F(\sum_i \alpha_i*x_i)||^2
174: With the constraint that \sum_i\alpha_i = 1
175: Where x_i is the solution from the ith subsolver.
176: */
177: static PetscErrorCode SNESCompositeApply_AdditiveOptimal(SNES snes,Vec X,Vec B,Vec F,PetscReal *fnorm)
178: {
179: PetscErrorCode ierr;
180: SNES_Composite *jac = (SNES_Composite*)snes->data;
181: SNES_CompositeLink next = jac->head;
182: Vec *Xes = jac->Xes,*Fes = jac->Fes;
183: PetscInt i,j;
184: PetscScalar tot,total,ftf;
185: PetscReal min_fnorm;
186: PetscInt min_i;
187: SNESConvergedReason reason;
190: if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
192: if (snes->normschedule == SNES_NORM_ALWAYS) {
193: next = jac->head;
194: SNESSetInitialFunction(next->snes,F);
195: while (next->next) {
196: next = next->next;
197: SNESSetInitialFunction(next->snes,F);
198: }
199: }
201: next = jac->head;
202: i = 0;
203: VecCopy(X,Xes[i]);
204: SNESSolve(next->snes,B,Xes[i]);
205: SNESGetConvergedReason(next->snes,&reason);
206: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
207: snes->reason = SNES_DIVERGED_INNER;
208: return(0);
209: }
210: while (next->next) {
211: i++;
212: next = next->next;
213: VecCopy(X,Xes[i]);
214: SNESSolve(next->snes,B,Xes[i]);
215: SNESGetConvergedReason(next->snes,&reason);
216: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
217: snes->reason = SNES_DIVERGED_INNER;
218: return(0);
219: }
220: }
222: /* all the solutions are collected; combine optimally */
223: for (i=0;i<jac->n;i++) {
224: for (j=0;j<i+1;j++) {
225: VecDotBegin(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
226: }
227: VecDotBegin(Fes[i],F,&jac->g[i]);
228: }
230: for (i=0;i<jac->n;i++) {
231: for (j=0;j<i+1;j++) {
232: VecDotEnd(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
233: if (i == j) jac->fnorms[i] = PetscSqrtReal(PetscRealPart(jac->h[i + j*jac->n]));
234: }
235: VecDotEnd(Fes[i],F,&jac->g[i]);
236: }
238: ftf = (*fnorm)*(*fnorm);
240: for (i=0; i<jac->n; i++) {
241: for (j=i+1;j<jac->n;j++) {
242: jac->h[i + j*jac->n] = jac->h[j + i*jac->n];
243: }
244: }
246: for (i=0; i<jac->n; i++) {
247: for (j=0; j<jac->n; j++) {
248: jac->h[i + j*jac->n] = jac->h[i + j*jac->n] - jac->g[j] - jac->g[i] + ftf;
249: }
250: jac->beta[i] = ftf - jac->g[i];
251: }
253: #if defined(PETSC_MISSING_LAPACK_GELSS)
254: SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"SNESCOMPOSITE with ADDITIVEOPTIMAL requires the LAPACK GELSS routine.");
255: #else
256: jac->info = 0;
257: jac->rcond = -1.;
258: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
259: #if defined(PETSC_USE_COMPLEX)
260: PetscStackCall("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));
261: #else
262: PetscStackCall("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));
263: #endif
264: PetscFPTrapPop();
265: if (jac->info < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"Bad argument to GELSS");
266: if (jac->info > 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD failed to converge");
267: #endif
268: tot = 0.;
269: total = 0.;
270: for (i=0; i<jac->n; i++) {
271: if (PetscIsInfOrNanScalar(jac->beta[i])) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD generated inconsistent output");
272: PetscInfo2(snes,"%d: %f\n",i,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);
279: VecNorm(F,NORM_2,fnorm);
281: /* take the minimum-normed candidate if it beats the combination by a factor of rtol or the combination has stagnated */
282: min_fnorm = jac->fnorms[0];
283: min_i = 0;
284: for (i=0; i<jac->n; i++) {
285: if (jac->fnorms[i] < min_fnorm) {
286: min_fnorm = jac->fnorms[i];
287: min_i = i;
288: }
289: }
291: /* stagnation or divergence restart to the solution of the solver that failed the least */
292: if (PetscRealPart(total) < jac->stol || min_fnorm*jac->rtol < *fnorm) {
293: VecCopy(jac->Xes[min_i],X);
294: VecCopy(jac->Fes[min_i],F);
295: *fnorm = min_fnorm;
296: }
297: return(0);
298: }
302: static PetscErrorCode SNESSetUp_Composite(SNES snes)
303: {
304: PetscErrorCode ierr;
305: DM dm;
306: SNES_Composite *jac = (SNES_Composite*)snes->data;
307: SNES_CompositeLink next = jac->head;
308: PetscInt n=0,i;
309: Vec F;
312: SNESGetDM(snes,&dm);
313: while (next) {
314: n++;
315: SNESSetDM(next->snes,dm);
316: SNESSetFromOptions(next->snes);
317: next = next->next;
318: }
319: jac->nsnes = n;
320: SNESGetFunction(snes,&F,NULL,NULL);
321: if (jac->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
322: VecDuplicateVecs(F,jac->nsnes,&jac->Xes);
323: PetscMalloc(sizeof(Vec)*n,&jac->Fes);
324: PetscMalloc(sizeof(PetscReal)*n,&jac->fnorms);
325: next = jac->head;
326: i = 0;
327: while (next) {
328: SNESGetFunction(next->snes,&F,NULL,NULL);
329: jac->Fes[i] = F;
330: PetscObjectReference((PetscObject)F);
331: next = next->next;
332: i++;
333: }
334: /* allocate the subspace direct solve area */
335: jac->nrhs = 1;
336: jac->lda = jac->nsnes;
337: jac->ldb = jac->nsnes;
338: jac->n = jac->nsnes;
340: PetscMalloc1(jac->n*jac->n,&jac->h);
341: PetscMalloc1(jac->n,&jac->beta);
342: PetscMalloc1(jac->n,&jac->s);
343: PetscMalloc1(jac->n,&jac->g);
344: jac->lwork = 12*jac->n;
345: #if PETSC_USE_COMPLEX
346: PetscMalloc(sizeof(PetscReal)*jac->lwork,&jac->rwork);
347: #endif
348: PetscMalloc(sizeof(PetscScalar)*jac->lwork,&jac->work);
349: }
351: return(0);
352: }
356: static PetscErrorCode SNESReset_Composite(SNES snes)
357: {
358: SNES_Composite *jac = (SNES_Composite*)snes->data;
359: PetscErrorCode ierr;
360: SNES_CompositeLink next = jac->head;
363: while (next) {
364: SNESReset(next->snes);
365: next = next->next;
366: }
367: VecDestroy(&jac->Xorig);
368: if (jac->Xes) {VecDestroyVecs(jac->nsnes,&jac->Xes);}
369: if (jac->Fes) {VecDestroyVecs(jac->nsnes,&jac->Fes);}
370: PetscFree(jac->fnorms);
371: PetscFree(jac->h);
372: PetscFree(jac->s);
373: PetscFree(jac->g);
374: PetscFree(jac->beta);
375: PetscFree(jac->work);
376: PetscFree(jac->rwork);
377: return(0);
378: }
382: static PetscErrorCode SNESDestroy_Composite(SNES snes)
383: {
384: SNES_Composite *jac = (SNES_Composite*)snes->data;
385: PetscErrorCode ierr;
386: SNES_CompositeLink next = jac->head,next_tmp;
389: SNESReset_Composite(snes);
390: while (next) {
391: SNESDestroy(&next->snes);
392: next_tmp = next;
393: next = next->next;
394: PetscFree(next_tmp);
395: }
396: PetscFree(snes->data);
397: return(0);
398: }
402: static PetscErrorCode SNESSetFromOptions_Composite(SNES snes)
403: {
404: SNES_Composite *jac = (SNES_Composite*)snes->data;
405: PetscErrorCode ierr;
406: PetscInt nmax = 8,i;
407: SNES_CompositeLink next;
408: char *sneses[8];
409: PetscReal dmps[8];
410: PetscBool flg;
413: PetscOptionsHead("Composite preconditioner options");
414: PetscOptionsEnum("-snes_composite_type","Type of composition","SNESCompositeSetType",SNESCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
415: if (flg) {
416: SNESCompositeSetType(snes,jac->type);
417: }
418: PetscOptionsStringArray("-snes_composite_sneses","List of composite solvers","SNESCompositeAddSNES",sneses,&nmax,&flg);
419: if (flg) {
420: for (i=0; i<nmax; i++) {
421: SNESCompositeAddSNES(snes,sneses[i]);
422: PetscFree(sneses[i]); /* deallocate string sneses[i], which is allocated in PetscOptionsStringArray() */
423: }
424: }
425: PetscOptionsRealArray("-snes_composite_damping","Damping of the additive composite solvers","SNESCompositeSetDamping",dmps,&nmax,&flg);
426: if (flg) {
427: for (i=0; i<nmax; i++) {
428: SNESCompositeSetDamping(snes,i,dmps[i]);
429: }
430: }
431: PetscOptionsReal("-snes_composite_stol","Step tolerance for restart on the additive composite solvers","",jac->stol,&jac->stol,NULL);
432: PetscOptionsReal("-snes_composite_rtol","Residual tolerance for the additive composite solvers","",jac->rtol,&jac->rtol,NULL);
433: PetscOptionsTail();
435: next = jac->head;
436: while (next) {
437: SNESSetFromOptions(next->snes);
438: next = next->next;
439: }
440: return(0);
441: }
445: static PetscErrorCode SNESView_Composite(SNES snes,PetscViewer viewer)
446: {
447: SNES_Composite *jac = (SNES_Composite*)snes->data;
448: PetscErrorCode ierr;
449: SNES_CompositeLink next = jac->head;
450: PetscBool iascii;
453: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
454: if (iascii) {
455: PetscViewerASCIIPrintf(viewer,"Composite SNES type - %s\n",SNESCompositeTypes[jac->type]);
456: PetscViewerASCIIPrintf(viewer,"SNESes on composite preconditioner follow\n");
457: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
458: }
459: if (iascii) {
460: PetscViewerASCIIPushTab(viewer);
461: }
462: while (next) {
463: SNESView(next->snes,viewer);
464: next = next->next;
465: }
466: if (iascii) {
467: PetscViewerASCIIPopTab(viewer);
468: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
469: }
470: return(0);
471: }
473: /* ------------------------------------------------------------------------------*/
477: static PetscErrorCode SNESCompositeSetType_Composite(SNES snes,SNESCompositeType type)
478: {
479: SNES_Composite *jac = (SNES_Composite*)snes->data;
482: jac->type = type;
483: return(0);
484: }
488: static PetscErrorCode SNESCompositeAddSNES_Composite(SNES snes,SNESType type)
489: {
490: SNES_Composite *jac;
491: SNES_CompositeLink next,ilink;
492: PetscErrorCode ierr;
493: PetscInt cnt = 0;
494: const char *prefix;
495: char newprefix[8];
496: DM dm;
499: PetscNewLog(snes,&ilink);
500: ilink->next = 0;
501: SNESCreate(PetscObjectComm((PetscObject)snes),&ilink->snes);
502: PetscLogObjectParent((PetscObject)snes,(PetscObject)ilink->snes);
503: SNESGetDM(snes,&dm);
504: SNESSetDM(ilink->snes,dm);
505: SNESSetTolerances(ilink->snes,ilink->snes->abstol,ilink->snes->rtol,ilink->snes->stol,1,ilink->snes->max_funcs);
506: jac = (SNES_Composite*)snes->data;
507: next = jac->head;
508: if (!next) {
509: jac->head = ilink;
510: ilink->previous = NULL;
511: } else {
512: cnt++;
513: while (next->next) {
514: next = next->next;
515: cnt++;
516: }
517: next->next = ilink;
518: ilink->previous = next;
519: }
520: SNESGetOptionsPrefix(snes,&prefix);
521: SNESSetOptionsPrefix(ilink->snes,prefix);
522: sprintf(newprefix,"sub_%d_",(int)cnt);
523: SNESAppendOptionsPrefix(ilink->snes,newprefix);
524: PetscObjectIncrementTabLevel((PetscObject)ilink->snes,(PetscObject)snes,1);
525: SNESSetType(ilink->snes,type);
526: SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY);
527: ilink->dmp = 1.0;
528: jac->nsnes++;
529: return(0);
530: }
534: static PetscErrorCode SNESCompositeGetSNES_Composite(SNES snes,PetscInt n,SNES *subsnes)
535: {
536: SNES_Composite *jac;
537: SNES_CompositeLink next;
538: PetscInt i;
541: jac = (SNES_Composite*)snes->data;
542: next = jac->head;
543: for (i=0; i<n; i++) {
544: if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
545: next = next->next;
546: }
547: *subsnes = next->snes;
548: return(0);
549: }
551: /* -------------------------------------------------------------------------------- */
554: /*@C
555: SNESCompositeSetType - Sets the type of composite preconditioner.
557: Logically Collective on SNES
559: Input Parameter:
560: + snes - the preconditioner context
561: - type - SNES_COMPOSITE_ADDITIVE (default), SNES_COMPOSITE_MULTIPLICATIVE
563: Options Database Key:
564: . -snes_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
566: Level: Developer
568: .keywords: SNES, set, type, composite preconditioner, additive, multiplicative
569: @*/
570: PetscErrorCode SNESCompositeSetType(SNES snes,SNESCompositeType type)
571: {
577: PetscTryMethod(snes,"SNESCompositeSetType_C",(SNES,SNESCompositeType),(snes,type));
578: return(0);
579: }
583: /*@C
584: SNESCompositeAddSNES - Adds another SNES to the composite SNES.
586: Collective on SNES
588: Input Parameters:
589: + snes - the preconditioner context
590: - type - the type of the new preconditioner
592: Level: Developer
594: .keywords: SNES, composite preconditioner, add
595: @*/
596: PetscErrorCode SNESCompositeAddSNES(SNES snes,SNESType type)
597: {
602: PetscTryMethod(snes,"SNESCompositeAddSNES_C",(SNES,SNESType),(snes,type));
603: return(0);
604: }
607: /*@
608: SNESCompositeGetSNES - Gets one of the SNES objects in the composite SNES.
610: Not Collective
612: Input Parameter:
613: + snes - the preconditioner context
614: - n - the number of the snes requested
616: Output Parameters:
617: . subsnes - the SNES requested
619: Level: Developer
621: .keywords: SNES, get, composite preconditioner, sub preconditioner
623: .seealso: SNESCompositeAddSNES()
624: @*/
625: PetscErrorCode SNESCompositeGetSNES(SNES snes,PetscInt n,SNES *subsnes)
626: {
632: PetscUseMethod(snes,"SNESCompositeGetSNES_C",(SNES,PetscInt,SNES*),(snes,n,subsnes));
633: return(0);
634: }
638: static PetscErrorCode SNESCompositeSetDamping_Composite(SNES snes,PetscInt n,PetscReal dmp)
639: {
640: SNES_Composite *jac;
641: SNES_CompositeLink next;
642: PetscInt i;
645: jac = (SNES_Composite*)snes->data;
646: next = jac->head;
647: for (i=0; i<n; i++) {
648: if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
649: next = next->next;
650: }
651: next->dmp = dmp;
652: return(0);
653: }
657: /*@
658: SNESCompositeSetDamping - Sets the damping of a subsolver when using additive composite SNES.
660: Not Collective
662: Input Parameter:
663: + snes - the preconditioner context
664: . n - the number of the snes requested
665: - dmp - the damping
667: Level: Developer
669: .keywords: SNES, get, composite preconditioner, sub preconditioner
671: .seealso: SNESCompositeAddSNES()
672: @*/
673: PetscErrorCode SNESCompositeSetDamping(SNES snes,PetscInt n,PetscReal dmp)
674: {
679: PetscUseMethod(snes,"SNESCompositeSetDamping_C",(SNES,PetscInt,PetscReal),(snes,n,dmp));
680: return(0);
681: }
685: PetscErrorCode SNESSolve_Composite(SNES snes)
686: {
687: Vec F;
688: Vec X;
689: Vec B;
690: PetscInt i;
691: PetscReal fnorm = 0.0;
693: SNESNormSchedule normtype;
694: SNES_Composite *comp = (SNES_Composite*)snes->data;
697: X = snes->vec_sol;
698: F = snes->vec_func;
699: B = snes->vec_rhs;
701: PetscObjectSAWsTakeAccess((PetscObject)snes);
702: snes->iter = 0;
703: snes->norm = 0.;
704: PetscObjectSAWsGrantAccess((PetscObject)snes);
705: snes->reason = SNES_CONVERGED_ITERATING;
706: SNESGetNormSchedule(snes, &normtype);
707: if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
708: if (!snes->vec_func_init_set) {
709: SNESComputeFunction(snes,X,F);
710: if (snes->domainerror) {
711: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
712: return(0);
713: }
714: } else snes->vec_func_init_set = PETSC_FALSE;
716: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
717: if (PetscIsInfOrNanReal(fnorm)) {
718: snes->reason = SNES_DIVERGED_FNORM_NAN;
719: return(0);
720: }
721: PetscObjectSAWsTakeAccess((PetscObject)snes);
722: snes->iter = 0;
723: snes->norm = fnorm;
724: PetscObjectSAWsGrantAccess((PetscObject)snes);
725: SNESLogConvergenceHistory(snes,snes->norm,0);
726: SNESMonitor(snes,0,snes->norm);
728: /* test convergence */
729: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
730: if (snes->reason) return(0);
731: } else {
732: PetscObjectSAWsGrantAccess((PetscObject)snes);
733: SNESLogConvergenceHistory(snes,snes->norm,0);
734: SNESMonitor(snes,0,snes->norm);
735: }
737: /* Call general purpose update function */
738: if (snes->ops->update) {
739: (*snes->ops->update)(snes, snes->iter);
740: }
742: for (i = 0; i < snes->max_its; i++) {
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 {
750: SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported SNESComposite type");
751: }
752: if (snes->reason < 0) break;
754: if ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY)) {
755: SNESComputeFunction(snes,X,F);
756: if (snes->domainerror) {
757: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
758: break;
759: }
760: VecNorm(F, NORM_2, &fnorm);
761: if (PetscIsInfOrNanReal(fnorm)) {
762: snes->reason = SNES_DIVERGED_FNORM_NAN;
763: break;
764: }
765: }
766: /* Monitor convergence */
767: PetscObjectSAWsTakeAccess((PetscObject)snes);
768: snes->iter = i+1;
769: snes->norm = fnorm;
770: PetscObjectSAWsGrantAccess((PetscObject)snes);
771: SNESLogConvergenceHistory(snes,snes->norm,0);
772: SNESMonitor(snes,snes->iter,snes->norm);
773: /* Test for convergence */
774: if (normtype == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);}
775: if (snes->reason) break;
776: /* Call general purpose update function */
777: if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
778: }
779: if (normtype == SNES_NORM_ALWAYS) {
780: if (i == snes->max_its) {
781: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
782: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
783: }
784: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
785: return(0);
786: }
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: Concepts: composing solvers
801: .seealso: SNESCreate(), SNESSetType(), SNESType (for list of available types), SNES,
802: SNESSHELL, SNESCompositeSetType(), SNESCompositeSpecialSetAlpha(), SNESCompositeAddSNES(),
803: SNESCompositeGetSNES()
805: M*/
809: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
810: {
812: SNES_Composite *jac;
815: PetscNewLog(snes,&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->data = (void*)jac;
825: jac->type = SNES_COMPOSITE_ADDITIVEOPTIMAL;
826: jac->Fes = NULL;
827: jac->Xes = NULL;
828: jac->fnorms = NULL;
829: jac->nsnes = 0;
830: jac->head = 0;
831: jac->stol = 0.1;
832: jac->rtol = 1.1;
834: jac->h = NULL;
835: jac->s = NULL;
836: jac->beta = NULL;
837: jac->work = NULL;
838: jac->rwork = NULL;
840: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetType_C",SNESCompositeSetType_Composite);
841: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeAddSNES_C",SNESCompositeAddSNES_Composite);
842: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeGetSNES_C",SNESCompositeGetSNES_Composite);
843: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetDamping_C",SNESCompositeSetDamping_Composite);
844: return(0);
845: }