Actual source code: glee.c
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: .seealso: TSGLEERegisterDestroy()
145: @*/
146: PetscErrorCode TSGLEERegisterAll(void)
147: {
151: if (TSGLEERegisterAllCalled) return(0);
152: TSGLEERegisterAllCalled = PETSC_TRUE;
154: {
155: #define GAMMA 0.5
156: /* y-eps form */
157: const PetscInt
158: p = 1,
159: s = 3,
160: r = 2;
161: const PetscReal
162: A[3][3] = {{1.0,0,0},{0,0.5,0},{0,0.5,0.5}},
163: B[2][3] = {{1.0,0,0},{-2.0,1.0,1.0}},
164: U[3][2] = {{1.0,0},{1.0,0.5},{1.0,0.5}},
165: V[2][2] = {{1,0},{0,1}},
166: S[2] = {1,0},
167: F[2] = {1,0},
168: Fembed[2] = {1,1-GAMMA},
169: Ferror[2] = {0,1},
170: Serror[2] = {1,0};
171: 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);
172: }
173: {
174: #undef GAMMA
175: #define GAMMA 0.0
176: /* y-eps form */
177: const PetscInt
178: p = 2,
179: s = 3,
180: r = 2;
181: const PetscReal
182: A[3][3] = {{0,0,0},{1,0,0},{0.25,0.25,0}},
183: 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}},
184: U[3][2] = {{1,0},{1,10},{1,-1}},
185: V[2][2] = {{1,0},{0,1}},
186: S[2] = {1,0},
187: F[2] = {1,0},
188: Fembed[2] = {1,1-GAMMA},
189: Ferror[2] = {0,1},
190: Serror[2] = {1,0};
191: 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);
192: }
193: {
194: #undef GAMMA
195: #define GAMMA 0.0
196: /* y-y~ form */
197: const PetscInt
198: p = 2,
199: s = 4,
200: r = 2;
201: const PetscReal
202: 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}},
203: 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}},
204: U[4][2] = {{0,1},{75.0/58.0,-17.0/58.0},{0,1},{0,1}},
205: V[2][2] = {{1,0},{0,1}},
206: S[2] = {1,1},
207: F[2] = {1,0},
208: Fembed[2] = {0,1},
209: Ferror[2] = {-1.0/(1.0-GAMMA),1.0/(1.0-GAMMA)},
210: Serror[2] = {1.0-GAMMA,1.0};
211: 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);
212: }
213: {
214: #undef GAMMA
215: #define GAMMA 0.0
216: /* y-y~ form */
217: const PetscInt
218: p = 2,
219: s = 5,
220: r = 2;
221: const PetscReal
222: A[5][5] = {{0,0,0,0,0},
223: {-0.94079244066783383269,0,0,0,0},
224: {0.64228187778301907108,0.10915356933958500042,0,0,0},
225: {-0.51764297742287450812,0.74414270351096040738,-0.71404164927824538121,0,0},
226: {-0.44696561556825969206,-0.76768425657590196518,0.20111608138142987881,0.93828186737840469796,0}},
227: B[2][5] = {{-0.029309178948150356153,-0.49671981884013874923,0.34275801517650053274,0.32941112623949194988,0.85385985637229662276},
228: {0.78133219686062535272,0.074238691892675897635,0.57957363498384957966,-0.24638502829674959968,-0.18875949544040123033}},
229: U[5][2] = {{0.16911424754448327735,0.83088575245551672265},
230: {0.53638465733199574340,0.46361534266800425660},
231: {0.39901579167169582526,0.60098420832830417474},
232: {0.87689005530618575480,0.12310994469381424520},
233: {0.99056100455550913009,0.0094389954444908699092}},
234: V[2][2] = {{1,0},{0,1}},
235: S[2] = {1,1},
236: F[2] = {1,0},
237: Fembed[2] = {0,1},
238: Ferror[2] = {-1.0/(1.0-GAMMA),1.0/(1.0-GAMMA)},
239: Serror[2] = {1.0-GAMMA,1.0};
240: 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);
241: }
242: {
243: #undef GAMMA
244: #define GAMMA 0.0
245: /* y-y~ form */
246: const PetscInt
247: p = 3,
248: s = 5,
249: r = 2;
250: const PetscReal
251: A[5][5] = {{0,0,0,0,0},
252: {- 2169604947363702313.0 / 24313474998937147335.0,0,0,0,0},
253: {46526746497697123895.0 / 94116917485856474137.0,-10297879244026594958.0 / 49199457603717988219.0,0,0,0},
254: {23364788935845982499.0 / 87425311444725389446.0,-79205144337496116638.0 / 148994349441340815519.0,40051189859317443782.0 / 36487615018004984309.0,0,0},
255: {42089522664062539205.0 / 124911313006412840286.0,-15074384760342762939.0 / 137927286865289746282.0,-62274678522253371016.0 / 125918573676298591413.0,13755475729852471739.0 / 79257927066651693390.0,0}},
256: B[2][5] = {{61546696837458703723.0 / 56982519523786160813.0,-55810892792806293355.0 / 206957624151308356511.0,24061048952676379087.0 / 158739347956038723465.0,3577972206874351339.0 / 7599733370677197135.0,-59449832954780563947.0 / 137360038685338563670.0},
257: {- 9738262186984159168.0 / 99299082461487742983.0,-32797097931948613195.0 / 61521565616362163366.0,42895514606418420631.0 / 71714201188501437336.0,22608567633166065068.0 / 55371917805607957003.0,94655809487476459565.0 / 151517167160302729021.0}},
258: U[5][2] = {{70820309139834661559.0 / 80863923579509469826.0,10043614439674808267.0 / 80863923579509469826.0},
259: {161694774978034105510.0 / 106187653640211060371.0,-55507121337823045139.0 / 106187653640211060371.0},
260: {78486094644566264568.0 / 88171030896733822981.0,9684936252167558413.0 / 88171030896733822981.0},
261: {65394922146334854435.0 / 84570853840405479554.0,19175931694070625119.0 / 84570853840405479554.0},
262: {8607282770183754108.0 / 108658046436496925911.0,100050763666313171803.0 / 108658046436496925911.0}},
263: V[2][2] = {{1,0},{0,1}},
264: S[2] = {1,1},
265: F[2] = {1,0},
266: Fembed[2] = {0,1},
267: Ferror[2] = {-1.0/(1.0-GAMMA),1.0/(1.0-GAMMA)},
268: Serror[2] = {1.0-GAMMA,1.0};
269: 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);
270: }
271: {
272: #undef GAMMA
273: #define GAMMA 0.25
274: /* y-eps form */
275: const PetscInt
276: p = 2,
277: s = 6,
278: r = 2;
279: const PetscReal
280: A[6][6] = {{0,0,0,0,0,0},
281: {1,0,0,0,0,0},
282: {0,0,0,0,0,0},
283: {0,0,0.5,0,0,0},
284: {0,0,0.25,0.25,0,0},
285: {0,0,0.25,0.25,0.5,0}},
286: B[2][6] = {{0.5,0.5,0,0,0,0},
287: {-2.0/3.0,-2.0/3.0,1.0/3.0,1.0/3.0,1.0/3.0,1.0/3.0}},
288: U[6][2] = {{1,0},{1,0},{1,0.75},{1,0.75},{1,0.75},{1,0.75}},
289: V[2][2] = {{1,0},{0,1}},
290: S[2] = {1,0},
291: F[2] = {1,0},
292: Fembed[2] = {1,1-GAMMA},
293: Ferror[2] = {0,1},
294: Serror[2] = {1,0};
295: 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);
296: }
297: {
298: #undef GAMMA
299: #define GAMMA 0.0
300: /* y-eps form */
301: const PetscInt
302: p = 3,
303: s = 8,
304: r = 2;
305: const PetscReal
306: A[8][8] = {{0,0,0,0,0,0,0,0},
307: {0.5,0,0,0,0,0,0,0},
308: {-1,2,0,0,0,0,0,0},
309: {1.0/6.0,2.0/3.0,1.0/6.0,0,0,0,0,0},
310: {0,0,0,0,0,0,0,0},
311: {-7.0/24.0,1.0/3.0,1.0/12.0,-1.0/8.0,0.5,0,0,0},
312: {7.0/6.0,-4.0/3.0,-1.0/3.0,0.5,-1.0,2.0,0,0},
313: {0,0,0,0,1.0/6.0,2.0/3.0,1.0/6.0,0}},
314: B[2][8] = {{1.0/6.0,2.0/3.0,1.0/6.0,0,0,0,0,0},
315: {-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}},
316: U[8][2] = {{1,0},{1,0},{1,0},{1,0},{1,1},{1,1},{1,1},{1,1}},
317: V[2][2] = {{1,0},{0,1}},
318: S[2] = {1,0},
319: F[2] = {1,0},
320: Fembed[2] = {1,1-GAMMA},
321: Ferror[2] = {0,1},
322: Serror[2] = {1,0};
323: 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);
324: }
325: {
326: #undef GAMMA
327: #define GAMMA 0.25
328: /* y-eps form */
329: const PetscInt
330: p = 2,
331: s = 9,
332: r = 2;
333: const PetscReal
334: A[9][9] = {{0,0,0,0,0,0,0,0,0},
335: {0.585786437626904966,0,0,0,0,0,0,0,0},
336: {0.149999999999999994,0.849999999999999978,0,0,0,0,0,0,0},
337: {0,0,0,0,0,0,0,0,0},
338: {0,0,0,0.292893218813452483,0,0,0,0,0},
339: {0,0,0,0.0749999999999999972,0.424999999999999989,0,0,0,0},
340: {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0,0,0},
341: {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0.292893218813452483,0,0},
342: {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0.0749999999999999972,0.424999999999999989,0}},
343: B[2][9] = {{0.353553390593273786,0.353553390593273786,0.292893218813452483,0,0,0,0,0,0},
344: {-0.471404520791031678,-0.471404520791031678,-0.390524291751269959,0.235702260395515839,0.235702260395515839,0.195262145875634979,0.235702260395515839,0.235702260395515839,0.195262145875634979}},
345: 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}},
346: V[2][2] = {{1,0},{0,1}},
347: S[2] = {1,0},
348: F[2] = {1,0},
349: Fembed[2] = {1,1-GAMMA},
350: Ferror[2] = {0,1},
351: Serror[2] = {1,0};
352: 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);
353: }
355: return(0);
356: }
358: /*@C
359: TSGLEERegisterDestroy - Frees the list of schemes that were registered by TSGLEERegister().
361: Not Collective
363: Level: advanced
365: .seealso: TSGLEERegister(), TSGLEERegisterAll()
366: @*/
367: PetscErrorCode TSGLEERegisterDestroy(void)
368: {
370: GLEETableauLink link;
373: while ((link = GLEETableauList)) {
374: GLEETableau t = &link->tab;
375: GLEETableauList = link->next;
376: PetscFree5(t->A,t->B,t->U,t->V,t->c);
377: PetscFree2(t->S,t->F);
378: PetscFree (t->Fembed);
379: PetscFree (t->Ferror);
380: PetscFree (t->Serror);
381: PetscFree (t->binterp);
382: PetscFree (t->name);
383: PetscFree (link);
384: }
385: TSGLEERegisterAllCalled = PETSC_FALSE;
386: return(0);
387: }
389: /*@C
390: TSGLEEInitializePackage - This function initializes everything in the TSGLEE package. It is called
391: from TSInitializePackage().
393: Level: developer
395: .seealso: PetscInitialize()
396: @*/
397: PetscErrorCode TSGLEEInitializePackage(void)
398: {
402: if (TSGLEEPackageInitialized) return(0);
403: TSGLEEPackageInitialized = PETSC_TRUE;
404: TSGLEERegisterAll();
405: PetscObjectComposedDataRegister(&explicit_stage_time_id);
406: PetscRegisterFinalize(TSGLEEFinalizePackage);
407: return(0);
408: }
410: /*@C
411: TSGLEEFinalizePackage - This function destroys everything in the TSGLEE package. It is
412: called from PetscFinalize().
414: Level: developer
416: .seealso: PetscFinalize()
417: @*/
418: PetscErrorCode TSGLEEFinalizePackage(void)
419: {
423: TSGLEEPackageInitialized = PETSC_FALSE;
424: TSGLEERegisterDestroy();
425: return(0);
426: }
428: /*@C
429: TSGLEERegister - register an GLEE scheme by providing the entries in the Butcher tableau
431: Not Collective, but the same schemes should be registered on all processes on which they will be used
433: Input Parameters:
434: + name - identifier for method
435: . order - order of method
436: . s - number of stages
437: . r - number of steps
438: . gamma - LTE ratio
439: . A - stage coefficients (dimension s*s, row-major)
440: . B - step completion coefficients (dimension r*s, row-major)
441: . U - method coefficients (dimension s*r, row-major)
442: . V - method coefficients (dimension r*r, row-major)
443: . S - starting coefficients
444: . F - finishing coefficients
445: . c - abscissa (dimension s; NULL to use row sums of A)
446: . Fembed - step completion coefficients for embedded method
447: . Ferror - error computation coefficients
448: . Serror - error initialization coefficients
449: . pinterp - order of interpolation (0 if unavailable)
450: - binterp - array of interpolation coefficients (NULL if unavailable)
452: Notes:
453: Several GLEE methods are provided, this function is only needed to create new methods.
455: Level: advanced
457: .seealso: TSGLEE
458: @*/
459: PetscErrorCode TSGLEERegister(TSGLEEType name,PetscInt order,PetscInt s, PetscInt r,
460: PetscReal gamma,
461: const PetscReal A[],const PetscReal B[],
462: const PetscReal U[],const PetscReal V[],
463: const PetscReal S[],const PetscReal F[],
464: const PetscReal c[],
465: const PetscReal Fembed[],const PetscReal Ferror[],
466: const PetscReal Serror[],
467: PetscInt pinterp, const PetscReal binterp[])
468: {
469: PetscErrorCode ierr;
470: GLEETableauLink link;
471: GLEETableau t;
472: PetscInt i,j;
475: TSGLEEInitializePackage();
476: PetscNew(&link);
477: t = &link->tab;
478: PetscStrallocpy(name,&t->name);
479: t->order = order;
480: t->s = s;
481: t->r = r;
482: t->gamma = gamma;
483: PetscMalloc5(s*s,&t->A,r*r,&t->V,s,&t->c,r*s,&t->B,s*r,&t->U);
484: PetscMalloc2(r,&t->S,r,&t->F);
485: PetscArraycpy(t->A,A,s*s);
486: PetscArraycpy(t->B,B,r*s);
487: PetscArraycpy(t->U,U,s*r);
488: PetscArraycpy(t->V,V,r*r);
489: PetscArraycpy(t->S,S,r);
490: PetscArraycpy(t->F,F,r);
491: if (c) {
492: PetscArraycpy(t->c,c,s);
493: } else {
494: for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j];
495: }
496: PetscMalloc1(r,&t->Fembed);
497: PetscMalloc1(r,&t->Ferror);
498: PetscMalloc1(r,&t->Serror);
499: PetscArraycpy(t->Fembed,Fembed,r);
500: PetscArraycpy(t->Ferror,Ferror,r);
501: PetscArraycpy(t->Serror,Serror,r);
502: t->pinterp = pinterp;
503: PetscMalloc1(s*pinterp,&t->binterp);
504: PetscArraycpy(t->binterp,binterp,s*pinterp);
506: link->next = GLEETableauList;
507: GLEETableauList = link;
508: return(0);
509: }
511: static PetscErrorCode TSEvaluateStep_GLEE(TS ts,PetscInt order,Vec X,PetscBool *done)
512: {
513: TS_GLEE *glee = (TS_GLEE*) ts->data;
514: GLEETableau tab = glee->tableau;
515: PetscReal h, *B = tab->B, *V = tab->V,
516: *F = tab->F,
517: *Fembed = tab->Fembed;
518: PetscInt s = tab->s, r = tab->r, i, j;
519: Vec *Y = glee->Y, *YdotStage = glee->YdotStage;
520: PetscScalar *ws = glee->swork, *wr = glee->rwork;
521: PetscErrorCode ierr;
525: switch (glee->status) {
526: case TS_STEP_INCOMPLETE:
527: case TS_STEP_PENDING:
528: h = ts->time_step; break;
529: case TS_STEP_COMPLETE:
530: h = ts->ptime - ts->ptime_prev; break;
531: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
532: }
534: if (order == tab->order) {
536: /* Note: Irrespective of whether status is TS_STEP_INCOMPLETE
537: or TS_STEP_COMPLETE, glee->X has the solution at the
538: beginning of the time step. So no need to roll-back.
539: */
540: if (glee->status == TS_STEP_INCOMPLETE) {
541: for (i=0; i<r; i++) {
542: VecZeroEntries(Y[i]);
543: for (j=0; j<r; j++) wr[j] = V[i*r+j];
544: VecMAXPY(Y[i],r,wr,glee->X);
545: for (j=0; j<s; j++) ws[j] = h*B[i*s+j];
546: VecMAXPY(Y[i],s,ws,YdotStage);
547: }
548: VecZeroEntries(X);
549: for (j=0; j<r; j++) wr[j] = F[j];
550: VecMAXPY(X,r,wr,Y);
551: } else {VecCopy(ts->vec_sol,X);}
552: return(0);
554: } else if (order == tab->order-1) {
556: /* Complete with the embedded method (Fembed) */
557: for (i=0; i<r; i++) {
558: VecZeroEntries(Y[i]);
559: for (j=0; j<r; j++) wr[j] = V[i*r+j];
560: VecMAXPY(Y[i],r,wr,glee->X);
561: for (j=0; j<s; j++) ws[j] = h*B[i*s+j];
562: VecMAXPY(Y[i],s,ws,YdotStage);
563: }
564: VecZeroEntries(X);
565: for (j=0; j<r; j++) wr[j] = Fembed[j];
566: VecMAXPY(X,r,wr,Y);
568: if (done) *done = PETSC_TRUE;
569: return(0);
570: }
571: if (done) *done = PETSC_FALSE;
572: else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"GLEE '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order);
573: return(0);
574: }
576: static PetscErrorCode TSStep_GLEE(TS ts)
577: {
578: TS_GLEE *glee = (TS_GLEE*)ts->data;
579: GLEETableau tab = glee->tableau;
580: const PetscInt s = tab->s, r = tab->r;
581: PetscReal *A = tab->A, *U = tab->U,
582: *F = tab->F,
583: *c = tab->c;
584: Vec *Y = glee->Y, *X = glee->X,
585: *YStage = glee->YStage,
586: *YdotStage = glee->YdotStage,
587: W = glee->W;
588: SNES snes;
589: PetscScalar *ws = glee->swork, *wr = glee->rwork;
590: TSAdapt adapt;
591: PetscInt i,j,reject,next_scheme,its,lits;
592: PetscReal next_time_step;
593: PetscReal t;
594: PetscBool accept;
595: PetscErrorCode ierr;
598: PetscCitationsRegister(citation,&cited);
600: for (i=0; i<r; i++) { VecCopy(Y[i],X[i]); }
602: TSGetSNES(ts,&snes);
603: next_time_step = ts->time_step;
604: t = ts->ptime;
605: accept = PETSC_TRUE;
606: glee->status = TS_STEP_INCOMPLETE;
608: for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) {
610: PetscReal h = ts->time_step;
611: TSPreStep(ts);
613: for (i=0; i<s; i++) {
615: glee->stage_time = t + h*c[i];
616: TSPreStage(ts,glee->stage_time);
618: if (A[i*s+i] == 0) { /* Explicit stage */
619: VecZeroEntries(YStage[i]);
620: for (j=0; j<r; j++) wr[j] = U[i*r+j];
621: VecMAXPY(YStage[i],r,wr,X);
622: for (j=0; j<i; j++) ws[j] = h*A[i*s+j];
623: VecMAXPY(YStage[i],i,ws,YdotStage);
624: } else { /* Implicit stage */
625: glee->scoeff = 1.0/A[i*s+i];
626: /* compute right-hand-side */
627: VecZeroEntries(W);
628: for (j=0; j<r; j++) wr[j] = U[i*r+j];
629: VecMAXPY(W,r,wr,X);
630: for (j=0; j<i; j++) ws[j] = h*A[i*s+j];
631: VecMAXPY(W,i,ws,YdotStage);
632: VecScale(W,glee->scoeff/h);
633: /* set initial guess */
634: VecCopy(i>0 ? YStage[i-1] : ts->vec_sol,YStage[i]);
635: /* solve for this stage */
636: SNESSolve(snes,W,YStage[i]);
637: SNESGetIterationNumber(snes,&its);
638: SNESGetLinearSolveIterations(snes,&lits);
639: ts->snes_its += its; ts->ksp_its += lits;
640: }
641: TSGetAdapt(ts,&adapt);
642: TSAdaptCheckStage(adapt,ts,glee->stage_time,YStage[i],&accept);
643: if (!accept) goto reject_step;
644: TSPostStage(ts,glee->stage_time,i,YStage);
645: TSComputeRHSFunction(ts,t+h*c[i],YStage[i],YdotStage[i]);
646: }
647: TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);
648: glee->status = TS_STEP_PENDING;
650: /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */
651: TSGetAdapt(ts,&adapt);
652: TSAdaptCandidatesClear(adapt);
653: TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);
654: TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);
655: if (accept) {
656: /* ignore next_scheme for now */
657: ts->ptime += ts->time_step;
658: ts->time_step = next_time_step;
659: glee->status = TS_STEP_COMPLETE;
660: /* compute and store the global error */
661: /* Note: this is not needed if TSAdaptGLEE is not used */
662: TSGetTimeError(ts,0,&(glee->yGErr));
663: PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);
664: break;
665: } else { /* Roll back the current step */
666: for (j=0; j<r; j++) wr[j] = F[j];
667: VecMAXPY(ts->vec_sol,r,wr,X);
668: ts->time_step = next_time_step;
669: glee->status = TS_STEP_INCOMPLETE;
670: }
671: reject_step: continue;
672: }
673: if (glee->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
674: return(0);
675: }
677: static PetscErrorCode TSInterpolate_GLEE(TS ts,PetscReal itime,Vec X)
678: {
679: TS_GLEE *glee = (TS_GLEE*)ts->data;
680: PetscInt s=glee->tableau->s, pinterp=glee->tableau->pinterp,i,j;
681: PetscReal h,tt,t;
682: PetscScalar *b;
683: const PetscReal *B = glee->tableau->binterp;
684: PetscErrorCode ierr;
687: if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSGLEE %s does not have an interpolation formula",glee->tableau->name);
688: switch (glee->status) {
689: case TS_STEP_INCOMPLETE:
690: case TS_STEP_PENDING:
691: h = ts->time_step;
692: t = (itime - ts->ptime)/h;
693: break;
694: case TS_STEP_COMPLETE:
695: h = ts->ptime - ts->ptime_prev;
696: t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
697: break;
698: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
699: }
700: PetscMalloc1(s,&b);
701: for (i=0; i<s; i++) b[i] = 0;
702: for (j=0,tt=t; j<pinterp; j++,tt*=t) {
703: for (i=0; i<s; i++) {
704: b[i] += h * B[i*pinterp+j] * tt;
705: }
706: }
707: VecCopy(glee->YStage[0],X);
708: VecMAXPY(X,s,b,glee->YdotStage);
709: PetscFree(b);
710: return(0);
711: }
713: /*------------------------------------------------------------*/
714: static PetscErrorCode TSReset_GLEE(TS ts)
715: {
716: TS_GLEE *glee = (TS_GLEE*)ts->data;
717: PetscInt s, r;
721: if (!glee->tableau) return(0);
722: s = glee->tableau->s;
723: r = glee->tableau->r;
724: VecDestroyVecs(r,&glee->Y);
725: VecDestroyVecs(r,&glee->X);
726: VecDestroyVecs(s,&glee->YStage);
727: VecDestroyVecs(s,&glee->YdotStage);
728: VecDestroy(&glee->Ydot);
729: VecDestroy(&glee->yGErr);
730: VecDestroy(&glee->W);
731: PetscFree2(glee->swork,glee->rwork);
732: return(0);
733: }
735: static PetscErrorCode TSGLEEGetVecs(TS ts,DM dm,Vec *Ydot)
736: {
737: TS_GLEE *glee = (TS_GLEE*)ts->data;
741: if (Ydot) {
742: if (dm && dm != ts->dm) {
743: DMGetNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);
744: } else *Ydot = glee->Ydot;
745: }
746: return(0);
747: }
750: static PetscErrorCode TSGLEERestoreVecs(TS ts,DM dm,Vec *Ydot)
751: {
755: if (Ydot) {
756: if (dm && dm != ts->dm) {
757: DMRestoreNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);
758: }
759: }
760: return(0);
761: }
763: /*
764: This defines the nonlinear equation that is to be solved with SNES
765: */
766: static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes,Vec X,Vec F,TS ts)
767: {
768: TS_GLEE *glee = (TS_GLEE*)ts->data;
769: DM dm,dmsave;
770: Vec Ydot;
771: PetscReal shift = glee->scoeff / ts->time_step;
775: SNESGetDM(snes,&dm);
776: TSGLEEGetVecs(ts,dm,&Ydot);
777: /* Set Ydot = shift*X */
778: VecCopy(X,Ydot);
779: VecScale(Ydot,shift);
780: dmsave = ts->dm;
781: ts->dm = dm;
783: TSComputeIFunction(ts,glee->stage_time,X,Ydot,F,PETSC_FALSE);
785: ts->dm = dmsave;
786: TSGLEERestoreVecs(ts,dm,&Ydot);
787: return(0);
788: }
790: static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes,Vec X,Mat A,Mat B,TS ts)
791: {
792: TS_GLEE *glee = (TS_GLEE*)ts->data;
793: DM dm,dmsave;
794: Vec Ydot;
795: PetscReal shift = glee->scoeff / ts->time_step;
799: SNESGetDM(snes,&dm);
800: TSGLEEGetVecs(ts,dm,&Ydot);
801: /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */
802: dmsave = ts->dm;
803: ts->dm = dm;
805: TSComputeIJacobian(ts,glee->stage_time,X,Ydot,shift,A,B,PETSC_FALSE);
807: ts->dm = dmsave;
808: TSGLEERestoreVecs(ts,dm,&Ydot);
809: return(0);
810: }
812: static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine,DM coarse,void *ctx)
813: {
815: return(0);
816: }
818: static PetscErrorCode DMRestrictHook_TSGLEE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
819: {
821: return(0);
822: }
825: static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm,DM subdm,void *ctx)
826: {
828: return(0);
829: }
831: static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
832: {
834: return(0);
835: }
837: static PetscErrorCode TSSetUp_GLEE(TS ts)
838: {
839: TS_GLEE *glee = (TS_GLEE*)ts->data;
840: GLEETableau tab;
841: PetscInt s,r;
843: DM dm;
846: if (!glee->tableau) {
847: TSGLEESetType(ts,TSGLEEDefaultType);
848: }
849: tab = glee->tableau;
850: s = tab->s;
851: r = tab->r;
852: VecDuplicateVecs(ts->vec_sol,r,&glee->Y);
853: VecDuplicateVecs(ts->vec_sol,r,&glee->X);
854: VecDuplicateVecs(ts->vec_sol,s,&glee->YStage);
855: VecDuplicateVecs(ts->vec_sol,s,&glee->YdotStage);
856: VecDuplicate(ts->vec_sol,&glee->Ydot);
857: VecDuplicate(ts->vec_sol,&glee->yGErr);
858: VecZeroEntries(glee->yGErr);
859: VecDuplicate(ts->vec_sol,&glee->W);
860: PetscMalloc2(s,&glee->swork,r,&glee->rwork);
861: TSGetDM(ts,&dm);
862: DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);
863: DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);
864: return(0);
865: }
867: PetscErrorCode TSStartingMethod_GLEE(TS ts)
868: {
869: TS_GLEE *glee = (TS_GLEE*)ts->data;
870: GLEETableau tab = glee->tableau;
871: PetscInt r=tab->r,i;
872: PetscReal *S=tab->S;
876: for (i=0; i<r; i++) {
877: VecZeroEntries(glee->Y[i]);
878: VecAXPY(glee->Y[i],S[i],ts->vec_sol);
879: }
881: return(0);
882: }
885: /*------------------------------------------------------------*/
887: static PetscErrorCode TSSetFromOptions_GLEE(PetscOptionItems *PetscOptionsObject,TS ts)
888: {
890: char gleetype[256];
893: PetscOptionsHead(PetscOptionsObject,"GLEE ODE solver options");
894: {
895: GLEETableauLink link;
896: PetscInt count,choice;
897: PetscBool flg;
898: const char **namelist;
900: PetscStrncpy(gleetype,TSGLEEDefaultType,sizeof(gleetype));
901: for (link=GLEETableauList,count=0; link; link=link->next,count++) ;
902: PetscMalloc1(count,(char***)&namelist);
903: for (link=GLEETableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
904: PetscOptionsEList("-ts_glee_type","Family of GLEE method","TSGLEESetType",(const char*const*)namelist,count,gleetype,&choice,&flg);
905: TSGLEESetType(ts,flg ? namelist[choice] : gleetype);
906: PetscFree(namelist);
907: }
908: PetscOptionsTail();
909: return(0);
910: }
912: static PetscErrorCode TSView_GLEE(TS ts,PetscViewer viewer)
913: {
914: TS_GLEE *glee = (TS_GLEE*)ts->data;
915: GLEETableau tab = glee->tableau;
916: PetscBool iascii;
920: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
921: if (iascii) {
922: TSGLEEType gleetype;
923: char buf[512];
924: TSGLEEGetType(ts,&gleetype);
925: PetscViewerASCIIPrintf(viewer," GLEE type %s\n",gleetype);
926: PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);
927: PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);
928: /* Note: print out r as well */
929: }
930: return(0);
931: }
933: static PetscErrorCode TSLoad_GLEE(TS ts,PetscViewer viewer)
934: {
936: SNES snes;
937: TSAdapt tsadapt;
940: TSGetAdapt(ts,&tsadapt);
941: TSAdaptLoad(tsadapt,viewer);
942: TSGetSNES(ts,&snes);
943: SNESLoad(snes,viewer);
944: /* function and Jacobian context for SNES when used with TS is always ts object */
945: SNESSetFunction(snes,NULL,NULL,ts);
946: SNESSetJacobian(snes,NULL,NULL,NULL,ts);
947: return(0);
948: }
950: /*@C
951: TSGLEESetType - Set the type of GLEE scheme
953: Logically collective
955: Input Parameter:
956: + ts - timestepping context
957: - gleetype - type of GLEE-scheme
959: Level: intermediate
961: .seealso: TSGLEEGetType(), TSGLEE
962: @*/
963: PetscErrorCode TSGLEESetType(TS ts,TSGLEEType gleetype)
964: {
970: PetscTryMethod(ts,"TSGLEESetType_C",(TS,TSGLEEType),(ts,gleetype));
971: return(0);
972: }
974: /*@C
975: TSGLEEGetType - Get the type of GLEE scheme
977: Logically collective
979: Input Parameter:
980: . ts - timestepping context
982: Output Parameter:
983: . gleetype - type of GLEE-scheme
985: Level: intermediate
987: .seealso: TSGLEESetType()
988: @*/
989: PetscErrorCode TSGLEEGetType(TS ts,TSGLEEType *gleetype)
990: {
995: PetscUseMethod(ts,"TSGLEEGetType_C",(TS,TSGLEEType*),(ts,gleetype));
996: return(0);
997: }
999: PetscErrorCode TSGLEEGetType_GLEE(TS ts,TSGLEEType *gleetype)
1000: {
1001: TS_GLEE *glee = (TS_GLEE*)ts->data;
1005: if (!glee->tableau) {
1006: TSGLEESetType(ts,TSGLEEDefaultType);
1007: }
1008: *gleetype = glee->tableau->name;
1009: return(0);
1010: }
1011: PetscErrorCode TSGLEESetType_GLEE(TS ts,TSGLEEType gleetype)
1012: {
1013: TS_GLEE *glee = (TS_GLEE*)ts->data;
1014: PetscErrorCode ierr;
1015: PetscBool match;
1016: GLEETableauLink link;
1019: if (glee->tableau) {
1020: PetscStrcmp(glee->tableau->name,gleetype,&match);
1021: if (match) return(0);
1022: }
1023: for (link = GLEETableauList; link; link=link->next) {
1024: PetscStrcmp(link->tab.name,gleetype,&match);
1025: if (match) {
1026: TSReset_GLEE(ts);
1027: glee->tableau = &link->tab;
1028: return(0);
1029: }
1030: }
1031: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",gleetype);
1032: }
1034: static PetscErrorCode TSGetStages_GLEE(TS ts,PetscInt *ns,Vec **Y)
1035: {
1036: TS_GLEE *glee = (TS_GLEE*)ts->data;
1039: if (ns) *ns = glee->tableau->s;
1040: if (Y) *Y = glee->YStage;
1041: return(0);
1042: }
1044: PetscErrorCode TSGetSolutionComponents_GLEE(TS ts,PetscInt *n,Vec *Y)
1045: {
1046: TS_GLEE *glee = (TS_GLEE*)ts->data;
1047: GLEETableau tab = glee->tableau;
1048: PetscErrorCode ierr;
1051: if (!Y) *n = tab->r;
1052: else {
1053: if ((*n >= 0) && (*n < tab->r)) {
1054: VecCopy(glee->Y[*n],*Y);
1055: } else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Second argument (%d) out of range[%d,%d].",*n,0,tab->r-1);
1056: }
1057: return(0);
1058: }
1060: PetscErrorCode TSGetAuxSolution_GLEE(TS ts,Vec *X)
1061: {
1062: TS_GLEE *glee = (TS_GLEE*)ts->data;
1063: GLEETableau tab = glee->tableau;
1064: PetscReal *F = tab->Fembed;
1065: PetscInt r = tab->r;
1066: Vec *Y = glee->Y;
1067: PetscScalar *wr = glee->rwork;
1068: PetscInt i;
1069: PetscErrorCode ierr;
1072: VecZeroEntries(*X);
1073: for (i=0; i<r; i++) wr[i] = F[i];
1074: VecMAXPY((*X),r,wr,Y);
1075: return(0);
1076: }
1078: PetscErrorCode TSGetTimeError_GLEE(TS ts,PetscInt n,Vec *X)
1079: {
1080: TS_GLEE *glee = (TS_GLEE*)ts->data;
1081: GLEETableau tab = glee->tableau;
1082: PetscReal *F = tab->Ferror;
1083: PetscInt r = tab->r;
1084: Vec *Y = glee->Y;
1085: PetscScalar *wr = glee->rwork;
1086: PetscInt i;
1087: PetscErrorCode ierr;
1090: VecZeroEntries(*X);
1091: if (n==0){
1092: for (i=0; i<r; i++) wr[i] = F[i];
1093: VecMAXPY((*X),r,wr,Y);
1094: } else if (n==-1) {
1095: *X=glee->yGErr;
1096: }
1097: return(0);
1098: }
1100: PetscErrorCode TSSetTimeError_GLEE(TS ts,Vec X)
1101: {
1102: TS_GLEE *glee = (TS_GLEE*)ts->data;
1103: GLEETableau tab = glee->tableau;
1104: PetscReal *S = tab->Serror;
1105: PetscInt r = tab->r,i;
1106: Vec *Y = glee->Y;
1107: PetscErrorCode ierr;
1110: if (r != 2) SETERRQ2(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSSetTimeError_GLEE not supported for '%s' with r=%D.",tab->name,tab->r);
1111: for (i=1; i<r; i++) {
1112: VecCopy(ts->vec_sol,Y[i]);
1113: VecAXPBY(Y[i],S[0],S[1],X);
1114: VecCopy(X,glee->yGErr);
1115: }
1116: return(0);
1117: }
1119: static PetscErrorCode TSDestroy_GLEE(TS ts)
1120: {
1124: TSReset_GLEE(ts);
1125: if (ts->dm) {
1126: DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);
1127: DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);
1128: }
1129: PetscFree(ts->data);
1130: PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",NULL);
1131: PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",NULL);
1132: return(0);
1133: }
1135: /* ------------------------------------------------------------ */
1136: /*MC
1137: TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes
1139: The user should provide the right hand side of the equation
1140: using TSSetRHSFunction().
1142: Notes:
1143: The default is TSGLEE35, it can be changed with TSGLEESetType() or -ts_glee_type
1145: Level: beginner
1147: .seealso: TSCreate(), TS, TSSetType(), TSGLEESetType(), TSGLEEGetType(),
1148: TSGLEE23, TTSGLEE24, TSGLEE35, TSGLEE25I, TSGLEEEXRK2A,
1149: TSGLEERK32G1, TSGLEERK285EX, TSGLEEType, TSGLEERegister()
1151: M*/
1152: PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts)
1153: {
1154: TS_GLEE *th;
1155: PetscErrorCode ierr;
1158: TSGLEEInitializePackage();
1160: ts->ops->reset = TSReset_GLEE;
1161: ts->ops->destroy = TSDestroy_GLEE;
1162: ts->ops->view = TSView_GLEE;
1163: ts->ops->load = TSLoad_GLEE;
1164: ts->ops->setup = TSSetUp_GLEE;
1165: ts->ops->step = TSStep_GLEE;
1166: ts->ops->interpolate = TSInterpolate_GLEE;
1167: ts->ops->evaluatestep = TSEvaluateStep_GLEE;
1168: ts->ops->setfromoptions = TSSetFromOptions_GLEE;
1169: ts->ops->getstages = TSGetStages_GLEE;
1170: ts->ops->snesfunction = SNESTSFormFunction_GLEE;
1171: ts->ops->snesjacobian = SNESTSFormJacobian_GLEE;
1172: ts->ops->getsolutioncomponents = TSGetSolutionComponents_GLEE;
1173: ts->ops->getauxsolution = TSGetAuxSolution_GLEE;
1174: ts->ops->gettimeerror = TSGetTimeError_GLEE;
1175: ts->ops->settimeerror = TSSetTimeError_GLEE;
1176: ts->ops->startingmethod = TSStartingMethod_GLEE;
1177: ts->default_adapt_type = TSADAPTGLEE;
1179: ts->usessnes = PETSC_TRUE;
1181: PetscNewLog(ts,&th);
1182: ts->data = (void*)th;
1184: PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",TSGLEEGetType_GLEE);
1185: PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",TSGLEESetType_GLEE);
1186: return(0);
1187: }