Actual source code: gl.c
petsc-3.3-p7 2013-05-11
2: #include <../src/ts/impls/implicit/gl/gl.h> /*I "petscts.h" I*/
3: #include <petscblaslapack.h>
5: static const char *TSGLErrorDirections[] = {"FORWARD","BACKWARD","TSGLErrorDirection","TSGLERROR_",0};
6: static PetscFList TSGLList;
7: static PetscFList TSGLAcceptList;
8: static PetscBool TSGLPackageInitialized;
9: static PetscBool TSGLRegisterAllCalled;
11: /* C++ does not promote int64_t to scalar or int32_t for std::pow() */
12: static PetscScalar Pow(PetscScalar b,PetscInt p)
13: {
14: return PetscPowScalar(b,(int)p);
15: }
17: /* This function is pure */
18: static PetscScalar Factorial(PetscInt n)
19: {
20: PetscInt i;
21: if (n < 12) { /* Can compute with 32-bit integers */
22: PetscInt f = 1;
23: for (i=2; i<=n; i++) f *= i;
24: return (PetscScalar)f;
25: } else {
26: PetscScalar f = 1.;
27: for (i=2; i<=n; i++) f *= (PetscScalar)i;
28: return f;
29: }
30: }
32: /* This function is pure */
33: static PetscScalar CPowF(PetscScalar c,PetscInt p)
34: {
35: return Pow(c,p)/Factorial(p);
36: }
40: static PetscErrorCode TSGLSchemeCreate(PetscInt p,PetscInt q,PetscInt r,PetscInt s,const PetscScalar *c,
41: const PetscScalar *a,const PetscScalar *b,const PetscScalar *u,const PetscScalar *v,TSGLScheme *inscheme)
42: {
43: TSGLScheme scheme;
44: PetscInt j;
49: if (p < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Scheme order must be positive");
50: if (r < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"At least one item must be carried between steps");
51: if (s < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"At least one stage is required");
53: *inscheme = 0;
54: PetscNew(struct _TSGLScheme,&scheme);
55: scheme->p = p;
56: scheme->q = q;
57: scheme->r = r;
58: scheme->s = s;
60: PetscMalloc5(s,PetscScalar,&scheme->c,s*s,PetscScalar,&scheme->a,r*s,PetscScalar,&scheme->b,r*s,PetscScalar,&scheme->u,r*r,PetscScalar,&scheme->v);
61: PetscMemcpy(scheme->c,c,s*sizeof(PetscScalar));
62: for (j=0; j<s*s; j++) scheme->a[j] = (PetscAbsScalar(a[j]) < 1e-12) ? 0 : a[j];
63: for (j=0; j<r*s; j++) scheme->b[j] = (PetscAbsScalar(b[j]) < 1e-12) ? 0 : b[j];
64: for (j=0; j<s*r; j++) scheme->u[j] = (PetscAbsScalar(u[j]) < 1e-12) ? 0 : u[j];
65: for (j=0; j<r*r; j++) scheme->v[j] = (PetscAbsScalar(v[j]) < 1e-12) ? 0 : v[j];
67: PetscMalloc6(r,PetscScalar,&scheme->alpha,r,PetscScalar,&scheme->beta,r,PetscScalar,&scheme->gamma,3*s,PetscScalar,&scheme->phi,3*r,PetscScalar,&scheme->psi,r,PetscScalar,&scheme->stage_error);
68: {
69: PetscInt i,j,k,ss=s+2;
70: PetscBLASInt m,n,one=1,*ipiv,lwork=4*((s+3)*3+3),info,ldb;
71: PetscReal rcond,*sing,*workreal;
72: PetscScalar *ImV,*H,*bmat,*workscalar,*c=scheme->c,*a=scheme->a,*b=scheme->b,*u=scheme->u,*v=scheme->v;
73: #if !defined(PETSC_MISSING_LAPACK_GELSS)
74: PetscBLASInt rank;
75: #endif
76: PetscMalloc7(PetscSqr(r),PetscScalar,&ImV,3*s,PetscScalar,&H,3*ss,PetscScalar,&bmat,lwork,PetscScalar,&workscalar,5*(3+r),PetscReal,&workreal,r+s,PetscReal,&sing,r+s,PetscBLASInt,&ipiv);
78: /* column-major input */
79: for (i=0; i<r-1; i++) {
80: for (j=0; j<r-1; j++) {
81: ImV[i+j*r] = 1.0*(i==j) - v[(i+1)*r+j+1];
82: }
83: }
84: /* Build right hand side for alpha (tp - glm.B(2:end,:)*(glm.c.^(p)./factorial(p))) */
85: for (i=1; i<r; i++) {
86: scheme->alpha[i] = 1./Factorial(p+1-i);
87: for (j=0; j<s; j++) scheme->alpha[i] -= b[i*s+j]*CPowF(c[j],p);
88: }
89: m = PetscBLASIntCast(r-1);
90: n = PetscBLASIntCast(r);
91: LAPACKgesv_(&m,&one,ImV,&n,ipiv,scheme->alpha+1,&n,&info);
92: if (info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GESV");
93: if (info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
95: /* Build right hand side for beta (tp1 - glm.B(2:end,:)*(glm.c.^(p+1)./factorial(p+1)) - e.alpha) */
96: for (i=1; i<r; i++) {
97: scheme->beta[i] = 1./Factorial(p+2-i) - scheme->alpha[i];
98: for (j=0; j<s; j++) scheme->beta[i] -= b[i*s+j]*CPowF(c[j],p+1);
99: }
100: LAPACKgetrs_("No transpose",&m,&one,ImV,&n,ipiv,scheme->beta+1,&n,&info);
101: if (info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GETRS");
102: if (info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Should not happen");
104: /* Build stage_error vector
105: xi = glm.c.^(p+1)/factorial(p+1) - glm.A*glm.c.^p/factorial(p) + glm.U(:,2:end)*e.alpha;
106: */
107: for (i=0; i<s; i++) {
108: scheme->stage_error[i] = CPowF(c[i],p+1);
109: for (j=0; j<s; j++) scheme->stage_error[i] -= a[i*s+j]*CPowF(c[j],p);
110: for (j=1; j<r; j++) scheme->stage_error[i] += u[i*r+j]*scheme->alpha[j];
111: }
113: /* alpha[0] (epsilon in B,J,W 2007)
114: epsilon = 1/factorial(p+1) - B(1,:)*c.^p/factorial(p) + V(1,2:end)*e.alpha;
115: */
116: scheme->alpha[0] = 1./Factorial(p+1);
117: for (j=0; j<s; j++) scheme->alpha[0] -= b[0*s+j]*CPowF(c[j],p);
118: for (j=1; j<r; j++) scheme->alpha[0] += v[0*r+j]*scheme->alpha[j];
120: /* right hand side for gamma (glm.B(2:end,:)*e.xi - e.epsilon*eye(s-1,1)) */
121: for (i=1; i<r; i++) {
122: scheme->gamma[i] = (i==1 ? -1. : 0)*scheme->alpha[0];
123: for (j=0; j<s; j++) scheme->gamma[i] += b[i*s+j]*scheme->stage_error[j];
124: }
125: LAPACKgetrs_("No transpose",&m,&one,ImV,&n,ipiv,scheme->gamma+1,&n,&info);
126: if (info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GETRS");
127: if (info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Should not happen");
129: /* beta[0] (rho in B,J,W 2007)
130: e.rho = 1/factorial(p+2) - glm.B(1,:)*glm.c.^(p+1)/factorial(p+1) ...
131: + glm.V(1,2:end)*e.beta;% - e.epsilon;
132: % Note: The paper (B,J,W 2007) includes the last term in their definition
133: * */
134: scheme->beta[0] = 1./Factorial(p+2);
135: for (j=0; j<s; j++) scheme->beta[0] -= b[0*s+j]*CPowF(c[j],p+1);
136: for (j=1; j<r; j++) scheme->beta[0] += v[0*r+j]*scheme->beta[j];
138: /* gamma[0] (sigma in B,J,W 2007)
139: * e.sigma = glm.B(1,:)*e.xi + glm.V(1,2:end)*e.gamma;
140: * */
141: scheme->gamma[0] = 0.0;
142: for (j=0; j<s; j++) scheme->gamma[0] += b[0*s+j]*scheme->stage_error[j];
143: for (j=1; j<r; j++) scheme->gamma[0] += v[0*s+j]*scheme->gamma[j];
145: /* Assemble H
146: * % Determine the error estimators phi
147: H = [[cpow(glm.c,p) + C*e.alpha] [cpow(glm.c,p+1) + C*e.beta] ...
148: [e.xi - C*(e.gamma + 0*e.epsilon*eye(s-1,1))]]';
149: % Paper has formula above without the 0, but that term must be left
150: % out to satisfy the conditions they propose and to make the
151: % example schemes work
152: e.H = H;
153: e.phi = (H \ [1 0 0;1 1 0;0 0 -1])';
154: e.psi = -e.phi*C;
155: * */
156: for (j=0; j<s; j++) {
157: H[0+j*3] = CPowF(c[j],p);
158: H[1+j*3] = CPowF(c[j],p+1);
159: H[2+j*3] = scheme->stage_error[j];
160: for (k=1; k<r; k++) {
161: H[0+j*3] += CPowF(c[j],k-1)*scheme->alpha[k];
162: H[1+j*3] += CPowF(c[j],k-1)*scheme->beta[k];
163: H[2+j*3] -= CPowF(c[j],k-1)*scheme->gamma[k];
164: }
165: }
166: bmat[0+0*ss] = 1.; bmat[0+1*ss] = 0.; bmat[0+2*ss] = 0.;
167: bmat[1+0*ss] = 1.; bmat[1+1*ss] = 1.; bmat[1+2*ss] = 0.;
168: bmat[2+0*ss] = 0.; bmat[2+1*ss] = 0.; bmat[2+2*ss] = -1.;
169: m = 3;
170: n = PetscBLASIntCast(s);
171: ldb = PetscBLASIntCast(ss);
172: rcond = 1e-12;
173: #if defined(PETSC_MISSING_LAPACK_GELSS)
174: /* ESSL does not have this routine */
175: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GELSS - Lapack routine is unavailable\nNot able to run GL time stepping.");
176: #else
177: #if defined(PETSC_USE_COMPLEX)
178: /* ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO ) */
179: LAPACKgelss_(&m,&n,&m,H,&m,bmat,&ldb,sing,&rcond,&rank,workscalar,&lwork,workreal,&info);
180: #else
181: /* DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO ) */
182: LAPACKgelss_(&m,&n,&m,H,&m,bmat,&ldb,sing,&rcond,&rank,workscalar,&lwork,&info);
183: #endif
184: if (info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS");
185: if (info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge");
186: #endif
188: for (j=0; j<3; j++) {
189: for (k=0; k<s; k++) {
190: scheme->phi[k+j*s] = bmat[k+j*ss];
191: }
192: }
194: /* the other part of the error estimator, psi in B,J,W 2007 */
195: scheme->psi[0*r+0] = 0.;
196: scheme->psi[1*r+0] = 0.;
197: scheme->psi[2*r+0] = 0.;
198: for (j=1; j<r; j++) {
199: scheme->psi[0*r+j] = 0.;
200: scheme->psi[1*r+j] = 0.;
201: scheme->psi[2*r+j] = 0.;
202: for (k=0; k<s; k++) {
203: scheme->psi[0*r+j] -= CPowF(c[k],j-1)*scheme->phi[0*s+k];
204: scheme->psi[1*r+j] -= CPowF(c[k],j-1)*scheme->phi[1*s+k];
205: scheme->psi[2*r+j] -= CPowF(c[k],j-1)*scheme->phi[2*s+k];
206: }
207: }
208: PetscFree7(ImV,H,bmat,workscalar,workreal,sing,ipiv);
209: }
210: /* Check which properties are satisfied */
211: scheme->stiffly_accurate = PETSC_TRUE;
212: if (scheme->c[s-1] != 1.) scheme->stiffly_accurate = PETSC_FALSE;
213: for (j=0; j<s; j++) if (a[(s-1)*s+j] != b[j]) scheme->stiffly_accurate = PETSC_FALSE;
214: for (j=0; j<r; j++) if (u[(s-1)*r+j] != v[j]) scheme->stiffly_accurate = PETSC_FALSE;
215: scheme->fsal = scheme->stiffly_accurate; /* FSAL is stronger */
216: for (j=0; j<s-1; j++) if (r>1 && b[1*s+j] != 0.) scheme->fsal = PETSC_FALSE;
217: if (b[1*s+r-1] != 1.) scheme->fsal = PETSC_FALSE;
218: for (j=0; j<r; j++) if (r>1 && v[1*r+j] != 0.) scheme->fsal = PETSC_FALSE;
220: *inscheme = scheme;
221: return(0);
222: }
226: static PetscErrorCode TSGLSchemeDestroy(TSGLScheme sc)
227: {
231: PetscFree5(sc->c,sc->a,sc->b,sc->u,sc->v);
232: PetscFree6(sc->alpha,sc->beta,sc->gamma,sc->phi,sc->psi,sc->stage_error);
233: PetscFree(sc);
234: return(0);
235: }
239: static PetscErrorCode TSGLDestroy_Default(TS_GL *gl)
240: {
242: PetscInt i;
245: for (i=0; i<gl->nschemes; i++) {
246: if (gl->schemes[i]) {TSGLSchemeDestroy(gl->schemes[i]);}
247: }
248: PetscFree(gl->schemes);
249: gl->nschemes = 0;
250: PetscMemzero(gl->type_name,sizeof(gl->type_name));
251: return(0);
252: }
256: static PetscErrorCode TSGLViewTable_Private(PetscViewer viewer,PetscInt m,PetscInt n,const PetscScalar a[],const char name[])
257: {
259: PetscBool iascii;
260: PetscInt i,j;
263: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
264: if (iascii) {
265: PetscViewerASCIIPrintf(viewer,"%30s = [",name);
266: for (i=0; i<m; i++) {
267: if (i) {PetscViewerASCIIPrintf(viewer,"%30s [","");}
268: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
269: for (j=0; j<n; j++) {
270: PetscViewerASCIIPrintf(viewer," %12.8g",PetscRealPart(a[i*n+j]));
271: }
272: PetscViewerASCIIPrintf(viewer,"]\n");
273: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
274: }
275: } else {
276: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
277: }
278: return(0);
279: }
284: static PetscErrorCode TSGLSchemeView(TSGLScheme sc,PetscBool view_details,PetscViewer viewer)
285: {
287: PetscBool iascii;
290: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
291: if (iascii) {
292: PetscViewerASCIIPrintf(viewer,"GL scheme p,q,r,s = %d,%d,%d,%d\n",sc->p,sc->q,sc->r,sc->s);
293: PetscViewerASCIIPushTab(viewer);
294: PetscViewerASCIIPrintf(viewer,"Stiffly accurate: %s, FSAL: %s\n",sc->stiffly_accurate?"yes":"no",sc->fsal?"yes":"no");
295: PetscViewerASCIIPrintf(viewer,"Leading error constants: %10.3e %10.3e %10.3e\n",
296: PetscRealPart(sc->alpha[0]),PetscRealPart(sc->beta[0]),PetscRealPart(sc->gamma[0]));
297: TSGLViewTable_Private(viewer,1,sc->s,sc->c,"Abscissas c");
298: if (view_details) {
299: TSGLViewTable_Private(viewer,sc->s,sc->s,sc->a,"A");
300: TSGLViewTable_Private(viewer,sc->r,sc->s,sc->b,"B");
301: TSGLViewTable_Private(viewer,sc->s,sc->r,sc->u,"U");
302: TSGLViewTable_Private(viewer,sc->r,sc->r,sc->v,"V");
304: TSGLViewTable_Private(viewer,3,sc->s,sc->phi,"Error estimate phi");
305: TSGLViewTable_Private(viewer,3,sc->r,sc->psi,"Error estimate psi");
306: TSGLViewTable_Private(viewer,1,sc->r,sc->alpha,"Modify alpha");
307: TSGLViewTable_Private(viewer,1,sc->r,sc->beta,"Modify beta");
308: TSGLViewTable_Private(viewer,1,sc->r,sc->gamma,"Modify gamma");
309: TSGLViewTable_Private(viewer,1,sc->s,sc->stage_error,"Stage error xi");
310: }
311: PetscViewerASCIIPopTab(viewer);
312: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
313: return(0);
314: }
318: static PetscErrorCode TSGLEstimateHigherMoments_Default(TSGLScheme sc,PetscReal h,Vec Ydot[],Vec Xold[],Vec hm[])
319: {
321: PetscInt i;
324: if (sc->r > 64 || sc->s > 64) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Ridiculous number of stages or items passed between stages");
325: /* build error vectors*/
326: for (i=0; i<3; i++) {
327: PetscScalar phih[64];
328: PetscInt j;
329: for (j=0; j<sc->s; j++) phih[j] = sc->phi[i*sc->s+j]*h;
330: VecZeroEntries(hm[i]);
331: VecMAXPY(hm[i],sc->s,phih,Ydot);
332: VecMAXPY(hm[i],sc->r,&sc->psi[i*sc->r],Xold);
333: }
334: return(0);
335: }
339: static PetscErrorCode TSGLCompleteStep_Rescale(TSGLScheme sc,PetscReal h,TSGLScheme next_sc,PetscReal next_h,Vec Ydot[],Vec Xold[],Vec X[])
340: {
342: PetscScalar brow[32],vrow[32];
343: PetscInt i,j,r,s;
346: /* Build the new solution from (X,Ydot) */
347: r = sc->r;
348: s = sc->s;
349: for (i=0; i<r; i++) {
350: VecZeroEntries(X[i]);
351: for (j=0; j<s; j++) brow[j] = h*sc->b[i*s+j];
352: VecMAXPY(X[i],s,brow,Ydot);
353: for (j=0; j<r; j++) vrow[j] = sc->v[i*r+j];
354: VecMAXPY(X[i],r,vrow,Xold);
355: }
356: return(0);
357: }
361: static PetscErrorCode TSGLCompleteStep_RescaleAndModify(TSGLScheme sc,PetscReal h,TSGLScheme next_sc,PetscReal next_h,Vec Ydot[],Vec Xold[],Vec X[])
362: {
364: PetscScalar brow[32],vrow[32];
365: PetscReal ratio;
366: PetscInt i,j,p,r,s;
369: /* Build the new solution from (X,Ydot) */
370: p = sc->p;
371: r = sc->r;
372: s = sc->s;
373: ratio = next_h/h;
374: for (i=0; i<r; i++) {
375: VecZeroEntries(X[i]);
376: for (j=0; j<s; j++) {
377: brow[j] = h*(Pow(ratio,i)*sc->b[i*s+j]
378: + (Pow(ratio,i) - Pow(ratio,p+1))*(+ sc->alpha[i]*sc->phi[0*s+j])
379: + (Pow(ratio,i) - Pow(ratio,p+2))*(+ sc->beta [i]*sc->phi[1*s+j]
380: + sc->gamma[i]*sc->phi[2*s+j]));
381: }
382: VecMAXPY(X[i],s,brow,Ydot);
383: for (j=0; j<r; j++) {
384: vrow[j] = (Pow(ratio,i)*sc->v[i*r+j]
385: + (Pow(ratio,i) - Pow(ratio,p+1))*(+ sc->alpha[i]*sc->psi[0*r+j])
386: + (Pow(ratio,i) - Pow(ratio,p+2))*(+ sc->beta [i]*sc->psi[1*r+j]
387: + sc->gamma[i]*sc->psi[2*r+j]));
388: }
389: VecMAXPY(X[i],r,vrow,Xold);
390: }
391: if (r < next_sc->r) {
392: if (r+1 != next_sc->r) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cannot accommodate jump in r greater than 1");
393: VecZeroEntries(X[r]);
394: for (j=0; j<s; j++) brow[j] = h*Pow(ratio,p+1)*sc->phi[0*s+j];
395: VecMAXPY(X[r],s,brow,Ydot);
396: for (j=0; j<r; j++) vrow[j] = Pow(ratio,p+1)*sc->psi[0*r+j];
397: VecMAXPY(X[r],r,vrow,Xold);
398: }
399: return(0);
400: }
402: EXTERN_C_BEGIN
405: PetscErrorCode TSGLCreate_IRKS(TS ts)
406: {
407: TS_GL *gl = (TS_GL*)ts->data;
411: gl->Destroy = TSGLDestroy_Default;
412: gl->EstimateHigherMoments = TSGLEstimateHigherMoments_Default;
413: gl->CompleteStep = TSGLCompleteStep_RescaleAndModify;
414: PetscMalloc(10*sizeof(TSGLScheme),&gl->schemes);
415: gl->nschemes = 0;
417: {
418: /* p=1,q=1, r=s=2, A- and L-stable with error estimates of order 2 and 3
419: * Listed in Butcher & Podhaisky 2006. On error estimation in general linear methods for stiff ODE.
420: * irks(0.3,0,[.3,1],[1],1)
421: * Note: can be made to have classical order (not stage order) 2 by replacing 0.3 with 1-sqrt(1/2)
422: * but doing so would sacrifice the error estimator.
423: */
424: const PetscScalar c[2] = {3./10., 1.}
425: ,a[2][2] = {{3./10., 0}, {7./10., 3./10.}}
426: ,b[2][2] = {{7./10., 3./10.}, {0,1}}
427: ,u[2][2] = {{1,0},{1,0}}
428: ,v[2][2] = {{1,0},{0,0}};
429: TSGLSchemeCreate(1,1,2,2,c,*a,*b,*u,*v,&gl->schemes[gl->nschemes++]);
430: }
432: {
433: /* p=q=2, r=s=3: irks(4/9,0,[1:3]/3,[0.33852],1) */
434: /* http://www.math.auckland.ac.nz/~hpod/atlas/i2a.html */
435: const PetscScalar c[3] = {1./3., 2./3., 1}
436: ,a[3][3] = {{4./9. ,0 , 0},
437: {1.03750643704090e+00 , 4./9., 0},
438: {7.67024779410304e-01 , -3.81140216918943e-01, 4./9.}}
439: ,b[3][3] = {{0.767024779410304, -0.381140216918943, 4./9.},
440: {0.000000000000000, 0.000000000000000, 1.000000000000000},
441: {-2.075048385225385, 0.621728385225383, 1.277197204924873}}
442: ,u[3][3] = {{1.0000000000000000, -0.1111111111111109, -0.0925925925925922},
443: {1.0000000000000000, -0.8152842148186744, -0.4199095530877056},
444: {1.0000000000000000, 0.1696709930641948, 0.0539741070314165}}
445: ,v[3][3] = {{1.0000000000000000, 0.1696709930641948, 0.0539741070314165},
446: {0.000000000000000, 0.000000000000000, 0.000000000000000},
447: {0.000000000000000, 0.176122795075129, 0.000000000000000}};
448: TSGLSchemeCreate(2,2,3,3,c,*a,*b,*u,*v,&gl->schemes[gl->nschemes++]);
449: }
450: {
451: /* p=q=3, r=s=4: irks(9/40,0,[1:4]/4,[0.3312 1.0050],[0.49541 1;1 0]) */
452: const PetscScalar c[4] = {0.25,0.5,0.75,1.0}
453: ,a[4][4] = {{9./40. , 0, 0, 0},
454: {2.11286958887701e-01 , 9./40. , 0, 0},
455: {9.46338294287584e-01 , -3.42942861246094e-01, 9./40. , 0},
456: {0.521490453970721 , -0.662474225622980, 0.490476425459734, 9./40. }}
457: ,b[4][4] = {{0.521490453970721 , -0.662474225622980, 0.490476425459734, 9./40. },
458: {0.000000000000000 , 0.000000000000000, 0.000000000000000, 1.000000000000000},
459: {-0.084677029310348 , 1.390757514776085, -1.568157386206001, 2.023192696767826},
460: {0.465383797936408 , 1.478273530625148, -1.930836081010182, 1.644872111193354}}
461: ,u[4][4] = {{1.00000000000000000 , 0.02500000000001035, -0.02499999999999053, -0.00442708333332865},
462: {1.00000000000000000 , 0.06371304111232945, -0.04032173972189845, -0.01389438413189452},
463: {1.00000000000000000 , -0.07839543304147778, 0.04738685705116663, 0.02032603595928376},
464: {1.00000000000000000 , 0.42550734619251651, 0.10800718022400080, -0.01726712647760034}}
465: ,v[4][4] = {{1.00000000000000000 , 0.42550734619251651, 0.10800718022400080, -0.01726712647760034},
466: {0.000000000000000 , 0.000000000000000, 0.000000000000000, 0.000000000000000},
467: {0.000000000000000 , -1.761115796027561, -0.521284157173780, 0.258249384305463},
468: {0.000000000000000 , -1.657693358744728, -1.052227765232394, 0.521284157173780}};
469: TSGLSchemeCreate(3,3,4,4,c,*a,*b,*u,*v,&gl->schemes[gl->nschemes++]);
470: }
471: {
472: /* p=q=4, r=s=5:
473: irks(3/11,0,[1:5]/5, [0.1715 -0.1238 0.6617],...
474: [ -0.0812 0.4079 1.0000
475: 1.0000 0 0
476: 0.8270 1.0000 0])
477: */
478: const PetscScalar c[5] = {0.2,0.4,0.6,0.8,1.0}
479: ,a[5][5] = {{2.72727272727352e-01 , 0.00000000000000e+00, 0.00000000000000e+00 , 0.00000000000000e+00 , 0.00000000000000e+00},
480: {-1.03980153733431e-01, 2.72727272727405e-01, 0.00000000000000e+00, 0.00000000000000e+00 , 0.00000000000000e+00},
481: {-1.58615400341492e+00, 7.44168951881122e-01, 2.72727272727309e-01, 0.00000000000000e+00 , 0.00000000000000e+00},
482: {-8.73658042865628e-01, 5.37884671894595e-01, -1.63298538799523e-01, 2.72727272726996e-01 , 0.00000000000000e+00},
483: {2.95489397443992e-01 , -1.18481693910097e+00 , -6.68029812659953e-01 , 1.00716687860943e+00 , 2.72727272727288e-01}}
484: ,b[5][5] = {{2.95489397443992e-01 , -1.18481693910097e+00 , -6.68029812659953e-01 , 1.00716687860943e+00 , 2.72727272727288e-01},
485: {0.00000000000000e+00 , 1.11022302462516e-16 , -2.22044604925031e-16 , 0.00000000000000e+00 , 1.00000000000000e+00},
486: {-4.05882503986005e+00, -4.00924006567769e+00, -1.38930610972481e+00, 4.45223930308488e+00 , 6.32331093108427e-01},
487: {8.35690179937017e+00 , -2.26640927349732e+00 , 6.86647884973826e+00 , -5.22595158025740e+00 , 4.50893068837431e+00},
488: {1.27656267027479e+01 , 2.80882153840821e+00 , 8.91173096522890e+00 , -1.07936444078906e+01 , 4.82534148988854e+00}}
489: ,u[5][5] = {{1.00000000000000e+00 , -7.27272727273551e-02 , -3.45454545454419e-02 , -4.12121212119565e-03 ,-2.96969696964014e-04},
490: {1.00000000000000e+00 , 2.31252881006154e-01 , -8.29487834416481e-03 , -9.07191207681020e-03 ,-1.70378403743473e-03},
491: {1.00000000000000e+00 , 1.16925777880663e+00 , 3.59268562942635e-02 , -4.09013451730615e-02 ,-1.02411119670164e-02},
492: {1.00000000000000e+00 , 1.02634463704356e+00 , 1.59375044913405e-01 , 1.89673015035370e-03 ,-4.89987231897569e-03},
493: {1.00000000000000e+00 , 1.27746320298021e+00 , 2.37186008132728e-01 , -8.28694373940065e-02 ,-5.34396510196430e-02}}
494: ,v[5][5] = {{1.00000000000000e+00 , 1.27746320298021e+00 , 2.37186008132728e-01 , -8.28694373940065e-02 ,-5.34396510196430e-02},
495: {0.00000000000000e+00 , -1.77635683940025e-15 , -1.99840144432528e-15 , -9.99200722162641e-16 ,-3.33066907387547e-16},
496: {0.00000000000000e+00 , 4.37280081906924e+00 , 5.49221645016377e-02 , -8.88913177394943e-02 , 1.12879077989154e-01},
497: {0.00000000000000e+00 , -1.22399504837280e+01 , -5.21287338448645e+00 , -8.03952325565291e-01 , 4.60298678047147e-01},
498: {0.00000000000000e+00 , -1.85178762883829e+01 , -5.21411849862624e+00 , -1.04283436528809e+00 , 7.49030161063651e-01}};
499: TSGLSchemeCreate(4,4,5,5,c,*a,*b,*u,*v,&gl->schemes[gl->nschemes++]);
500: }
501: {
502: /* p=q=5, r=s=6;
503: irks(1/3,0,[1:6]/6,...
504: [-0.0489 0.4228 -0.8814 0.9021],...
505: [-0.3474 -0.6617 0.6294 0.2129
506: 0.0044 -0.4256 -0.1427 -0.8936
507: -0.8267 0.4821 0.1371 -0.2557
508: -0.4426 -0.3855 -0.7514 0.3014])
509: */
510: const PetscScalar c[6] = {1./6, 2./6, 3./6, 4./6, 5./6, 1.}
511: ,a[6][6] = {{ 3.33333333333940e-01, 0 , 0 , 0 , 0 , 0 },
512: { -8.64423857333350e-02, 3.33333333332888e-01, 0 , 0 , 0 , 0 },
513: { -2.16850174258252e+00, -2.23619072028839e+00, 3.33333333335204e-01, 0 , 0 , 0 },
514: { -4.73160970138997e+00, -3.89265344629268e+00, -2.76318716520933e-01, 3.33333333335759e-01, 0 , 0 },
515: { -6.75187540297338e+00, -7.90756533769377e+00, 7.90245051802259e-01, -4.48352364517632e-01, 3.33333333328483e-01, 0 },
516: { -4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00, -2.23969848002481e+00, 6.62807710124007e-01, 3.33333333335440e-01}}
517: ,b[6][6] = {{ -4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00, -2.23969848002481e+00, 6.62807710124007e-01, 3.33333333335440e-01},
518: { -8.88178419700125e-16, 4.44089209850063e-16, -1.54737334057131e-15, -8.88178419700125e-16, 0.00000000000000e+00, 1.00000000000001e+00},
519: { -2.87780425770651e+01, -1.13520448264971e+01, 2.62002318943161e+01, 2.56943874812797e+01, -3.06702268304488e+01, 6.68067773510103e+00},
520: { 5.47971245256474e+01, 6.80366875868284e+01, -6.50952588861999e+01, -8.28643975339097e+01, 8.17416943896414e+01, -1.17819043489036e+01},
521: { -2.33332114788869e+02, 6.12942539462634e+01, -4.91850135865944e+01, 1.82716844135480e+02, -1.29788173979395e+02, 3.09968095651099e+01},
522: { -1.72049132343751e+02, 8.60194713593999e+00, 7.98154219170200e-01, 1.50371386053218e+02, -1.18515423962066e+02, 2.50898277784663e+01}}
523: ,u[6][6] = {{ 1.00000000000000e+00, -1.66666666666870e-01, -4.16666666664335e-02, -3.85802469124815e-03, -2.25051440302250e-04, -9.64506172339142e-06},
524: { 1.00000000000000e+00, 8.64423857327162e-02, -4.11484912671353e-02, -1.11450903217645e-02, -1.47651050487126e-03, -1.34395070766826e-04},
525: { 1.00000000000000e+00, 4.57135912953434e+00, 1.06514719719137e+00, 1.33517564218007e-01, 1.11365952968659e-02, 6.12382756769504e-04},
526: { 1.00000000000000e+00, 9.23391519753404e+00, 2.22431212392095e+00, 2.91823807741891e-01, 2.52058456411084e-02, 1.22800542949647e-03},
527: { 1.00000000000000e+00, 1.48175480533865e+01, 3.73439117461835e+00, 5.14648336541804e-01, 4.76430038853402e-02, 2.56798515502156e-03},
528: { 1.00000000000000e+00, 1.50512347758335e+01, 4.10099701165164e+00, 5.66039141003603e-01, 3.91213893800891e-02, -2.99136269067853e-03}}
529: ,v[6][6] = {{ 1.00000000000000e+00, 1.50512347758335e+01, 4.10099701165164e+00, 5.66039141003603e-01, 3.91213893800891e-02, -2.99136269067853e-03},
530: { 0.00000000000000e+00, -4.88498130835069e-15, -6.43929354282591e-15, -3.55271367880050e-15, -1.22124532708767e-15, -3.12250225675825e-16},
531: { 0.00000000000000e+00, 1.22250171233141e+01, -1.77150760606169e+00, 3.54516769879390e-01, 6.22298845883398e-01, 2.31647447450276e-01},
532: { 0.00000000000000e+00, -4.48339457331040e+01, -3.57363126641880e-01, 5.18750173123425e-01, 6.55727990241799e-02, 1.63175368287079e-01},
533: { 0.00000000000000e+00, 1.37297394708005e+02, -1.60145272991317e+00, -5.05319555199441e+00, 1.55328940390990e-01, 9.16629423682464e-01},
534: { 0.00000000000000e+00, 1.05703241119022e+02, -1.16610260983038e+00, -2.99767252773859e+00, -1.13472315553890e-01, 1.09742849254729e+00}};
535: TSGLSchemeCreate(5,5,6,6,c,*a,*b,*u,*v,&gl->schemes[gl->nschemes++]);
536: }
537: return(0);
538: }
539: EXTERN_C_END
543: /*@C
544: TSGLSetType - sets the class of general linear method to use for time-stepping
546: Collective on TS
548: Input Parameters:
549: + ts - the TS context
550: - type - a method
552: Options Database Key:
553: . -ts_gl_type <type> - sets the method, use -help for a list of available method (e.g. irks)
555: Notes:
556: See "petsc/include/petscts.h" for available methods (for instance)
557: . TSGL_IRKS - Diagonally implicit methods with inherent Runge-Kutta stability (for stiff problems)
559: Normally, it is best to use the TSSetFromOptions() command and
560: then set the TSGL type from the options database rather than by using
561: this routine. Using the options database provides the user with
562: maximum flexibility in evaluating the many different solvers.
563: The TSGLSetType() routine is provided for those situations where it
564: is necessary to set the timestepping solver independently of the
565: command line or options database. This might be the case, for example,
566: when the choice of solver changes during the execution of the
567: program, and the user's application is taking responsibility for
568: choosing the appropriate method.
570: Level: intermediate
572: .keywords: TS, TSGL, set, type
573: @*/
574: PetscErrorCode TSGLSetType(TS ts,const TSGLType type)
575: {
581: PetscTryMethod(ts,"TSGLSetType_C",(TS,const TSGLType),(ts,type));
582: return(0);
583: }
587: /*@C
588: TSGLSetAcceptType - sets the acceptance test
590: Time integrators that need to control error must have the option to reject a time step based on local error
591: estimates. This function allows different schemes to be set.
593: Logically Collective on TS
595: Input Parameters:
596: + ts - the TS context
597: - type - the type
599: Options Database Key:
600: . -ts_gl_accept_type <type> - sets the method used to determine whether to accept or reject a step
602: Level: intermediate
604: .seealso: TS, TSGL, TSGLAcceptRegisterDynamic(), TSGLAdapt, set type
605: @*/
606: PetscErrorCode TSGLSetAcceptType(TS ts,const TSGLAcceptType type)
607: {
613: PetscTryMethod(ts,"TSGLSetAcceptType_C",(TS,const TSGLAcceptType),(ts,type));
614: return(0);
615: }
619: /*@C
620: TSGLGetAdapt - gets the TSGLAdapt object from the TS
622: Not Collective
624: Input Parameter:
625: . ts - the TS context
627: Output Parameter:
628: . adapt - the TSGLAdapt context
630: Notes:
631: This allows the user set options on the TSGLAdapt object. Usually it is better to do this using the options
632: database, so this function is rarely needed.
634: Level: advanced
636: .seealso: TSGLAdapt, TSGLAdaptRegisterDynamic()
637: @*/
638: PetscErrorCode TSGLGetAdapt(TS ts,TSGLAdapt *adapt)
639: {
645: PetscUseMethod(ts,"TSGLGetAdapt_C",(TS,TSGLAdapt*),(ts,adapt));
646: return(0);
647: }
649: EXTERN_C_BEGIN
652: PetscErrorCode TSGLAccept_Always(TS ts,PetscReal tleft,PetscReal h,const PetscReal enorms[],PetscBool *accept)
653: {
655: *accept = PETSC_TRUE;
656: return(0);
657: }
658: EXTERN_C_END
662: static PetscErrorCode TSGLUpdateWRMS(TS ts)
663: {
664: TS_GL *gl = (TS_GL*)ts->data;
666: PetscScalar *x,*w;
667: PetscInt n,i;
670: VecGetArray(gl->X[0],&x);
671: VecGetArray(gl->W,&w);
672: VecGetLocalSize(gl->W,&n);
673: for (i=0; i<n; i++) {
674: w[i] = 1./(gl->wrms_atol + gl->wrms_rtol*PetscAbsScalar(x[i]));
675: }
676: VecRestoreArray(gl->X[0],&x);
677: VecRestoreArray(gl->W,&w);
678: return(0);
679: }
683: static PetscErrorCode TSGLVecNormWRMS(TS ts,Vec X,PetscReal *nrm)
684: {
685: TS_GL *gl = (TS_GL*)ts->data;
687: PetscScalar *x,*w;
688: PetscReal sum = 0.0,gsum;
689: PetscInt n,N,i;
692: VecGetArray(X,&x);
693: VecGetArray(gl->W,&w);
694: VecGetLocalSize(gl->W,&n);
695: for (i=0; i<n; i++) {
696: sum += PetscAbsScalar(PetscSqr(x[i]*w[i]));
697: }
698: VecRestoreArray(X,&x);
699: VecRestoreArray(gl->W,&w);
700: MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,((PetscObject)ts)->comm);
701: VecGetSize(gl->W,&N);
702: *nrm = PetscSqrtReal(gsum/(1.*N));
703: return(0);
704: }
706: EXTERN_C_BEGIN
709: PetscErrorCode TSGLSetType_GL(TS ts,const TSGLType type)
710: {
711: PetscErrorCode ierr,(*r)(TS);
712: PetscBool same;
713: TS_GL *gl = (TS_GL*)ts->data;
716: if (gl->type_name[0]) {
717: PetscStrcmp(gl->type_name,type,&same);
718: if (same) return(0);
719: (*gl->Destroy)(gl);
720: }
722: PetscFListFind(TSGLList,((PetscObject)ts)->comm,type,PETSC_TRUE,(PetscVoidStarFunction)&r);
723: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGL type \"%s\" given",type);
724: (*r)(ts);
725: PetscStrcpy(gl->type_name,type);
726: return(0);
727: }
731: PetscErrorCode TSGLSetAcceptType_GL(TS ts,const TSGLAcceptType type)
732: {
734: TSGLAcceptFunction r;
735: TS_GL *gl = (TS_GL*)ts->data;
738: PetscFListFind(TSGLAcceptList,((PetscObject)ts)->comm,type,PETSC_TRUE,(PetscVoidStarFunction)&r);
739: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLAccept type \"%s\" given",type);
740: gl->Accept = r;
741: PetscStrncpy(gl->accept_name,type,sizeof(gl->accept_name));
742: return(0);
743: }
747: PetscErrorCode TSGLGetAdapt_GL(TS ts,TSGLAdapt *adapt)
748: {
750: TS_GL *gl = (TS_GL*)ts->data;
753: if (!gl->adapt) {
754: TSGLAdaptCreate(((PetscObject)ts)->comm,&gl->adapt);
755: PetscObjectIncrementTabLevel((PetscObject)gl->adapt,(PetscObject)ts,1);
756: PetscLogObjectParent(ts,gl->adapt);
757: }
758: *adapt = gl->adapt;
759: return(0);
760: }
761: EXTERN_C_END
765: static PetscErrorCode TSGLChooseNextScheme(TS ts,PetscReal h,const PetscReal hmnorm[],PetscInt *next_scheme,PetscReal *next_h,PetscBool *finish)
766: {
768: TS_GL *gl = (TS_GL*)ts->data;
769: PetscInt i,n,cur_p,cur,next_sc,candidates[64],orders[64];
770: PetscReal errors[64],costs[64],tleft;
773: cur = -1;
774: cur_p = gl->schemes[gl->current_scheme]->p;
775: tleft = ts->max_time - (ts->ptime + ts->time_step);
776: for (i=0,n=0; i<gl->nschemes; i++) {
777: TSGLScheme sc = gl->schemes[i];
778: if (sc->p < gl->min_order || gl->max_order < sc->p) continue;
779: if (sc->p == cur_p - 1) {
780: errors[n] = PetscAbsScalar(sc->alpha[0])*hmnorm[0];
781: } else if (sc->p == cur_p) {
782: errors[n] = PetscAbsScalar(sc->alpha[0])*hmnorm[1];
783: } else if (sc->p == cur_p+1) {
784: errors[n] = PetscAbsScalar(sc->alpha[0])*(hmnorm[2]+hmnorm[3]);
785: } else continue;
786: candidates[n] = i;
787: orders[n] = PetscMin(sc->p,sc->q); /* order of global truncation error */
788: costs[n] = sc->s; /* estimate the cost as the number of stages */
789: if (i == gl->current_scheme) cur = n;
790: n++;
791: }
792: if (cur < 0 || gl->nschemes <= cur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Current scheme not found in scheme list");
793: TSGLAdaptChoose(gl->adapt,n,orders,errors,costs,cur,h,tleft,&next_sc,next_h,finish);
794: *next_scheme = candidates[next_sc];
795: PetscInfo7(ts,"Adapt chose scheme %d (%d,%d,%d,%d) with step size %6.2e, finish=%d\n",*next_scheme,gl->schemes[*next_scheme]->p,gl->schemes[*next_scheme]->q,gl->schemes[*next_scheme]->r,gl->schemes[*next_scheme]->s,*next_h,*finish);
796: return(0);
797: }
801: static PetscErrorCode TSGLGetMaxSizes(TS ts,PetscInt *max_r,PetscInt *max_s)
802: {
803: TS_GL *gl = (TS_GL*)ts->data;
806: *max_r = gl->schemes[gl->nschemes-1]->r;
807: *max_s = gl->schemes[gl->nschemes-1]->s;
808: return(0);
809: }
813: static PetscErrorCode TSSolve_GL(TS ts)
814: {
815: TS_GL *gl = (TS_GL*)ts->data;
816: PetscInt i,k,its,lits,max_r,max_s;
817: PetscBool final_step,finish;
818: SNESConvergedReason snesreason;
819: PetscErrorCode ierr;
822: TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);
824: TSGLGetMaxSizes(ts,&max_r,&max_s);
825: VecCopy(ts->vec_sol,gl->X[0]);
826: for (i=1; i<max_r; i++) {
827: VecZeroEntries(gl->X[i]);
828: }
829: TSGLUpdateWRMS(ts);
831: if (0) {
832: /* Find consistent initial data for DAE */
833: gl->stage_time = ts->ptime + ts->time_step;
834: gl->shift = 1./ts->time_step;
835: gl->stage = 0;
836: VecCopy(ts->vec_sol,gl->Z);
837: VecCopy(ts->vec_sol,gl->Y);
838: SNESSolve(ts->snes,PETSC_NULL,gl->Y);
839: SNESGetIterationNumber(ts->snes,&its);
840: SNESGetLinearSolveIterations(ts->snes,&lits);
841: SNESGetConvergedReason(ts->snes,&snesreason);
842: ts->snes_its += its; ts->ksp_its += lits;
843: if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) {
844: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
845: PetscInfo2(ts,"Step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
846: return(0);
847: }
848: }
850: if (gl->current_scheme < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"A starting scheme has not been provided");
852: for (k=0,final_step=PETSC_FALSE,finish=PETSC_FALSE; k<ts->max_steps && !finish; k++) {
853: PetscInt j,r,s,next_scheme = 0,rejections;
854: PetscReal h,hmnorm[4],enorm[3],next_h;
855: PetscBool accept;
856: const PetscScalar *c,*a,*u;
857: Vec *X,*Ydot,Y;
858: TSGLScheme scheme = gl->schemes[gl->current_scheme];
860: r = scheme->r; s = scheme->s;
861: c = scheme->c;
862: a = scheme->a; u = scheme->u;
863: h = ts->time_step;
864: X = gl->X; Ydot = gl->Ydot; Y = gl->Y;
866: if (ts->ptime > ts->max_time) break;
868: /*
869: We only call PreStep at the start of each STEP, not each STAGE. This is because it is
870: possible to fail (have to restart a step) after multiple stages.
871: */
872: TSPreStep(ts);
874: rejections = 0;
875: while (1) {
876: for (i=0; i<s; i++) {
877: PetscScalar shift = gl->shift = 1./PetscRealPart(h*a[i*s+i]);
878: gl->stage = i;
879: gl->stage_time = ts->ptime + PetscRealPart(c[i])*h;
881: /*
882: * Stage equation: Y = h A Y' + U X
883: * We assume that A is lower-triangular so that we can solve the stages (Y,Y') sequentially
884: * Build the affine vector z_i = -[1/(h a_ii)](h sum_j a_ij y'_j + sum_j u_ij x_j)
885: * Then y'_i = z + 1/(h a_ii) y_i
886: */
887: VecZeroEntries(gl->Z);
888: for (j=0; j<r; j++) {
889: VecAXPY(gl->Z,-shift*u[i*r+j],X[j]);
890: }
891: for (j=0; j<i; j++) {
892: VecAXPY(gl->Z,-shift*h*a[i*s+j],Ydot[j]);
893: }
894: /* Note: Z is used within function evaluation, Ydot = Z + shift*Y */
896: /* Compute an estimate of Y to start Newton iteration */
897: if (gl->extrapolate) {
898: if (i==0) {
899: /* Linear extrapolation on the first stage */
900: VecWAXPY(Y,c[i]*h,X[1],X[0]);
901: } else {
902: /* Linear extrapolation from the last stage */
903: VecAXPY(Y,(c[i]-c[i-1])*h,Ydot[i-1]);
904: }
905: } else if (i==0) { /* Directly use solution from the last step, otherwise reuse the last stage (do nothing) */
906: VecCopy(X[0],Y);
907: }
909: /* Solve this stage (Ydot[i] is computed during function evaluation) */
910: SNESSolve(ts->snes,PETSC_NULL,Y);
911: SNESGetIterationNumber(ts->snes,&its);
912: SNESGetLinearSolveIterations(ts->snes,&lits);
913: SNESGetConvergedReason(ts->snes,&snesreason);
914: ts->snes_its += its; ts->ksp_its += lits;
915: if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) {
916: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
917: PetscInfo2(ts,"Step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
918: return(0);
919: }
920: }
922: gl->stage_time = ts->ptime + ts->time_step;
924: (*gl->EstimateHigherMoments)(scheme,h,Ydot,gl->X,gl->himom);
925: /* hmnorm[i] = h^{p+i}x^{(p+i)} with i=0,1,2; hmnorm[3] = h^{p+2}(dx'/dx) x^{(p+1)} */
926: for (i=0; i<3; i++) {
927: TSGLVecNormWRMS(ts,gl->himom[i],&hmnorm[i+1]);
928: }
929: enorm[0] = PetscRealPart(scheme->alpha[0])*hmnorm[1];
930: enorm[1] = PetscRealPart(scheme->beta[0]) *hmnorm[2];
931: enorm[2] = PetscRealPart(scheme->gamma[0])*hmnorm[3];
932: (*gl->Accept)(ts,ts->max_time-gl->stage_time,h,enorm,&accept);
933: if (accept) goto accepted;
934: rejections++;
935: PetscInfo3(ts,"Step %D (t=%g) not accepted, rejections=%D\n",k,gl->stage_time,rejections);
936: if (rejections > gl->max_step_rejections) break;
937: /*
938: There are lots of reasons why a step might be rejected, including solvers not converging and other factors that
939: TSGLChooseNextScheme does not support. Additionally, the error estimates may be very screwed up, so I'm not
940: convinced that it's safe to just compute a new error estimate using the same interface as the current adaptor
941: (the adaptor interface probably has to change). Here we make an arbitrary and naive choice. This assumes that
942: steps were written in Nordsieck form. The "correct" method would be to re-complete the previous time step with
943: the correct "next" step size. It is unclear to me whether the present ad-hoc method of rescaling X is stable.
944: */
945: h *= 0.5;
946: for (i=1; i<scheme->r; i++) {
947: VecScale(X[i],Pow(0.5,i));
948: }
949: }
950: SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Time step %D (t=%g) not accepted after %D failures\n",k,gl->stage_time,rejections);
952: accepted:
953: /* This term is not error, but it *would* be the leading term for a lower order method */
954: TSGLVecNormWRMS(ts,gl->X[scheme->r-1],&hmnorm[0]);
955: /* Correct scaling so that these are equivalent to norms of the Nordsieck vectors */
957: PetscInfo4(ts,"Last moment norm %10.2e, estimated error norms %10.2e %10.2e %10.2e\n",hmnorm[0],enorm[0],enorm[1],enorm[2]);
958: if (!final_step) {
959: TSGLChooseNextScheme(ts,h,hmnorm,&next_scheme,&next_h,&final_step);
960: } else {
961: /* Dummy values to complete the current step in a consistent manner */
962: next_scheme = gl->current_scheme;
963: next_h = h;
964: finish = PETSC_TRUE;
965: }
967: X = gl->Xold;
968: gl->Xold = gl->X;
969: gl->X = X;
970: (*gl->CompleteStep)(scheme,h,gl->schemes[next_scheme],next_h,Ydot,gl->Xold,gl->X);
972: TSGLUpdateWRMS(ts);
974: /* Post the solution for the user, we could avoid this copy with a small bit of cleverness */
975: VecCopy(gl->X[0],ts->vec_sol);
976: ts->ptime += h;
977: ts->steps++;
979: TSPostStep(ts);
980: TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);
982: gl->current_scheme = next_scheme;
983: ts->time_step = next_h;
984: }
985: return(0);
986: }
988: /*------------------------------------------------------------*/
992: static PetscErrorCode TSReset_GL(TS ts)
993: {
994: TS_GL *gl = (TS_GL*)ts->data;
995: PetscInt max_r,max_s;
996: PetscErrorCode ierr;
999: if (gl->setupcalled) {
1000: TSGLGetMaxSizes(ts,&max_r,&max_s);
1001: VecDestroyVecs(max_r,&gl->Xold);
1002: VecDestroyVecs(max_r,&gl->X);
1003: VecDestroyVecs(max_s,&gl->Ydot);
1004: VecDestroyVecs(3,&gl->himom);
1005: VecDestroy(&gl->W);
1006: VecDestroy(&gl->Y);
1007: VecDestroy(&gl->Z);
1008: }
1009: gl->setupcalled = PETSC_FALSE;
1010: return(0);
1011: }
1015: static PetscErrorCode TSDestroy_GL(TS ts)
1016: {
1017: TS_GL *gl = (TS_GL*)ts->data;
1018: PetscErrorCode ierr;
1021: TSReset_GL(ts);
1022: if (gl->adapt) {TSGLAdaptDestroy(&gl->adapt);}
1023: if (gl->Destroy) {(*gl->Destroy)(gl);}
1024: PetscFree(ts->data);
1025: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSGLSetType_C", "",PETSC_NULL);
1026: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSGLSetAcceptType_C","",PETSC_NULL);
1027: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSGLGetAdapt_C", "",PETSC_NULL);
1028: return(0);
1029: }
1031: /*
1032: This defines the nonlinear equation that is to be solved with SNES
1033: g(x) = f(t,x,z+shift*x) = 0
1034: */
1037: static PetscErrorCode SNESTSFormFunction_GL(SNES snes,Vec x,Vec f,TS ts)
1038: {
1039: TS_GL *gl = (TS_GL*)ts->data;
1040: PetscErrorCode ierr;
1043: VecWAXPY(gl->Ydot[gl->stage],gl->shift,x,gl->Z);
1044: TSComputeIFunction(ts,gl->stage_time,x,gl->Ydot[gl->stage],f,PETSC_FALSE);
1045: return(0);
1046: }
1050: static PetscErrorCode SNESTSFormJacobian_GL(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *str,TS ts)
1051: {
1052: TS_GL *gl = (TS_GL*)ts->data;
1053: PetscErrorCode ierr;
1056: /* gl->Xdot will have already been computed in SNESTSFormFunction_GL */
1057: TSComputeIJacobian(ts,gl->stage_time,x,gl->Ydot[gl->stage],gl->shift,A,B,str,PETSC_FALSE);
1058: return(0);
1059: }
1064: static PetscErrorCode TSSetUp_GL(TS ts)
1065: {
1066: TS_GL *gl = (TS_GL*)ts->data;
1067: PetscInt max_r,max_s;
1068: PetscErrorCode ierr;
1071: gl->setupcalled = PETSC_TRUE;
1072: TSGLGetMaxSizes(ts,&max_r,&max_s);
1073: VecDuplicateVecs(ts->vec_sol,max_r,&gl->X);
1074: VecDuplicateVecs(ts->vec_sol,max_r,&gl->Xold);
1075: VecDuplicateVecs(ts->vec_sol,max_s,&gl->Ydot);
1076: VecDuplicateVecs(ts->vec_sol,3,&gl->himom);
1077: VecDuplicate(ts->vec_sol,&gl->W);
1078: VecDuplicate(ts->vec_sol,&gl->Y);
1079: VecDuplicate(ts->vec_sol,&gl->Z);
1081: /* Default acceptance tests and adaptivity */
1082: if (!gl->Accept) {TSGLSetAcceptType(ts,TSGLACCEPT_ALWAYS);}
1083: if (!gl->adapt) {TSGLGetAdapt(ts,&gl->adapt);}
1085: if (gl->current_scheme < 0) {
1086: PetscInt i;
1087: for (i=0; ; i++) {
1088: if (gl->schemes[i]->p == gl->start_order) break;
1089: if (i+1 == gl->nschemes) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No schemes available with requested start order %d",i);
1090: }
1091: gl->current_scheme = i;
1092: }
1093: return(0);
1094: }
1095: /*------------------------------------------------------------*/
1099: static PetscErrorCode TSSetFromOptions_GL(TS ts)
1100: {
1101: TS_GL *gl = (TS_GL*)ts->data;
1102: char tname[256] = TSGL_IRKS,completef[256] = "rescale-and-modify";
1106: PetscOptionsHead("General Linear ODE solver options");
1107: {
1108: PetscBool flg;
1109: PetscOptionsList("-ts_gl_type","Type of GL method","TSGLSetType",TSGLList,gl->type_name[0]?gl->type_name:tname,tname,sizeof(tname),&flg);
1110: if (flg || !gl->type_name[0]) {
1111: TSGLSetType(ts,tname);
1112: }
1113: PetscOptionsInt("-ts_gl_max_step_rejections","Maximum number of times to attempt a step","None",gl->max_step_rejections,&gl->max_step_rejections,PETSC_NULL);
1114: PetscOptionsInt("-ts_gl_max_order","Maximum order to try","TSGLSetMaxOrder",gl->max_order,&gl->max_order,PETSC_NULL);
1115: PetscOptionsInt("-ts_gl_min_order","Minimum order to try","TSGLSetMinOrder",gl->min_order,&gl->min_order,PETSC_NULL);
1116: PetscOptionsInt("-ts_gl_start_order","Initial order to try","TSGLSetMinOrder",gl->start_order,&gl->start_order,PETSC_NULL);
1117: PetscOptionsEnum("-ts_gl_error_direction","Which direction to look when estimating error","TSGLSetErrorDirection",TSGLErrorDirections,(PetscEnum)gl->error_direction,(PetscEnum*)&gl->error_direction,PETSC_NULL);
1118: PetscOptionsBool("-ts_gl_extrapolate","Extrapolate stage solution from previous solution (sometimes unstable)","TSGLSetExtrapolate",gl->extrapolate,&gl->extrapolate,PETSC_NULL);
1119: PetscOptionsReal("-ts_gl_atol","Absolute tolerance","TSGLSetTolerances",gl->wrms_atol,&gl->wrms_atol,PETSC_NULL);
1120: PetscOptionsReal("-ts_gl_rtol","Relative tolerance","TSGLSetTolerances",gl->wrms_rtol,&gl->wrms_rtol,PETSC_NULL);
1121: PetscOptionsString("-ts_gl_complete","Method to use for completing the step","none",completef,completef,sizeof completef,&flg);
1122: if (flg) {
1123: PetscBool match1,match2;
1124: PetscStrcmp(completef,"rescale",&match1);
1125: PetscStrcmp(completef,"rescale-and-modify",&match2);
1126: if (match1) gl->CompleteStep = TSGLCompleteStep_Rescale;
1127: else if (match2) gl->CompleteStep = TSGLCompleteStep_RescaleAndModify;
1128: else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"%s",completef);
1129: }
1130: {
1131: char type[256] = TSGLACCEPT_ALWAYS;
1132: PetscOptionsList("-ts_gl_accept_type","Method to use for determining whether to accept a step","TSGLSetAcceptType",TSGLAcceptList,gl->accept_name[0]?gl->accept_name:type,type,sizeof type,&flg);
1133: if (flg || !gl->accept_name[0]) {
1134: TSGLSetAcceptType(ts,type);
1135: }
1136: }
1137: SNESSetFromOptions(ts->snes);
1138: {
1139: TSGLAdapt adapt;
1140: TSGLGetAdapt(ts,&adapt);
1141: TSGLAdaptSetFromOptions(adapt);
1142: }
1143: }
1144: PetscOptionsTail();
1145: return(0);
1146: }
1150: static PetscErrorCode TSView_GL(TS ts,PetscViewer viewer)
1151: {
1152: TS_GL *gl = (TS_GL*)ts->data;
1153: PetscInt i;
1154: PetscBool iascii,details;
1155: PetscErrorCode ierr;
1158: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1159: if (iascii) {
1160: PetscViewerASCIIPrintf(viewer," min order %D, max order %D, current order %D\n",gl->min_order,gl->max_order,gl->schemes[gl->current_scheme]->p);
1161: PetscViewerASCIIPrintf(viewer," Error estimation: %s\n",TSGLErrorDirections[gl->error_direction]);
1162: PetscViewerASCIIPrintf(viewer," Extrapolation: %s\n",gl->extrapolate?"yes":"no");
1163: PetscViewerASCIIPrintf(viewer," Acceptance test: %s\n",gl->accept_name[0]?gl->accept_name:"(not yet set)");
1164: PetscViewerASCIIPushTab(viewer);
1165: TSGLAdaptView(gl->adapt,viewer);
1166: PetscViewerASCIIPopTab(viewer);
1167: PetscViewerASCIIPrintf(viewer," type: %s\n",gl->type_name[0]?gl->type_name:"(not yet set)");
1168: PetscViewerASCIIPrintf(viewer,"Schemes within family (%d):\n",gl->nschemes);
1169: details = PETSC_FALSE;
1170: PetscOptionsGetBool(((PetscObject)ts)->prefix,"-ts_gl_view_detailed",&details,PETSC_NULL);
1171: PetscViewerASCIIPushTab(viewer);
1172: for (i=0; i<gl->nschemes; i++) {
1173: TSGLSchemeView(gl->schemes[i],details,viewer);
1174: }
1175: if (gl->View) {
1176: (*gl->View)(gl,viewer);
1177: }
1178: PetscViewerASCIIPopTab(viewer);
1179: }
1180: SNESView(ts->snes,viewer);
1181: return(0);
1182: }
1186: /*@C
1187: TSGLRegister - see TSGLRegisterDynamic()
1189: Level: advanced
1190: @*/
1191: PetscErrorCode TSGLRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(TS))
1192: {
1194: char fullname[PETSC_MAX_PATH_LEN];
1197: PetscFListConcat(path,name,fullname);
1198: PetscFListAdd(&TSGLList,sname,fullname,(void(*)(void))function);
1199: return(0);
1200: }
1204: /*@C
1205: TSGLAcceptRegister - see TSGLAcceptRegisterDynamic()
1207: Level: advanced
1208: @*/
1209: PetscErrorCode TSGLAcceptRegister(const char sname[],const char path[],const char name[],TSGLAcceptFunction function)
1210: {
1212: char fullname[PETSC_MAX_PATH_LEN];
1215: PetscFListConcat(path,name,fullname);
1216: PetscFListAdd(&TSGLAcceptList,sname,fullname,(void(*)(void))function);
1217: return(0);
1218: }
1222: /*@C
1223: TSGLRegisterAll - Registers all of the general linear methods in TSGL
1225: Not Collective
1227: Level: advanced
1229: .keywords: TS, TSGL, register, all
1231: .seealso: TSGLRegisterDestroy()
1232: @*/
1233: PetscErrorCode TSGLRegisterAll(const char path[])
1234: {
1238: if (TSGLRegisterAllCalled) return(0);
1239: TSGLRegisterAllCalled = PETSC_TRUE;
1241: TSGLRegisterDynamic(TSGL_IRKS,path,"TSGLCreate_IRKS",TSGLCreate_IRKS);
1242: TSGLAcceptRegisterDynamic(TSGLACCEPT_ALWAYS,path,"TSGLAccept_Always",TSGLAccept_Always);
1243: return(0);
1244: }
1248: /*@C
1249: TSGLRegisterDestroy - Frees the list of schemes that were registered by TSGLRegister()/TSGLRegisterDynamic().
1251: Not Collective
1253: Level: advanced
1255: .keywords: TSGL, register, destroy
1256: .seealso: TSGLRegister(), TSGLRegisterAll(), TSGLRegisterDynamic()
1257: @*/
1258: PetscErrorCode TSGLRegisterDestroy(void)
1259: {
1263: PetscFListDestroy(&TSGLList);
1264: TSGLRegisterAllCalled = PETSC_FALSE;
1265: return(0);
1266: }
1271: /*@C
1272: TSGLInitializePackage - This function initializes everything in the TSGL package. It is called
1273: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_GL()
1274: when using static libraries.
1276: Input Parameter:
1277: path - The dynamic library path, or PETSC_NULL
1279: Level: developer
1281: .keywords: TS, TSGL, initialize, package
1282: .seealso: PetscInitialize()
1283: @*/
1284: PetscErrorCode TSGLInitializePackage(const char path[])
1285: {
1289: if (TSGLPackageInitialized) return(0);
1290: TSGLPackageInitialized = PETSC_TRUE;
1291: TSGLRegisterAll(PETSC_NULL);
1292: PetscRegisterFinalize(TSGLFinalizePackage);
1293: return(0);
1294: }
1298: /*@C
1299: TSGLFinalizePackage - This function destroys everything in the TSGL package. It is
1300: called from PetscFinalize().
1302: Level: developer
1304: .keywords: Petsc, destroy, package
1305: .seealso: PetscFinalize()
1306: @*/
1307: PetscErrorCode TSGLFinalizePackage(void)
1308: {
1310: TSGLPackageInitialized = PETSC_FALSE;
1311: TSGLRegisterAllCalled = PETSC_FALSE;
1312: TSGLList = PETSC_NULL;
1313: TSGLAcceptList = PETSC_NULL;
1314: return(0);
1315: }
1317: /* ------------------------------------------------------------ */
1318: /*MC
1319: TSGL - DAE solver using implicit General Linear methods
1321: These methods contain Runge-Kutta and multistep schemes as special cases. These special cases have some fundamental
1322: limitations. For example, diagonally implicit Runge-Kutta cannot have stage order greater than 1 which limits their
1323: applicability to very stiff systems. Meanwhile, multistep methods cannot be A-stable for order greater than 2 and BDF
1324: are not 0-stable for order greater than 6. GL methods can be A- and L-stable with arbitrarily high stage order and
1325: reliable error estimates for both 1 and 2 orders higher to facilitate adaptive step sizes and adaptive order schemes.
1326: All this is possible while preserving a singly diagonally implicit structure.
1328: Options database keys:
1329: + -ts_gl_type <type> - the class of general linear method (irks)
1330: . -ts_gl_rtol <tol> - relative error
1331: . -ts_gl_atol <tol> - absolute error
1332: . -ts_gl_min_order <p> - minimum order method to consider (default=1)
1333: . -ts_gl_max_order <p> - maximum order method to consider (default=3)
1334: . -ts_gl_start_order <p> - order of starting method (default=1)
1335: . -ts_gl_complete <method> - method to use for completing the step (rescale-and-modify or rescale)
1336: - -ts_adapt_type <method> - adaptive controller to use (none step both)
1338: Notes:
1339: This integrator can be applied to DAE.
1341: Diagonally implicit general linear (DIGL) methods are a generalization of diagonally implicit Runge-Kutta (DIRK).
1342: They are represented by the tableau
1344: .vb
1345: A | U
1346: -------
1347: B | V
1348: .ve
1350: combined with a vector c of abscissa. "Diagonally implicit" means that A is lower triangular.
1351: A step of the general method reads
1353: .vb
1354: [ Y ] = [A U] [ Y' ]
1355: [X^k] = [B V] [X^{k-1}]
1356: .ve
1358: where Y is the multivector of stage values, Y' is the multivector of stage derivatives, X^k is the Nordsieck vector of
1359: the solution at step k. The Nordsieck vector consists of the first r moments of the solution, given by
1361: .vb
1362: X = [x_0,x_1,...,x_{r-1}] = [x, h x', h^2 x'', ..., h^{r-1} x^{(r-1)} ]
1363: .ve
1365: If A is lower triangular, we can solve the stages (Y,Y') sequentially
1367: .vb
1368: y_i = h sum_{j=0}^{s-1} (a_ij y'_j) + sum_{j=0}^{r-1} u_ij x_j, i=0,...,{s-1}
1369: .ve
1371: and then construct the pieces to carry to the next step
1373: .vb
1374: xx_i = h sum_{j=0}^{s-1} b_ij y'_j + sum_{j=0}^{r-1} v_ij x_j, i=0,...,{r-1}
1375: .ve
1377: Note that when the equations are cast in implicit form, we are using the stage equation to define y'_i
1378: in terms of y_i and known stuff (y_j for j<i and x_j for all j).
1381: Error estimation
1383: At present, the most attractive GL methods for stiff problems are singly diagonally implicit schemes which posses
1384: Inherent Runge-Kutta Stability (IRKS). These methods have r=s, the number of items passed between steps is equal to
1385: the number of stages. The order and stage-order are one less than the number of stages. We use the error estimates
1386: in the 2007 paper which provide the following estimates
1388: .vb
1389: h^{p+1} X^{(p+1)} = phi_0^T Y' + [0 psi_0^T] Xold
1390: h^{p+2} X^{(p+2)} = phi_1^T Y' + [0 psi_1^T] Xold
1391: h^{p+2} (dx'/dx) X^{(p+1)} = phi_2^T Y' + [0 psi_2^T] Xold
1392: .ve
1394: These estimates are accurate to O(h^{p+3}).
1396: Changing the step size
1398: We use the generalized "rescale and modify" scheme, see equation (4.5) of the 2007 paper.
1400: Level: beginner
1402: References:
1403: John Butcher and Z. Jackieweicz and W. Wright, On error propagation in general linear methods for
1404: ordinary differential equations, Journal of Complexity, Vol 23 (4-6), 2007.
1406: John Butcher, Numerical methods for ordinary differential equations, second edition, Wiley, 2009.
1408: .seealso: TSCreate(), TS, TSSetType()
1410: M*/
1411: EXTERN_C_BEGIN
1414: PetscErrorCode TSCreate_GL(TS ts)
1415: {
1416: TS_GL *gl;
1420: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
1421: TSGLInitializePackage(PETSC_NULL);
1422: #endif
1424: PetscNewLog(ts,TS_GL,&gl);
1425: ts->data = (void*)gl;
1427: ts->ops->reset = TSReset_GL;
1428: ts->ops->destroy = TSDestroy_GL;
1429: ts->ops->view = TSView_GL;
1430: ts->ops->setup = TSSetUp_GL;
1431: ts->ops->solve = TSSolve_GL;
1432: ts->ops->setfromoptions = TSSetFromOptions_GL;
1433: ts->ops->snesfunction = SNESTSFormFunction_GL;
1434: ts->ops->snesjacobian = SNESTSFormJacobian_GL;
1436: gl->max_step_rejections = 1;
1437: gl->min_order = 1;
1438: gl->max_order = 3;
1439: gl->start_order = 1;
1440: gl->current_scheme = -1;
1441: gl->extrapolate = PETSC_FALSE;
1443: gl->wrms_atol = 1e-8;
1444: gl->wrms_rtol = 1e-5;
1446: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSGLSetType_C", "TSGLSetType_GL", &TSGLSetType_GL);
1447: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSGLSetAcceptType_C","TSGLSetAcceptType_GL",&TSGLSetAcceptType_GL);
1448: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSGLGetAdapt_C", "TSGLGetAdapt_GL", &TSGLGetAdapt_GL);
1449: return(0);
1450: }
1451: EXTERN_C_END