Actual source code: snescomposite.c
petsc-3.12.5 2020-03-29
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: @*/
592: PetscErrorCode SNESCompositeSetType(SNES snes,SNESCompositeType type)
593: {
599: PetscTryMethod(snes,"SNESCompositeSetType_C",(SNES,SNESCompositeType),(snes,type));
600: return(0);
601: }
603: /*@C
604: SNESCompositeAddSNES - Adds another SNES to the composite SNES.
606: Collective on SNES
608: Input Parameters:
609: + snes - the preconditioner context
610: - type - the type of the new preconditioner
612: Level: Developer
614: @*/
615: PetscErrorCode SNESCompositeAddSNES(SNES snes,SNESType type)
616: {
621: PetscTryMethod(snes,"SNESCompositeAddSNES_C",(SNES,SNESType),(snes,type));
622: return(0);
623: }
624: /*@
625: SNESCompositeGetSNES - Gets one of the SNES objects in the composite SNES.
627: Not Collective
629: Input Parameter:
630: + snes - the preconditioner context
631: - n - the number of the snes requested
633: Output Parameters:
634: . subsnes - the SNES requested
636: Level: Developer
638: .seealso: SNESCompositeAddSNES()
639: @*/
640: PetscErrorCode SNESCompositeGetSNES(SNES snes,PetscInt n,SNES *subsnes)
641: {
647: PetscUseMethod(snes,"SNESCompositeGetSNES_C",(SNES,PetscInt,SNES*),(snes,n,subsnes));
648: return(0);
649: }
651: /*@
652: SNESCompositeGetNumber - Get the number of subsolvers in the composite SNES.
654: Logically Collective on SNES
656: Input Parameter:
657: snes - the preconditioner context
659: Output Parameter:
660: n - the number of subsolvers
662: Level: Developer
664: @*/
665: PetscErrorCode SNESCompositeGetNumber(SNES snes,PetscInt *n)
666: {
667: SNES_Composite *jac;
668: SNES_CompositeLink next;
671: jac = (SNES_Composite*)snes->data;
672: next = jac->head;
674: *n = 0;
675: while (next) {
676: *n = *n + 1;
677: next = next->next;
678: }
679: return(0);
680: }
682: static PetscErrorCode SNESCompositeSetDamping_Composite(SNES snes,PetscInt n,PetscReal dmp)
683: {
684: SNES_Composite *jac;
685: SNES_CompositeLink next;
686: PetscInt i;
689: jac = (SNES_Composite*)snes->data;
690: next = jac->head;
691: for (i=0; i<n; i++) {
692: if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
693: next = next->next;
694: }
695: next->dmp = dmp;
696: return(0);
697: }
699: /*@
700: SNESCompositeSetDamping - Sets the damping of a subsolver when using additive composite SNES.
702: Not Collective
704: Input Parameter:
705: + snes - the preconditioner context
706: . n - the number of the snes requested
707: - dmp - the damping
709: Level: Developer
711: .seealso: SNESCompositeAddSNES()
712: @*/
713: PetscErrorCode SNESCompositeSetDamping(SNES snes,PetscInt n,PetscReal dmp)
714: {
719: PetscUseMethod(snes,"SNESCompositeSetDamping_C",(SNES,PetscInt,PetscReal),(snes,n,dmp));
720: return(0);
721: }
723: static PetscErrorCode SNESSolve_Composite(SNES snes)
724: {
725: Vec F;
726: Vec X;
727: Vec B;
728: PetscInt i;
729: PetscReal fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
730: PetscErrorCode ierr;
731: SNESNormSchedule normtype;
732: SNES_Composite *comp = (SNES_Composite*)snes->data;
735: X = snes->vec_sol;
736: F = snes->vec_func;
737: B = snes->vec_rhs;
739: PetscObjectSAWsTakeAccess((PetscObject)snes);
740: snes->iter = 0;
741: snes->norm = 0.;
742: comp->innerFailures = 0;
743: PetscObjectSAWsGrantAccess((PetscObject)snes);
744: SNESSetWorkVecs(snes, 1);
745: snes->reason = SNES_CONVERGED_ITERATING;
746: SNESGetNormSchedule(snes, &normtype);
747: if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
748: if (!snes->vec_func_init_set) {
749: SNESComputeFunction(snes,X,F);
750: } else snes->vec_func_init_set = PETSC_FALSE;
752: if (snes->xl && snes->xu) {
753: SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
754: } else {
755: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
756: }
757: SNESCheckFunctionNorm(snes,fnorm);
758: PetscObjectSAWsTakeAccess((PetscObject)snes);
759: snes->iter = 0;
760: snes->norm = fnorm;
761: PetscObjectSAWsGrantAccess((PetscObject)snes);
762: SNESLogConvergenceHistory(snes,snes->norm,0);
763: SNESMonitor(snes,0,snes->norm);
765: /* test convergence */
766: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
767: if (snes->reason) return(0);
768: } else {
769: PetscObjectSAWsGrantAccess((PetscObject)snes);
770: SNESLogConvergenceHistory(snes,snes->norm,0);
771: SNESMonitor(snes,0,snes->norm);
772: }
774: /* Call general purpose update function */
775: if (snes->ops->update) {
776: (*snes->ops->update)(snes, snes->iter);
777: }
779: for (i = 0; i < snes->max_its; i++) {
780: /* Copy the state before modification by application of the composite solver;
781: we will subtract the new state after application */
782: VecCopy(X, snes->work[0]);
784: if (comp->type == SNES_COMPOSITE_ADDITIVE) {
785: SNESCompositeApply_Additive(snes,X,B,F,&fnorm);
786: } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
787: SNESCompositeApply_Multiplicative(snes,X,B,F,&fnorm);
788: } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
789: SNESCompositeApply_AdditiveOptimal(snes,X,B,F,&fnorm);
790: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported SNESComposite type");
791: if (snes->reason < 0) break;
793: /* Compute the solution update for convergence testing */
794: VecAXPY(snes->work[0], -1.0, X);
795: VecScale(snes->work[0], -1.0);
797: if ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY)) {
798: SNESComputeFunction(snes,X,F);
800: if (snes->xl && snes->xu) {
801: VecNormBegin(X, NORM_2, &xnorm);
802: VecNormBegin(snes->work[0], NORM_2, &snorm);
803: SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
804: VecNormEnd(X, NORM_2, &xnorm);
805: VecNormEnd(snes->work[0], NORM_2, &snorm);
806: } else {
807: VecNormBegin(F, NORM_2, &fnorm);
808: VecNormBegin(X, NORM_2, &xnorm);
809: VecNormBegin(snes->work[0], NORM_2, &snorm);
811: VecNormEnd(F, NORM_2, &fnorm);
812: VecNormEnd(X, NORM_2, &xnorm);
813: VecNormEnd(snes->work[0], NORM_2, &snorm);
814: }
815: SNESCheckFunctionNorm(snes,fnorm);
816: } else if (normtype == SNES_NORM_ALWAYS) {
817: VecNormBegin(X, NORM_2, &xnorm);
818: VecNormBegin(snes->work[0], NORM_2, &snorm);
819: VecNormEnd(X, NORM_2, &xnorm);
820: VecNormEnd(snes->work[0], NORM_2, &snorm);
821: }
822: /* Monitor convergence */
823: PetscObjectSAWsTakeAccess((PetscObject)snes);
824: snes->iter = i+1;
825: snes->norm = fnorm;
826: snes->xnorm = xnorm;
827: snes->ynorm = snorm;
828: PetscObjectSAWsGrantAccess((PetscObject)snes);
829: SNESLogConvergenceHistory(snes,snes->norm,0);
830: SNESMonitor(snes,snes->iter,snes->norm);
831: /* Test for convergence */
832: if (normtype == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,xnorm,snorm,fnorm,&snes->reason,snes->cnvP);}
833: if (snes->reason) break;
834: /* Call general purpose update function */
835: if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
836: }
837: if (normtype == SNES_NORM_ALWAYS) {
838: if (i == snes->max_its) {
839: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
840: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
841: }
842: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
843: return(0);
844: }
846: /* -------------------------------------------------------------------------------------------*/
848: /*MC
849: SNESCOMPOSITE - Build a preconditioner by composing together several nonlinear solvers
851: Options Database Keys:
852: + -snes_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
853: - -snes_composite_sneses - <snes0,snes1,...> list of SNESes to compose
855: Level: intermediate
857: .seealso: SNESCreate(), SNESSetType(), SNESType (for list of available types), SNES,
858: SNESSHELL, SNESCompositeSetType(), SNESCompositeSpecialSetAlpha(), SNESCompositeAddSNES(),
859: SNESCompositeGetSNES()
861: References:
862: . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
863: SIAM Review, 57(4), 2015
865: M*/
867: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
868: {
870: SNES_Composite *jac;
873: PetscNewLog(snes,&jac);
875: snes->ops->solve = SNESSolve_Composite;
876: snes->ops->setup = SNESSetUp_Composite;
877: snes->ops->reset = SNESReset_Composite;
878: snes->ops->destroy = SNESDestroy_Composite;
879: snes->ops->setfromoptions = SNESSetFromOptions_Composite;
880: snes->ops->view = SNESView_Composite;
882: snes->usesksp = PETSC_FALSE;
884: snes->alwayscomputesfinalresidual = PETSC_FALSE;
886: snes->data = (void*)jac;
887: jac->type = SNES_COMPOSITE_ADDITIVEOPTIMAL;
888: jac->Fes = NULL;
889: jac->Xes = NULL;
890: jac->fnorms = NULL;
891: jac->nsnes = 0;
892: jac->head = 0;
893: jac->stol = 0.1;
894: jac->rtol = 1.1;
896: jac->h = NULL;
897: jac->s = NULL;
898: jac->beta = NULL;
899: jac->work = NULL;
900: jac->rwork = NULL;
902: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetType_C",SNESCompositeSetType_Composite);
903: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeAddSNES_C",SNESCompositeAddSNES_Composite);
904: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeGetSNES_C",SNESCompositeGetSNES_Composite);
905: PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetDamping_C",SNESCompositeSetDamping_Composite);
906: return(0);
907: }