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