Actual source code: snescomposite.c
petsc-3.8.4 2018-03-24
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",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;
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: #if defined(PETSC_MISSING_LAPACK_GELSS)
272: SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"SNESCOMPOSITE with ADDITIVEOPTIMAL requires the LAPACK GELSS routine.");
273: #else
274: jac->info = 0;
275: jac->rcond = -1.;
276: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
277: #if defined(PETSC_USE_COMPLEX)
278: 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));
279: #else
280: 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));
281: #endif
282: PetscFPTrapPop();
283: if (jac->info < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"Bad argument to GELSS");
284: if (jac->info > 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD failed to converge");
285: #endif
286: tot = 0.;
287: total = 0.;
288: for (i=0; i<jac->n; i++) {
289: if (snes->errorifnotconverged && PetscIsInfOrNanScalar(jac->beta[i])) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD generated inconsistent output");
290: PetscInfo2(snes,"%D: %g\n",i,(double)PetscRealPart(jac->beta[i]));
291: tot += jac->beta[i];
292: total += PetscAbsScalar(jac->beta[i]);
293: }
294: VecScale(X,(1. - tot));
295: VecMAXPY(X,jac->n,jac->beta,Xes);
296: SNESComputeFunction(snes,X,F);
298: if (snes->xl && snes->xu) {
299: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
300: } else {
301: VecNorm(F, NORM_2, fnorm);
302: }
304: /* take the minimum-normed candidate if it beats the combination by a factor of rtol or the combination has stagnated */
305: min_fnorm = jac->fnorms[0];
306: min_i = 0;
307: for (i=0; i<jac->n; i++) {
308: if (jac->fnorms[i] < min_fnorm) {
309: min_fnorm = jac->fnorms[i];
310: min_i = i;
311: }
312: }
314: /* stagnation or divergence restart to the solution of the solver that failed the least */
315: if (PetscRealPart(total) < jac->stol || min_fnorm*jac->rtol < *fnorm) {
316: VecCopy(jac->Xes[min_i],X);
317: VecCopy(jac->Fes[min_i],F);
318: *fnorm = min_fnorm;
319: }
320: return(0);
321: }
323: static PetscErrorCode SNESSetUp_Composite(SNES snes)
324: {
325: PetscErrorCode ierr;
326: DM dm;
327: SNES_Composite *jac = (SNES_Composite*)snes->data;
328: SNES_CompositeLink next = jac->head;
329: PetscInt n=0,i;
330: Vec F;
333: SNESGetDM(snes,&dm);
335: if (snes->ops->computevariablebounds) {
336: /* SNESVI only ever calls computevariablebounds once, so calling it once here is justified */
337: if (!snes->xl) {VecDuplicate(snes->vec_sol,&snes->xl);}
338: if (!snes->xu) {VecDuplicate(snes->vec_sol,&snes->xu);}
339: (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);
340: }
342: while (next) {
343: n++;
344: SNESSetDM(next->snes,dm);
345: SNESSetApplicationContext(next->snes, snes->user);
346: if (snes->xl && snes->xu) {
347: if (snes->ops->computevariablebounds) {
348: SNESVISetComputeVariableBounds(next->snes, snes->ops->computevariablebounds);
349: } else {
350: SNESVISetVariableBounds(next->snes,snes->xl,snes->xu);
351: }
352: }
354: next = next->next;
355: }
356: jac->nsnes = n;
357: SNESGetFunction(snes,&F,NULL,NULL);
358: if (jac->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
359: VecDuplicateVecs(F,jac->nsnes,&jac->Xes);
360: PetscMalloc1(n,&jac->Fes);
361: PetscMalloc1(n,&jac->fnorms);
362: next = jac->head;
363: i = 0;
364: while (next) {
365: SNESGetFunction(next->snes,&F,NULL,NULL);
366: jac->Fes[i] = F;
367: PetscObjectReference((PetscObject)F);
368: next = next->next;
369: i++;
370: }
371: /* allocate the subspace direct solve area */
372: jac->nrhs = 1;
373: jac->lda = jac->nsnes;
374: jac->ldb = jac->nsnes;
375: jac->n = jac->nsnes;
377: PetscMalloc1(jac->n*jac->n,&jac->h);
378: PetscMalloc1(jac->n,&jac->beta);
379: PetscMalloc1(jac->n,&jac->s);
380: PetscMalloc1(jac->n,&jac->g);
381: jac->lwork = 12*jac->n;
382: #if defined(PETSC_USE_COMPLEX)
383: PetscMalloc1(jac->lwork,&jac->rwork);
384: #endif
385: PetscMalloc1(jac->lwork,&jac->work);
386: }
388: return(0);
389: }
391: static PetscErrorCode SNESReset_Composite(SNES snes)
392: {
393: SNES_Composite *jac = (SNES_Composite*)snes->data;
394: PetscErrorCode ierr;
395: SNES_CompositeLink next = jac->head;
398: while (next) {
399: SNESReset(next->snes);
400: next = next->next;
401: }
402: VecDestroy(&jac->Xorig);
403: if (jac->Xes) {VecDestroyVecs(jac->nsnes,&jac->Xes);}
404: if (jac->Fes) {VecDestroyVecs(jac->nsnes,&jac->Fes);}
405: PetscFree(jac->fnorms);
406: PetscFree(jac->h);
407: PetscFree(jac->s);
408: PetscFree(jac->g);
409: PetscFree(jac->beta);
410: PetscFree(jac->work);
411: PetscFree(jac->rwork);
412: return(0);
413: }
415: static PetscErrorCode SNESDestroy_Composite(SNES snes)
416: {
417: SNES_Composite *jac = (SNES_Composite*)snes->data;
418: PetscErrorCode ierr;
419: SNES_CompositeLink next = jac->head,next_tmp;
422: SNESReset_Composite(snes);
423: while (next) {
424: SNESDestroy(&next->snes);
425: next_tmp = next;
426: next = next->next;
427: PetscFree(next_tmp);
428: }
429: PetscFree(snes->data);
430: return(0);
431: }
433: static PetscErrorCode SNESSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,SNES snes)
434: {
435: SNES_Composite *jac = (SNES_Composite*)snes->data;
436: PetscErrorCode ierr;
437: PetscInt nmax = 8,i;
438: SNES_CompositeLink next;
439: char *sneses[8];
440: PetscReal dmps[8];
441: PetscBool flg;
444: PetscOptionsHead(PetscOptionsObject,"Composite preconditioner options");
445: PetscOptionsEnum("-snes_composite_type","Type of composition","SNESCompositeSetType",SNESCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
446: if (flg) {
447: SNESCompositeSetType(snes,jac->type);
448: }
449: PetscOptionsStringArray("-snes_composite_sneses","List of composite solvers","SNESCompositeAddSNES",sneses,&nmax,&flg);
450: if (flg) {
451: for (i=0; i<nmax; i++) {
452: SNESCompositeAddSNES(snes,sneses[i]);
453: PetscFree(sneses[i]); /* deallocate string sneses[i], which is allocated in PetscOptionsStringArray() */
454: }
455: }
456: PetscOptionsRealArray("-snes_composite_damping","Damping of the additive composite solvers","SNESCompositeSetDamping",dmps,&nmax,&flg);
457: if (flg) {
458: for (i=0; i<nmax; i++) {
459: SNESCompositeSetDamping(snes,i,dmps[i]);
460: }
461: }
462: PetscOptionsReal("-snes_composite_stol","Step tolerance for restart on the additive composite solvers","",jac->stol,&jac->stol,NULL);
463: PetscOptionsReal("-snes_composite_rtol","Residual tolerance for the additive composite solvers","",jac->rtol,&jac->rtol,NULL);
464: PetscOptionsTail();
466: next = jac->head;
467: while (next) {
468: SNESSetFromOptions(next->snes);
469: next = next->next;
470: }
471: return(0);
472: }
474: static PetscErrorCode SNESView_Composite(SNES snes,PetscViewer viewer)
475: {
476: SNES_Composite *jac = (SNES_Composite*)snes->data;
477: PetscErrorCode ierr;
478: SNES_CompositeLink next = jac->head;
479: PetscBool iascii;
482: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
483: if (iascii) {
484: PetscViewerASCIIPrintf(viewer," type - %s\n",SNESCompositeTypes[jac->type]);
485: PetscViewerASCIIPrintf(viewer," SNESes on composite preconditioner follow\n");
486: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
487: }
488: if (iascii) {
489: PetscViewerASCIIPushTab(viewer);
490: }
491: while (next) {
492: SNESView(next->snes,viewer);
493: next = next->next;
494: }
495: if (iascii) {
496: PetscViewerASCIIPopTab(viewer);
497: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
498: }
499: return(0);
500: }
502: /* ------------------------------------------------------------------------------*/
504: static PetscErrorCode SNESCompositeSetType_Composite(SNES snes,SNESCompositeType type)
505: {
506: SNES_Composite *jac = (SNES_Composite*)snes->data;
509: jac->type = type;
510: return(0);
511: }
513: static PetscErrorCode SNESCompositeAddSNES_Composite(SNES snes,SNESType type)
514: {
515: SNES_Composite *jac;
516: SNES_CompositeLink next,ilink;
517: PetscErrorCode ierr;
518: PetscInt cnt = 0;
519: const char *prefix;
520: char newprefix[20];
521: DM dm;
524: PetscNewLog(snes,&ilink);
525: ilink->next = 0;
526: SNESCreate(PetscObjectComm((PetscObject)snes),&ilink->snes);
527: PetscLogObjectParent((PetscObject)snes,(PetscObject)ilink->snes);
528: SNESGetDM(snes,&dm);
529: SNESSetDM(ilink->snes,dm);
530: SNESSetTolerances(ilink->snes,ilink->snes->abstol,ilink->snes->rtol,ilink->snes->stol,1,ilink->snes->max_funcs);
531: PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)ilink->snes);
532: jac = (SNES_Composite*)snes->data;
533: next = jac->head;
534: if (!next) {
535: jac->head = ilink;
536: ilink->previous = NULL;
537: } else {
538: cnt++;
539: while (next->next) {
540: next = next->next;
541: cnt++;
542: }
543: next->next = ilink;
544: ilink->previous = next;
545: }
546: SNESGetOptionsPrefix(snes,&prefix);
547: SNESSetOptionsPrefix(ilink->snes,prefix);
548: sprintf(newprefix,"sub_%d_",(int)cnt);
549: SNESAppendOptionsPrefix(ilink->snes,newprefix);
550: PetscObjectIncrementTabLevel((PetscObject)ilink->snes,(PetscObject)snes,1);
551: SNESSetType(ilink->snes,type);
552: SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY);
554: ilink->dmp = 1.0;
555: jac->nsnes++;
556: return(0);
557: }
559: static PetscErrorCode SNESCompositeGetSNES_Composite(SNES snes,PetscInt n,SNES *subsnes)
560: {
561: SNES_Composite *jac;
562: SNES_CompositeLink next;
563: PetscInt i;
566: jac = (SNES_Composite*)snes->data;
567: next = jac->head;
568: for (i=0; i<n; i++) {
569: if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
570: next = next->next;
571: }
572: *subsnes = next->snes;
573: return(0);
574: }
576: /* -------------------------------------------------------------------------------- */
577: /*@C
578: SNESCompositeSetType - Sets the type of composite preconditioner.
580: Logically Collective on SNES
582: Input Parameter:
583: + snes - the preconditioner context
584: - type - SNES_COMPOSITE_ADDITIVE (default), SNES_COMPOSITE_MULTIPLICATIVE
586: Options Database Key:
587: . -snes_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
589: Level: Developer
591: .keywords: SNES, set, type, composite preconditioner, additive, multiplicative
592: @*/
593: PetscErrorCode SNESCompositeSetType(SNES snes,SNESCompositeType type)
594: {
600: PetscTryMethod(snes,"SNESCompositeSetType_C",(SNES,SNESCompositeType),(snes,type));
601: return(0);
602: }
604: /*@C
605: SNESCompositeAddSNES - Adds another SNES to the composite SNES.
607: Collective on SNES
609: Input Parameters:
610: + snes - the preconditioner context
611: - type - the type of the new preconditioner
613: Level: Developer
615: .keywords: SNES, composite preconditioner, add
616: @*/
617: PetscErrorCode SNESCompositeAddSNES(SNES snes,SNESType type)
618: {
623: PetscTryMethod(snes,"SNESCompositeAddSNES_C",(SNES,SNESType),(snes,type));
624: return(0);
625: }
626: /*@
627: SNESCompositeGetSNES - Gets one of the SNES objects in the composite SNES.
629: Not Collective
631: Input Parameter:
632: + snes - the preconditioner context
633: - n - the number of the snes requested
635: Output Parameters:
636: . subsnes - the SNES requested
638: Level: Developer
640: .keywords: SNES, get, composite preconditioner, sub preconditioner
642: .seealso: SNESCompositeAddSNES()
643: @*/
644: PetscErrorCode SNESCompositeGetSNES(SNES snes,PetscInt n,SNES *subsnes)
645: {
651: PetscUseMethod(snes,"SNESCompositeGetSNES_C",(SNES,PetscInt,SNES*),(snes,n,subsnes));
652: return(0);
653: }
655: /*@
656: SNESCompositeGetNumber - Get the number of subsolvers in the composite SNES.
658: Logically Collective on SNES
660: Input Parameter:
661: snes - the preconditioner context
663: Output Parameter:
664: n - the number of subsolvers
666: Level: Developer
668: .keywords: SNES, composite preconditioner
669: @*/
670: PetscErrorCode SNESCompositeGetNumber(SNES snes,PetscInt *n)
671: {
672: SNES_Composite *jac;
673: SNES_CompositeLink next;
676: jac = (SNES_Composite*)snes->data;
677: next = jac->head;
679: *n = 0;
680: while (next) {
681: *n = *n + 1;
682: next = next->next;
683: }
684: return(0);
685: }
687: static PetscErrorCode SNESCompositeSetDamping_Composite(SNES snes,PetscInt n,PetscReal dmp)
688: {
689: SNES_Composite *jac;
690: SNES_CompositeLink next;
691: PetscInt i;
694: jac = (SNES_Composite*)snes->data;
695: next = jac->head;
696: for (i=0; i<n; i++) {
697: if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
698: next = next->next;
699: }
700: next->dmp = dmp;
701: return(0);
702: }
704: /*@
705: SNESCompositeSetDamping - Sets the damping of a subsolver when using additive composite SNES.
707: Not Collective
709: Input Parameter:
710: + snes - the preconditioner context
711: . n - the number of the snes requested
712: - dmp - the damping
714: Level: Developer
716: .keywords: SNES, get, composite preconditioner, sub preconditioner
718: .seealso: SNESCompositeAddSNES()
719: @*/
720: PetscErrorCode SNESCompositeSetDamping(SNES snes,PetscInt n,PetscReal dmp)
721: {
726: PetscUseMethod(snes,"SNESCompositeSetDamping_C",(SNES,PetscInt,PetscReal),(snes,n,dmp));
727: return(0);
728: }
730: static PetscErrorCode SNESSolve_Composite(SNES snes)
731: {
732: Vec F;
733: Vec X;
734: Vec B;
735: PetscInt i;
736: PetscReal fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
737: PetscErrorCode ierr;
738: SNESNormSchedule normtype;
739: SNES_Composite *comp = (SNES_Composite*)snes->data;
742: X = snes->vec_sol;
743: F = snes->vec_func;
744: B = snes->vec_rhs;
746: PetscObjectSAWsTakeAccess((PetscObject)snes);
747: snes->iter = 0;
748: snes->norm = 0.;
749: comp->innerFailures = 0;
750: PetscObjectSAWsGrantAccess((PetscObject)snes);
751: SNESSetWorkVecs(snes, 1);
752: snes->reason = SNES_CONVERGED_ITERATING;
753: SNESGetNormSchedule(snes, &normtype);
754: if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
755: if (!snes->vec_func_init_set) {
756: SNESComputeFunction(snes,X,F);
757: } else snes->vec_func_init_set = PETSC_FALSE;
759: if (snes->xl && snes->xu) {
760: SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
761: } else {
762: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
763: }
764: SNESCheckFunctionNorm(snes,fnorm);
765: PetscObjectSAWsTakeAccess((PetscObject)snes);
766: snes->iter = 0;
767: snes->norm = fnorm;
768: PetscObjectSAWsGrantAccess((PetscObject)snes);
769: SNESLogConvergenceHistory(snes,snes->norm,0);
770: SNESMonitor(snes,0,snes->norm);
772: /* test convergence */
773: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
774: if (snes->reason) return(0);
775: } else {
776: PetscObjectSAWsGrantAccess((PetscObject)snes);
777: SNESLogConvergenceHistory(snes,snes->norm,0);
778: SNESMonitor(snes,0,snes->norm);
779: }
781: /* Call general purpose update function */
782: if (snes->ops->update) {
783: (*snes->ops->update)(snes, snes->iter);
784: }
786: for (i = 0; i < snes->max_its; i++) {
787: /* Copy the state before modification by application of the composite solver;
788: we will subtract the new state after application */
789: VecCopy(X, snes->work[0]);
791: if (comp->type == SNES_COMPOSITE_ADDITIVE) {
792: SNESCompositeApply_Additive(snes,X,B,F,&fnorm);
793: } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
794: SNESCompositeApply_Multiplicative(snes,X,B,F,&fnorm);
795: } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
796: SNESCompositeApply_AdditiveOptimal(snes,X,B,F,&fnorm);
797: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported SNESComposite type");
798: if (snes->reason < 0) break;
800: /* Compute the solution update for convergence testing */
801: VecAXPY(snes->work[0], -1.0, X);
802: VecScale(snes->work[0], -1.0);
804: if ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY)) {
805: SNESComputeFunction(snes,X,F);
807: if (snes->xl && snes->xu) {
808: VecNormBegin(X, NORM_2, &xnorm);
809: VecNormBegin(snes->work[0], NORM_2, &snorm);
810: SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
811: VecNormEnd(X, NORM_2, &xnorm);
812: VecNormEnd(snes->work[0], NORM_2, &snorm);
813: } else {
814: VecNormBegin(F, NORM_2, &fnorm);
815: VecNormBegin(X, NORM_2, &xnorm);
816: VecNormBegin(snes->work[0], NORM_2, &snorm);
818: VecNormEnd(F, NORM_2, &fnorm);
819: VecNormEnd(X, NORM_2, &xnorm);
820: VecNormEnd(snes->work[0], NORM_2, &snorm);
821: }
822: SNESCheckFunctionNorm(snes,fnorm);
823: } else if (normtype == SNES_NORM_ALWAYS) {
824: VecNormBegin(X, NORM_2, &xnorm);
825: VecNormBegin(snes->work[0], NORM_2, &snorm);
826: VecNormEnd(X, NORM_2, &xnorm);
827: VecNormEnd(snes->work[0], NORM_2, &snorm);
828: }
829: /* Monitor convergence */
830: PetscObjectSAWsTakeAccess((PetscObject)snes);
831: snes->iter = i+1;
832: snes->norm = fnorm;
833: PetscObjectSAWsGrantAccess((PetscObject)snes);
834: SNESLogConvergenceHistory(snes,snes->norm,0);
835: SNESMonitor(snes,snes->iter,snes->norm);
836: /* Test for convergence */
837: if (normtype == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,xnorm,snorm,fnorm,&snes->reason,snes->cnvP);}
838: if (snes->reason) break;
839: /* Call general purpose update function */
840: if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
841: }
842: if (normtype == SNES_NORM_ALWAYS) {
843: if (i == snes->max_its) {
844: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
845: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
846: }
847: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
848: return(0);
849: }
851: /* -------------------------------------------------------------------------------------------*/
853: /*MC
854: SNESCOMPOSITE - Build a preconditioner by composing together several nonlinear solvers
856: Options Database Keys:
857: + -snes_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
858: - -snes_composite_sneses - <snes0,snes1,...> list of SNESes to compose
860: Level: intermediate
862: Concepts: composing solvers
864: .seealso: SNESCreate(), SNESSetType(), SNESType (for list of available types), SNES,
865: SNESSHELL, SNESCompositeSetType(), SNESCompositeSpecialSetAlpha(), SNESCompositeAddSNES(),
866: SNESCompositeGetSNES()
868: References:
869: . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
870: SIAM Review, 57(4), 2015
872: M*/
874: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
875: {
877: SNES_Composite *jac;
880: PetscNewLog(snes,&jac);
882: snes->ops->solve = SNESSolve_Composite;
883: snes->ops->setup = SNESSetUp_Composite;
884: snes->ops->reset = SNESReset_Composite;
885: snes->ops->destroy = SNESDestroy_Composite;
886: snes->ops->setfromoptions = SNESSetFromOptions_Composite;
887: snes->ops->view = SNESView_Composite;
889: snes->alwayscomputesfinalresidual = PETSC_FALSE;
891: snes->data = (void*)jac;
892: jac->type = SNES_COMPOSITE_ADDITIVEOPTIMAL;
893: jac->Fes = NULL;
894: jac->Xes = NULL;
895: jac->fnorms = NULL;
896: jac->nsnes = 0;
897: jac->head = 0;
898: jac->stol = 0.1;
899: jac->rtol = 1.1;
901: jac->h = NULL;
902: jac->s = NULL;
903: jac->beta = NULL;
904: jac->work = NULL;
905: jac->rwork = NULL;
907: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetType_C",SNESCompositeSetType_Composite);
908: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeAddSNES_C",SNESCompositeAddSNES_Composite);
909: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeGetSNES_C",SNESCompositeGetSNES_Composite);
910: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetDamping_C",SNESCompositeSetDamping_Composite);
911: return(0);
912: }