Actual source code: glee.c
petsc-3.9.4 2018-09-11
1: /*
2: Code for time stepping with the General Linear with Error Estimation method
5: Notes:
6: The general system is written as
8: Udot = F(t,U)
10: */
11: #include <petsc/private/tsimpl.h>
12: #include <petscdm.h>
14: static PetscBool cited = PETSC_FALSE;
15: static const char citation[] =
16: "@ARTICLE{Constantinescu_TR2016b,\n"
17: " author = {Constantinescu, E.M.},\n"
18: " title = {Estimating Global Errors in Time Stepping},\n"
19: " journal = {ArXiv e-prints},\n"
20: " year = 2016,\n"
21: " adsurl = {http://adsabs.harvard.edu/abs/2015arXiv150305166C}\n}\n";
23: static TSGLEEType TSGLEEDefaultType = TSGLEE35;
24: static PetscBool TSGLEERegisterAllCalled;
25: static PetscBool TSGLEEPackageInitialized;
26: static PetscInt explicit_stage_time_id;
28: typedef struct _GLEETableau *GLEETableau;
29: struct _GLEETableau {
30: char *name;
31: PetscInt order; /* Classical approximation order of the method i*/
32: PetscInt s; /* Number of stages */
33: PetscInt r; /* Number of steps */
34: PetscReal gamma; /* LTE ratio */
35: PetscReal *A,*B,*U,*V,*S,*F,*c; /* Tableau */
36: PetscReal *Fembed; /* Embedded final method coefficients */
37: PetscReal *Ferror; /* Coefficients for computing error */
38: PetscReal *Serror; /* Coefficients for initializing the error */
39: PetscInt pinterp; /* Interpolation order */
40: PetscReal *binterp; /* Interpolation coefficients */
41: PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */
42: };
43: typedef struct _GLEETableauLink *GLEETableauLink;
44: struct _GLEETableauLink {
45: struct _GLEETableau tab;
46: GLEETableauLink next;
47: };
48: static GLEETableauLink GLEETableauList;
50: typedef struct {
51: GLEETableau tableau;
52: Vec *Y; /* Solution vector (along with auxiliary solution y~ or eps) */
53: Vec *X; /* Temporary solution vector */
54: Vec *YStage; /* Stage values */
55: Vec *YdotStage; /* Stage right hand side */
56: Vec W; /* Right-hand-side for implicit stage solve */
57: Vec Ydot; /* Work vector holding Ydot during residual evaluation */
58: Vec yGErr; /* Vector holding the global error after a step is completed */
59: PetscScalar *swork; /* Scalar work (size of the number of stages)*/
60: PetscScalar *rwork; /* Scalar work (size of the number of steps)*/
61: PetscReal scoeff; /* shift = scoeff/dt */
62: PetscReal stage_time;
63: TSStepStatus status;
64: } TS_GLEE;
66: /*MC
67: TSGLEE23 - Second order three stage GLEE method
69: This method has three stages.
70: s = 3, r = 2
72: Level: advanced
74: .seealso: TSGLEE
75: M*/
76: /*MC
77: TSGLEE24 - Second order four stage GLEE method
79: This method has four stages.
80: s = 4, r = 2
82: Level: advanced
84: .seealso: TSGLEE
85: M*/
86: /*MC
87: TSGLEE25i - Second order five stage GLEE method
89: This method has five stages.
90: s = 5, r = 2
92: Level: advanced
94: .seealso: TSGLEE
95: M*/
96: /*MC
97: TSGLEE35 - Third order five stage GLEE method
99: This method has five stages.
100: s = 5, r = 2
102: Level: advanced
104: .seealso: TSGLEE
105: M*/
106: /*MC
107: TSGLEEEXRK2A - Second order six stage GLEE method
109: This method has six stages.
110: s = 6, r = 2
112: Level: advanced
114: .seealso: TSGLEE
115: M*/
116: /*MC
117: TSGLEERK32G1 - Third order eight stage GLEE method
119: This method has eight stages.
120: s = 8, r = 2
122: Level: advanced
124: .seealso: TSGLEE
125: M*/
126: /*MC
127: TSGLEERK285EX - Second order nine stage GLEE method
129: This method has nine stages.
130: s = 9, r = 2
132: Level: advanced
134: .seealso: TSGLEE
135: M*/
137: /*@C
138: TSGLEERegisterAll - Registers all of the General Linear with Error Estimation methods in TSGLEE
140: Not Collective, but should be called by all processes which will need the schemes to be registered
142: Level: advanced
144: .keywords: TS, TSGLEE, register, all
146: .seealso: TSGLEERegisterDestroy()
147: @*/
148: PetscErrorCode TSGLEERegisterAll(void)
149: {
153: if (TSGLEERegisterAllCalled) return(0);
154: TSGLEERegisterAllCalled = PETSC_TRUE;
156: {
157: #define GAMMA 0.5
158: /* y-eps form */
159: const PetscInt
160: p = 1,
161: s = 3,
162: r = 2;
163: const PetscReal
164: A[3][3] = {{1.0,0,0},{0,0.5,0},{0,0.5,0.5}},
165: B[2][3] = {{1.0,0,0},{-2.0,1.0,1.0}},
166: U[3][2] = {{1.0,0},{1.0,0.5},{1.0,0.5}},
167: V[2][2] = {{1,0},{0,1}},
168: S[2] = {1,0},
169: F[2] = {1,0},
170: Fembed[2] = {1,1-GAMMA},
171: Ferror[2] = {0,1},
172: Serror[2] = {1,0};
173: TSGLEERegister(TSGLEEi1,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
174: }
175: {
176: #undef GAMMA
177: #define GAMMA 0.0
178: /* y-eps form */
179: const PetscInt
180: p = 2,
181: s = 3,
182: r = 2;
183: const PetscReal
184: A[3][3] = {{0,0,0},{1,0,0},{0.25,0.25,0}},
185: B[2][3] = {{1.0/12.0,1.0/12.0,5.0/6.0},{1.0/12.0,1.0/12.0,-1.0/6.0}},
186: U[3][2] = {{1,0},{1,10},{1,-1}},
187: V[2][2] = {{1,0},{0,1}},
188: S[2] = {1,0},
189: F[2] = {1,0},
190: Fembed[2] = {1,1-GAMMA},
191: Ferror[2] = {0,1},
192: Serror[2] = {1,0};
193: TSGLEERegister(TSGLEE23,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
194: }
195: {
196: #undef GAMMA
197: #define GAMMA 0.0
198: /* y-y~ form */
199: const PetscInt
200: p = 2,
201: s = 4,
202: r = 2;
203: const PetscReal
204: A[4][4] = {{0,0,0,0},{0.75,0,0,0},{0.25,29.0/60.0,0,0},{-21.0/44.0,145.0/44.0,-20.0/11.0,0}},
205: B[2][4] = {{109.0/275.0,58.0/75.0,-37.0/110.0,1.0/6.0},{3.0/11.0,0,75.0/88.0,-1.0/8.0}},
206: U[4][2] = {{0,1},{75.0/58.0,-17.0/58.0},{0,1},{0,1}},
207: V[2][2] = {{1,0},{0,1}},
208: S[2] = {1,1},
209: F[2] = {1,0},
210: Fembed[2] = {0,1},
211: Ferror[2] = {-1.0/(1.0-GAMMA),1.0/(1.0-GAMMA)},
212: Serror[2] = {1.0-GAMMA,1.0};
213: TSGLEERegister(TSGLEE24,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
214: }
215: {
216: #undef GAMMA
217: #define GAMMA 0.0
218: /* y-y~ form */
219: const PetscInt
220: p = 2,
221: s = 5,
222: r = 2;
223: const PetscReal
224: A[5][5] = {{0,0,0,0,0},
225: {-0.94079244066783383269,0,0,0,0},
226: {0.64228187778301907108,0.10915356933958500042,0,0,0},
227: {-0.51764297742287450812,0.74414270351096040738,-0.71404164927824538121,0,0},
228: {-0.44696561556825969206,-0.76768425657590196518,0.20111608138142987881,0.93828186737840469796,0}},
229: B[2][5] = {{-0.029309178948150356153,-0.49671981884013874923,0.34275801517650053274,0.32941112623949194988,0.85385985637229662276},
230: {0.78133219686062535272,0.074238691892675897635,0.57957363498384957966,-0.24638502829674959968,-0.18875949544040123033}},
231: U[5][2] = {{0.16911424754448327735,0.83088575245551672265},
232: {0.53638465733199574340,0.46361534266800425660},
233: {0.39901579167169582526,0.60098420832830417474},
234: {0.87689005530618575480,0.12310994469381424520},
235: {0.99056100455550913009,0.0094389954444908699092}},
236: V[2][2] = {{1,0},{0,1}},
237: S[2] = {1,1},
238: F[2] = {1,0},
239: Fembed[2] = {0,1},
240: Ferror[2] = {-1.0/(1.0-GAMMA),1.0/(1.0-GAMMA)},
241: Serror[2] = {1.0-GAMMA,1.0};
242: TSGLEERegister(TSGLEE25I,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
243: }
244: {
245: #undef GAMMA
246: #define GAMMA 0.0
247: /* y-y~ form */
248: const PetscInt
249: p = 3,
250: s = 5,
251: r = 2;
252: const PetscReal
253: A[5][5] = {{0,0,0,0,0},
254: {- 2169604947363702313.0 / 24313474998937147335.0,0,0,0,0},
255: {46526746497697123895.0 / 94116917485856474137.0,-10297879244026594958.0 / 49199457603717988219.0,0,0,0},
256: {23364788935845982499.0 / 87425311444725389446.0,-79205144337496116638.0 / 148994349441340815519.0,40051189859317443782.0 / 36487615018004984309.0,0,0},
257: {42089522664062539205.0 / 124911313006412840286.0,-15074384760342762939.0 / 137927286865289746282.0,-62274678522253371016.0 / 125918573676298591413.0,13755475729852471739.0 / 79257927066651693390.0,0}},
258: B[2][5] = {{61546696837458703723.0 / 56982519523786160813.0,-55810892792806293355.0 / 206957624151308356511.0,24061048952676379087.0 / 158739347956038723465.0,3577972206874351339.0 / 7599733370677197135.0,-59449832954780563947.0 / 137360038685338563670.0},
259: {- 9738262186984159168.0 / 99299082461487742983.0,-32797097931948613195.0 / 61521565616362163366.0,42895514606418420631.0 / 71714201188501437336.0,22608567633166065068.0 / 55371917805607957003.0,94655809487476459565.0 / 151517167160302729021.0}},
260: U[5][2] = {{70820309139834661559.0 / 80863923579509469826.0,10043614439674808267.0 / 80863923579509469826.0},
261: {161694774978034105510.0 / 106187653640211060371.0,-55507121337823045139.0 / 106187653640211060371.0},
262: {78486094644566264568.0 / 88171030896733822981.0,9684936252167558413.0 / 88171030896733822981.0},
263: {65394922146334854435.0 / 84570853840405479554.0,19175931694070625119.0 / 84570853840405479554.0},
264: {8607282770183754108.0 / 108658046436496925911.0,100050763666313171803.0 / 108658046436496925911.0}},
265: V[2][2] = {{1,0},{0,1}},
266: S[2] = {1,1},
267: F[2] = {1,0},
268: Fembed[2] = {0,1},
269: Ferror[2] = {-1.0/(1.0-GAMMA),1.0/(1.0-GAMMA)},
270: Serror[2] = {1.0-GAMMA,1.0};
271: TSGLEERegister(TSGLEE35,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
272: }
273: {
274: #undef GAMMA
275: #define GAMMA 0.25
276: /* y-eps form */
277: const PetscInt
278: p = 2,
279: s = 6,
280: r = 2;
281: const PetscReal
282: A[6][6] = {{0,0,0,0,0,0},
283: {1,0,0,0,0,0},
284: {0,0,0,0,0,0},
285: {0,0,0.5,0,0,0},
286: {0,0,0.25,0.25,0,0},
287: {0,0,0.25,0.25,0.5,0}},
288: B[2][6] = {{0.5,0.5,0,0,0,0},
289: {-2.0/3.0,-2.0/3.0,1.0/3.0,1.0/3.0,1.0/3.0,1.0/3.0}},
290: U[6][2] = {{1,0},{1,0},{1,0.75},{1,0.75},{1,0.75},{1,0.75}},
291: V[2][2] = {{1,0},{0,1}},
292: S[2] = {1,0},
293: F[2] = {1,0},
294: Fembed[2] = {1,1-GAMMA},
295: Ferror[2] = {0,1},
296: Serror[2] = {1,0};
297: TSGLEERegister(TSGLEEEXRK2A,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
298: }
299: {
300: #undef GAMMA
301: #define GAMMA 0.0
302: /* y-eps form */
303: const PetscInt
304: p = 3,
305: s = 8,
306: r = 2;
307: const PetscReal
308: A[8][8] = {{0,0,0,0,0,0,0,0},
309: {0.5,0,0,0,0,0,0,0},
310: {-1,2,0,0,0,0,0,0},
311: {1.0/6.0,2.0/3.0,1.0/6.0,0,0,0,0,0},
312: {0,0,0,0,0,0,0,0},
313: {-7.0/24.0,1.0/3.0,1.0/12.0,-1.0/8.0,0.5,0,0,0},
314: {7.0/6.0,-4.0/3.0,-1.0/3.0,0.5,-1.0,2.0,0,0},
315: {0,0,0,0,1.0/6.0,2.0/3.0,1.0/6.0,0}},
316: B[2][8] = {{1.0/6.0,2.0/3.0,1.0/6.0,0,0,0,0,0},
317: {-1.0/6.0,-2.0/3.0,-1.0/6.0,0,1.0/6.0,2.0/3.0,1.0/6.0,0}},
318: U[8][2] = {{1,0},{1,0},{1,0},{1,0},{1,1},{1,1},{1,1},{1,1}},
319: V[2][2] = {{1,0},{0,1}},
320: S[2] = {1,0},
321: F[2] = {1,0},
322: Fembed[2] = {1,1-GAMMA},
323: Ferror[2] = {0,1},
324: Serror[2] = {1,0};
325: TSGLEERegister(TSGLEERK32G1,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
326: }
327: {
328: #undef GAMMA
329: #define GAMMA 0.25
330: /* y-eps form */
331: const PetscInt
332: p = 2,
333: s = 9,
334: r = 2;
335: const PetscReal
336: A[9][9] = {{0,0,0,0,0,0,0,0,0},
337: {0.585786437626904966,0,0,0,0,0,0,0,0},
338: {0.149999999999999994,0.849999999999999978,0,0,0,0,0,0,0},
339: {0,0,0,0,0,0,0,0,0},
340: {0,0,0,0.292893218813452483,0,0,0,0,0},
341: {0,0,0,0.0749999999999999972,0.424999999999999989,0,0,0,0},
342: {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0,0,0},
343: {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0.292893218813452483,0,0},
344: {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0.0749999999999999972,0.424999999999999989,0}},
345: B[2][9] = {{0.353553390593273786,0.353553390593273786,0.292893218813452483,0,0,0,0,0,0},
346: {-0.471404520791031678,-0.471404520791031678,-0.390524291751269959,0.235702260395515839,0.235702260395515839,0.195262145875634979,0.235702260395515839,0.235702260395515839,0.195262145875634979}},
347: U[9][2] = {{1,0},{1,0},{1,0},{1,0.75},{1,0.75},{1,0.75},{1,0.75},{1,0.75},{1,0.75}},
348: V[2][2] = {{1,0},{0,1}},
349: S[2] = {1,0},
350: F[2] = {1,0},
351: Fembed[2] = {1,1-GAMMA},
352: Ferror[2] = {0,1},
353: Serror[2] = {1,0};
354: TSGLEERegister(TSGLEERK285EX,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
355: }
357: return(0);
358: }
360: /*@C
361: TSGLEERegisterDestroy - Frees the list of schemes that were registered by TSGLEERegister().
363: Not Collective
365: Level: advanced
367: .keywords: TSGLEE, register, destroy
368: .seealso: TSGLEERegister(), TSGLEERegisterAll()
369: @*/
370: PetscErrorCode TSGLEERegisterDestroy(void)
371: {
373: GLEETableauLink link;
376: while ((link = GLEETableauList)) {
377: GLEETableau t = &link->tab;
378: GLEETableauList = link->next;
379: PetscFree5(t->A,t->B,t->U,t->V,t->c);
380: PetscFree2(t->S,t->F);
381: PetscFree (t->Fembed);
382: PetscFree (t->Ferror);
383: PetscFree (t->Serror);
384: PetscFree (t->binterp);
385: PetscFree (t->name);
386: PetscFree (link);
387: }
388: TSGLEERegisterAllCalled = PETSC_FALSE;
389: return(0);
390: }
392: /*@C
393: TSGLEEInitializePackage - This function initializes everything in the TSGLEE package. It is called
394: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_GLEE()
395: when using static libraries.
397: Level: developer
399: .keywords: TS, TSGLEE, initialize, package
400: .seealso: PetscInitialize()
401: @*/
402: PetscErrorCode TSGLEEInitializePackage(void)
403: {
407: if (TSGLEEPackageInitialized) return(0);
408: TSGLEEPackageInitialized = PETSC_TRUE;
409: TSGLEERegisterAll();
410: PetscObjectComposedDataRegister(&explicit_stage_time_id);
411: PetscRegisterFinalize(TSGLEEFinalizePackage);
412: return(0);
413: }
415: /*@C
416: TSGLEEFinalizePackage - This function destroys everything in the TSGLEE package. It is
417: called from PetscFinalize().
419: Level: developer
421: .keywords: Petsc, destroy, package
422: .seealso: PetscFinalize()
423: @*/
424: PetscErrorCode TSGLEEFinalizePackage(void)
425: {
429: TSGLEEPackageInitialized = PETSC_FALSE;
430: TSGLEERegisterDestroy();
431: return(0);
432: }
434: /*@C
435: TSGLEERegister - register an GLEE scheme by providing the entries in the Butcher tableau
437: Not Collective, but the same schemes should be registered on all processes on which they will be used
439: Input Parameters:
440: + name - identifier for method
441: . order - order of method
442: . s - number of stages
443: . r - number of steps
444: . gamma - LTE ratio
445: . A - stage coefficients (dimension s*s, row-major)
446: . B - step completion coefficients (dimension r*s, row-major)
447: . U - method coefficients (dimension s*r, row-major)
448: . V - method coefficients (dimension r*r, row-major)
449: . S - starting coefficients
450: . F - finishing coefficients
451: . c - abscissa (dimension s; NULL to use row sums of A)
452: . Fembed - step completion coefficients for embedded method
453: . Ferror - error computation coefficients
454: . Serror - error initialization coefficients
455: . pinterp - order of interpolation (0 if unavailable)
456: - binterp - array of interpolation coefficients (NULL if unavailable)
458: Notes:
459: Several GLEE methods are provided, this function is only needed to create new methods.
461: Level: advanced
463: .keywords: TS, register
465: .seealso: TSGLEE
466: @*/
467: PetscErrorCode TSGLEERegister(TSGLEEType name,PetscInt order,PetscInt s, PetscInt r,
468: PetscReal gamma,
469: const PetscReal A[],const PetscReal B[],
470: const PetscReal U[],const PetscReal V[],
471: const PetscReal S[],const PetscReal F[],
472: const PetscReal c[],
473: const PetscReal Fembed[],const PetscReal Ferror[],
474: const PetscReal Serror[],
475: PetscInt pinterp, const PetscReal binterp[])
476: {
477: PetscErrorCode ierr;
478: GLEETableauLink link;
479: GLEETableau t;
480: PetscInt i,j;
483: PetscMalloc(sizeof(*link),&link);
484: PetscMemzero(link,sizeof(*link));
485: t = &link->tab;
486: PetscStrallocpy(name,&t->name);
487: t->order = order;
488: t->s = s;
489: t->r = r;
490: t->gamma = gamma;
491: PetscMalloc5(s*s,&t->A,r*r,&t->V,s,&t->c,r*s,&t->B,s*r,&t->U);
492: PetscMalloc2(r,&t->S,r,&t->F);
493: PetscMemcpy(t->A,A,s*s*sizeof(A[0]));
494: PetscMemcpy(t->B,B,r*s*sizeof(B[0]));
495: PetscMemcpy(t->U,U,s*r*sizeof(U[0]));
496: PetscMemcpy(t->V,V,r*r*sizeof(V[0]));
497: PetscMemcpy(t->S,S,r *sizeof(S[0]));
498: PetscMemcpy(t->F,F,r *sizeof(F[0]));
499: if (c) {
500: PetscMemcpy(t->c,c,s*sizeof(c[0]));
501: } else {
502: for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j];
503: }
504: PetscMalloc1(r,&t->Fembed);
505: PetscMalloc1(r,&t->Ferror);
506: PetscMalloc1(r,&t->Serror);
507: PetscMemcpy(t->Fembed,Fembed,r*sizeof(Fembed[0]));
508: PetscMemcpy(t->Ferror,Ferror,r*sizeof(Ferror[0]));
509: PetscMemcpy(t->Serror,Serror,r*sizeof(Serror[0]));
510: t->pinterp = pinterp;
511: PetscMalloc1(s*pinterp,&t->binterp);
512: PetscMemcpy(t->binterp,binterp,s*pinterp*sizeof(binterp[0]));
514: link->next = GLEETableauList;
515: GLEETableauList = link;
516: return(0);
517: }
519: static PetscErrorCode TSEvaluateStep_GLEE(TS ts,PetscInt order,Vec X,PetscBool *done)
520: {
521: TS_GLEE *glee = (TS_GLEE*) ts->data;
522: GLEETableau tab = glee->tableau;
523: PetscReal h, *B = tab->B, *V = tab->V,
524: *F = tab->F,
525: *Fembed = tab->Fembed;
526: PetscInt s = tab->s, r = tab->r, i, j;
527: Vec *Y = glee->Y, *YdotStage = glee->YdotStage;
528: PetscScalar *ws = glee->swork, *wr = glee->rwork;
529: PetscErrorCode ierr;
533: switch (glee->status) {
534: case TS_STEP_INCOMPLETE:
535: case TS_STEP_PENDING:
536: h = ts->time_step; break;
537: case TS_STEP_COMPLETE:
538: h = ts->ptime - ts->ptime_prev; break;
539: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
540: }
542: if (order == tab->order) {
544: /* Note: Irrespective of whether status is TS_STEP_INCOMPLETE
545: or TS_STEP_COMPLETE, glee->X has the solution at the
546: beginning of the time step. So no need to roll-back.
547: */
548: if (glee->status == TS_STEP_INCOMPLETE) {
549: for (i=0; i<r; i++) {
550: VecZeroEntries(Y[i]);
551: for (j=0; j<r; j++) wr[j] = V[i*r+j];
552: VecMAXPY(Y[i],r,wr,glee->X);
553: for (j=0; j<s; j++) ws[j] = h*B[i*s+j];
554: VecMAXPY(Y[i],s,ws,YdotStage);
555: }
556: VecZeroEntries(X);
557: for (j=0; j<r; j++) wr[j] = F[j];
558: VecMAXPY(X,r,wr,Y);
559: } else {VecCopy(ts->vec_sol,X);}
560: return(0);
562: } else if (order == tab->order-1) {
564: /* Complete with the embedded method (Fembed) */
565: for (i=0; i<r; i++) {
566: VecZeroEntries(Y[i]);
567: for (j=0; j<r; j++) wr[j] = V[i*r+j];
568: VecMAXPY(Y[i],r,wr,glee->X);
569: for (j=0; j<s; j++) ws[j] = h*B[i*s+j];
570: VecMAXPY(Y[i],s,ws,YdotStage);
571: }
572: VecZeroEntries(X);
573: for (j=0; j<r; j++) wr[j] = Fembed[j];
574: VecMAXPY(X,r,wr,Y);
576: if (done) *done = PETSC_TRUE;
577: return(0);
578: }
579: if (done) *done = PETSC_FALSE;
580: else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"GLEE '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order);
581: return(0);
582: }
584: static PetscErrorCode TSStep_GLEE(TS ts)
585: {
586: TS_GLEE *glee = (TS_GLEE*)ts->data;
587: GLEETableau tab = glee->tableau;
588: const PetscInt s = tab->s, r = tab->r;
589: PetscReal *A = tab->A, *U = tab->U,
590: *F = tab->F,
591: *c = tab->c;
592: Vec *Y = glee->Y, *X = glee->X,
593: *YStage = glee->YStage,
594: *YdotStage = glee->YdotStage,
595: W = glee->W;
596: SNES snes;
597: PetscScalar *ws = glee->swork, *wr = glee->rwork;
598: TSAdapt adapt;
599: PetscInt i,j,reject,next_scheme,its,lits;
600: PetscReal next_time_step;
601: PetscReal t;
602: PetscBool accept;
603: PetscErrorCode ierr;
606: PetscCitationsRegister(citation,&cited);
608: for (i=0; i<r; i++) { VecCopy(Y[i],X[i]); }
610: TSGetSNES(ts,&snes);
611: next_time_step = ts->time_step;
612: t = ts->ptime;
613: accept = PETSC_TRUE;
614: glee->status = TS_STEP_INCOMPLETE;
616: for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) {
618: PetscReal h = ts->time_step;
619: TSPreStep(ts);
621: for (i=0; i<s; i++) {
623: glee->stage_time = t + h*c[i];
624: TSPreStage(ts,glee->stage_time);
626: if (A[i*s+i] == 0) { /* Explicit stage */
627: VecZeroEntries(YStage[i]);
628: for (j=0; j<r; j++) wr[j] = U[i*r+j];
629: VecMAXPY(YStage[i],r,wr,X);
630: for (j=0; j<i; j++) ws[j] = h*A[i*s+j];
631: VecMAXPY(YStage[i],i,ws,YdotStage);
632: } else { /* Implicit stage */
633: glee->scoeff = 1.0/A[i*s+i];
634: /* compute right-hand-side */
635: VecZeroEntries(W);
636: for (j=0; j<r; j++) wr[j] = U[i*r+j];
637: VecMAXPY(W,r,wr,X);
638: for (j=0; j<i; j++) ws[j] = h*A[i*s+j];
639: VecMAXPY(W,i,ws,YdotStage);
640: VecScale(W,glee->scoeff/h);
641: /* set initial guess */
642: VecCopy(i>0 ? YStage[i-1] : ts->vec_sol,YStage[i]);
643: /* solve for this stage */
644: SNESSolve(snes,W,YStage[i]);
645: SNESGetIterationNumber(snes,&its);
646: SNESGetLinearSolveIterations(snes,&lits);
647: ts->snes_its += its; ts->ksp_its += lits;
648: }
649: TSGetAdapt(ts,&adapt);
650: TSAdaptCheckStage(adapt,ts,glee->stage_time,YStage[i],&accept);
651: if (!accept) goto reject_step;
652: TSPostStage(ts,glee->stage_time,i,YStage);
653: TSComputeRHSFunction(ts,t+h*c[i],YStage[i],YdotStage[i]);
654: }
655: TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);
656: glee->status = TS_STEP_PENDING;
658: /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */
659: TSGetAdapt(ts,&adapt);
660: TSAdaptCandidatesClear(adapt);
661: TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);
662: TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);
663: if (accept) {
664: /* ignore next_scheme for now */
665: ts->ptime += ts->time_step;
666: ts->time_step = next_time_step;
667: glee->status = TS_STEP_COMPLETE;
668: /* compute and store the global error */
669: /* Note: this is not needed if TSAdaptGLEE is not used */
670: TSGetTimeError(ts,0,&(glee->yGErr));
671: PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);
672: break;
673: } else { /* Roll back the current step */
674: for (j=0; j<r; j++) wr[j] = F[j];
675: VecMAXPY(ts->vec_sol,r,wr,X);
676: ts->time_step = next_time_step;
677: glee->status = TS_STEP_INCOMPLETE;
678: }
679: reject_step: continue;
680: }
681: if (glee->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
682: return(0);
683: }
685: static PetscErrorCode TSInterpolate_GLEE(TS ts,PetscReal itime,Vec X)
686: {
687: TS_GLEE *glee = (TS_GLEE*)ts->data;
688: PetscInt s=glee->tableau->s, pinterp=glee->tableau->pinterp,i,j;
689: PetscReal h,tt,t;
690: PetscScalar *b;
691: const PetscReal *B = glee->tableau->binterp;
692: PetscErrorCode ierr;
695: if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSGLEE %s does not have an interpolation formula",glee->tableau->name);
696: switch (glee->status) {
697: case TS_STEP_INCOMPLETE:
698: case TS_STEP_PENDING:
699: h = ts->time_step;
700: t = (itime - ts->ptime)/h;
701: break;
702: case TS_STEP_COMPLETE:
703: h = ts->ptime - ts->ptime_prev;
704: t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
705: break;
706: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
707: }
708: PetscMalloc1(s,&b);
709: for (i=0; i<s; i++) b[i] = 0;
710: for (j=0,tt=t; j<pinterp; j++,tt*=t) {
711: for (i=0; i<s; i++) {
712: b[i] += h * B[i*pinterp+j] * tt;
713: }
714: }
715: VecCopy(glee->YStage[0],X);
716: VecMAXPY(X,s,b,glee->YdotStage);
717: PetscFree(b);
718: return(0);
719: }
721: /*------------------------------------------------------------*/
722: static PetscErrorCode TSReset_GLEE(TS ts)
723: {
724: TS_GLEE *glee = (TS_GLEE*)ts->data;
725: PetscInt s, r;
729: if (!glee->tableau) return(0);
730: s = glee->tableau->s;
731: r = glee->tableau->r;
732: VecDestroyVecs(r,&glee->Y);
733: VecDestroyVecs(r,&glee->X);
734: VecDestroyVecs(s,&glee->YStage);
735: VecDestroyVecs(s,&glee->YdotStage);
736: VecDestroy(&glee->Ydot);
737: VecDestroy(&glee->yGErr);
738: VecDestroy(&glee->W);
739: PetscFree2(glee->swork,glee->rwork);
740: return(0);
741: }
743: static PetscErrorCode TSGLEEGetVecs(TS ts,DM dm,Vec *Ydot)
744: {
745: TS_GLEE *glee = (TS_GLEE*)ts->data;
749: if (Ydot) {
750: if (dm && dm != ts->dm) {
751: DMGetNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);
752: } else *Ydot = glee->Ydot;
753: }
754: return(0);
755: }
758: static PetscErrorCode TSGLEERestoreVecs(TS ts,DM dm,Vec *Ydot)
759: {
763: if (Ydot) {
764: if (dm && dm != ts->dm) {
765: DMRestoreNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);
766: }
767: }
768: return(0);
769: }
771: /*
772: This defines the nonlinear equation that is to be solved with SNES
773: */
774: static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes,Vec X,Vec F,TS ts)
775: {
776: TS_GLEE *glee = (TS_GLEE*)ts->data;
777: DM dm,dmsave;
778: Vec Ydot;
779: PetscReal shift = glee->scoeff / ts->time_step;
783: SNESGetDM(snes,&dm);
784: TSGLEEGetVecs(ts,dm,&Ydot);
785: /* Set Ydot = shift*X */
786: VecCopy(X,Ydot);
787: VecScale(Ydot,shift);
788: dmsave = ts->dm;
789: ts->dm = dm;
791: TSComputeIFunction(ts,glee->stage_time,X,Ydot,F,PETSC_FALSE);
793: ts->dm = dmsave;
794: TSGLEERestoreVecs(ts,dm,&Ydot);
795: return(0);
796: }
798: static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes,Vec X,Mat A,Mat B,TS ts)
799: {
800: TS_GLEE *glee = (TS_GLEE*)ts->data;
801: DM dm,dmsave;
802: Vec Ydot;
803: PetscReal shift = glee->scoeff / ts->time_step;
807: SNESGetDM(snes,&dm);
808: TSGLEEGetVecs(ts,dm,&Ydot);
809: /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */
810: dmsave = ts->dm;
811: ts->dm = dm;
813: TSComputeIJacobian(ts,glee->stage_time,X,Ydot,shift,A,B,PETSC_FALSE);
815: ts->dm = dmsave;
816: TSGLEERestoreVecs(ts,dm,&Ydot);
817: return(0);
818: }
820: static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine,DM coarse,void *ctx)
821: {
823: return(0);
824: }
826: static PetscErrorCode DMRestrictHook_TSGLEE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
827: {
829: return(0);
830: }
833: static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm,DM subdm,void *ctx)
834: {
836: return(0);
837: }
839: static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
840: {
842: return(0);
843: }
845: static PetscErrorCode TSSetUp_GLEE(TS ts)
846: {
847: TS_GLEE *glee = (TS_GLEE*)ts->data;
848: GLEETableau tab;
849: PetscInt s,r;
851: DM dm;
854: if (!glee->tableau) {
855: TSGLEESetType(ts,TSGLEEDefaultType);
856: }
857: tab = glee->tableau;
858: s = tab->s;
859: r = tab->r;
860: VecDuplicateVecs(ts->vec_sol,r,&glee->Y);
861: VecDuplicateVecs(ts->vec_sol,r,&glee->X);
862: VecDuplicateVecs(ts->vec_sol,s,&glee->YStage);
863: VecDuplicateVecs(ts->vec_sol,s,&glee->YdotStage);
864: VecDuplicate(ts->vec_sol,&glee->Ydot);
865: VecDuplicate(ts->vec_sol,&glee->yGErr);
866: VecZeroEntries(glee->yGErr);
867: VecDuplicate(ts->vec_sol,&glee->W);
868: PetscMalloc2(s,&glee->swork,r,&glee->rwork);
869: TSGetDM(ts,&dm);
870: DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);
871: DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);
872: return(0);
873: }
875: PetscErrorCode TSStartingMethod_GLEE(TS ts)
876: {
877: TS_GLEE *glee = (TS_GLEE*)ts->data;
878: GLEETableau tab = glee->tableau;
879: PetscInt r=tab->r,i;
880: PetscReal *S=tab->S;
884: for (i=0; i<r; i++) {
885: VecZeroEntries(glee->Y[i]);
886: VecAXPY(glee->Y[i],S[i],ts->vec_sol);
887: }
889: return(0);
890: }
893: /*------------------------------------------------------------*/
895: static PetscErrorCode TSSetFromOptions_GLEE(PetscOptionItems *PetscOptionsObject,TS ts)
896: {
898: char gleetype[256];
901: PetscOptionsHead(PetscOptionsObject,"GLEE ODE solver options");
902: {
903: GLEETableauLink link;
904: PetscInt count,choice;
905: PetscBool flg;
906: const char **namelist;
908: PetscStrncpy(gleetype,TSGLEEDefaultType,sizeof(gleetype));
909: for (link=GLEETableauList,count=0; link; link=link->next,count++) ;
910: PetscMalloc1(count,(char***)&namelist);
911: for (link=GLEETableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
912: PetscOptionsEList("-ts_glee_type","Family of GLEE method","TSGLEESetType",(const char*const*)namelist,count,gleetype,&choice,&flg);
913: TSGLEESetType(ts,flg ? namelist[choice] : gleetype);
914: PetscFree(namelist);
915: }
916: PetscOptionsTail();
917: return(0);
918: }
920: static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[])
921: {
923: PetscInt i;
924: size_t left,count;
925: char *p;
928: for (i=0,p=buf,left=len; i<n; i++) {
929: PetscSNPrintfCount(p,left,fmt,&count,x[i]);
930: if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer");
931: left -= count;
932: p += count;
933: *p++ = ' ';
934: }
935: p[i ? 0 : -1] = 0;
936: return(0);
937: }
939: static PetscErrorCode TSView_GLEE(TS ts,PetscViewer viewer)
940: {
941: TS_GLEE *glee = (TS_GLEE*)ts->data;
942: GLEETableau tab = glee->tableau;
943: PetscBool iascii;
947: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
948: if (iascii) {
949: TSGLEEType gleetype;
950: char buf[512];
951: TSGLEEGetType(ts,&gleetype);
952: PetscViewerASCIIPrintf(viewer," GLEE type %s\n",gleetype);
953: PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);
954: PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);
955: /* Note: print out r as well */
956: }
957: return(0);
958: }
960: static PetscErrorCode TSLoad_GLEE(TS ts,PetscViewer viewer)
961: {
963: SNES snes;
964: TSAdapt tsadapt;
967: TSGetAdapt(ts,&tsadapt);
968: TSAdaptLoad(tsadapt,viewer);
969: TSGetSNES(ts,&snes);
970: SNESLoad(snes,viewer);
971: /* function and Jacobian context for SNES when used with TS is always ts object */
972: SNESSetFunction(snes,NULL,NULL,ts);
973: SNESSetJacobian(snes,NULL,NULL,NULL,ts);
974: return(0);
975: }
977: /*@C
978: TSGLEESetType - Set the type of GLEE scheme
980: Logically collective
982: Input Parameter:
983: + ts - timestepping context
984: - gleetype - type of GLEE-scheme
986: Level: intermediate
988: .seealso: TSGLEEGetType(), TSGLEE
989: @*/
990: PetscErrorCode TSGLEESetType(TS ts,TSGLEEType gleetype)
991: {
997: PetscTryMethod(ts,"TSGLEESetType_C",(TS,TSGLEEType),(ts,gleetype));
998: return(0);
999: }
1001: /*@C
1002: TSGLEEGetType - Get the type of GLEE scheme
1004: Logically collective
1006: Input Parameter:
1007: . ts - timestepping context
1009: Output Parameter:
1010: . gleetype - type of GLEE-scheme
1012: Level: intermediate
1014: .seealso: TSGLEESetType()
1015: @*/
1016: PetscErrorCode TSGLEEGetType(TS ts,TSGLEEType *gleetype)
1017: {
1022: PetscUseMethod(ts,"TSGLEEGetType_C",(TS,TSGLEEType*),(ts,gleetype));
1023: return(0);
1024: }
1026: PetscErrorCode TSGLEEGetType_GLEE(TS ts,TSGLEEType *gleetype)
1027: {
1028: TS_GLEE *glee = (TS_GLEE*)ts->data;
1032: if (!glee->tableau) {
1033: TSGLEESetType(ts,TSGLEEDefaultType);
1034: }
1035: *gleetype = glee->tableau->name;
1036: return(0);
1037: }
1038: PetscErrorCode TSGLEESetType_GLEE(TS ts,TSGLEEType gleetype)
1039: {
1040: TS_GLEE *glee = (TS_GLEE*)ts->data;
1041: PetscErrorCode ierr;
1042: PetscBool match;
1043: GLEETableauLink link;
1046: if (glee->tableau) {
1047: PetscStrcmp(glee->tableau->name,gleetype,&match);
1048: if (match) return(0);
1049: }
1050: for (link = GLEETableauList; link; link=link->next) {
1051: PetscStrcmp(link->tab.name,gleetype,&match);
1052: if (match) {
1053: TSReset_GLEE(ts);
1054: glee->tableau = &link->tab;
1055: return(0);
1056: }
1057: }
1058: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",gleetype);
1059: return(0);
1060: }
1062: static PetscErrorCode TSGetStages_GLEE(TS ts,PetscInt *ns,Vec **Y)
1063: {
1064: TS_GLEE *glee = (TS_GLEE*)ts->data;
1067: *ns = glee->tableau->s;
1068: if(Y) *Y = glee->YStage;
1069: return(0);
1070: }
1072: PetscErrorCode TSGetSolutionComponents_GLEE(TS ts,PetscInt *n,Vec *Y)
1073: {
1074: TS_GLEE *glee = (TS_GLEE*)ts->data;
1075: GLEETableau tab = glee->tableau;
1076: PetscErrorCode ierr;
1079: if (!Y) *n = tab->r;
1080: else {
1081: if ((*n >= 0) && (*n < tab->r)) {
1082: VecCopy(glee->Y[*n],*Y);
1083: } else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Second argument (%d) out of range[%d,%d].",*n,0,tab->r-1);
1084: }
1085: return(0);
1086: }
1088: PetscErrorCode TSGetAuxSolution_GLEE(TS ts,Vec *X)
1089: {
1090: TS_GLEE *glee = (TS_GLEE*)ts->data;
1091: GLEETableau tab = glee->tableau;
1092: PetscReal *F = tab->Fembed;
1093: PetscInt r = tab->r;
1094: Vec *Y = glee->Y;
1095: PetscScalar *wr = glee->rwork;
1096: PetscInt i;
1097: PetscErrorCode ierr;
1100: VecZeroEntries(*X);
1101: for (i=0; i<r; i++) wr[i] = F[i];
1102: VecMAXPY((*X),r,wr,Y);
1103: return(0);
1104: }
1106: PetscErrorCode TSGetTimeError_GLEE(TS ts,PetscInt n,Vec *X)
1107: {
1108: TS_GLEE *glee = (TS_GLEE*)ts->data;
1109: GLEETableau tab = glee->tableau;
1110: PetscReal *F = tab->Ferror;
1111: PetscInt r = tab->r;
1112: Vec *Y = glee->Y;
1113: PetscScalar *wr = glee->rwork;
1114: PetscInt i;
1115: PetscErrorCode ierr;
1118: VecZeroEntries(*X);
1119: if(n==0){
1120: for (i=0; i<r; i++) wr[i] = F[i];
1121: VecMAXPY((*X),r,wr,Y);
1122: } else if(n==-1) {
1123: *X=glee->yGErr;
1124: }
1125: return(0);
1126: }
1128: PetscErrorCode TSSetTimeError_GLEE(TS ts,Vec X)
1129: {
1130: TS_GLEE *glee = (TS_GLEE*)ts->data;
1131: GLEETableau tab = glee->tableau;
1132: PetscReal *S = tab->Serror;
1133: PetscInt r = tab->r,i;
1134: Vec *Y = glee->Y;
1135: PetscErrorCode ierr;
1138: if (r != 2) SETERRQ2(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSSetTimeError_GLEE not supported for '%s' with r=%D.",tab->name,tab->r);
1139: for (i=1; i<r; i++) {
1140: VecCopy(ts->vec_sol,Y[i]);
1141: VecAXPBY(Y[i],S[0],S[1],X);
1142: VecCopy(X,glee->yGErr);
1143: }
1144: return(0);
1145: }
1147: static PetscErrorCode TSDestroy_GLEE(TS ts)
1148: {
1152: TSReset_GLEE(ts);
1153: if (ts->dm) {
1154: DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);
1155: DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);
1156: }
1157: PetscFree(ts->data);
1158: PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",NULL);
1159: PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",NULL);
1160: return(0);
1161: }
1163: /* ------------------------------------------------------------ */
1164: /*MC
1165: TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes
1167: The user should provide the right hand side of the equation
1168: using TSSetRHSFunction().
1170: Notes:
1171: The default is TSGLEE35, it can be changed with TSGLEESetType() or -ts_glee_type
1173: Level: beginner
1175: .seealso: TSCreate(), TS, TSSetType(), TSGLEESetType(), TSGLEEGetType(),
1176: TSGLEE23, TTSGLEE24, TSGLEE35, TSGLEE25I, TSGLEEEXRK2A,
1177: TSGLEERK32G1, TSGLEERK285EX, TSGLEEType, TSGLEERegister()
1179: M*/
1180: PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts)
1181: {
1182: TS_GLEE *th;
1183: PetscErrorCode ierr;
1186: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
1187: TSGLEEInitializePackage();
1188: #endif
1190: ts->ops->reset = TSReset_GLEE;
1191: ts->ops->destroy = TSDestroy_GLEE;
1192: ts->ops->view = TSView_GLEE;
1193: ts->ops->load = TSLoad_GLEE;
1194: ts->ops->setup = TSSetUp_GLEE;
1195: ts->ops->step = TSStep_GLEE;
1196: ts->ops->interpolate = TSInterpolate_GLEE;
1197: ts->ops->evaluatestep = TSEvaluateStep_GLEE;
1198: ts->ops->setfromoptions = TSSetFromOptions_GLEE;
1199: ts->ops->getstages = TSGetStages_GLEE;
1200: ts->ops->snesfunction = SNESTSFormFunction_GLEE;
1201: ts->ops->snesjacobian = SNESTSFormJacobian_GLEE;
1202: ts->ops->getsolutioncomponents = TSGetSolutionComponents_GLEE;
1203: ts->ops->getauxsolution = TSGetAuxSolution_GLEE;
1204: ts->ops->gettimeerror = TSGetTimeError_GLEE;
1205: ts->ops->settimeerror = TSSetTimeError_GLEE;
1206: ts->ops->startingmethod = TSStartingMethod_GLEE;
1207: ts->default_adapt_type = TSADAPTGLEE;
1209: ts->usessnes = PETSC_TRUE;
1211: PetscNewLog(ts,&th);
1212: ts->data = (void*)th;
1214: PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",TSGLEEGetType_GLEE);
1215: PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",TSGLEESetType_GLEE);
1216: return(0);
1217: }