Actual source code: snescomposite.c
petsc-3.7.5 2017-01-01
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;
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;
50: static PetscErrorCode SNESCompositeApply_Multiplicative(SNES snes,Vec X,Vec B,Vec F,PetscReal *fnorm)
51: {
52: PetscErrorCode ierr;
53: SNES_Composite *jac = (SNES_Composite*)snes->data;
54: SNES_CompositeLink next = jac->head;
55: Vec FSub;
56: SNESConvergedReason reason;
59: if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
60: if (snes->normschedule == SNES_NORM_ALWAYS) {
61: SNESSetInitialFunction(next->snes,F);
62: }
63: SNESSolve(next->snes,B,X);
64: SNESGetConvergedReason(next->snes,&reason);
65: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
66: jac->innerFailures++;
67: if (jac->innerFailures >= snes->maxFailures) {
68: snes->reason = SNES_DIVERGED_INNER;
69: return(0);
70: }
71: }
73: while (next->next) {
74: /* only copy the function over in the case where the functions correspond */
75: if (next->snes->pcside == PC_RIGHT && next->snes->normschedule != SNES_NORM_NONE) {
76: SNESGetFunction(next->snes,&FSub,NULL,NULL);
77: next = next->next;
78: SNESSetInitialFunction(next->snes,FSub);
79: } else {
80: next = next->next;
81: }
82: SNESSolve(next->snes,B,X);
83: SNESGetConvergedReason(next->snes,&reason);
84: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
85: jac->innerFailures++;
86: if (jac->innerFailures >= snes->maxFailures) {
87: snes->reason = SNES_DIVERGED_INNER;
88: return(0);
89: }
90: }
91: }
92: if (next->snes->pcside == PC_RIGHT) {
93: SNESGetFunction(next->snes,&FSub,NULL,NULL);
94: VecCopy(FSub,F);
95: if (fnorm) {
96: if (snes->xl && snes->xu) {
97: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
98: } else {
99: VecNorm(F, NORM_2, fnorm);
100: }
101: SNESCheckFunctionNorm(snes,*fnorm);
102: }
103: } else if (snes->normschedule == SNES_NORM_ALWAYS) {
104: SNESComputeFunction(snes,X,F);
105: if (fnorm) {
106: if (snes->xl && snes->xu) {
107: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
108: } else {
109: VecNorm(F, NORM_2, fnorm);
110: }
111: SNESCheckFunctionNorm(snes,*fnorm);
112: }
113: }
114: return(0);
115: }
119: static PetscErrorCode SNESCompositeApply_Additive(SNES snes,Vec X,Vec B,Vec F,PetscReal *fnorm)
120: {
121: PetscErrorCode ierr;
122: SNES_Composite *jac = (SNES_Composite*)snes->data;
123: SNES_CompositeLink next = jac->head;
124: Vec Y,Xorig;
125: SNESConvergedReason reason;
128: Y = snes->vec_sol_update;
129: if (!jac->Xorig) {VecDuplicate(X,&jac->Xorig);}
130: Xorig = jac->Xorig;
131: VecCopy(X,Xorig);
132: if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
133: if (snes->normschedule == SNES_NORM_ALWAYS) {
134: SNESSetInitialFunction(next->snes,F);
135: while (next->next) {
136: next = next->next;
137: SNESSetInitialFunction(next->snes,F);
138: }
139: }
140: next = jac->head;
141: VecCopy(Xorig,Y);
142: SNESSolve(next->snes,B,Y);
143: SNESGetConvergedReason(next->snes,&reason);
144: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
145: jac->innerFailures++;
146: if (jac->innerFailures >= snes->maxFailures) {
147: snes->reason = SNES_DIVERGED_INNER;
148: return(0);
149: }
150: }
151: VecAXPY(Y,-1.0,Xorig);
152: VecAXPY(X,next->dmp,Y);
153: while (next->next) {
154: next = next->next;
155: VecCopy(Xorig,Y);
156: SNESSolve(next->snes,B,Y);
157: SNESGetConvergedReason(next->snes,&reason);
158: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
159: jac->innerFailures++;
160: if (jac->innerFailures >= snes->maxFailures) {
161: snes->reason = SNES_DIVERGED_INNER;
162: return(0);
163: }
164: }
165: VecAXPY(Y,-1.0,Xorig);
166: VecAXPY(X,next->dmp,Y);
167: }
168: if (snes->normschedule == SNES_NORM_ALWAYS) {
169: SNESComputeFunction(snes,X,F);
170: if (fnorm) {
171: if (snes->xl && snes->xu) {
172: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
173: } else {
174: VecNorm(F, NORM_2, fnorm);
175: }
176: SNESCheckFunctionNorm(snes,*fnorm);
177: }
178: }
179: return(0);
180: }
184: /* approximately solve the overdetermined system:
186: 2*F(x_i)\cdot F(\x_j)\alpha_i = 0
187: \alpha_i = 1
189: Which minimizes the L2 norm of the linearization of:
190: ||F(\sum_i \alpha_i*x_i)||^2
192: With the constraint that \sum_i\alpha_i = 1
193: Where x_i is the solution from the ith subsolver.
194: */
195: static PetscErrorCode SNESCompositeApply_AdditiveOptimal(SNES snes,Vec X,Vec B,Vec F,PetscReal *fnorm)
196: {
197: PetscErrorCode ierr;
198: SNES_Composite *jac = (SNES_Composite*)snes->data;
199: SNES_CompositeLink next = jac->head;
200: Vec *Xes = jac->Xes,*Fes = jac->Fes;
201: PetscInt i,j;
202: PetscScalar tot,total,ftf;
203: PetscReal min_fnorm;
204: PetscInt min_i;
205: SNESConvergedReason reason;
208: if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
210: if (snes->normschedule == SNES_NORM_ALWAYS) {
211: next = jac->head;
212: SNESSetInitialFunction(next->snes,F);
213: while (next->next) {
214: next = next->next;
215: SNESSetInitialFunction(next->snes,F);
216: }
217: }
219: next = jac->head;
220: i = 0;
221: VecCopy(X,Xes[i]);
222: SNESSolve(next->snes,B,Xes[i]);
223: SNESGetConvergedReason(next->snes,&reason);
224: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
225: jac->innerFailures++;
226: if (jac->innerFailures >= snes->maxFailures) {
227: snes->reason = SNES_DIVERGED_INNER;
228: return(0);
229: }
230: }
231: while (next->next) {
232: i++;
233: next = next->next;
234: VecCopy(X,Xes[i]);
235: SNESSolve(next->snes,B,Xes[i]);
236: SNESGetConvergedReason(next->snes,&reason);
237: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
238: jac->innerFailures++;
239: if (jac->innerFailures >= snes->maxFailures) {
240: snes->reason = SNES_DIVERGED_INNER;
241: return(0);
242: }
243: }
244: }
246: /* all the solutions are collected; combine optimally */
247: for (i=0;i<jac->n;i++) {
248: for (j=0;j<i+1;j++) {
249: VecDotBegin(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
250: }
251: VecDotBegin(Fes[i],F,&jac->g[i]);
252: }
254: for (i=0;i<jac->n;i++) {
255: for (j=0;j<i+1;j++) {
256: VecDotEnd(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
257: if (i == j) jac->fnorms[i] = PetscSqrtReal(PetscRealPart(jac->h[i + j*jac->n]));
258: }
259: VecDotEnd(Fes[i],F,&jac->g[i]);
260: }
262: ftf = (*fnorm)*(*fnorm);
264: for (i=0; i<jac->n; i++) {
265: for (j=i+1;j<jac->n;j++) {
266: jac->h[i + j*jac->n] = jac->h[j + i*jac->n];
267: }
268: }
270: for (i=0; i<jac->n; i++) {
271: for (j=0; j<jac->n; j++) {
272: jac->h[i + j*jac->n] = jac->h[i + j*jac->n] - jac->g[j] - jac->g[i] + ftf;
273: }
274: jac->beta[i] = ftf - jac->g[i];
275: }
277: #if defined(PETSC_MISSING_LAPACK_GELSS)
278: SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"SNESCOMPOSITE with ADDITIVEOPTIMAL requires the LAPACK GELSS routine.");
279: #else
280: jac->info = 0;
281: jac->rcond = -1.;
282: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
283: #if defined(PETSC_USE_COMPLEX)
284: 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));
285: #else
286: 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));
287: #endif
288: PetscFPTrapPop();
289: if (jac->info < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"Bad argument to GELSS");
290: if (jac->info > 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD failed to converge");
291: #endif
292: tot = 0.;
293: total = 0.;
294: for (i=0; i<jac->n; i++) {
295: if (snes->errorifnotconverged && PetscIsInfOrNanScalar(jac->beta[i])) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD generated inconsistent output");
296: PetscInfo2(snes,"%D: %g\n",i,(double)PetscRealPart(jac->beta[i]));
297: tot += jac->beta[i];
298: total += PetscAbsScalar(jac->beta[i]);
299: }
300: VecScale(X,(1. - tot));
301: VecMAXPY(X,jac->n,jac->beta,Xes);
302: SNESComputeFunction(snes,X,F);
304: if (snes->xl && snes->xu) {
305: SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
306: } else {
307: VecNorm(F, NORM_2, fnorm);
308: }
310: /* take the minimum-normed candidate if it beats the combination by a factor of rtol or the combination has stagnated */
311: min_fnorm = jac->fnorms[0];
312: min_i = 0;
313: for (i=0; i<jac->n; i++) {
314: if (jac->fnorms[i] < min_fnorm) {
315: min_fnorm = jac->fnorms[i];
316: min_i = i;
317: }
318: }
320: /* stagnation or divergence restart to the solution of the solver that failed the least */
321: if (PetscRealPart(total) < jac->stol || min_fnorm*jac->rtol < *fnorm) {
322: VecCopy(jac->Xes[min_i],X);
323: VecCopy(jac->Fes[min_i],F);
324: *fnorm = min_fnorm;
325: }
326: return(0);
327: }
331: static PetscErrorCode SNESSetUp_Composite(SNES snes)
332: {
333: PetscErrorCode ierr;
334: DM dm;
335: SNES_Composite *jac = (SNES_Composite*)snes->data;
336: SNES_CompositeLink next = jac->head;
337: PetscInt n=0,i;
338: Vec F;
341: SNESGetDM(snes,&dm);
343: if (snes->ops->computevariablebounds) {
344: /* SNESVI only ever calls computevariablebounds once, so calling it once here is justified */
345: if (!snes->xl) {VecDuplicate(snes->vec_sol,&snes->xl);}
346: if (!snes->xu) {VecDuplicate(snes->vec_sol,&snes->xu);}
347: (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);
348: }
350: while (next) {
351: n++;
352: SNESSetDM(next->snes,dm);
353: SNESSetApplicationContext(next->snes, snes->user);
354: if (snes->xl && snes->xu) {
355: if (snes->ops->computevariablebounds) {
356: SNESVISetComputeVariableBounds(next->snes, snes->ops->computevariablebounds);
357: } else {
358: SNESVISetVariableBounds(next->snes,snes->xl,snes->xu);
359: }
360: }
362: next = next->next;
363: }
364: jac->nsnes = n;
365: SNESGetFunction(snes,&F,NULL,NULL);
366: if (jac->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
367: VecDuplicateVecs(F,jac->nsnes,&jac->Xes);
368: PetscMalloc1(n,&jac->Fes);
369: PetscMalloc1(n,&jac->fnorms);
370: next = jac->head;
371: i = 0;
372: while (next) {
373: SNESGetFunction(next->snes,&F,NULL,NULL);
374: jac->Fes[i] = F;
375: PetscObjectReference((PetscObject)F);
376: next = next->next;
377: i++;
378: }
379: /* allocate the subspace direct solve area */
380: jac->nrhs = 1;
381: jac->lda = jac->nsnes;
382: jac->ldb = jac->nsnes;
383: jac->n = jac->nsnes;
385: PetscMalloc1(jac->n*jac->n,&jac->h);
386: PetscMalloc1(jac->n,&jac->beta);
387: PetscMalloc1(jac->n,&jac->s);
388: PetscMalloc1(jac->n,&jac->g);
389: jac->lwork = 12*jac->n;
390: #if PETSC_USE_COMPLEX
391: PetscMalloc1(jac->lwork,&jac->rwork);
392: #endif
393: PetscMalloc1(jac->lwork,&jac->work);
394: }
396: return(0);
397: }
401: static PetscErrorCode SNESReset_Composite(SNES snes)
402: {
403: SNES_Composite *jac = (SNES_Composite*)snes->data;
404: PetscErrorCode ierr;
405: SNES_CompositeLink next = jac->head;
408: while (next) {
409: SNESReset(next->snes);
410: next = next->next;
411: }
412: VecDestroy(&jac->Xorig);
413: if (jac->Xes) {VecDestroyVecs(jac->nsnes,&jac->Xes);}
414: if (jac->Fes) {VecDestroyVecs(jac->nsnes,&jac->Fes);}
415: PetscFree(jac->fnorms);
416: PetscFree(jac->h);
417: PetscFree(jac->s);
418: PetscFree(jac->g);
419: PetscFree(jac->beta);
420: PetscFree(jac->work);
421: PetscFree(jac->rwork);
422: return(0);
423: }
427: static PetscErrorCode SNESDestroy_Composite(SNES snes)
428: {
429: SNES_Composite *jac = (SNES_Composite*)snes->data;
430: PetscErrorCode ierr;
431: SNES_CompositeLink next = jac->head,next_tmp;
434: SNESReset_Composite(snes);
435: while (next) {
436: SNESDestroy(&next->snes);
437: next_tmp = next;
438: next = next->next;
439: PetscFree(next_tmp);
440: }
441: PetscFree(snes->data);
442: return(0);
443: }
447: static PetscErrorCode SNESSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,SNES snes)
448: {
449: SNES_Composite *jac = (SNES_Composite*)snes->data;
450: PetscErrorCode ierr;
451: PetscInt nmax = 8,i;
452: SNES_CompositeLink next;
453: char *sneses[8];
454: PetscReal dmps[8];
455: PetscBool flg;
458: PetscOptionsHead(PetscOptionsObject,"Composite preconditioner options");
459: PetscOptionsEnum("-snes_composite_type","Type of composition","SNESCompositeSetType",SNESCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
460: if (flg) {
461: SNESCompositeSetType(snes,jac->type);
462: }
463: PetscOptionsStringArray("-snes_composite_sneses","List of composite solvers","SNESCompositeAddSNES",sneses,&nmax,&flg);
464: if (flg) {
465: for (i=0; i<nmax; i++) {
466: SNESCompositeAddSNES(snes,sneses[i]);
467: PetscFree(sneses[i]); /* deallocate string sneses[i], which is allocated in PetscOptionsStringArray() */
468: }
469: }
470: PetscOptionsRealArray("-snes_composite_damping","Damping of the additive composite solvers","SNESCompositeSetDamping",dmps,&nmax,&flg);
471: if (flg) {
472: for (i=0; i<nmax; i++) {
473: SNESCompositeSetDamping(snes,i,dmps[i]);
474: }
475: }
476: PetscOptionsReal("-snes_composite_stol","Step tolerance for restart on the additive composite solvers","",jac->stol,&jac->stol,NULL);
477: PetscOptionsReal("-snes_composite_rtol","Residual tolerance for the additive composite solvers","",jac->rtol,&jac->rtol,NULL);
478: PetscOptionsTail();
480: next = jac->head;
481: while (next) {
482: SNESSetFromOptions(next->snes);
483: next = next->next;
484: }
485: return(0);
486: }
490: static PetscErrorCode SNESView_Composite(SNES snes,PetscViewer viewer)
491: {
492: SNES_Composite *jac = (SNES_Composite*)snes->data;
493: PetscErrorCode ierr;
494: SNES_CompositeLink next = jac->head;
495: PetscBool iascii;
498: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
499: if (iascii) {
500: PetscViewerASCIIPrintf(viewer,"Composite SNES type - %s\n",SNESCompositeTypes[jac->type]);
501: PetscViewerASCIIPrintf(viewer,"SNESes on composite preconditioner follow\n");
502: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
503: }
504: if (iascii) {
505: PetscViewerASCIIPushTab(viewer);
506: }
507: while (next) {
508: SNESView(next->snes,viewer);
509: next = next->next;
510: }
511: if (iascii) {
512: PetscViewerASCIIPopTab(viewer);
513: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
514: }
515: return(0);
516: }
518: /* ------------------------------------------------------------------------------*/
522: static PetscErrorCode SNESCompositeSetType_Composite(SNES snes,SNESCompositeType type)
523: {
524: SNES_Composite *jac = (SNES_Composite*)snes->data;
527: jac->type = type;
528: return(0);
529: }
533: static PetscErrorCode SNESCompositeAddSNES_Composite(SNES snes,SNESType type)
534: {
535: SNES_Composite *jac;
536: SNES_CompositeLink next,ilink;
537: PetscErrorCode ierr;
538: PetscInt cnt = 0;
539: const char *prefix;
540: char newprefix[8];
541: DM dm;
544: PetscNewLog(snes,&ilink);
545: ilink->next = 0;
546: SNESCreate(PetscObjectComm((PetscObject)snes),&ilink->snes);
547: PetscLogObjectParent((PetscObject)snes,(PetscObject)ilink->snes);
548: SNESGetDM(snes,&dm);
549: SNESSetDM(ilink->snes,dm);
550: SNESSetTolerances(ilink->snes,ilink->snes->abstol,ilink->snes->rtol,ilink->snes->stol,1,ilink->snes->max_funcs);
551: jac = (SNES_Composite*)snes->data;
552: next = jac->head;
553: if (!next) {
554: jac->head = ilink;
555: ilink->previous = NULL;
556: } else {
557: cnt++;
558: while (next->next) {
559: next = next->next;
560: cnt++;
561: }
562: next->next = ilink;
563: ilink->previous = next;
564: }
565: SNESGetOptionsPrefix(snes,&prefix);
566: SNESSetOptionsPrefix(ilink->snes,prefix);
567: sprintf(newprefix,"sub_%d_",(int)cnt);
568: SNESAppendOptionsPrefix(ilink->snes,newprefix);
569: PetscObjectIncrementTabLevel((PetscObject)ilink->snes,(PetscObject)snes,1);
570: SNESSetType(ilink->snes,type);
571: SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY);
573: ilink->dmp = 1.0;
574: jac->nsnes++;
575: return(0);
576: }
580: static PetscErrorCode SNESCompositeGetSNES_Composite(SNES snes,PetscInt n,SNES *subsnes)
581: {
582: SNES_Composite *jac;
583: SNES_CompositeLink next;
584: PetscInt i;
587: jac = (SNES_Composite*)snes->data;
588: next = jac->head;
589: for (i=0; i<n; i++) {
590: if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
591: next = next->next;
592: }
593: *subsnes = next->snes;
594: return(0);
595: }
597: /* -------------------------------------------------------------------------------- */
600: /*@C
601: SNESCompositeSetType - Sets the type of composite preconditioner.
603: Logically Collective on SNES
605: Input Parameter:
606: + snes - the preconditioner context
607: - type - SNES_COMPOSITE_ADDITIVE (default), SNES_COMPOSITE_MULTIPLICATIVE
609: Options Database Key:
610: . -snes_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
612: Level: Developer
614: .keywords: SNES, set, type, composite preconditioner, additive, multiplicative
615: @*/
616: PetscErrorCode SNESCompositeSetType(SNES snes,SNESCompositeType type)
617: {
623: PetscTryMethod(snes,"SNESCompositeSetType_C",(SNES,SNESCompositeType),(snes,type));
624: return(0);
625: }
629: /*@C
630: SNESCompositeAddSNES - Adds another SNES to the composite SNES.
632: Collective on SNES
634: Input Parameters:
635: + snes - the preconditioner context
636: - type - the type of the new preconditioner
638: Level: Developer
640: .keywords: SNES, composite preconditioner, add
641: @*/
642: PetscErrorCode SNESCompositeAddSNES(SNES snes,SNESType type)
643: {
648: PetscTryMethod(snes,"SNESCompositeAddSNES_C",(SNES,SNESType),(snes,type));
649: return(0);
650: }
653: /*@
654: SNESCompositeGetSNES - Gets one of the SNES objects in the composite SNES.
656: Not Collective
658: Input Parameter:
659: + snes - the preconditioner context
660: - n - the number of the snes requested
662: Output Parameters:
663: . subsnes - the SNES requested
665: Level: Developer
667: .keywords: SNES, get, composite preconditioner, sub preconditioner
669: .seealso: SNESCompositeAddSNES()
670: @*/
671: PetscErrorCode SNESCompositeGetSNES(SNES snes,PetscInt n,SNES *subsnes)
672: {
678: PetscUseMethod(snes,"SNESCompositeGetSNES_C",(SNES,PetscInt,SNES*),(snes,n,subsnes));
679: return(0);
680: }
684: /*@
685: SNESCompositeGetNumber - Get the number of subsolvers in the composite SNES.
687: Logically Collective on SNES
689: Input Parameter:
690: snes - the preconditioner context
692: Output Parameter:
693: n - the number of subsolvers
695: Level: Developer
697: .keywords: SNES, composite preconditioner
698: @*/
699: PetscErrorCode SNESCompositeGetNumber(SNES snes,PetscInt *n)
700: {
701: SNES_Composite *jac;
702: SNES_CompositeLink next;
705: jac = (SNES_Composite*)snes->data;
706: next = jac->head;
708: *n = 0;
709: while (next) {
710: *n = *n + 1;
711: next = next->next;
712: }
713: return(0);
714: }
718: static PetscErrorCode SNESCompositeSetDamping_Composite(SNES snes,PetscInt n,PetscReal dmp)
719: {
720: SNES_Composite *jac;
721: SNES_CompositeLink next;
722: PetscInt i;
725: jac = (SNES_Composite*)snes->data;
726: next = jac->head;
727: for (i=0; i<n; i++) {
728: if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
729: next = next->next;
730: }
731: next->dmp = dmp;
732: return(0);
733: }
737: /*@
738: SNESCompositeSetDamping - Sets the damping of a subsolver when using additive composite SNES.
740: Not Collective
742: Input Parameter:
743: + snes - the preconditioner context
744: . n - the number of the snes requested
745: - dmp - the damping
747: Level: Developer
749: .keywords: SNES, get, composite preconditioner, sub preconditioner
751: .seealso: SNESCompositeAddSNES()
752: @*/
753: PetscErrorCode SNESCompositeSetDamping(SNES snes,PetscInt n,PetscReal dmp)
754: {
759: PetscUseMethod(snes,"SNESCompositeSetDamping_C",(SNES,PetscInt,PetscReal),(snes,n,dmp));
760: return(0);
761: }
765: PetscErrorCode SNESSolve_Composite(SNES snes)
766: {
767: Vec F;
768: Vec X;
769: Vec B;
770: PetscInt i;
771: PetscReal fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
773: SNESNormSchedule normtype;
774: SNES_Composite *comp = (SNES_Composite*)snes->data;
777: X = snes->vec_sol;
778: F = snes->vec_func;
779: B = snes->vec_rhs;
781: PetscObjectSAWsTakeAccess((PetscObject)snes);
782: snes->iter = 0;
783: snes->norm = 0.;
784: comp->innerFailures = 0;
785: PetscObjectSAWsGrantAccess((PetscObject)snes);
786: SNESSetWorkVecs(snes, 1);
787: snes->reason = SNES_CONVERGED_ITERATING;
788: SNESGetNormSchedule(snes, &normtype);
789: if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
790: if (!snes->vec_func_init_set) {
791: SNESComputeFunction(snes,X,F);
792: } else snes->vec_func_init_set = PETSC_FALSE;
794: if (snes->xl && snes->xu) {
795: SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
796: } else {
797: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
798: }
799: SNESCheckFunctionNorm(snes,fnorm);
800: PetscObjectSAWsTakeAccess((PetscObject)snes);
801: snes->iter = 0;
802: snes->norm = fnorm;
803: PetscObjectSAWsGrantAccess((PetscObject)snes);
804: SNESLogConvergenceHistory(snes,snes->norm,0);
805: SNESMonitor(snes,0,snes->norm);
807: /* test convergence */
808: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
809: if (snes->reason) return(0);
810: } else {
811: PetscObjectSAWsGrantAccess((PetscObject)snes);
812: SNESLogConvergenceHistory(snes,snes->norm,0);
813: SNESMonitor(snes,0,snes->norm);
814: }
816: /* Call general purpose update function */
817: if (snes->ops->update) {
818: (*snes->ops->update)(snes, snes->iter);
819: }
821: for (i = 0; i < snes->max_its; i++) {
822: /* Copy the state before modification by application of the composite solver;
823: we will subtract the new state after application */
824: VecCopy(X, snes->work[0]);
826: if (comp->type == SNES_COMPOSITE_ADDITIVE) {
827: SNESCompositeApply_Additive(snes,X,B,F,&fnorm);
828: } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
829: SNESCompositeApply_Multiplicative(snes,X,B,F,&fnorm);
830: } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
831: SNESCompositeApply_AdditiveOptimal(snes,X,B,F,&fnorm);
832: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported SNESComposite type");
833: if (snes->reason < 0) break;
835: /* Compute the solution update for convergence testing */
836: VecAXPY(snes->work[0], -1.0, X);
837: VecScale(snes->work[0], -1.0);
839: if ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY)) {
840: SNESComputeFunction(snes,X,F);
842: if (snes->xl && snes->xu) {
843: VecNormBegin(X, NORM_2, &xnorm);
844: VecNormBegin(snes->work[0], NORM_2, &snorm);
845: SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
846: VecNormEnd(X, NORM_2, &xnorm);
847: VecNormEnd(snes->work[0], NORM_2, &snorm);
848: } else {
849: VecNormBegin(F, NORM_2, &fnorm);
850: VecNormBegin(X, NORM_2, &xnorm);
851: VecNormBegin(snes->work[0], NORM_2, &snorm);
853: VecNormEnd(F, NORM_2, &fnorm);
854: VecNormEnd(X, NORM_2, &xnorm);
855: VecNormEnd(snes->work[0], NORM_2, &snorm);
856: }
857: SNESCheckFunctionNorm(snes,fnorm);
858: } else if (normtype == SNES_NORM_ALWAYS) {
859: VecNormBegin(X, NORM_2, &xnorm);
860: VecNormBegin(snes->work[0], NORM_2, &snorm);
861: VecNormEnd(X, NORM_2, &xnorm);
862: VecNormEnd(snes->work[0], NORM_2, &snorm);
863: }
864: /* Monitor convergence */
865: PetscObjectSAWsTakeAccess((PetscObject)snes);
866: snes->iter = i+1;
867: snes->norm = fnorm;
868: PetscObjectSAWsGrantAccess((PetscObject)snes);
869: SNESLogConvergenceHistory(snes,snes->norm,0);
870: SNESMonitor(snes,snes->iter,snes->norm);
871: /* Test for convergence */
872: if (normtype == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,xnorm,snorm,fnorm,&snes->reason,snes->cnvP);}
873: if (snes->reason) break;
874: /* Call general purpose update function */
875: if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
876: }
877: if (normtype == SNES_NORM_ALWAYS) {
878: if (i == snes->max_its) {
879: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
880: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
881: }
882: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
883: return(0);
884: }
886: /* -------------------------------------------------------------------------------------------*/
888: /*MC
889: SNESCOMPOSITE - Build a preconditioner by composing together several nonlinear solvers
891: Options Database Keys:
892: + -snes_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
893: - -snes_composite_sneses - <snes0,snes1,...> list of SNESes to compose
895: Level: intermediate
897: Concepts: composing solvers
899: .seealso: SNESCreate(), SNESSetType(), SNESType (for list of available types), SNES,
900: SNESSHELL, SNESCompositeSetType(), SNESCompositeSpecialSetAlpha(), SNESCompositeAddSNES(),
901: SNESCompositeGetSNES()
903: References:
904: . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
905: SIAM Review, 57(4), 2015
907: M*/
911: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
912: {
914: SNES_Composite *jac;
917: PetscNewLog(snes,&jac);
919: snes->ops->solve = SNESSolve_Composite;
920: snes->ops->setup = SNESSetUp_Composite;
921: snes->ops->reset = SNESReset_Composite;
922: snes->ops->destroy = SNESDestroy_Composite;
923: snes->ops->setfromoptions = SNESSetFromOptions_Composite;
924: snes->ops->view = SNESView_Composite;
926: snes->data = (void*)jac;
927: jac->type = SNES_COMPOSITE_ADDITIVEOPTIMAL;
928: jac->Fes = NULL;
929: jac->Xes = NULL;
930: jac->fnorms = NULL;
931: jac->nsnes = 0;
932: jac->head = 0;
933: jac->stol = 0.1;
934: jac->rtol = 1.1;
936: jac->h = NULL;
937: jac->s = NULL;
938: jac->beta = NULL;
939: jac->work = NULL;
940: jac->rwork = NULL;
942: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetType_C",SNESCompositeSetType_Composite);
943: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeAddSNES_C",SNESCompositeAddSNES_Composite);
944: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeGetSNES_C",SNESCompositeGetSNES_Composite);
945: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetDamping_C",SNESCompositeSetDamping_Composite);
946: return(0);
947: }