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: PetscErrorCode ierr;
51: SNES_Composite *jac = (SNES_Composite*)snes->data;
52: SNES_CompositeLink next = jac->head;
53: Vec FSub;
54: SNESConvergedReason reason;
57: if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
58: if (snes->normschedule == SNES_NORM_ALWAYS) {
59: SNESSetInitialFunction(next->snes,F);
60: }
61: SNESSolve(next->snes,B,X);
62: SNESGetConvergedReason(next->snes,&reason);
63: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
64: jac->innerFailures++;
65: if (jac->innerFailures >= snes->maxFailures) {
66: snes->reason = SNES_DIVERGED_INNER;
67: return(0);
68: }
69: }
71: while (next->next) {
72: /* only copy the function over in the case where the functions correspond */
73: if (next->snes->npcside== PC_RIGHT && next->snes->normschedule != SNES_NORM_NONE) {
74: SNESGetFunction(next->snes,&FSub,NULL,NULL);
75: next = next->next;
76: SNESSetInitialFunction(next->snes,FSub);
77: } else {
78: next = next->next;
79: }
80: SNESSolve(next->snes,B,X);
81: SNESGetConvergedReason(next->snes,&reason);
82: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
83: jac->innerFailures++;
84: if (jac->innerFailures >= snes->maxFailures) {
85: snes->reason = SNES_DIVERGED_INNER;
86: return(0);
87: }
88: }
89: }
90: if (next->snes->npcside== PC_RIGHT) {
91: SNESGetFunction(next->snes,&FSub,NULL,NULL);
92: VecCopy(FSub,F);
93: if (fnorm) {
94: if (snes->xl && snes->xu) {
95: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
96: } else {
97: VecNorm(F, NORM_2, fnorm);
98: }
99: SNESCheckFunctionNorm(snes,*fnorm);
100: }
101: } else if (snes->normschedule == SNES_NORM_ALWAYS) {
102: SNESComputeFunction(snes,X,F);
103: if (fnorm) {
104: if (snes->xl && snes->xu) {
105: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
106: } else {
107: VecNorm(F, NORM_2, fnorm);
108: }
109: SNESCheckFunctionNorm(snes,*fnorm);
110: }
111: }
112: return(0);
113: }
115: static PetscErrorCode SNESCompositeApply_Additive(SNES snes,Vec X,Vec B,Vec F,PetscReal *fnorm)
116: {
117: PetscErrorCode ierr;
118: SNES_Composite *jac = (SNES_Composite*)snes->data;
119: SNES_CompositeLink next = jac->head;
120: Vec Y,Xorig;
121: SNESConvergedReason reason;
124: Y = snes->vec_sol_update;
125: if (!jac->Xorig) {VecDuplicate(X,&jac->Xorig);}
126: Xorig = jac->Xorig;
127: VecCopy(X,Xorig);
128: if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
129: if (snes->normschedule == SNES_NORM_ALWAYS) {
130: SNESSetInitialFunction(next->snes,F);
131: while (next->next) {
132: next = next->next;
133: SNESSetInitialFunction(next->snes,F);
134: }
135: }
136: next = jac->head;
137: VecCopy(Xorig,Y);
138: SNESSolve(next->snes,B,Y);
139: SNESGetConvergedReason(next->snes,&reason);
140: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
141: jac->innerFailures++;
142: if (jac->innerFailures >= snes->maxFailures) {
143: snes->reason = SNES_DIVERGED_INNER;
144: return(0);
145: }
146: }
147: VecAXPY(Y,-1.0,Xorig);
148: VecAXPY(X,next->dmp,Y);
149: while (next->next) {
150: next = next->next;
151: VecCopy(Xorig,Y);
152: SNESSolve(next->snes,B,Y);
153: SNESGetConvergedReason(next->snes,&reason);
154: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
155: jac->innerFailures++;
156: if (jac->innerFailures >= snes->maxFailures) {
157: snes->reason = SNES_DIVERGED_INNER;
158: return(0);
159: }
160: }
161: VecAXPY(Y,-1.0,Xorig);
162: VecAXPY(X,next->dmp,Y);
163: }
164: if (snes->normschedule == SNES_NORM_ALWAYS) {
165: SNESComputeFunction(snes,X,F);
166: if (fnorm) {
167: if (snes->xl && snes->xu) {
168: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
169: } else {
170: VecNorm(F, NORM_2, fnorm);
171: }
172: SNESCheckFunctionNorm(snes,*fnorm);
173: }
174: }
175: return(0);
176: }
178: /* approximately solve the overdetermined system:
180: 2*F(x_i)\cdot F(\x_j)\alpha_i = 0
181: \alpha_i = 1
183: Which minimizes the L2 norm of the linearization of:
184: ||F(\sum_i \alpha_i*x_i)||^2
186: With the constraint that \sum_i\alpha_i = 1
187: Where x_i is the solution from the ith subsolver.
188: */
189: static PetscErrorCode SNESCompositeApply_AdditiveOptimal(SNES snes,Vec X,Vec B,Vec F,PetscReal *fnorm)
190: {
191: PetscErrorCode ierr;
192: SNES_Composite *jac = (SNES_Composite*)snes->data;
193: SNES_CompositeLink next = jac->head;
194: Vec *Xes = jac->Xes,*Fes = jac->Fes;
195: PetscInt i,j;
196: PetscScalar tot,total,ftf;
197: PetscReal min_fnorm;
198: PetscInt min_i;
199: SNESConvergedReason reason;
202: if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
204: if (snes->normschedule == SNES_NORM_ALWAYS) {
205: next = jac->head;
206: SNESSetInitialFunction(next->snes,F);
207: while (next->next) {
208: next = next->next;
209: SNESSetInitialFunction(next->snes,F);
210: }
211: }
213: next = jac->head;
214: i = 0;
215: VecCopy(X,Xes[i]);
216: SNESSolve(next->snes,B,Xes[i]);
217: SNESGetConvergedReason(next->snes,&reason);
218: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
219: jac->innerFailures++;
220: if (jac->innerFailures >= snes->maxFailures) {
221: snes->reason = SNES_DIVERGED_INNER;
222: return(0);
223: }
224: }
225: while (next->next) {
226: i++;
227: next = next->next;
228: VecCopy(X,Xes[i]);
229: SNESSolve(next->snes,B,Xes[i]);
230: SNESGetConvergedReason(next->snes,&reason);
231: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
232: jac->innerFailures++;
233: if (jac->innerFailures >= snes->maxFailures) {
234: snes->reason = SNES_DIVERGED_INNER;
235: return(0);
236: }
237: }
238: }
240: /* all the solutions are collected; combine optimally */
241: for (i=0;i<jac->n;i++) {
242: for (j=0;j<i+1;j++) {
243: VecDotBegin(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
244: }
245: VecDotBegin(Fes[i],F,&jac->g[i]);
246: }
248: for (i=0;i<jac->n;i++) {
249: for (j=0;j<i+1;j++) {
250: VecDotEnd(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
251: if (i == j) jac->fnorms[i] = PetscSqrtReal(PetscRealPart(jac->h[i + j*jac->n]));
252: }
253: VecDotEnd(Fes[i],F,&jac->g[i]);
254: }
256: ftf = (*fnorm)*(*fnorm);
258: for (i=0; i<jac->n; i++) {
259: for (j=i+1;j<jac->n;j++) {
260: jac->h[i + j*jac->n] = jac->h[j + i*jac->n];
261: }
262: }
264: for (i=0; i<jac->n; i++) {
265: for (j=0; j<jac->n; j++) {
266: jac->h[i + j*jac->n] = jac->h[i + j*jac->n] - jac->g[j] - jac->g[i] + ftf;
267: }
268: jac->beta[i] = ftf - jac->g[i];
269: }
271: jac->info = 0;
272: jac->rcond = -1.;
273: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
274: #if defined(PETSC_USE_COMPLEX)
275: 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));
276: #else
277: 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));
278: #endif
279: PetscFPTrapPop();
280: if (jac->info < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"Bad argument to GELSS");
281: if (jac->info > 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD failed to converge");
282: tot = 0.;
283: total = 0.;
284: for (i=0; i<jac->n; i++) {
285: if (snes->errorifnotconverged && PetscIsInfOrNanScalar(jac->beta[i])) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD generated inconsistent output");
286: PetscInfo2(snes,"%D: %g\n",i,(double)PetscRealPart(jac->beta[i]));
287: tot += jac->beta[i];
288: total += PetscAbsScalar(jac->beta[i]);
289: }
290: VecScale(X,(1. - tot));
291: VecMAXPY(X,jac->n,jac->beta,Xes);
292: SNESComputeFunction(snes,X,F);
294: if (snes->xl && snes->xu) {
295: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
296: } else {
297: VecNorm(F, NORM_2, fnorm);
298: }
300: /* take the minimum-normed candidate if it beats the combination by a factor of rtol or the combination has stagnated */
301: min_fnorm = jac->fnorms[0];
302: min_i = 0;
303: for (i=0; i<jac->n; i++) {
304: if (jac->fnorms[i] < min_fnorm) {
305: min_fnorm = jac->fnorms[i];
306: min_i = i;
307: }
308: }
310: /* stagnation or divergence restart to the solution of the solver that failed the least */
311: if (PetscRealPart(total) < jac->stol || min_fnorm*jac->rtol < *fnorm) {
312: VecCopy(jac->Xes[min_i],X);
313: VecCopy(jac->Fes[min_i],F);
314: *fnorm = min_fnorm;
315: }
316: return(0);
317: }
319: static PetscErrorCode SNESSetUp_Composite(SNES snes)
320: {
321: PetscErrorCode ierr;
322: DM dm;
323: SNES_Composite *jac = (SNES_Composite*)snes->data;
324: SNES_CompositeLink next = jac->head;
325: PetscInt n=0,i;
326: Vec F;
329: SNESGetDM(snes,&dm);
331: if (snes->ops->computevariablebounds) {
332: /* SNESVI only ever calls computevariablebounds once, so calling it once here is justified */
333: if (!snes->xl) {VecDuplicate(snes->vec_sol,&snes->xl);}
334: if (!snes->xu) {VecDuplicate(snes->vec_sol,&snes->xu);}
335: (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);
336: }
338: while (next) {
339: n++;
340: SNESSetDM(next->snes,dm);
341: SNESSetApplicationContext(next->snes, snes->user);
342: if (snes->xl && snes->xu) {
343: if (snes->ops->computevariablebounds) {
344: SNESVISetComputeVariableBounds(next->snes, snes->ops->computevariablebounds);
345: } else {
346: SNESVISetVariableBounds(next->snes,snes->xl,snes->xu);
347: }
348: }
350: next = next->next;
351: }
352: jac->nsnes = n;
353: SNESGetFunction(snes,&F,NULL,NULL);
354: if (jac->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
355: VecDuplicateVecs(F,jac->nsnes,&jac->Xes);
356: PetscMalloc1(n,&jac->Fes);
357: PetscMalloc1(n,&jac->fnorms);
358: next = jac->head;
359: i = 0;
360: while (next) {
361: SNESGetFunction(next->snes,&F,NULL,NULL);
362: jac->Fes[i] = F;
363: PetscObjectReference((PetscObject)F);
364: next = next->next;
365: i++;
366: }
367: /* allocate the subspace direct solve area */
368: jac->nrhs = 1;
369: jac->lda = jac->nsnes;
370: jac->ldb = jac->nsnes;
371: jac->n = jac->nsnes;
373: PetscMalloc1(jac->n*jac->n,&jac->h);
374: PetscMalloc1(jac->n,&jac->beta);
375: PetscMalloc1(jac->n,&jac->s);
376: PetscMalloc1(jac->n,&jac->g);
377: jac->lwork = 12*jac->n;
378: #if defined(PETSC_USE_COMPLEX)
379: PetscMalloc1(jac->lwork,&jac->rwork);
380: #endif
381: PetscMalloc1(jac->lwork,&jac->work);
382: }
384: return(0);
385: }
387: static PetscErrorCode SNESReset_Composite(SNES snes)
388: {
389: SNES_Composite *jac = (SNES_Composite*)snes->data;
390: PetscErrorCode ierr;
391: SNES_CompositeLink next = jac->head;
394: while (next) {
395: SNESReset(next->snes);
396: next = next->next;
397: }
398: VecDestroy(&jac->Xorig);
399: if (jac->Xes) {VecDestroyVecs(jac->nsnes,&jac->Xes);}
400: if (jac->Fes) {VecDestroyVecs(jac->nsnes,&jac->Fes);}
401: PetscFree(jac->fnorms);
402: PetscFree(jac->h);
403: PetscFree(jac->s);
404: PetscFree(jac->g);
405: PetscFree(jac->beta);
406: PetscFree(jac->work);
407: PetscFree(jac->rwork);
408: return(0);
409: }
411: static PetscErrorCode SNESDestroy_Composite(SNES snes)
412: {
413: SNES_Composite *jac = (SNES_Composite*)snes->data;
414: PetscErrorCode ierr;
415: SNES_CompositeLink next = jac->head,next_tmp;
418: SNESReset_Composite(snes);
419: while (next) {
420: SNESDestroy(&next->snes);
421: next_tmp = next;
422: next = next->next;
423: PetscFree(next_tmp);
424: }
425: PetscFree(snes->data);
426: return(0);
427: }
429: static PetscErrorCode SNESSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,SNES snes)
430: {
431: SNES_Composite *jac = (SNES_Composite*)snes->data;
432: PetscErrorCode ierr;
433: PetscInt nmax = 8,i;
434: SNES_CompositeLink next;
435: char *sneses[8];
436: PetscReal dmps[8];
437: PetscBool flg;
440: PetscOptionsHead(PetscOptionsObject,"Composite preconditioner options");
441: PetscOptionsEnum("-snes_composite_type","Type of composition","SNESCompositeSetType",SNESCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
442: if (flg) {
443: SNESCompositeSetType(snes,jac->type);
444: }
445: PetscOptionsStringArray("-snes_composite_sneses","List of composite solvers","SNESCompositeAddSNES",sneses,&nmax,&flg);
446: if (flg) {
447: for (i=0; i<nmax; i++) {
448: SNESCompositeAddSNES(snes,sneses[i]);
449: PetscFree(sneses[i]); /* deallocate string sneses[i], which is allocated in PetscOptionsStringArray() */
450: }
451: }
452: PetscOptionsRealArray("-snes_composite_damping","Damping of the additive composite solvers","SNESCompositeSetDamping",dmps,&nmax,&flg);
453: if (flg) {
454: for (i=0; i<nmax; i++) {
455: SNESCompositeSetDamping(snes,i,dmps[i]);
456: }
457: }
458: PetscOptionsReal("-snes_composite_stol","Step tolerance for restart on the additive composite solvers","",jac->stol,&jac->stol,NULL);
459: PetscOptionsReal("-snes_composite_rtol","Residual tolerance for the additive composite solvers","",jac->rtol,&jac->rtol,NULL);
460: PetscOptionsTail();
462: next = jac->head;
463: while (next) {
464: SNESSetFromOptions(next->snes);
465: next = next->next;
466: }
467: return(0);
468: }
470: static PetscErrorCode SNESView_Composite(SNES snes,PetscViewer viewer)
471: {
472: SNES_Composite *jac = (SNES_Composite*)snes->data;
473: PetscErrorCode ierr;
474: SNES_CompositeLink next = jac->head;
475: PetscBool iascii;
478: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
479: if (iascii) {
480: PetscViewerASCIIPrintf(viewer," type - %s\n",SNESCompositeTypes[jac->type]);
481: PetscViewerASCIIPrintf(viewer," SNESes on composite preconditioner follow\n");
482: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
483: }
484: if (iascii) {
485: PetscViewerASCIIPushTab(viewer);
486: }
487: while (next) {
488: SNESView(next->snes,viewer);
489: next = next->next;
490: }
491: if (iascii) {
492: PetscViewerASCIIPopTab(viewer);
493: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
494: }
495: return(0);
496: }
498: /* ------------------------------------------------------------------------------*/
500: static PetscErrorCode SNESCompositeSetType_Composite(SNES snes,SNESCompositeType type)
501: {
502: SNES_Composite *jac = (SNES_Composite*)snes->data;
505: jac->type = type;
506: return(0);
507: }
509: static PetscErrorCode SNESCompositeAddSNES_Composite(SNES snes,SNESType type)
510: {
511: SNES_Composite *jac;
512: SNES_CompositeLink next,ilink;
513: PetscErrorCode ierr;
514: PetscInt cnt = 0;
515: const char *prefix;
516: char newprefix[20];
517: DM dm;
520: PetscNewLog(snes,&ilink);
521: ilink->next = NULL;
522: SNESCreate(PetscObjectComm((PetscObject)snes),&ilink->snes);
523: PetscLogObjectParent((PetscObject)snes,(PetscObject)ilink->snes);
524: SNESGetDM(snes,&dm);
525: SNESSetDM(ilink->snes,dm);
526: SNESSetTolerances(ilink->snes,snes->abstol,snes->rtol,snes->stol,1,snes->max_funcs);
527: PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)ilink->snes);
528: jac = (SNES_Composite*)snes->data;
529: next = jac->head;
530: if (!next) {
531: jac->head = ilink;
532: ilink->previous = NULL;
533: } else {
534: cnt++;
535: while (next->next) {
536: next = next->next;
537: cnt++;
538: }
539: next->next = ilink;
540: ilink->previous = next;
541: }
542: SNESGetOptionsPrefix(snes,&prefix);
543: SNESSetOptionsPrefix(ilink->snes,prefix);
544: PetscSNPrintf(newprefix,sizeof(newprefix),"sub_%d_",(int)cnt);
545: SNESAppendOptionsPrefix(ilink->snes,newprefix);
546: PetscObjectIncrementTabLevel((PetscObject)ilink->snes,(PetscObject)snes,1);
547: SNESSetType(ilink->snes,type);
548: SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY);
550: ilink->dmp = 1.0;
551: jac->nsnes++;
552: return(0);
553: }
555: static PetscErrorCode SNESCompositeGetSNES_Composite(SNES snes,PetscInt n,SNES *subsnes)
556: {
557: SNES_Composite *jac;
558: SNES_CompositeLink next;
559: PetscInt i;
562: jac = (SNES_Composite*)snes->data;
563: next = jac->head;
564: for (i=0; i<n; i++) {
565: if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
566: next = next->next;
567: }
568: *subsnes = next->snes;
569: return(0);
570: }
572: /* -------------------------------------------------------------------------------- */
573: /*@C
574: SNESCompositeSetType - Sets the type of composite preconditioner.
576: Logically Collective on SNES
578: Input Parameters:
579: + snes - the preconditioner context
580: - type - SNES_COMPOSITE_ADDITIVE (default), SNES_COMPOSITE_MULTIPLICATIVE
582: Options Database Key:
583: . -snes_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
585: Level: Developer
587: @*/
588: PetscErrorCode SNESCompositeSetType(SNES snes,SNESCompositeType type)
589: {
595: PetscTryMethod(snes,"SNESCompositeSetType_C",(SNES,SNESCompositeType),(snes,type));
596: return(0);
597: }
599: /*@C
600: SNESCompositeAddSNES - Adds another SNES to the composite SNES.
602: Collective on SNES
604: Input Parameters:
605: + snes - the preconditioner context
606: - type - the type of the new preconditioner
608: Level: Developer
610: @*/
611: PetscErrorCode SNESCompositeAddSNES(SNES snes,SNESType type)
612: {
617: PetscTryMethod(snes,"SNESCompositeAddSNES_C",(SNES,SNESType),(snes,type));
618: return(0);
619: }
621: /*@
622: SNESCompositeGetSNES - Gets one of the SNES objects in the composite SNES.
624: Not Collective
626: Input Parameters:
627: + snes - the preconditioner context
628: - n - the number of the snes requested
630: Output Parameters:
631: . subsnes - the SNES requested
633: Level: Developer
635: .seealso: SNESCompositeAddSNES()
636: @*/
637: PetscErrorCode SNESCompositeGetSNES(SNES snes,PetscInt n,SNES *subsnes)
638: {
644: PetscUseMethod(snes,"SNESCompositeGetSNES_C",(SNES,PetscInt,SNES*),(snes,n,subsnes));
645: return(0);
646: }
648: /*@
649: SNESCompositeGetNumber - Get the number of subsolvers in the composite SNES.
651: Logically Collective on SNES
653: Input Parameter:
654: snes - the preconditioner context
656: Output Parameter:
657: n - the number of subsolvers
659: Level: Developer
661: @*/
662: PetscErrorCode SNESCompositeGetNumber(SNES snes,PetscInt *n)
663: {
664: SNES_Composite *jac;
665: SNES_CompositeLink next;
668: jac = (SNES_Composite*)snes->data;
669: next = jac->head;
671: *n = 0;
672: while (next) {
673: *n = *n + 1;
674: next = next->next;
675: }
676: return(0);
677: }
679: static PetscErrorCode SNESCompositeSetDamping_Composite(SNES snes,PetscInt n,PetscReal dmp)
680: {
681: SNES_Composite *jac;
682: SNES_CompositeLink next;
683: PetscInt i;
686: jac = (SNES_Composite*)snes->data;
687: next = jac->head;
688: for (i=0; i<n; i++) {
689: if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
690: next = next->next;
691: }
692: next->dmp = dmp;
693: return(0);
694: }
696: /*@
697: SNESCompositeSetDamping - Sets the damping of a subsolver when using additive composite SNES.
699: Not Collective
701: Input Parameters:
702: + snes - the preconditioner context
703: . n - the number of the snes requested
704: - dmp - the damping
706: Level: Developer
708: .seealso: SNESCompositeAddSNES()
709: @*/
710: PetscErrorCode SNESCompositeSetDamping(SNES snes,PetscInt n,PetscReal dmp)
711: {
716: PetscUseMethod(snes,"SNESCompositeSetDamping_C",(SNES,PetscInt,PetscReal),(snes,n,dmp));
717: return(0);
718: }
720: static PetscErrorCode SNESSolve_Composite(SNES snes)
721: {
722: Vec F,X,B,Y;
723: PetscInt i;
724: PetscReal fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
725: PetscErrorCode ierr;
726: SNESNormSchedule normtype;
727: SNES_Composite *comp = (SNES_Composite*)snes->data;
730: X = snes->vec_sol;
731: F = snes->vec_func;
732: B = snes->vec_rhs;
733: Y = snes->vec_sol_update;
735: PetscObjectSAWsTakeAccess((PetscObject)snes);
736: snes->iter = 0;
737: snes->norm = 0.;
738: comp->innerFailures = 0;
739: PetscObjectSAWsGrantAccess((PetscObject)snes);
740: snes->reason = SNES_CONVERGED_ITERATING;
741: SNESGetNormSchedule(snes, &normtype);
742: if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
743: if (!snes->vec_func_init_set) {
744: SNESComputeFunction(snes,X,F);
745: } else snes->vec_func_init_set = PETSC_FALSE;
747: if (snes->xl && snes->xu) {
748: SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
749: } else {
750: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
751: }
752: SNESCheckFunctionNorm(snes,fnorm);
753: PetscObjectSAWsTakeAccess((PetscObject)snes);
754: snes->iter = 0;
755: snes->norm = fnorm;
756: PetscObjectSAWsGrantAccess((PetscObject)snes);
757: SNESLogConvergenceHistory(snes,snes->norm,0);
758: SNESMonitor(snes,0,snes->norm);
760: /* test convergence */
761: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
762: if (snes->reason) return(0);
763: } else {
764: PetscObjectSAWsGrantAccess((PetscObject)snes);
765: SNESLogConvergenceHistory(snes,snes->norm,0);
766: SNESMonitor(snes,0,snes->norm);
767: }
769: for (i = 0; i < snes->max_its; i++) {
770: /* Call general purpose update function */
771: if (snes->ops->update) {
772: (*snes->ops->update)(snes, snes->iter);
773: }
775: /* Copy the state before modification by application of the composite solver;
776: we will subtract the new state after application */
777: VecCopy(X, Y);
779: if (comp->type == SNES_COMPOSITE_ADDITIVE) {
780: SNESCompositeApply_Additive(snes,X,B,F,&fnorm);
781: } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
782: SNESCompositeApply_Multiplicative(snes,X,B,F,&fnorm);
783: } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
784: SNESCompositeApply_AdditiveOptimal(snes,X,B,F,&fnorm);
785: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported SNESComposite type");
786: if (snes->reason < 0) break;
788: /* Compute the solution update for convergence testing */
789: VecAYPX(Y, -1.0, X);
791: if ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY)) {
792: SNESComputeFunction(snes,X,F);
794: if (snes->xl && snes->xu) {
795: VecNormBegin(X, NORM_2, &xnorm);
796: VecNormBegin(Y, NORM_2, &snorm);
797: SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
798: VecNormEnd(X, NORM_2, &xnorm);
799: VecNormEnd(Y, NORM_2, &snorm);
800: } else {
801: VecNormBegin(F, NORM_2, &fnorm);
802: VecNormBegin(X, NORM_2, &xnorm);
803: VecNormBegin(Y, NORM_2, &snorm);
805: VecNormEnd(F, NORM_2, &fnorm);
806: VecNormEnd(X, NORM_2, &xnorm);
807: VecNormEnd(Y, NORM_2, &snorm);
808: }
809: SNESCheckFunctionNorm(snes,fnorm);
810: } else if (normtype == SNES_NORM_ALWAYS) {
811: VecNormBegin(X, NORM_2, &xnorm);
812: VecNormBegin(Y, NORM_2, &snorm);
813: VecNormEnd(X, NORM_2, &xnorm);
814: VecNormEnd(Y, NORM_2, &snorm);
815: }
816: /* Monitor convergence */
817: PetscObjectSAWsTakeAccess((PetscObject)snes);
818: snes->iter = i+1;
819: snes->norm = fnorm;
820: snes->xnorm = xnorm;
821: snes->ynorm = snorm;
822: PetscObjectSAWsGrantAccess((PetscObject)snes);
823: SNESLogConvergenceHistory(snes,snes->norm,0);
824: SNESMonitor(snes,snes->iter,snes->norm);
825: /* Test for convergence */
826: if (normtype == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,xnorm,snorm,fnorm,&snes->reason,snes->cnvP);}
827: if (snes->reason) break;
828: }
829: if (normtype == SNES_NORM_ALWAYS) {
830: if (i == snes->max_its) {
831: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
832: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
833: }
834: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
835: return(0);
836: }
838: /* -------------------------------------------------------------------------------------------*/
840: /*MC
841: SNESCOMPOSITE - Build a preconditioner by composing together several nonlinear solvers
843: Options Database Keys:
844: + -snes_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
845: - -snes_composite_sneses - <snes0,snes1,...> list of SNESes to compose
847: Level: intermediate
849: .seealso: SNESCreate(), SNESSetType(), SNESType (for list of available types), SNES,
850: SNESSHELL, SNESCompositeSetType(), SNESCompositeSpecialSetAlpha(), SNESCompositeAddSNES(),
851: SNESCompositeGetSNES()
853: References:
854: . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
855: SIAM Review, 57(4), 2015
857: M*/
859: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
860: {
862: SNES_Composite *jac;
865: PetscNewLog(snes,&jac);
867: snes->ops->solve = SNESSolve_Composite;
868: snes->ops->setup = SNESSetUp_Composite;
869: snes->ops->reset = SNESReset_Composite;
870: snes->ops->destroy = SNESDestroy_Composite;
871: snes->ops->setfromoptions = SNESSetFromOptions_Composite;
872: snes->ops->view = SNESView_Composite;
874: snes->usesksp = PETSC_FALSE;
876: snes->alwayscomputesfinalresidual = PETSC_FALSE;
878: snes->data = (void*)jac;
879: jac->type = SNES_COMPOSITE_ADDITIVEOPTIMAL;
880: jac->Fes = NULL;
881: jac->Xes = NULL;
882: jac->fnorms = NULL;
883: jac->nsnes = 0;
884: jac->head = NULL;
885: jac->stol = 0.1;
886: jac->rtol = 1.1;
888: jac->h = NULL;
889: jac->s = NULL;
890: jac->beta = NULL;
891: jac->work = NULL;
892: jac->rwork = NULL;
894: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetType_C",SNESCompositeSetType_Composite);
895: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeAddSNES_C",SNESCompositeAddSNES_Composite);
896: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeGetSNES_C",SNESCompositeGetSNES_Composite);
897: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetDamping_C",SNESCompositeSetDamping_Composite);
898: return(0);
899: }