Actual source code: glle.c
petsc-3.9.4 2018-09-11
2: #include <../src/ts/impls/implicit/glle/glle.h>
3: #include <petscdm.h>
4: #include <petscblaslapack.h>
6: static const char *TSGLLEErrorDirections[] = {"FORWARD","BACKWARD","TSGLLEErrorDirection","TSGLLEERROR_",0};
7: static PetscFunctionList TSGLLEList;
8: static PetscFunctionList TSGLLEAcceptList;
9: static PetscBool TSGLLEPackageInitialized;
10: static PetscBool TSGLLERegisterAllCalled;
12: /* This function is pure */
13: static PetscScalar Factorial(PetscInt n)
14: {
15: PetscInt i;
16: if (n < 12) { /* Can compute with 32-bit integers */
17: PetscInt f = 1;
18: for (i=2; i<=n; i++) f *= i;
19: return (PetscScalar)f;
20: } else {
21: PetscScalar f = 1.;
22: for (i=2; i<=n; i++) f *= (PetscScalar)i;
23: return f;
24: }
25: }
27: /* This function is pure */
28: static PetscScalar CPowF(PetscScalar c,PetscInt p)
29: {
30: return PetscPowRealInt(PetscRealPart(c),p)/Factorial(p);
31: }
33: static PetscErrorCode TSGLLEGetVecs(TS ts,DM dm,Vec *Z,Vec *Ydotstage)
34: {
35: TS_GLLE *gl = (TS_GLLE*)ts->data;
39: if (Z) {
40: if (dm && dm != ts->dm) {
41: DMGetNamedGlobalVector(dm,"TSGLLE_Z",Z);
42: } else *Z = gl->Z;
43: }
44: if (Ydotstage) {
45: if (dm && dm != ts->dm) {
46: DMGetNamedGlobalVector(dm,"TSGLLE_Ydot",Ydotstage);
47: } else *Ydotstage = gl->Ydot[gl->stage];
48: }
49: return(0);
50: }
53: static PetscErrorCode TSGLLERestoreVecs(TS ts,DM dm,Vec *Z,Vec *Ydotstage)
54: {
58: if (Z) {
59: if (dm && dm != ts->dm) {
60: DMRestoreNamedGlobalVector(dm,"TSGLLE_Z",Z);
61: }
62: }
63: if (Ydotstage) {
65: if (dm && dm != ts->dm) {
66: DMRestoreNamedGlobalVector(dm,"TSGLLE_Ydot",Ydotstage);
67: }
68: }
69: return(0);
70: }
72: static PetscErrorCode DMCoarsenHook_TSGLLE(DM fine,DM coarse,void *ctx)
73: {
75: return(0);
76: }
78: static PetscErrorCode DMRestrictHook_TSGLLE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
79: {
80: TS ts = (TS)ctx;
82: Vec Ydot,Ydot_c;
85: TSGLLEGetVecs(ts,fine,NULL,&Ydot);
86: TSGLLEGetVecs(ts,coarse,NULL,&Ydot_c);
87: MatRestrict(restrct,Ydot,Ydot_c);
88: VecPointwiseMult(Ydot_c,rscale,Ydot_c);
89: TSGLLERestoreVecs(ts,fine,NULL,&Ydot);
90: TSGLLERestoreVecs(ts,coarse,NULL,&Ydot_c);
91: return(0);
92: }
94: static PetscErrorCode DMSubDomainHook_TSGLLE(DM dm,DM subdm,void *ctx)
95: {
97: return(0);
98: }
100: static PetscErrorCode DMSubDomainRestrictHook_TSGLLE(DM dm,VecScatter gscat, VecScatter lscat,DM subdm,void *ctx)
101: {
102: TS ts = (TS)ctx;
104: Vec Ydot,Ydot_s;
107: TSGLLEGetVecs(ts,dm,NULL,&Ydot);
108: TSGLLEGetVecs(ts,subdm,NULL,&Ydot_s);
110: VecScatterBegin(gscat,Ydot,Ydot_s,INSERT_VALUES,SCATTER_FORWARD);
111: VecScatterEnd(gscat,Ydot,Ydot_s,INSERT_VALUES,SCATTER_FORWARD);
113: TSGLLERestoreVecs(ts,dm,NULL,&Ydot);
114: TSGLLERestoreVecs(ts,subdm,NULL,&Ydot_s);
115: return(0);
116: }
118: static PetscErrorCode TSGLLESchemeCreate(PetscInt p,PetscInt q,PetscInt r,PetscInt s,const PetscScalar *c,
119: const PetscScalar *a,const PetscScalar *b,const PetscScalar *u,const PetscScalar *v,TSGLLEScheme *inscheme)
120: {
121: TSGLLEScheme scheme;
122: PetscInt j;
126: if (p < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Scheme order must be positive");
127: if (r < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"At least one item must be carried between steps");
128: if (s < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"At least one stage is required");
130: *inscheme = 0;
131: PetscNew(&scheme);
132: scheme->p = p;
133: scheme->q = q;
134: scheme->r = r;
135: scheme->s = s;
137: PetscMalloc5(s,&scheme->c,s*s,&scheme->a,r*s,&scheme->b,r*s,&scheme->u,r*r,&scheme->v);
138: PetscMemcpy(scheme->c,c,s*sizeof(PetscScalar));
139: for (j=0; j<s*s; j++) scheme->a[j] = (PetscAbsScalar(a[j]) < 1e-12) ? 0 : a[j];
140: for (j=0; j<r*s; j++) scheme->b[j] = (PetscAbsScalar(b[j]) < 1e-12) ? 0 : b[j];
141: for (j=0; j<s*r; j++) scheme->u[j] = (PetscAbsScalar(u[j]) < 1e-12) ? 0 : u[j];
142: for (j=0; j<r*r; j++) scheme->v[j] = (PetscAbsScalar(v[j]) < 1e-12) ? 0 : v[j];
144: PetscMalloc6(r,&scheme->alpha,r,&scheme->beta,r,&scheme->gamma,3*s,&scheme->phi,3*r,&scheme->psi,r,&scheme->stage_error);
145: {
146: PetscInt i,j,k,ss=s+2;
147: PetscBLASInt m,n,one=1,*ipiv,lwork=4*((s+3)*3+3),info,ldb;
148: PetscReal rcond,*sing,*workreal;
149: PetscScalar *ImV,*H,*bmat,*workscalar,*c=scheme->c,*a=scheme->a,*b=scheme->b,*u=scheme->u,*v=scheme->v;
150: #if !defined(PETSC_MISSING_LAPACK_GELSS)
151: PetscBLASInt rank;
152: #endif
153: PetscMalloc7(PetscSqr(r),&ImV,3*s,&H,3*ss,&bmat,lwork,&workscalar,5*(3+r),&workreal,r+s,&sing,r+s,&ipiv);
155: /* column-major input */
156: for (i=0; i<r-1; i++) {
157: for (j=0; j<r-1; j++) ImV[i+j*r] = 1.0*(i==j) - v[(i+1)*r+j+1];
158: }
159: /* Build right hand side for alpha (tp - glm.B(2:end,:)*(glm.c.^(p)./factorial(p))) */
160: for (i=1; i<r; i++) {
161: scheme->alpha[i] = 1./Factorial(p+1-i);
162: for (j=0; j<s; j++) scheme->alpha[i] -= b[i*s+j]*CPowF(c[j],p);
163: }
164: PetscBLASIntCast(r-1,&m);
165: PetscBLASIntCast(r,&n);
166: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&m,&one,ImV,&n,ipiv,scheme->alpha+1,&n,&info));
167: if (info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GESV");
168: if (info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
170: /* Build right hand side for beta (tp1 - glm.B(2:end,:)*(glm.c.^(p+1)./factorial(p+1)) - e.alpha) */
171: for (i=1; i<r; i++) {
172: scheme->beta[i] = 1./Factorial(p+2-i) - scheme->alpha[i];
173: for (j=0; j<s; j++) scheme->beta[i] -= b[i*s+j]*CPowF(c[j],p+1);
174: }
175: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("No transpose",&m,&one,ImV,&n,ipiv,scheme->beta+1,&n,&info));
176: if (info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GETRS");
177: if (info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Should not happen");
179: /* Build stage_error vector
180: xi = glm.c.^(p+1)/factorial(p+1) - glm.A*glm.c.^p/factorial(p) + glm.U(:,2:end)*e.alpha;
181: */
182: for (i=0; i<s; i++) {
183: scheme->stage_error[i] = CPowF(c[i],p+1);
184: for (j=0; j<s; j++) scheme->stage_error[i] -= a[i*s+j]*CPowF(c[j],p);
185: for (j=1; j<r; j++) scheme->stage_error[i] += u[i*r+j]*scheme->alpha[j];
186: }
188: /* alpha[0] (epsilon in B,J,W 2007)
189: epsilon = 1/factorial(p+1) - B(1,:)*c.^p/factorial(p) + V(1,2:end)*e.alpha;
190: */
191: scheme->alpha[0] = 1./Factorial(p+1);
192: for (j=0; j<s; j++) scheme->alpha[0] -= b[0*s+j]*CPowF(c[j],p);
193: for (j=1; j<r; j++) scheme->alpha[0] += v[0*r+j]*scheme->alpha[j];
195: /* right hand side for gamma (glm.B(2:end,:)*e.xi - e.epsilon*eye(s-1,1)) */
196: for (i=1; i<r; i++) {
197: scheme->gamma[i] = (i==1 ? -1. : 0)*scheme->alpha[0];
198: for (j=0; j<s; j++) scheme->gamma[i] += b[i*s+j]*scheme->stage_error[j];
199: }
200: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("No transpose",&m,&one,ImV,&n,ipiv,scheme->gamma+1,&n,&info));
201: if (info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GETRS");
202: if (info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Should not happen");
204: /* beta[0] (rho in B,J,W 2007)
205: e.rho = 1/factorial(p+2) - glm.B(1,:)*glm.c.^(p+1)/factorial(p+1) ...
206: + glm.V(1,2:end)*e.beta;% - e.epsilon;
207: % Note: The paper (B,J,W 2007) includes the last term in their definition
208: * */
209: scheme->beta[0] = 1./Factorial(p+2);
210: for (j=0; j<s; j++) scheme->beta[0] -= b[0*s+j]*CPowF(c[j],p+1);
211: for (j=1; j<r; j++) scheme->beta[0] += v[0*r+j]*scheme->beta[j];
213: /* gamma[0] (sigma in B,J,W 2007)
214: * e.sigma = glm.B(1,:)*e.xi + glm.V(1,2:end)*e.gamma;
215: * */
216: scheme->gamma[0] = 0.0;
217: for (j=0; j<s; j++) scheme->gamma[0] += b[0*s+j]*scheme->stage_error[j];
218: for (j=1; j<r; j++) scheme->gamma[0] += v[0*s+j]*scheme->gamma[j];
220: /* Assemble H
221: * % Determine the error estimators phi
222: H = [[cpow(glm.c,p) + C*e.alpha] [cpow(glm.c,p+1) + C*e.beta] ...
223: [e.xi - C*(e.gamma + 0*e.epsilon*eye(s-1,1))]]';
224: % Paper has formula above without the 0, but that term must be left
225: % out to satisfy the conditions they propose and to make the
226: % example schemes work
227: e.H = H;
228: e.phi = (H \ [1 0 0;1 1 0;0 0 -1])';
229: e.psi = -e.phi*C;
230: * */
231: for (j=0; j<s; j++) {
232: H[0+j*3] = CPowF(c[j],p);
233: H[1+j*3] = CPowF(c[j],p+1);
234: H[2+j*3] = scheme->stage_error[j];
235: for (k=1; k<r; k++) {
236: H[0+j*3] += CPowF(c[j],k-1)*scheme->alpha[k];
237: H[1+j*3] += CPowF(c[j],k-1)*scheme->beta[k];
238: H[2+j*3] -= CPowF(c[j],k-1)*scheme->gamma[k];
239: }
240: }
241: bmat[0+0*ss] = 1.; bmat[0+1*ss] = 0.; bmat[0+2*ss] = 0.;
242: bmat[1+0*ss] = 1.; bmat[1+1*ss] = 1.; bmat[1+2*ss] = 0.;
243: bmat[2+0*ss] = 0.; bmat[2+1*ss] = 0.; bmat[2+2*ss] = -1.;
244: m = 3;
245: PetscBLASIntCast(s,&n);
246: PetscBLASIntCast(ss,&ldb);
247: rcond = 1e-12;
248: #if defined(PETSC_MISSING_LAPACK_GELSS)
249: /* ESSL does not have this routine */
250: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GELSS - Lapack routine is unavailable\nNot able to run GL time stepping.");
251: #else
252: #if defined(PETSC_USE_COMPLEX)
253: /* ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO) */
254: PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&m,&n,&m,H,&m,bmat,&ldb,sing,&rcond,&rank,workscalar,&lwork,workreal,&info));
255: #else
256: /* DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO) */
257: PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&m,&n,&m,H,&m,bmat,&ldb,sing,&rcond,&rank,workscalar,&lwork,&info));
258: #endif
259: if (info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS");
260: if (info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge");
261: #endif
263: for (j=0; j<3; j++) {
264: for (k=0; k<s; k++) scheme->phi[k+j*s] = bmat[k+j*ss];
265: }
267: /* the other part of the error estimator, psi in B,J,W 2007 */
268: scheme->psi[0*r+0] = 0.;
269: scheme->psi[1*r+0] = 0.;
270: scheme->psi[2*r+0] = 0.;
271: for (j=1; j<r; j++) {
272: scheme->psi[0*r+j] = 0.;
273: scheme->psi[1*r+j] = 0.;
274: scheme->psi[2*r+j] = 0.;
275: for (k=0; k<s; k++) {
276: scheme->psi[0*r+j] -= CPowF(c[k],j-1)*scheme->phi[0*s+k];
277: scheme->psi[1*r+j] -= CPowF(c[k],j-1)*scheme->phi[1*s+k];
278: scheme->psi[2*r+j] -= CPowF(c[k],j-1)*scheme->phi[2*s+k];
279: }
280: }
281: PetscFree7(ImV,H,bmat,workscalar,workreal,sing,ipiv);
282: }
283: /* Check which properties are satisfied */
284: scheme->stiffly_accurate = PETSC_TRUE;
285: if (scheme->c[s-1] != 1.) scheme->stiffly_accurate = PETSC_FALSE;
286: for (j=0; j<s; j++) if (a[(s-1)*s+j] != b[j]) scheme->stiffly_accurate = PETSC_FALSE;
287: for (j=0; j<r; j++) if (u[(s-1)*r+j] != v[j]) scheme->stiffly_accurate = PETSC_FALSE;
288: scheme->fsal = scheme->stiffly_accurate; /* FSAL is stronger */
289: for (j=0; j<s-1; j++) if (r>1 && b[1*s+j] != 0.) scheme->fsal = PETSC_FALSE;
290: if (b[1*s+r-1] != 1.) scheme->fsal = PETSC_FALSE;
291: for (j=0; j<r; j++) if (r>1 && v[1*r+j] != 0.) scheme->fsal = PETSC_FALSE;
293: *inscheme = scheme;
294: return(0);
295: }
297: static PetscErrorCode TSGLLESchemeDestroy(TSGLLEScheme sc)
298: {
302: PetscFree5(sc->c,sc->a,sc->b,sc->u,sc->v);
303: PetscFree6(sc->alpha,sc->beta,sc->gamma,sc->phi,sc->psi,sc->stage_error);
304: PetscFree(sc);
305: return(0);
306: }
308: static PetscErrorCode TSGLLEDestroy_Default(TS_GLLE *gl)
309: {
311: PetscInt i;
314: for (i=0; i<gl->nschemes; i++) {
315: if (gl->schemes[i]) {TSGLLESchemeDestroy(gl->schemes[i]);}
316: }
317: PetscFree(gl->schemes);
318: gl->nschemes = 0;
319: PetscMemzero(gl->type_name,sizeof(gl->type_name));
320: return(0);
321: }
323: static PetscErrorCode TSGLLEViewTable_Private(PetscViewer viewer,PetscInt m,PetscInt n,const PetscScalar a[],const char name[])
324: {
326: PetscBool iascii;
327: PetscInt i,j;
330: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
331: if (iascii) {
332: PetscViewerASCIIPrintf(viewer,"%30s = [",name);
333: for (i=0; i<m; i++) {
334: if (i) {PetscViewerASCIIPrintf(viewer,"%30s [","");}
335: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
336: for (j=0; j<n; j++) {
337: PetscViewerASCIIPrintf(viewer," %12.8g",PetscRealPart(a[i*n+j]));
338: }
339: PetscViewerASCIIPrintf(viewer,"]\n");
340: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
341: }
342: }
343: return(0);
344: }
347: static PetscErrorCode TSGLLESchemeView(TSGLLEScheme sc,PetscBool view_details,PetscViewer viewer)
348: {
350: PetscBool iascii;
353: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
354: if (iascii) {
355: PetscViewerASCIIPrintf(viewer,"GL scheme p,q,r,s = %d,%d,%d,%d\n",sc->p,sc->q,sc->r,sc->s);
356: PetscViewerASCIIPushTab(viewer);
357: PetscViewerASCIIPrintf(viewer,"Stiffly accurate: %s, FSAL: %s\n",sc->stiffly_accurate ? "yes" : "no",sc->fsal ? "yes" : "no");
358: PetscViewerASCIIPrintf(viewer,"Leading error constants: %10.3e %10.3e %10.3e\n",
359: PetscRealPart(sc->alpha[0]),PetscRealPart(sc->beta[0]),PetscRealPart(sc->gamma[0]));
360: TSGLLEViewTable_Private(viewer,1,sc->s,sc->c,"Abscissas c");
361: if (view_details) {
362: TSGLLEViewTable_Private(viewer,sc->s,sc->s,sc->a,"A");
363: TSGLLEViewTable_Private(viewer,sc->r,sc->s,sc->b,"B");
364: TSGLLEViewTable_Private(viewer,sc->s,sc->r,sc->u,"U");
365: TSGLLEViewTable_Private(viewer,sc->r,sc->r,sc->v,"V");
367: TSGLLEViewTable_Private(viewer,3,sc->s,sc->phi,"Error estimate phi");
368: TSGLLEViewTable_Private(viewer,3,sc->r,sc->psi,"Error estimate psi");
369: TSGLLEViewTable_Private(viewer,1,sc->r,sc->alpha,"Modify alpha");
370: TSGLLEViewTable_Private(viewer,1,sc->r,sc->beta,"Modify beta");
371: TSGLLEViewTable_Private(viewer,1,sc->r,sc->gamma,"Modify gamma");
372: TSGLLEViewTable_Private(viewer,1,sc->s,sc->stage_error,"Stage error xi");
373: }
374: PetscViewerASCIIPopTab(viewer);
375: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
376: return(0);
377: }
379: static PetscErrorCode TSGLLEEstimateHigherMoments_Default(TSGLLEScheme sc,PetscReal h,Vec Ydot[],Vec Xold[],Vec hm[])
380: {
382: PetscInt i;
385: if (sc->r > 64 || sc->s > 64) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Ridiculous number of stages or items passed between stages");
386: /* build error vectors*/
387: for (i=0; i<3; i++) {
388: PetscScalar phih[64];
389: PetscInt j;
390: for (j=0; j<sc->s; j++) phih[j] = sc->phi[i*sc->s+j]*h;
391: VecZeroEntries(hm[i]);
392: VecMAXPY(hm[i],sc->s,phih,Ydot);
393: VecMAXPY(hm[i],sc->r,&sc->psi[i*sc->r],Xold);
394: }
395: return(0);
396: }
398: static PetscErrorCode TSGLLECompleteStep_Rescale(TSGLLEScheme sc,PetscReal h,TSGLLEScheme next_sc,PetscReal next_h,Vec Ydot[],Vec Xold[],Vec X[])
399: {
401: PetscScalar brow[32],vrow[32];
402: PetscInt i,j,r,s;
405: /* Build the new solution from (X,Ydot) */
406: r = sc->r;
407: s = sc->s;
408: for (i=0; i<r; i++) {
409: VecZeroEntries(X[i]);
410: for (j=0; j<s; j++) brow[j] = h*sc->b[i*s+j];
411: VecMAXPY(X[i],s,brow,Ydot);
412: for (j=0; j<r; j++) vrow[j] = sc->v[i*r+j];
413: VecMAXPY(X[i],r,vrow,Xold);
414: }
415: return(0);
416: }
418: static PetscErrorCode TSGLLECompleteStep_RescaleAndModify(TSGLLEScheme sc,PetscReal h,TSGLLEScheme next_sc,PetscReal next_h,Vec Ydot[],Vec Xold[],Vec X[])
419: {
421: PetscScalar brow[32],vrow[32];
422: PetscReal ratio;
423: PetscInt i,j,p,r,s;
426: /* Build the new solution from (X,Ydot) */
427: p = sc->p;
428: r = sc->r;
429: s = sc->s;
430: ratio = next_h/h;
431: for (i=0; i<r; i++) {
432: VecZeroEntries(X[i]);
433: for (j=0; j<s; j++) {
434: brow[j] = h*(PetscPowRealInt(ratio,i)*sc->b[i*s+j]
435: + (PetscPowRealInt(ratio,i) - PetscPowRealInt(ratio,p+1))*(+ sc->alpha[i]*sc->phi[0*s+j])
436: + (PetscPowRealInt(ratio,i) - PetscPowRealInt(ratio,p+2))*(+ sc->beta [i]*sc->phi[1*s+j]
437: + sc->gamma[i]*sc->phi[2*s+j]));
438: }
439: VecMAXPY(X[i],s,brow,Ydot);
440: for (j=0; j<r; j++) {
441: vrow[j] = (PetscPowRealInt(ratio,i)*sc->v[i*r+j]
442: + (PetscPowRealInt(ratio,i) - PetscPowRealInt(ratio,p+1))*(+ sc->alpha[i]*sc->psi[0*r+j])
443: + (PetscPowRealInt(ratio,i) - PetscPowRealInt(ratio,p+2))*(+ sc->beta [i]*sc->psi[1*r+j]
444: + sc->gamma[i]*sc->psi[2*r+j]));
445: }
446: VecMAXPY(X[i],r,vrow,Xold);
447: }
448: if (r < next_sc->r) {
449: if (r+1 != next_sc->r) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cannot accommodate jump in r greater than 1");
450: VecZeroEntries(X[r]);
451: for (j=0; j<s; j++) brow[j] = h*PetscPowRealInt(ratio,p+1)*sc->phi[0*s+j];
452: VecMAXPY(X[r],s,brow,Ydot);
453: for (j=0; j<r; j++) vrow[j] = PetscPowRealInt(ratio,p+1)*sc->psi[0*r+j];
454: VecMAXPY(X[r],r,vrow,Xold);
455: }
456: return(0);
457: }
459: static PetscErrorCode TSGLLECreate_IRKS(TS ts)
460: {
461: TS_GLLE *gl = (TS_GLLE*)ts->data;
465: gl->Destroy = TSGLLEDestroy_Default;
466: gl->EstimateHigherMoments = TSGLLEEstimateHigherMoments_Default;
467: gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify;
468: PetscMalloc1(10,&gl->schemes);
469: gl->nschemes = 0;
471: {
472: /* p=1,q=1, r=s=2, A- and L-stable with error estimates of order 2 and 3
473: * Listed in Butcher & Podhaisky 2006. On error estimation in general linear methods for stiff ODE.
474: * irks(0.3,0,[.3,1],[1],1)
475: * Note: can be made to have classical order (not stage order) 2 by replacing 0.3 with 1-sqrt(1/2)
476: * but doing so would sacrifice the error estimator.
477: */
478: const PetscScalar c[2] = {3./10., 1.};
479: const PetscScalar a[2][2] = {{3./10., 0}, {7./10., 3./10.}};
480: const PetscScalar b[2][2] = {{7./10., 3./10.}, {0,1}};
481: const PetscScalar u[2][2] = {{1,0},{1,0}};
482: const PetscScalar v[2][2] = {{1,0},{0,0}};
483: TSGLLESchemeCreate(1,1,2,2,c,*a,*b,*u,*v,&gl->schemes[gl->nschemes++]);
484: }
486: {
487: /* p=q=2, r=s=3: irks(4/9,0,[1:3]/3,[0.33852],1) */
488: /* http://www.math.auckland.ac.nz/~hpod/atlas/i2a.html */
489: const PetscScalar c[3] = {1./3., 2./3., 1}
490: ,a[3][3] = {{4./9. ,0 , 0},
491: {1.03750643704090e+00 , 4./9., 0},
492: {7.67024779410304e-01 , -3.81140216918943e-01, 4./9.}}
493: ,b[3][3] = {{0.767024779410304, -0.381140216918943, 4./9.},
494: {0.000000000000000, 0.000000000000000, 1.000000000000000},
495: {-2.075048385225385, 0.621728385225383, 1.277197204924873}}
496: ,u[3][3] = {{1.0000000000000000, -0.1111111111111109, -0.0925925925925922},
497: {1.0000000000000000, -0.8152842148186744, -0.4199095530877056},
498: {1.0000000000000000, 0.1696709930641948, 0.0539741070314165}}
499: ,v[3][3] = {{1.0000000000000000, 0.1696709930641948, 0.0539741070314165},
500: {0.000000000000000, 0.000000000000000, 0.000000000000000},
501: {0.000000000000000, 0.176122795075129, 0.000000000000000}};
502: TSGLLESchemeCreate(2,2,3,3,c,*a,*b,*u,*v,&gl->schemes[gl->nschemes++]);
503: }
504: {
505: /* p=q=3, r=s=4: irks(9/40,0,[1:4]/4,[0.3312 1.0050],[0.49541 1;1 0]) */
506: const PetscScalar c[4] = {0.25,0.5,0.75,1.0}
507: ,a[4][4] = {{9./40. , 0, 0, 0},
508: {2.11286958887701e-01 , 9./40. , 0, 0},
509: {9.46338294287584e-01 , -3.42942861246094e-01, 9./40. , 0},
510: {0.521490453970721 , -0.662474225622980, 0.490476425459734, 9./40. }}
511: ,b[4][4] = {{0.521490453970721 , -0.662474225622980, 0.490476425459734, 9./40. },
512: {0.000000000000000 , 0.000000000000000, 0.000000000000000, 1.000000000000000},
513: {-0.084677029310348 , 1.390757514776085, -1.568157386206001, 2.023192696767826},
514: {0.465383797936408 , 1.478273530625148, -1.930836081010182, 1.644872111193354}}
515: ,u[4][4] = {{1.00000000000000000 , 0.02500000000001035, -0.02499999999999053, -0.00442708333332865},
516: {1.00000000000000000 , 0.06371304111232945, -0.04032173972189845, -0.01389438413189452},
517: {1.00000000000000000 , -0.07839543304147778, 0.04738685705116663, 0.02032603595928376},
518: {1.00000000000000000 , 0.42550734619251651, 0.10800718022400080, -0.01726712647760034}}
519: ,v[4][4] = {{1.00000000000000000 , 0.42550734619251651, 0.10800718022400080, -0.01726712647760034},
520: {0.000000000000000 , 0.000000000000000, 0.000000000000000, 0.000000000000000},
521: {0.000000000000000 , -1.761115796027561, -0.521284157173780, 0.258249384305463},
522: {0.000000000000000 , -1.657693358744728, -1.052227765232394, 0.521284157173780}};
523: TSGLLESchemeCreate(3,3,4,4,c,*a,*b,*u,*v,&gl->schemes[gl->nschemes++]);
524: }
525: {
526: /* p=q=4, r=s=5:
527: irks(3/11,0,[1:5]/5, [0.1715 -0.1238 0.6617],...
528: [ -0.0812 0.4079 1.0000
529: 1.0000 0 0
530: 0.8270 1.0000 0])
531: */
532: const PetscScalar c[5] = {0.2,0.4,0.6,0.8,1.0}
533: ,a[5][5] = {{2.72727272727352e-01 , 0.00000000000000e+00, 0.00000000000000e+00 , 0.00000000000000e+00 , 0.00000000000000e+00},
534: {-1.03980153733431e-01, 2.72727272727405e-01, 0.00000000000000e+00, 0.00000000000000e+00 , 0.00000000000000e+00},
535: {-1.58615400341492e+00, 7.44168951881122e-01, 2.72727272727309e-01, 0.00000000000000e+00 , 0.00000000000000e+00},
536: {-8.73658042865628e-01, 5.37884671894595e-01, -1.63298538799523e-01, 2.72727272726996e-01 , 0.00000000000000e+00},
537: {2.95489397443992e-01 , -1.18481693910097e+00 , -6.68029812659953e-01 , 1.00716687860943e+00 , 2.72727272727288e-01}}
538: ,b[5][5] = {{2.95489397443992e-01 , -1.18481693910097e+00 , -6.68029812659953e-01 , 1.00716687860943e+00 , 2.72727272727288e-01},
539: {0.00000000000000e+00 , 1.11022302462516e-16 , -2.22044604925031e-16 , 0.00000000000000e+00 , 1.00000000000000e+00},
540: {-4.05882503986005e+00, -4.00924006567769e+00, -1.38930610972481e+00, 4.45223930308488e+00 , 6.32331093108427e-01},
541: {8.35690179937017e+00 , -2.26640927349732e+00 , 6.86647884973826e+00 , -5.22595158025740e+00 , 4.50893068837431e+00},
542: {1.27656267027479e+01 , 2.80882153840821e+00 , 8.91173096522890e+00 , -1.07936444078906e+01 , 4.82534148988854e+00}}
543: ,u[5][5] = {{1.00000000000000e+00 , -7.27272727273551e-02 , -3.45454545454419e-02 , -4.12121212119565e-03 ,-2.96969696964014e-04},
544: {1.00000000000000e+00 , 2.31252881006154e-01 , -8.29487834416481e-03 , -9.07191207681020e-03 ,-1.70378403743473e-03},
545: {1.00000000000000e+00 , 1.16925777880663e+00 , 3.59268562942635e-02 , -4.09013451730615e-02 ,-1.02411119670164e-02},
546: {1.00000000000000e+00 , 1.02634463704356e+00 , 1.59375044913405e-01 , 1.89673015035370e-03 ,-4.89987231897569e-03},
547: {1.00000000000000e+00 , 1.27746320298021e+00 , 2.37186008132728e-01 , -8.28694373940065e-02 ,-5.34396510196430e-02}}
548: ,v[5][5] = {{1.00000000000000e+00 , 1.27746320298021e+00 , 2.37186008132728e-01 , -8.28694373940065e-02 ,-5.34396510196430e-02},
549: {0.00000000000000e+00 , -1.77635683940025e-15 , -1.99840144432528e-15 , -9.99200722162641e-16 ,-3.33066907387547e-16},
550: {0.00000000000000e+00 , 4.37280081906924e+00 , 5.49221645016377e-02 , -8.88913177394943e-02 , 1.12879077989154e-01},
551: {0.00000000000000e+00 , -1.22399504837280e+01 , -5.21287338448645e+00 , -8.03952325565291e-01 , 4.60298678047147e-01},
552: {0.00000000000000e+00 , -1.85178762883829e+01 , -5.21411849862624e+00 , -1.04283436528809e+00 , 7.49030161063651e-01}};
553: TSGLLESchemeCreate(4,4,5,5,c,*a,*b,*u,*v,&gl->schemes[gl->nschemes++]);
554: }
555: {
556: /* p=q=5, r=s=6;
557: irks(1/3,0,[1:6]/6,...
558: [-0.0489 0.4228 -0.8814 0.9021],...
559: [-0.3474 -0.6617 0.6294 0.2129
560: 0.0044 -0.4256 -0.1427 -0.8936
561: -0.8267 0.4821 0.1371 -0.2557
562: -0.4426 -0.3855 -0.7514 0.3014])
563: */
564: const PetscScalar c[6] = {1./6, 2./6, 3./6, 4./6, 5./6, 1.}
565: ,a[6][6] = {{ 3.33333333333940e-01, 0 , 0 , 0 , 0 , 0 },
566: { -8.64423857333350e-02, 3.33333333332888e-01, 0 , 0 , 0 , 0 },
567: { -2.16850174258252e+00, -2.23619072028839e+00, 3.33333333335204e-01, 0 , 0 , 0 },
568: { -4.73160970138997e+00, -3.89265344629268e+00, -2.76318716520933e-01, 3.33333333335759e-01, 0 , 0 },
569: { -6.75187540297338e+00, -7.90756533769377e+00, 7.90245051802259e-01, -4.48352364517632e-01, 3.33333333328483e-01, 0 },
570: { -4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00, -2.23969848002481e+00, 6.62807710124007e-01, 3.33333333335440e-01}}
571: ,b[6][6] = {{ -4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00, -2.23969848002481e+00, 6.62807710124007e-01, 3.33333333335440e-01},
572: { -8.88178419700125e-16, 4.44089209850063e-16, -1.54737334057131e-15, -8.88178419700125e-16, 0.00000000000000e+00, 1.00000000000001e+00},
573: { -2.87780425770651e+01, -1.13520448264971e+01, 2.62002318943161e+01, 2.56943874812797e+01, -3.06702268304488e+01, 6.68067773510103e+00},
574: { 5.47971245256474e+01, 6.80366875868284e+01, -6.50952588861999e+01, -8.28643975339097e+01, 8.17416943896414e+01, -1.17819043489036e+01},
575: { -2.33332114788869e+02, 6.12942539462634e+01, -4.91850135865944e+01, 1.82716844135480e+02, -1.29788173979395e+02, 3.09968095651099e+01},
576: { -1.72049132343751e+02, 8.60194713593999e+00, 7.98154219170200e-01, 1.50371386053218e+02, -1.18515423962066e+02, 2.50898277784663e+01}}
577: ,u[6][6] = {{ 1.00000000000000e+00, -1.66666666666870e-01, -4.16666666664335e-02, -3.85802469124815e-03, -2.25051440302250e-04, -9.64506172339142e-06},
578: { 1.00000000000000e+00, 8.64423857327162e-02, -4.11484912671353e-02, -1.11450903217645e-02, -1.47651050487126e-03, -1.34395070766826e-04},
579: { 1.00000000000000e+00, 4.57135912953434e+00, 1.06514719719137e+00, 1.33517564218007e-01, 1.11365952968659e-02, 6.12382756769504e-04},
580: { 1.00000000000000e+00, 9.23391519753404e+00, 2.22431212392095e+00, 2.91823807741891e-01, 2.52058456411084e-02, 1.22800542949647e-03},
581: { 1.00000000000000e+00, 1.48175480533865e+01, 3.73439117461835e+00, 5.14648336541804e-01, 4.76430038853402e-02, 2.56798515502156e-03},
582: { 1.00000000000000e+00, 1.50512347758335e+01, 4.10099701165164e+00, 5.66039141003603e-01, 3.91213893800891e-02, -2.99136269067853e-03}}
583: ,v[6][6] = {{ 1.00000000000000e+00, 1.50512347758335e+01, 4.10099701165164e+00, 5.66039141003603e-01, 3.91213893800891e-02, -2.99136269067853e-03},
584: { 0.00000000000000e+00, -4.88498130835069e-15, -6.43929354282591e-15, -3.55271367880050e-15, -1.22124532708767e-15, -3.12250225675825e-16},
585: { 0.00000000000000e+00, 1.22250171233141e+01, -1.77150760606169e+00, 3.54516769879390e-01, 6.22298845883398e-01, 2.31647447450276e-01},
586: { 0.00000000000000e+00, -4.48339457331040e+01, -3.57363126641880e-01, 5.18750173123425e-01, 6.55727990241799e-02, 1.63175368287079e-01},
587: { 0.00000000000000e+00, 1.37297394708005e+02, -1.60145272991317e+00, -5.05319555199441e+00, 1.55328940390990e-01, 9.16629423682464e-01},
588: { 0.00000000000000e+00, 1.05703241119022e+02, -1.16610260983038e+00, -2.99767252773859e+00, -1.13472315553890e-01, 1.09742849254729e+00}};
589: TSGLLESchemeCreate(5,5,6,6,c,*a,*b,*u,*v,&gl->schemes[gl->nschemes++]);
590: }
591: return(0);
592: }
594: /*@C
595: TSGLLESetType - sets the class of general linear method to use for time-stepping
597: Collective on TS
599: Input Parameters:
600: + ts - the TS context
601: - type - a method
603: Options Database Key:
604: . -ts_gl_type <type> - sets the method, use -help for a list of available method (e.g. irks)
606: Notes:
607: See "petsc/include/petscts.h" for available methods (for instance)
608: . TSGLLE_IRKS - Diagonally implicit methods with inherent Runge-Kutta stability (for stiff problems)
610: Normally, it is best to use the TSSetFromOptions() command and
611: then set the TSGLLE type from the options database rather than by using
612: this routine. Using the options database provides the user with
613: maximum flexibility in evaluating the many different solvers.
614: The TSGLLESetType() routine is provided for those situations where it
615: is necessary to set the timestepping solver independently of the
616: command line or options database. This might be the case, for example,
617: when the choice of solver changes during the execution of the
618: program, and the user's application is taking responsibility for
619: choosing the appropriate method.
621: Level: intermediate
623: .keywords: TS, TSGLLE, set, type
624: @*/
625: PetscErrorCode TSGLLESetType(TS ts,TSGLLEType type)
626: {
632: PetscTryMethod(ts,"TSGLLESetType_C",(TS,TSGLLEType),(ts,type));
633: return(0);
634: }
636: /*@C
637: TSGLLESetAcceptType - sets the acceptance test
639: Time integrators that need to control error must have the option to reject a time step based on local error
640: estimates. This function allows different schemes to be set.
642: Logically Collective on TS
644: Input Parameters:
645: + ts - the TS context
646: - type - the type
648: Options Database Key:
649: . -ts_gl_accept_type <type> - sets the method used to determine whether to accept or reject a step
651: Level: intermediate
653: .seealso: TS, TSGLLE, TSGLLEAcceptRegister(), TSGLLEAdapt, set type
654: @*/
655: PetscErrorCode TSGLLESetAcceptType(TS ts,TSGLLEAcceptType type)
656: {
662: PetscTryMethod(ts,"TSGLLESetAcceptType_C",(TS,TSGLLEAcceptType),(ts,type));
663: return(0);
664: }
666: /*@C
667: TSGLLEGetAdapt - gets the TSGLLEAdapt object from the TS
669: Not Collective
671: Input Parameter:
672: . ts - the TS context
674: Output Parameter:
675: . adapt - the TSGLLEAdapt context
677: Notes:
678: This allows the user set options on the TSGLLEAdapt object. Usually it is better to do this using the options
679: database, so this function is rarely needed.
681: Level: advanced
683: .seealso: TSGLLEAdapt, TSGLLEAdaptRegister()
684: @*/
685: PetscErrorCode TSGLLEGetAdapt(TS ts,TSGLLEAdapt *adapt)
686: {
692: PetscUseMethod(ts,"TSGLLEGetAdapt_C",(TS,TSGLLEAdapt*),(ts,adapt));
693: return(0);
694: }
696: static PetscErrorCode TSGLLEAccept_Always(TS ts,PetscReal tleft,PetscReal h,const PetscReal enorms[],PetscBool *accept)
697: {
699: *accept = PETSC_TRUE;
700: return(0);
701: }
703: static PetscErrorCode TSGLLEUpdateWRMS(TS ts)
704: {
705: TS_GLLE *gl = (TS_GLLE*)ts->data;
707: PetscScalar *x,*w;
708: PetscInt n,i;
711: VecGetArray(gl->X[0],&x);
712: VecGetArray(gl->W,&w);
713: VecGetLocalSize(gl->W,&n);
714: for (i=0; i<n; i++) w[i] = 1./(gl->wrms_atol + gl->wrms_rtol*PetscAbsScalar(x[i]));
715: VecRestoreArray(gl->X[0],&x);
716: VecRestoreArray(gl->W,&w);
717: return(0);
718: }
720: static PetscErrorCode TSGLLEVecNormWRMS(TS ts,Vec X,PetscReal *nrm)
721: {
722: TS_GLLE *gl = (TS_GLLE*)ts->data;
724: PetscScalar *x,*w;
725: PetscReal sum = 0.0,gsum;
726: PetscInt n,N,i;
729: VecGetArray(X,&x);
730: VecGetArray(gl->W,&w);
731: VecGetLocalSize(gl->W,&n);
732: for (i=0; i<n; i++) sum += PetscAbsScalar(PetscSqr(x[i]*w[i]));
733: VecRestoreArray(X,&x);
734: VecRestoreArray(gl->W,&w);
735: MPIU_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));
736: VecGetSize(gl->W,&N);
737: *nrm = PetscSqrtReal(gsum/(1.*N));
738: return(0);
739: }
741: static PetscErrorCode TSGLLESetType_GLLE(TS ts,TSGLLEType type)
742: {
743: PetscErrorCode ierr,(*r)(TS);
744: PetscBool same;
745: TS_GLLE *gl = (TS_GLLE*)ts->data;
748: if (gl->type_name[0]) {
749: PetscStrcmp(gl->type_name,type,&same);
750: if (same) return(0);
751: (*gl->Destroy)(gl);
752: }
754: PetscFunctionListFind(TSGLLEList,type,&r);
755: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLLE type \"%s\" given",type);
756: (*r)(ts);
757: PetscStrcpy(gl->type_name,type);
758: return(0);
759: }
761: static PetscErrorCode TSGLLESetAcceptType_GLLE(TS ts,TSGLLEAcceptType type)
762: {
763: PetscErrorCode ierr;
764: TSGLLEAcceptFunction r;
765: TS_GLLE *gl = (TS_GLLE*)ts->data;
768: PetscFunctionListFind(TSGLLEAcceptList,type,&r);
769: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLLEAccept type \"%s\" given",type);
770: gl->Accept = r;
771: PetscStrncpy(gl->accept_name,type,sizeof(gl->accept_name));
772: return(0);
773: }
775: static PetscErrorCode TSGLLEGetAdapt_GLLE(TS ts,TSGLLEAdapt *adapt)
776: {
778: TS_GLLE *gl = (TS_GLLE*)ts->data;
781: if (!gl->adapt) {
782: TSGLLEAdaptCreate(PetscObjectComm((PetscObject)ts),&gl->adapt);
783: PetscObjectIncrementTabLevel((PetscObject)gl->adapt,(PetscObject)ts,1);
784: PetscLogObjectParent((PetscObject)ts,(PetscObject)gl->adapt);
785: }
786: *adapt = gl->adapt;
787: return(0);
788: }
790: static PetscErrorCode TSGLLEChooseNextScheme(TS ts,PetscReal h,const PetscReal hmnorm[],PetscInt *next_scheme,PetscReal *next_h,PetscBool *finish)
791: {
793: TS_GLLE *gl = (TS_GLLE*)ts->data;
794: PetscInt i,n,cur_p,cur,next_sc,candidates[64],orders[64];
795: PetscReal errors[64],costs[64],tleft;
798: cur = -1;
799: cur_p = gl->schemes[gl->current_scheme]->p;
800: tleft = ts->max_time - (ts->ptime + ts->time_step);
801: for (i=0,n=0; i<gl->nschemes; i++) {
802: TSGLLEScheme sc = gl->schemes[i];
803: if (sc->p < gl->min_order || gl->max_order < sc->p) continue;
804: if (sc->p == cur_p - 1) errors[n] = PetscAbsScalar(sc->alpha[0])*hmnorm[0];
805: else if (sc->p == cur_p) errors[n] = PetscAbsScalar(sc->alpha[0])*hmnorm[1];
806: else if (sc->p == cur_p+1) errors[n] = PetscAbsScalar(sc->alpha[0])*(hmnorm[2]+hmnorm[3]);
807: else continue;
808: candidates[n] = i;
809: orders[n] = PetscMin(sc->p,sc->q); /* order of global truncation error */
810: costs[n] = sc->s; /* estimate the cost as the number of stages */
811: if (i == gl->current_scheme) cur = n;
812: n++;
813: }
814: if (cur < 0 || gl->nschemes <= cur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Current scheme not found in scheme list");
815: TSGLLEAdaptChoose(gl->adapt,n,orders,errors,costs,cur,h,tleft,&next_sc,next_h,finish);
816: *next_scheme = candidates[next_sc];
817: 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);
818: return(0);
819: }
821: static PetscErrorCode TSGLLEGetMaxSizes(TS ts,PetscInt *max_r,PetscInt *max_s)
822: {
823: TS_GLLE *gl = (TS_GLLE*)ts->data;
826: *max_r = gl->schemes[gl->nschemes-1]->r;
827: *max_s = gl->schemes[gl->nschemes-1]->s;
828: return(0);
829: }
831: static PetscErrorCode TSSolve_GLLE(TS ts)
832: {
833: TS_GLLE *gl = (TS_GLLE*)ts->data;
834: PetscInt i,k,its,lits,max_r,max_s;
835: PetscBool final_step,finish;
836: SNESConvergedReason snesreason;
837: PetscErrorCode ierr;
840: TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);
842: TSGLLEGetMaxSizes(ts,&max_r,&max_s);
843: VecCopy(ts->vec_sol,gl->X[0]);
844: for (i=1; i<max_r; i++) {
845: VecZeroEntries(gl->X[i]);
846: }
847: TSGLLEUpdateWRMS(ts);
849: if (0) {
850: /* Find consistent initial data for DAE */
851: gl->stage_time = ts->ptime + ts->time_step;
852: gl->scoeff = 1.;
853: gl->stage = 0;
855: VecCopy(ts->vec_sol,gl->Z);
856: VecCopy(ts->vec_sol,gl->Y);
857: SNESSolve(ts->snes,NULL,gl->Y);
858: SNESGetIterationNumber(ts->snes,&its);
859: SNESGetLinearSolveIterations(ts->snes,&lits);
860: SNESGetConvergedReason(ts->snes,&snesreason);
862: ts->snes_its += its; ts->ksp_its += lits;
863: if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) {
864: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
865: PetscInfo2(ts,"Step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
866: return(0);
867: }
868: }
870: if (gl->current_scheme < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"A starting scheme has not been provided");
872: for (k=0,final_step=PETSC_FALSE,finish=PETSC_FALSE; k<ts->max_steps && !finish; k++) {
873: PetscInt j,r,s,next_scheme = 0,rejections;
874: PetscReal h,hmnorm[4],enorm[3],next_h;
875: PetscBool accept;
876: const PetscScalar *c,*a,*u;
877: Vec *X,*Ydot,Y;
878: TSGLLEScheme scheme = gl->schemes[gl->current_scheme];
880: r = scheme->r; s = scheme->s;
881: c = scheme->c;
882: a = scheme->a; u = scheme->u;
883: h = ts->time_step;
884: X = gl->X; Ydot = gl->Ydot; Y = gl->Y;
886: if (ts->ptime > ts->max_time) break;
888: /*
889: We only call PreStep at the start of each STEP, not each STAGE. This is because it is
890: possible to fail (have to restart a step) after multiple stages.
891: */
892: TSPreStep(ts);
894: rejections = 0;
895: while (1) {
896: for (i=0; i<s; i++) {
897: PetscScalar shift;
898: gl->scoeff = 1./PetscRealPart(a[i*s+i]);
899: shift = gl->scoeff/ts->time_step;
900: gl->stage = i;
901: gl->stage_time = ts->ptime + PetscRealPart(c[i])*h;
903: /*
904: * Stage equation: Y = h A Y' + U X
905: * We assume that A is lower-triangular so that we can solve the stages (Y,Y') sequentially
906: * Build the affine vector z_i = -[1/(h a_ii)](h sum_j a_ij y'_j + sum_j u_ij x_j)
907: * Then y'_i = z + 1/(h a_ii) y_i
908: */
909: VecZeroEntries(gl->Z);
910: for (j=0; j<r; j++) {
911: VecAXPY(gl->Z,-shift*u[i*r+j],X[j]);
912: }
913: for (j=0; j<i; j++) {
914: VecAXPY(gl->Z,-shift*h*a[i*s+j],Ydot[j]);
915: }
916: /* Note: Z is used within function evaluation, Ydot = Z + shift*Y */
918: /* Compute an estimate of Y to start Newton iteration */
919: if (gl->extrapolate) {
920: if (i==0) {
921: /* Linear extrapolation on the first stage */
922: VecWAXPY(Y,c[i]*h,X[1],X[0]);
923: } else {
924: /* Linear extrapolation from the last stage */
925: VecAXPY(Y,(c[i]-c[i-1])*h,Ydot[i-1]);
926: }
927: } else if (i==0) { /* Directly use solution from the last step, otherwise reuse the last stage (do nothing) */
928: VecCopy(X[0],Y);
929: }
931: /* Solve this stage (Ydot[i] is computed during function evaluation) */
932: SNESSolve(ts->snes,NULL,Y);
933: SNESGetIterationNumber(ts->snes,&its);
934: SNESGetLinearSolveIterations(ts->snes,&lits);
935: SNESGetConvergedReason(ts->snes,&snesreason);
936: ts->snes_its += its; ts->ksp_its += lits;
937: if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) {
938: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
939: PetscInfo2(ts,"Step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
940: return(0);
941: }
942: }
944: gl->stage_time = ts->ptime + ts->time_step;
946: (*gl->EstimateHigherMoments)(scheme,h,Ydot,gl->X,gl->himom);
947: /* hmnorm[i] = h^{p+i}x^{(p+i)} with i=0,1,2; hmnorm[3] = h^{p+2}(dx'/dx) x^{(p+1)} */
948: for (i=0; i<3; i++) {
949: TSGLLEVecNormWRMS(ts,gl->himom[i],&hmnorm[i+1]);
950: }
951: enorm[0] = PetscRealPart(scheme->alpha[0])*hmnorm[1];
952: enorm[1] = PetscRealPart(scheme->beta[0]) *hmnorm[2];
953: enorm[2] = PetscRealPart(scheme->gamma[0])*hmnorm[3];
954: (*gl->Accept)(ts,ts->max_time-gl->stage_time,h,enorm,&accept);
955: if (accept) goto accepted;
956: rejections++;
957: PetscInfo3(ts,"Step %D (t=%g) not accepted, rejections=%D\n",k,gl->stage_time,rejections);
958: if (rejections > gl->max_step_rejections) break;
959: /*
960: There are lots of reasons why a step might be rejected, including solvers not converging and other factors that
961: TSGLLEChooseNextScheme does not support. Additionally, the error estimates may be very screwed up, so I'm not
962: convinced that it's safe to just compute a new error estimate using the same interface as the current adaptor
963: (the adaptor interface probably has to change). Here we make an arbitrary and naive choice. This assumes that
964: steps were written in Nordsieck form. The "correct" method would be to re-complete the previous time step with
965: the correct "next" step size. It is unclear to me whether the present ad-hoc method of rescaling X is stable.
966: */
967: h *= 0.5;
968: for (i=1; i<scheme->r; i++) {
969: VecScale(X[i],PetscPowRealInt(0.5,i));
970: }
971: }
972: SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Time step %D (t=%g) not accepted after %D failures\n",k,gl->stage_time,rejections);
974: accepted:
975: /* This term is not error, but it *would* be the leading term for a lower order method */
976: TSGLLEVecNormWRMS(ts,gl->X[scheme->r-1],&hmnorm[0]);
977: /* Correct scaling so that these are equivalent to norms of the Nordsieck vectors */
979: 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]);
980: if (!final_step) {
981: TSGLLEChooseNextScheme(ts,h,hmnorm,&next_scheme,&next_h,&final_step);
982: } else {
983: /* Dummy values to complete the current step in a consistent manner */
984: next_scheme = gl->current_scheme;
985: next_h = h;
986: finish = PETSC_TRUE;
987: }
989: X = gl->Xold;
990: gl->Xold = gl->X;
991: gl->X = X;
992: (*gl->CompleteStep)(scheme,h,gl->schemes[next_scheme],next_h,Ydot,gl->Xold,gl->X);
994: TSGLLEUpdateWRMS(ts);
996: /* Post the solution for the user, we could avoid this copy with a small bit of cleverness */
997: VecCopy(gl->X[0],ts->vec_sol);
998: ts->ptime += h;
999: ts->steps++;
1001: TSPostEvaluate(ts);
1002: TSPostStep(ts);
1003: TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);
1005: gl->current_scheme = next_scheme;
1006: ts->time_step = next_h;
1007: }
1008: return(0);
1009: }
1011: /*------------------------------------------------------------*/
1013: static PetscErrorCode TSReset_GLLE(TS ts)
1014: {
1015: TS_GLLE *gl = (TS_GLLE*)ts->data;
1016: PetscInt max_r,max_s;
1020: if (gl->setupcalled) {
1021: TSGLLEGetMaxSizes(ts,&max_r,&max_s);
1022: VecDestroyVecs(max_r,&gl->Xold);
1023: VecDestroyVecs(max_r,&gl->X);
1024: VecDestroyVecs(max_s,&gl->Ydot);
1025: VecDestroyVecs(3,&gl->himom);
1026: VecDestroy(&gl->W);
1027: VecDestroy(&gl->Y);
1028: VecDestroy(&gl->Z);
1029: }
1030: gl->setupcalled = PETSC_FALSE;
1031: return(0);
1032: }
1034: static PetscErrorCode TSDestroy_GLLE(TS ts)
1035: {
1036: TS_GLLE *gl = (TS_GLLE*)ts->data;
1040: TSReset_GLLE(ts);
1041: if (ts->dm) {
1042: DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSGLLE,DMRestrictHook_TSGLLE,ts);
1043: DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSGLLE,DMSubDomainRestrictHook_TSGLLE,ts);
1044: }
1045: if (gl->adapt) {TSGLLEAdaptDestroy(&gl->adapt);}
1046: if (gl->Destroy) {(*gl->Destroy)(gl);}
1047: PetscFree(ts->data);
1048: PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetType_C",NULL);
1049: PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetAcceptType_C",NULL);
1050: PetscObjectComposeFunction((PetscObject)ts,"TSGLLEGetAdapt_C",NULL);
1051: return(0);
1052: }
1054: /*
1055: This defines the nonlinear equation that is to be solved with SNES
1056: g(x) = f(t,x,z+shift*x) = 0
1057: */
1058: static PetscErrorCode SNESTSFormFunction_GLLE(SNES snes,Vec x,Vec f,TS ts)
1059: {
1060: TS_GLLE *gl = (TS_GLLE*)ts->data;
1062: Vec Z,Ydot;
1063: DM dm,dmsave;
1066: SNESGetDM(snes,&dm);
1067: TSGLLEGetVecs(ts,dm,&Z,&Ydot);
1068: VecWAXPY(Ydot,gl->scoeff/ts->time_step,x,Z);
1069: dmsave = ts->dm;
1070: ts->dm = dm;
1071: TSComputeIFunction(ts,gl->stage_time,x,Ydot,f,PETSC_FALSE);
1072: ts->dm = dmsave;
1073: TSGLLERestoreVecs(ts,dm,&Z,&Ydot);
1074: return(0);
1075: }
1077: static PetscErrorCode SNESTSFormJacobian_GLLE(SNES snes,Vec x,Mat A,Mat B,TS ts)
1078: {
1079: TS_GLLE *gl = (TS_GLLE*)ts->data;
1081: Vec Z,Ydot;
1082: DM dm,dmsave;
1085: SNESGetDM(snes,&dm);
1086: TSGLLEGetVecs(ts,dm,&Z,&Ydot);
1087: dmsave = ts->dm;
1088: ts->dm = dm;
1089: /* gl->Xdot will have already been computed in SNESTSFormFunction_GLLE */
1090: TSComputeIJacobian(ts,gl->stage_time,x,gl->Ydot[gl->stage],gl->scoeff/ts->time_step,A,B,PETSC_FALSE);
1091: ts->dm = dmsave;
1092: TSGLLERestoreVecs(ts,dm,&Z,&Ydot);
1093: return(0);
1094: }
1097: static PetscErrorCode TSSetUp_GLLE(TS ts)
1098: {
1099: TS_GLLE *gl = (TS_GLLE*)ts->data;
1100: PetscInt max_r,max_s;
1102: DM dm;
1105: gl->setupcalled = PETSC_TRUE;
1106: TSGLLEGetMaxSizes(ts,&max_r,&max_s);
1107: VecDuplicateVecs(ts->vec_sol,max_r,&gl->X);
1108: VecDuplicateVecs(ts->vec_sol,max_r,&gl->Xold);
1109: VecDuplicateVecs(ts->vec_sol,max_s,&gl->Ydot);
1110: VecDuplicateVecs(ts->vec_sol,3,&gl->himom);
1111: VecDuplicate(ts->vec_sol,&gl->W);
1112: VecDuplicate(ts->vec_sol,&gl->Y);
1113: VecDuplicate(ts->vec_sol,&gl->Z);
1115: /* Default acceptance tests and adaptivity */
1116: if (!gl->Accept) {TSGLLESetAcceptType(ts,TSGLLEACCEPT_ALWAYS);}
1117: if (!gl->adapt) {TSGLLEGetAdapt(ts,&gl->adapt);}
1119: if (gl->current_scheme < 0) {
1120: PetscInt i;
1121: for (i=0;; i++) {
1122: if (gl->schemes[i]->p == gl->start_order) break;
1123: if (i+1 == gl->nschemes) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No schemes available with requested start order %d",i);
1124: }
1125: gl->current_scheme = i;
1126: }
1127: TSGetDM(ts,&dm);
1128: DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLLE,DMRestrictHook_TSGLLE,ts);
1129: DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLLE,DMSubDomainRestrictHook_TSGLLE,ts);
1130: return(0);
1131: }
1132: /*------------------------------------------------------------*/
1134: static PetscErrorCode TSSetFromOptions_GLLE(PetscOptionItems *PetscOptionsObject,TS ts)
1135: {
1136: TS_GLLE *gl = (TS_GLLE*)ts->data;
1137: char tname[256] = TSGLLE_IRKS,completef[256] = "rescale-and-modify";
1141: PetscOptionsHead(PetscOptionsObject,"General Linear ODE solver options");
1142: {
1143: PetscBool flg;
1144: PetscOptionsFList("-ts_gl_type","Type of GL method","TSGLLESetType",TSGLLEList,gl->type_name[0] ? gl->type_name : tname,tname,sizeof(tname),&flg);
1145: if (flg || !gl->type_name[0]) {
1146: TSGLLESetType(ts,tname);
1147: }
1148: PetscOptionsInt("-ts_gl_max_step_rejections","Maximum number of times to attempt a step","None",gl->max_step_rejections,&gl->max_step_rejections,NULL);
1149: PetscOptionsInt("-ts_gl_max_order","Maximum order to try","TSGLLESetMaxOrder",gl->max_order,&gl->max_order,NULL);
1150: PetscOptionsInt("-ts_gl_min_order","Minimum order to try","TSGLLESetMinOrder",gl->min_order,&gl->min_order,NULL);
1151: PetscOptionsInt("-ts_gl_start_order","Initial order to try","TSGLLESetMinOrder",gl->start_order,&gl->start_order,NULL);
1152: PetscOptionsEnum("-ts_gl_error_direction","Which direction to look when estimating error","TSGLLESetErrorDirection",TSGLLEErrorDirections,(PetscEnum)gl->error_direction,(PetscEnum*)&gl->error_direction,NULL);
1153: PetscOptionsBool("-ts_gl_extrapolate","Extrapolate stage solution from previous solution (sometimes unstable)","TSGLLESetExtrapolate",gl->extrapolate,&gl->extrapolate,NULL);
1154: PetscOptionsReal("-ts_gl_atol","Absolute tolerance","TSGLLESetTolerances",gl->wrms_atol,&gl->wrms_atol,NULL);
1155: PetscOptionsReal("-ts_gl_rtol","Relative tolerance","TSGLLESetTolerances",gl->wrms_rtol,&gl->wrms_rtol,NULL);
1156: PetscOptionsString("-ts_gl_complete","Method to use for completing the step","none",completef,completef,sizeof(completef),&flg);
1157: if (flg) {
1158: PetscBool match1,match2;
1159: PetscStrcmp(completef,"rescale",&match1);
1160: PetscStrcmp(completef,"rescale-and-modify",&match2);
1161: if (match1) gl->CompleteStep = TSGLLECompleteStep_Rescale;
1162: else if (match2) gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify;
1163: else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"%s",completef);
1164: }
1165: {
1166: char type[256] = TSGLLEACCEPT_ALWAYS;
1167: PetscOptionsFList("-ts_gl_accept_type","Method to use for determining whether to accept a step","TSGLLESetAcceptType",TSGLLEAcceptList,gl->accept_name[0] ? gl->accept_name : type,type,sizeof(type),&flg);
1168: if (flg || !gl->accept_name[0]) {
1169: TSGLLESetAcceptType(ts,type);
1170: }
1171: }
1172: {
1173: TSGLLEAdapt adapt;
1174: TSGLLEGetAdapt(ts,&adapt);
1175: TSGLLEAdaptSetFromOptions(PetscOptionsObject,adapt);
1176: }
1177: }
1178: PetscOptionsTail();
1179: return(0);
1180: }
1182: static PetscErrorCode TSView_GLLE(TS ts,PetscViewer viewer)
1183: {
1184: TS_GLLE *gl = (TS_GLLE*)ts->data;
1185: PetscInt i;
1186: PetscBool iascii,details;
1190: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1191: if (iascii) {
1192: PetscViewerASCIIPrintf(viewer," min order %D, max order %D, current order %D\n",gl->min_order,gl->max_order,gl->schemes[gl->current_scheme]->p);
1193: PetscViewerASCIIPrintf(viewer," Error estimation: %s\n",TSGLLEErrorDirections[gl->error_direction]);
1194: PetscViewerASCIIPrintf(viewer," Extrapolation: %s\n",gl->extrapolate ? "yes" : "no");
1195: PetscViewerASCIIPrintf(viewer," Acceptance test: %s\n",gl->accept_name[0] ? gl->accept_name : "(not yet set)");
1196: PetscViewerASCIIPushTab(viewer);
1197: TSGLLEAdaptView(gl->adapt,viewer);
1198: PetscViewerASCIIPopTab(viewer);
1199: PetscViewerASCIIPrintf(viewer," type: %s\n",gl->type_name[0] ? gl->type_name : "(not yet set)");
1200: PetscViewerASCIIPrintf(viewer,"Schemes within family (%d):\n",gl->nschemes);
1201: details = PETSC_FALSE;
1202: PetscOptionsGetBool(((PetscObject)ts)->options,((PetscObject)ts)->prefix,"-ts_gl_view_detailed",&details,NULL);
1203: PetscViewerASCIIPushTab(viewer);
1204: for (i=0; i<gl->nschemes; i++) {
1205: TSGLLESchemeView(gl->schemes[i],details,viewer);
1206: }
1207: if (gl->View) {
1208: (*gl->View)(gl,viewer);
1209: }
1210: PetscViewerASCIIPopTab(viewer);
1211: }
1212: return(0);
1213: }
1215: /*@C
1216: TSGLLERegister - adds a TSGLLE implementation
1218: Not Collective
1220: Input Parameters:
1221: + name_scheme - name of user-defined general linear scheme
1222: - routine_create - routine to create method context
1224: Notes:
1225: TSGLLERegister() may be called multiple times to add several user-defined families.
1227: Sample usage:
1228: .vb
1229: TSGLLERegister("my_scheme",MySchemeCreate);
1230: .ve
1232: Then, your scheme can be chosen with the procedural interface via
1233: $ TSGLLESetType(ts,"my_scheme")
1234: or at runtime via the option
1235: $ -ts_gl_type my_scheme
1237: Level: advanced
1239: .keywords: TSGLLE, register
1241: .seealso: TSGLLERegisterAll()
1242: @*/
1243: PetscErrorCode TSGLLERegister(const char sname[],PetscErrorCode (*function)(TS))
1244: {
1248: PetscFunctionListAdd(&TSGLLEList,sname,function);
1249: return(0);
1250: }
1252: /*@C
1253: TSGLLEAcceptRegister - adds a TSGLLE acceptance scheme
1255: Not Collective
1257: Input Parameters:
1258: + name_scheme - name of user-defined acceptance scheme
1259: - routine_create - routine to create method context
1261: Notes:
1262: TSGLLEAcceptRegister() may be called multiple times to add several user-defined families.
1264: Sample usage:
1265: .vb
1266: TSGLLEAcceptRegister("my_scheme",MySchemeCreate);
1267: .ve
1269: Then, your scheme can be chosen with the procedural interface via
1270: $ TSGLLESetAcceptType(ts,"my_scheme")
1271: or at runtime via the option
1272: $ -ts_gl_accept_type my_scheme
1274: Level: advanced
1276: .keywords: TSGLLE, TSGLLEAcceptType, register
1278: .seealso: TSGLLERegisterAll()
1279: @*/
1280: PetscErrorCode TSGLLEAcceptRegister(const char sname[],TSGLLEAcceptFunction function)
1281: {
1285: PetscFunctionListAdd(&TSGLLEAcceptList,sname,function);
1286: return(0);
1287: }
1289: /*@C
1290: TSGLLERegisterAll - Registers all of the general linear methods in TSGLLE
1292: Not Collective
1294: Level: advanced
1296: .keywords: TS, TSGLLE, register, all
1298: .seealso: TSGLLERegisterDestroy()
1299: @*/
1300: PetscErrorCode TSGLLERegisterAll(void)
1301: {
1305: if (TSGLLERegisterAllCalled) return(0);
1306: TSGLLERegisterAllCalled = PETSC_TRUE;
1308: TSGLLERegister(TSGLLE_IRKS, TSGLLECreate_IRKS);
1309: TSGLLEAcceptRegister(TSGLLEACCEPT_ALWAYS,TSGLLEAccept_Always);
1310: return(0);
1311: }
1313: /*@C
1314: TSGLLEInitializePackage - This function initializes everything in the TSGLLE package. It is called
1315: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_GLLE()
1316: when using static libraries.
1318: Level: developer
1320: .keywords: TS, TSGLLE, initialize, package
1321: .seealso: PetscInitialize()
1322: @*/
1323: PetscErrorCode TSGLLEInitializePackage(void)
1324: {
1328: if (TSGLLEPackageInitialized) return(0);
1329: TSGLLEPackageInitialized = PETSC_TRUE;
1330: TSGLLERegisterAll();
1331: PetscRegisterFinalize(TSGLLEFinalizePackage);
1332: return(0);
1333: }
1335: /*@C
1336: TSGLLEFinalizePackage - This function destroys everything in the TSGLLE package. It is
1337: called from PetscFinalize().
1339: Level: developer
1341: .keywords: Petsc, destroy, package
1342: .seealso: PetscFinalize()
1343: @*/
1344: PetscErrorCode TSGLLEFinalizePackage(void)
1345: {
1349: PetscFunctionListDestroy(&TSGLLEList);
1350: PetscFunctionListDestroy(&TSGLLEAcceptList);
1351: TSGLLEPackageInitialized = PETSC_FALSE;
1352: TSGLLERegisterAllCalled = PETSC_FALSE;
1353: return(0);
1354: }
1356: /* ------------------------------------------------------------ */
1357: /*MC
1358: TSGLLE - DAE solver using implicit General Linear methods
1360: These methods contain Runge-Kutta and multistep schemes as special cases. These special cases have some fundamental
1361: limitations. For example, diagonally implicit Runge-Kutta cannot have stage order greater than 1 which limits their
1362: applicability to very stiff systems. Meanwhile, multistep methods cannot be A-stable for order greater than 2 and BDF
1363: are not 0-stable for order greater than 6. GL methods can be A- and L-stable with arbitrarily high stage order and
1364: reliable error estimates for both 1 and 2 orders higher to facilitate adaptive step sizes and adaptive order schemes.
1365: All this is possible while preserving a singly diagonally implicit structure.
1367: Options database keys:
1368: + -ts_gl_type <type> - the class of general linear method (irks)
1369: . -ts_gl_rtol <tol> - relative error
1370: . -ts_gl_atol <tol> - absolute error
1371: . -ts_gl_min_order <p> - minimum order method to consider (default=1)
1372: . -ts_gl_max_order <p> - maximum order method to consider (default=3)
1373: . -ts_gl_start_order <p> - order of starting method (default=1)
1374: . -ts_gl_complete <method> - method to use for completing the step (rescale-and-modify or rescale)
1375: - -ts_adapt_type <method> - adaptive controller to use (none step both)
1377: Notes:
1378: This integrator can be applied to DAE.
1380: Diagonally implicit general linear (DIGL) methods are a generalization of diagonally implicit Runge-Kutta (DIRK).
1381: They are represented by the tableau
1383: .vb
1384: A | U
1385: -------
1386: B | V
1387: .ve
1389: combined with a vector c of abscissa. "Diagonally implicit" means that A is lower triangular.
1390: A step of the general method reads
1392: .vb
1393: [ Y ] = [A U] [ Y' ]
1394: [X^k] = [B V] [X^{k-1}]
1395: .ve
1397: where Y is the multivector of stage values, Y' is the multivector of stage derivatives, X^k is the Nordsieck vector of
1398: the solution at step k. The Nordsieck vector consists of the first r moments of the solution, given by
1400: .vb
1401: X = [x_0,x_1,...,x_{r-1}] = [x, h x', h^2 x'', ..., h^{r-1} x^{(r-1)} ]
1402: .ve
1404: If A is lower triangular, we can solve the stages (Y,Y') sequentially
1406: .vb
1407: 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}
1408: .ve
1410: and then construct the pieces to carry to the next step
1412: .vb
1413: 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}
1414: .ve
1416: Note that when the equations are cast in implicit form, we are using the stage equation to define y'_i
1417: in terms of y_i and known stuff (y_j for j<i and x_j for all j).
1420: Error estimation
1422: At present, the most attractive GL methods for stiff problems are singly diagonally implicit schemes which posses
1423: Inherent Runge-Kutta Stability (IRKS). These methods have r=s, the number of items passed between steps is equal to
1424: the number of stages. The order and stage-order are one less than the number of stages. We use the error estimates
1425: in the 2007 paper which provide the following estimates
1427: .vb
1428: h^{p+1} X^{(p+1)} = phi_0^T Y' + [0 psi_0^T] Xold
1429: h^{p+2} X^{(p+2)} = phi_1^T Y' + [0 psi_1^T] Xold
1430: h^{p+2} (dx'/dx) X^{(p+1)} = phi_2^T Y' + [0 psi_2^T] Xold
1431: .ve
1433: These estimates are accurate to O(h^{p+3}).
1435: Changing the step size
1437: We use the generalized "rescale and modify" scheme, see equation (4.5) of the 2007 paper.
1439: Level: beginner
1441: References:
1442: + 1. - John Butcher and Z. Jackieweicz and W. Wright, On error propagation in general linear methods for
1443: ordinary differential equations, Journal of Complexity, Vol 23, 2007.
1444: - 2. - John Butcher, Numerical methods for ordinary differential equations, second edition, Wiley, 2009.
1446: .seealso: TSCreate(), TS, TSSetType()
1448: M*/
1449: PETSC_EXTERN PetscErrorCode TSCreate_GLLE(TS ts)
1450: {
1451: TS_GLLE *gl;
1455: TSGLLEInitializePackage();
1457: PetscNewLog(ts,&gl);
1458: ts->data = (void*)gl;
1460: ts->ops->reset = TSReset_GLLE;
1461: ts->ops->destroy = TSDestroy_GLLE;
1462: ts->ops->view = TSView_GLLE;
1463: ts->ops->setup = TSSetUp_GLLE;
1464: ts->ops->solve = TSSolve_GLLE;
1465: ts->ops->setfromoptions = TSSetFromOptions_GLLE;
1466: ts->ops->snesfunction = SNESTSFormFunction_GLLE;
1467: ts->ops->snesjacobian = SNESTSFormJacobian_GLLE;
1469: ts->usessnes = PETSC_TRUE;
1472: gl->max_step_rejections = 1;
1473: gl->min_order = 1;
1474: gl->max_order = 3;
1475: gl->start_order = 1;
1476: gl->current_scheme = -1;
1477: gl->extrapolate = PETSC_FALSE;
1479: gl->wrms_atol = 1e-8;
1480: gl->wrms_rtol = 1e-5;
1482: PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetType_C", &TSGLLESetType_GLLE);
1483: PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetAcceptType_C",&TSGLLESetAcceptType_GLLE);
1484: PetscObjectComposeFunction((PetscObject)ts,"TSGLLEGetAdapt_C", &TSGLLEGetAdapt_GLLE);
1485: return(0);
1486: }