Actual source code: rk.c
petsc-3.8.4 2018-03-24
1: /*
2: Code for time stepping with the Runge-Kutta method
4: Notes:
5: The general system is written as
7: Udot = F(t,U)
9: */
10: #include <petsc/private/tsimpl.h>
11: #include <petscdm.h>
13: static TSRKType TSRKDefault = TSRK3BS;
14: static PetscBool TSRKRegisterAllCalled;
15: static PetscBool TSRKPackageInitialized;
17: typedef struct _RKTableau *RKTableau;
18: struct _RKTableau {
19: char *name;
20: PetscInt order; /* Classical approximation order of the method i */
21: PetscInt s; /* Number of stages */
22: PetscInt p; /* Interpolation order */
23: PetscBool FSAL; /* flag to indicate if tableau is FSAL */
24: PetscReal *A,*b,*c; /* Tableau */
25: PetscReal *bembed; /* Embedded formula of order one less (order-1) */
26: PetscReal *binterp; /* Dense output formula */
27: PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */
28: };
29: typedef struct _RKTableauLink *RKTableauLink;
30: struct _RKTableauLink {
31: struct _RKTableau tab;
32: RKTableauLink next;
33: };
34: static RKTableauLink RKTableauList;
36: typedef struct {
37: RKTableau tableau;
38: Vec *Y; /* States computed during the step */
39: Vec *YdotRHS; /* Function evaluations for the non-stiff part */
40: Vec *VecDeltaLam; /* Increment of the adjoint sensitivity w.r.t IC at stage */
41: Vec *VecDeltaMu; /* Increment of the adjoint sensitivity w.r.t P at stage */
42: Vec *VecSensiTemp; /* Vector to be timed with Jacobian transpose */
43: Vec VecCostIntegral0; /* backup for roll-backs due to events */
44: PetscScalar *work; /* Scalar work */
45: PetscReal stage_time;
46: TSStepStatus status;
47: PetscReal ptime;
48: PetscReal time_step;
49: } TS_RK;
51: /*MC
52: TSRK1FE - First order forward Euler scheme.
54: This method has one stage.
56: Options database:
57: . -ts_rk_type 1fe
59: Level: advanced
61: .seealso: TSRK, TSRKType, TSRKSetType()
62: M*/
63: /*MC
64: TSRK2A - Second order RK scheme.
66: This method has two stages.
68: Options database:
69: . -ts_rk_type 2a
71: Level: advanced
73: .seealso: TSRK, TSRKType, TSRKSetType()
74: M*/
75: /*MC
76: TSRK3 - Third order RK scheme.
78: This method has three stages.
80: Options database:
81: . -ts_rk_type 3
83: Level: advanced
85: .seealso: TSRK, TSRKType, TSRKSetType()
86: M*/
87: /*MC
88: TSRK3BS - Third order RK scheme of Bogacki-Shampine with 2nd order embedded method.
90: This method has four stages with the First Same As Last (FSAL) property.
92: Options database:
93: . -ts_rk_type 3bs
95: Level: advanced
97: References: https://doi.org/10.1016/0893-9659(89)90079-7
99: .seealso: TSRK, TSRKType, TSRKSetType()
100: M*/
101: /*MC
102: TSRK4 - Fourth order RK scheme.
104: This is the classical Runge-Kutta method with four stages.
106: Options database:
107: . -ts_rk_type 4
109: Level: advanced
111: .seealso: TSRK, TSRKType, TSRKSetType()
112: M*/
113: /*MC
114: TSRK5F - Fifth order Fehlberg RK scheme with a 4th order embedded method.
116: This method has six stages.
118: Options database:
119: . -ts_rk_type 5f
121: Level: advanced
123: .seealso: TSRK, TSRKType, TSRKSetType()
124: M*/
125: /*MC
126: TSRK5DP - Fifth order Dormand-Prince RK scheme with the 4th order embedded method.
128: This method has seven stages with the First Same As Last (FSAL) property.
130: Options database:
131: . -ts_rk_type 5dp
133: Level: advanced
135: References: https://doi.org/10.1016/0771-050X(80)90013-3
137: .seealso: TSRK, TSRKType, TSRKSetType()
138: M*/
139: /*MC
140: TSRK5BS - Fifth order Bogacki-Shampine RK scheme with 4th order embedded method.
142: This method has eight stages with the First Same As Last (FSAL) property.
144: Options database:
145: . -ts_rk_type 5bs
147: Level: advanced
149: References: https://doi.org/10.1016/0898-1221(96)00141-1
151: .seealso: TSRK, TSRKType, TSRKSetType()
152: M*/
154: /*@C
155: TSRKRegisterAll - Registers all of the Runge-Kutta explicit methods in TSRK
157: Not Collective, but should be called by all processes which will need the schemes to be registered
159: Level: advanced
161: .keywords: TS, TSRK, register, all
163: .seealso: TSRKRegisterDestroy()
164: @*/
165: PetscErrorCode TSRKRegisterAll(void)
166: {
170: if (TSRKRegisterAllCalled) return(0);
171: TSRKRegisterAllCalled = PETSC_TRUE;
173: #define RC PetscRealConstant
174: {
175: const PetscReal
176: A[1][1] = {{0}},
177: b[1] = {RC(1.0)};
178: TSRKRegister(TSRK1FE,1,1,&A[0][0],b,NULL,NULL,0,NULL);
179: }
180: {
181: const PetscReal
182: A[2][2] = {{0,0},
183: {RC(1.0),0}},
184: b[2] = {RC(0.5),RC(0.5)},
185: bembed[2] = {RC(1.0),0};
186: TSRKRegister(TSRK2A,2,2,&A[0][0],b,NULL,bembed,0,NULL);
187: }
188: {
189: const PetscReal
190: A[3][3] = {{0,0,0},
191: {RC(2.0)/RC(3.0),0,0},
192: {RC(-1.0)/RC(3.0),RC(1.0),0}},
193: b[3] = {RC(0.25),RC(0.5),RC(0.25)};
194: TSRKRegister(TSRK3,3,3,&A[0][0],b,NULL,NULL,0,NULL);
195: }
196: {
197: const PetscReal
198: A[4][4] = {{0,0,0,0},
199: {RC(1.0)/RC(2.0),0,0,0},
200: {0,RC(3.0)/RC(4.0),0,0},
201: {RC(2.0)/RC(9.0),RC(1.0)/RC(3.0),RC(4.0)/RC(9.0),0}},
202: b[4] = {RC(2.0)/RC(9.0),RC(1.0)/RC(3.0),RC(4.0)/RC(9.0),0},
203: bembed[4] = {RC(7.0)/RC(24.0),RC(1.0)/RC(4.0),RC(1.0)/RC(3.0),RC(1.0)/RC(8.0)};
204: TSRKRegister(TSRK3BS,3,4,&A[0][0],b,NULL,bembed,0,NULL);
205: }
206: {
207: const PetscReal
208: A[4][4] = {{0,0,0,0},
209: {RC(0.5),0,0,0},
210: {0,RC(0.5),0,0},
211: {0,0,RC(1.0),0}},
212: b[4] = {RC(1.0)/RC(6.0),RC(1.0)/RC(3.0),RC(1.0)/RC(3.0),RC(1.0)/RC(6.0)};
213: TSRKRegister(TSRK4,4,4,&A[0][0],b,NULL,NULL,0,NULL);
214: }
215: {
216: const PetscReal
217: A[6][6] = {{0,0,0,0,0,0},
218: {RC(0.25),0,0,0,0,0},
219: {RC(3.0)/RC(32.0),RC(9.0)/RC(32.0),0,0,0,0},
220: {RC(1932.0)/RC(2197.0),RC(-7200.0)/RC(2197.0),RC(7296.0)/RC(2197.0),0,0,0},
221: {RC(439.0)/RC(216.0),RC(-8.0),RC(3680.0)/RC(513.0),RC(-845.0)/RC(4104.0),0,0},
222: {RC(-8.0)/RC(27.0),RC(2.0),RC(-3544.0)/RC(2565.0),RC(1859.0)/RC(4104.0),RC(-11.0)/RC(40.0),0}},
223: b[6] = {RC(16.0)/RC(135.0),0,RC(6656.0)/RC(12825.0),RC(28561.0)/RC(56430.0),RC(-9.0)/RC(50.0),RC(2.0)/RC(55.0)},
224: bembed[6] = {RC(25.0)/RC(216.0),0,RC(1408.0)/RC(2565.0),RC(2197.0)/RC(4104.0),RC(-1.0)/RC(5.0),0};
225: TSRKRegister(TSRK5F,5,6,&A[0][0],b,NULL,bembed,0,NULL);
226: }
227: {
228: const PetscReal
229: A[7][7] = {{0,0,0,0,0,0,0},
230: {RC(1.0)/RC(5.0),0,0,0,0,0,0},
231: {RC(3.0)/RC(40.0),RC(9.0)/RC(40.0),0,0,0,0,0},
232: {RC(44.0)/RC(45.0),RC(-56.0)/RC(15.0),RC(32.0)/RC(9.0),0,0,0,0},
233: {RC(19372.0)/RC(6561.0),RC(-25360.0)/RC(2187.0),RC(64448.0)/RC(6561.0),RC(-212.0)/RC(729.0),0,0,0},
234: {RC(9017.0)/RC(3168.0),RC(-355.0)/RC(33.0),RC(46732.0)/RC(5247.0),RC(49.0)/RC(176.0),RC(-5103.0)/RC(18656.0),0,0},
235: {RC(35.0)/RC(384.0),0,RC(500.0)/RC(1113.0),RC(125.0)/RC(192.0),RC(-2187.0)/RC(6784.0),RC(11.0)/RC(84.0),0}},
236: b[7] = {RC(35.0)/RC(384.0),0,RC(500.0)/RC(1113.0),RC(125.0)/RC(192.0),RC(-2187.0)/RC(6784.0),RC(11.0)/RC(84.0),0},
237: bembed[7] = {RC(5179.0)/RC(57600.0),0,RC(7571.0)/RC(16695.0),RC(393.0)/RC(640.0),RC(-92097.0)/RC(339200.0),RC(187.0)/RC(2100.0),RC(1.0)/RC(40.0)};
238: TSRKRegister(TSRK5DP,5,7,&A[0][0],b,NULL,bembed,0,NULL);
239: }
240: {
241: const PetscReal
242: A[8][8] = {{0,0,0,0,0,0,0,0},
243: {RC(1.0)/RC(6.0),0,0,0,0,0,0,0},
244: {RC(2.0)/RC(27.0),RC(4.0)/RC(27.0),0,0,0,0,0,0},
245: {RC(183.0)/RC(1372.0),RC(-162.0)/RC(343.0),RC(1053.0)/RC(1372.0),0,0,0,0,0},
246: {RC(68.0)/RC(297.0),RC(-4.0)/RC(11.0),RC(42.0)/RC(143.0),RC(1960.0)/RC(3861.0),0,0,0,0},
247: {RC(597.0)/RC(22528.0),RC(81.0)/RC(352.0),RC(63099.0)/RC(585728.0),RC(58653.0)/RC(366080.0),RC(4617.0)/RC(20480.0),0,0,0},
248: {RC(174197.0)/RC(959244.0),RC(-30942.0)/RC(79937.0),RC(8152137.0)/RC(19744439.0),RC(666106.0)/RC(1039181.0),RC(-29421.0)/RC(29068.0),RC(482048.0)/RC(414219.0),0,0},
249: {RC(587.0)/RC(8064.0),0,RC(4440339.0)/RC(15491840.0),RC(24353.0)/RC(124800.0),RC(387.0)/RC(44800.0),RC(2152.0)/RC(5985.0),RC(7267.0)/RC(94080.0),0}},
250: b[8] = {RC(587.0)/RC(8064.0),0,RC(4440339.0)/RC(15491840.0),RC(24353.0)/RC(124800.0),RC(387.0)/RC(44800.0),RC(2152.0)/RC(5985.0),RC(7267.0)/RC(94080.0),0},
251: bembed[8] = {RC(2479.0)/RC(34992.0),0,RC(123.0)/RC(416.0),RC(612941.0)/RC(3411720.0),RC(43.0)/RC(1440.0),RC(2272.0)/RC(6561.0),RC(79937.0)/RC(1113912.0),RC(3293.0)/RC(556956.0)};
252: TSRKRegister(TSRK5BS,5,8,&A[0][0],b,NULL,bembed,0,NULL);
253: }
254: #undef RC
255: return(0);
256: }
258: /*@C
259: TSRKRegisterDestroy - Frees the list of schemes that were registered by TSRKRegister().
261: Not Collective
263: Level: advanced
265: .keywords: TSRK, register, destroy
266: .seealso: TSRKRegister(), TSRKRegisterAll()
267: @*/
268: PetscErrorCode TSRKRegisterDestroy(void)
269: {
271: RKTableauLink link;
274: while ((link = RKTableauList)) {
275: RKTableau t = &link->tab;
276: RKTableauList = link->next;
277: PetscFree3(t->A,t->b,t->c);
278: PetscFree (t->bembed);
279: PetscFree (t->binterp);
280: PetscFree (t->name);
281: PetscFree (link);
282: }
283: TSRKRegisterAllCalled = PETSC_FALSE;
284: return(0);
285: }
287: /*@C
288: TSRKInitializePackage - This function initializes everything in the TSRK package. It is called
289: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RK()
290: when using static libraries.
292: Level: developer
294: .keywords: TS, TSRK, initialize, package
295: .seealso: PetscInitialize()
296: @*/
297: PetscErrorCode TSRKInitializePackage(void)
298: {
302: if (TSRKPackageInitialized) return(0);
303: TSRKPackageInitialized = PETSC_TRUE;
304: TSRKRegisterAll();
305: PetscRegisterFinalize(TSRKFinalizePackage);
306: return(0);
307: }
309: /*@C
310: TSRKFinalizePackage - This function destroys everything in the TSRK package. It is
311: called from PetscFinalize().
313: Level: developer
315: .keywords: Petsc, destroy, package
316: .seealso: PetscFinalize()
317: @*/
318: PetscErrorCode TSRKFinalizePackage(void)
319: {
323: TSRKPackageInitialized = PETSC_FALSE;
324: TSRKRegisterDestroy();
325: return(0);
326: }
328: /*@C
329: TSRKRegister - register an RK scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation
331: Not Collective, but the same schemes should be registered on all processes on which they will be used
333: Input Parameters:
334: + name - identifier for method
335: . order - approximation order of method
336: . s - number of stages, this is the dimension of the matrices below
337: . A - stage coefficients (dimension s*s, row-major)
338: . b - step completion table (dimension s; NULL to use last row of A)
339: . c - abscissa (dimension s; NULL to use row sums of A)
340: . bembed - completion table for embedded method (dimension s; NULL if not available)
341: . p - Order of the interpolation scheme, equal to the number of columns of binterp
342: - binterp - Coefficients of the interpolation formula (dimension s*p; NULL to reuse b with p=1)
344: Notes:
345: Several RK methods are provided, this function is only needed to create new methods.
347: Level: advanced
349: .keywords: TS, register
351: .seealso: TSRK
352: @*/
353: PetscErrorCode TSRKRegister(TSRKType name,PetscInt order,PetscInt s,
354: const PetscReal A[],const PetscReal b[],const PetscReal c[],
355: const PetscReal bembed[],PetscInt p,const PetscReal binterp[])
356: {
357: PetscErrorCode ierr;
358: RKTableauLink link;
359: RKTableau t;
360: PetscInt i,j;
370: PetscNew(&link);
371: t = &link->tab;
373: PetscStrallocpy(name,&t->name);
374: t->order = order;
375: t->s = s;
376: PetscMalloc3(s*s,&t->A,s,&t->b,s,&t->c);
377: PetscMemcpy(t->A,A,s*s*sizeof(A[0]));
378: if (b) { PetscMemcpy(t->b,b,s*sizeof(b[0])); }
379: else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i];
380: if (c) { PetscMemcpy(t->c,c,s*sizeof(c[0])); }
381: else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j];
382: t->FSAL = PETSC_TRUE;
383: for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE;
385: if (bembed) {
386: PetscMalloc1(s,&t->bembed);
387: PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));
388: }
390: if (!binterp) { p = 1; binterp = t->b; }
391: t->p = p;
392: PetscMalloc1(s*p,&t->binterp);
393: PetscMemcpy(t->binterp,binterp,s*p*sizeof(binterp[0]));
395: link->next = RKTableauList;
396: RKTableauList = link;
397: return(0);
398: }
400: /*
401: The step completion formula is
403: x1 = x0 + h b^T YdotRHS
405: This function can be called before or after ts->vec_sol has been updated.
406: Suppose we have a completion formula (b) and an embedded formula (be) of different order.
407: We can write
409: x1e = x0 + h be^T YdotRHS
410: = x1 - h b^T YdotRHS + h be^T YdotRHS
411: = x1 + h (be - b)^T YdotRHS
413: so we can evaluate the method with different order even after the step has been optimistically completed.
414: */
415: static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done)
416: {
417: TS_RK *rk = (TS_RK*)ts->data;
418: RKTableau tab = rk->tableau;
419: PetscScalar *w = rk->work;
420: PetscReal h;
421: PetscInt s = tab->s,j;
425: switch (rk->status) {
426: case TS_STEP_INCOMPLETE:
427: case TS_STEP_PENDING:
428: h = ts->time_step; break;
429: case TS_STEP_COMPLETE:
430: h = ts->ptime - ts->ptime_prev; break;
431: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
432: }
433: if (order == tab->order) {
434: if (rk->status == TS_STEP_INCOMPLETE) {
435: VecCopy(ts->vec_sol,X);
436: for (j=0; j<s; j++) w[j] = h*tab->b[j];
437: VecMAXPY(X,s,w,rk->YdotRHS);
438: } else {VecCopy(ts->vec_sol,X);}
439: return(0);
440: } else if (order == tab->order-1) {
441: if (!tab->bembed) goto unavailable;
442: if (rk->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (be) */
443: VecCopy(ts->vec_sol,X);
444: for (j=0; j<s; j++) w[j] = h*tab->bembed[j];
445: VecMAXPY(X,s,w,rk->YdotRHS);
446: } else { /* Rollback and re-complete using (be-b) */
447: VecCopy(ts->vec_sol,X);
448: for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]);
449: VecMAXPY(X,s,w,rk->YdotRHS);
450: if (ts->vec_costintegral && ts->costintegralfwd) {
451: VecCopy(rk->VecCostIntegral0,ts->vec_costintegral);
452: }
453: }
454: if (done) *done = PETSC_TRUE;
455: return(0);
456: }
457: unavailable:
458: if (done) *done = PETSC_FALSE;
459: else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"RK '%s' of order %D cannot evaluate step at order %D. Consider using -ts_adapt_type none or a different method that has an embedded estimate.",tab->name,tab->order,order);
460: return(0);
461: }
463: static PetscErrorCode TSForwardCostIntegral_RK(TS ts)
464: {
465: TS_RK *rk = (TS_RK*)ts->data;
466: RKTableau tab = rk->tableau;
467: const PetscInt s = tab->s;
468: const PetscReal *b = tab->b,*c = tab->c;
469: Vec *Y = rk->Y;
470: PetscInt i;
471: PetscErrorCode ierr;
474: /* backup cost integral */
475: VecCopy(ts->vec_costintegral,rk->VecCostIntegral0);
476: for (i=s-1; i>=0; i--) {
477: /* Evolve ts->vec_costintegral to compute integrals */
478: TSComputeCostIntegrand(ts,rk->ptime+rk->time_step*(1.0-c[i]),Y[i],ts->vec_costintegrand);
479: VecAXPY(ts->vec_costintegral,rk->time_step*b[i],ts->vec_costintegrand);
480: }
481: return(0);
482: }
484: static PetscErrorCode TSAdjointCostIntegral_RK(TS ts)
485: {
486: TS_RK *rk = (TS_RK*)ts->data;
487: RKTableau tab = rk->tableau;
488: const PetscInt s = tab->s;
489: const PetscReal *b = tab->b,*c = tab->c;
490: Vec *Y = rk->Y;
491: PetscInt i;
492: PetscErrorCode ierr;
495: for (i=s-1; i>=0; i--) {
496: /* Evolve ts->vec_costintegral to compute integrals */
497: TSComputeCostIntegrand(ts,ts->ptime-ts->time_step*(1.0-c[i]),Y[i],ts->vec_costintegrand);
498: VecAXPY(ts->vec_costintegral,-ts->time_step*b[i],ts->vec_costintegrand);
499: }
500: return(0);
501: }
503: static PetscErrorCode TSRollBack_RK(TS ts)
504: {
505: TS_RK *rk = (TS_RK*)ts->data;
506: RKTableau tab = rk->tableau;
507: const PetscInt s = tab->s;
508: const PetscReal *b = tab->b;
509: PetscScalar *w = rk->work;
510: Vec *YdotRHS = rk->YdotRHS;
511: PetscInt j;
512: PetscReal h;
513: PetscErrorCode ierr;
516: switch (rk->status) {
517: case TS_STEP_INCOMPLETE:
518: case TS_STEP_PENDING:
519: h = ts->time_step; break;
520: case TS_STEP_COMPLETE:
521: h = ts->ptime - ts->ptime_prev; break;
522: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
523: }
524: for (j=0; j<s; j++) w[j] = -h*b[j];
525: VecMAXPY(ts->vec_sol,s,w,YdotRHS);
526: return(0);
527: }
529: static PetscErrorCode TSStep_RK(TS ts)
530: {
531: TS_RK *rk = (TS_RK*)ts->data;
532: RKTableau tab = rk->tableau;
533: const PetscInt s = tab->s;
534: const PetscReal *A = tab->A,*c = tab->c;
535: PetscScalar *w = rk->work;
536: Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS;
537: PetscBool FSAL = tab->FSAL;
538: TSAdapt adapt;
539: PetscInt i,j;
540: PetscInt rejections = 0;
541: PetscBool stageok,accept = PETSC_TRUE;
542: PetscReal next_time_step = ts->time_step;
543: PetscErrorCode ierr;
546: if (ts->steprollback || ts->steprestart) FSAL = PETSC_FALSE;
547: if (FSAL) { VecCopy(YdotRHS[s-1],YdotRHS[0]); }
549: rk->status = TS_STEP_INCOMPLETE;
550: while (!ts->reason && rk->status != TS_STEP_COMPLETE) {
551: PetscReal t = ts->ptime;
552: PetscReal h = ts->time_step;
553: for (i=0; i<s; i++) {
554: rk->stage_time = t + h*c[i];
555: TSPreStage(ts,rk->stage_time);
556: VecCopy(ts->vec_sol,Y[i]);
557: for (j=0; j<i; j++) w[j] = h*A[i*s+j];
558: VecMAXPY(Y[i],i,w,YdotRHS);
559: TSPostStage(ts,rk->stage_time,i,Y);
560: TSGetAdapt(ts,&adapt);
561: TSAdaptCheckStage(adapt,ts,rk->stage_time,Y[i],&stageok);
562: if (!stageok) goto reject_step;
563: if (FSAL && !i) continue;
564: TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);
565: }
567: rk->status = TS_STEP_INCOMPLETE;
568: TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);
569: rk->status = TS_STEP_PENDING;
570: TSGetAdapt(ts,&adapt);
571: TSAdaptCandidatesClear(adapt);
572: TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);
573: TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
574: rk->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
575: if (!accept) { /* Roll back the current step */
576: TSRollBack_RK(ts);
577: ts->time_step = next_time_step;
578: goto reject_step;
579: }
581: if (ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation*/
582: rk->ptime = ts->ptime;
583: rk->time_step = ts->time_step;
584: }
586: ts->ptime += ts->time_step;
587: ts->time_step = next_time_step;
588: break;
590: reject_step:
591: ts->reject++; accept = PETSC_FALSE;
592: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
593: ts->reason = TS_DIVERGED_STEP_REJECTED;
594: PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
595: }
596: }
597: return(0);
598: }
600: static PetscErrorCode TSAdjointSetUp_RK(TS ts)
601: {
602: TS_RK *rk = (TS_RK*)ts->data;
603: RKTableau tab = rk->tableau;
604: PetscInt s = tab->s;
608: if (ts->adjointsetupcalled++) return(0);
609: VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecDeltaLam);
610: if(ts->vecs_sensip) {
611: VecDuplicateVecs(ts->vecs_sensip[0],s*ts->numcost,&rk->VecDeltaMu);
612: }
613: VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&rk->VecSensiTemp);
614: return(0);
615: }
617: static PetscErrorCode TSAdjointStep_RK(TS ts)
618: {
619: TS_RK *rk = (TS_RK*)ts->data;
620: RKTableau tab = rk->tableau;
621: const PetscInt s = tab->s;
622: const PetscReal *A = tab->A,*b = tab->b,*c = tab->c;
623: PetscScalar *w = rk->work;
624: Vec *Y = rk->Y,*VecDeltaLam = rk->VecDeltaLam,*VecDeltaMu = rk->VecDeltaMu,*VecSensiTemp = rk->VecSensiTemp;
625: PetscInt i,j,nadj;
626: PetscReal t = ts->ptime;
627: PetscErrorCode ierr;
628: PetscReal h = ts->time_step;
629: Mat J,Jp;
632: rk->status = TS_STEP_INCOMPLETE;
633: for (i=s-1; i>=0; i--) {
634: rk->stage_time = t + h*(1.0-c[i]);
635: for (nadj=0; nadj<ts->numcost; nadj++) {
636: VecCopy(ts->vecs_sensi[nadj],VecSensiTemp[nadj]);
637: VecScale(VecSensiTemp[nadj],-h*b[i]);
638: for (j=i+1; j<s; j++) {
639: VecAXPY(VecSensiTemp[nadj],-h*A[j*s+i],VecDeltaLam[nadj*s+j]);
640: }
641: }
642: /* Stage values of lambda */
643: TSGetRHSJacobian(ts,&J,&Jp,NULL,NULL);
644: TSComputeRHSJacobian(ts,rk->stage_time,Y[i],J,Jp);
645: if (ts->vec_costintegral) {
646: TSAdjointComputeDRDYFunction(ts,rk->stage_time,Y[i],ts->vecs_drdy);
647: }
648: for (nadj=0; nadj<ts->numcost; nadj++) {
649: MatMultTranspose(J,VecSensiTemp[nadj],VecDeltaLam[nadj*s+i]);
650: if (ts->vec_costintegral) {
651: VecAXPY(VecDeltaLam[nadj*s+i],-h*b[i],ts->vecs_drdy[nadj]);
652: }
653: }
655: /* Stage values of mu */
656: if(ts->vecs_sensip) {
657: TSAdjointComputeRHSJacobian(ts,rk->stage_time,Y[i],ts->Jacp);
658: if (ts->vec_costintegral) {
659: TSAdjointComputeDRDPFunction(ts,rk->stage_time,Y[i],ts->vecs_drdp);
660: }
662: for (nadj=0; nadj<ts->numcost; nadj++) {
663: MatMultTranspose(ts->Jacp,VecSensiTemp[nadj],VecDeltaMu[nadj*s+i]);
664: if (ts->vec_costintegral) {
665: VecAXPY(VecDeltaMu[nadj*s+i],-h*b[i],ts->vecs_drdp[nadj]);
666: }
667: }
668: }
669: }
671: for (j=0; j<s; j++) w[j] = 1.0;
672: for (nadj=0; nadj<ts->numcost; nadj++) {
673: VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecDeltaLam[nadj*s]);
674: if(ts->vecs_sensip) {
675: VecMAXPY(ts->vecs_sensip[nadj],s,w,&VecDeltaMu[nadj*s]);
676: }
677: }
678: rk->status = TS_STEP_COMPLETE;
679: return(0);
680: }
682: static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X)
683: {
684: TS_RK *rk = (TS_RK*)ts->data;
685: PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j;
686: PetscReal h;
687: PetscReal tt,t;
688: PetscScalar *b;
689: const PetscReal *B = rk->tableau->binterp;
690: PetscErrorCode ierr;
693: if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name);
695: switch (rk->status) {
696: case TS_STEP_INCOMPLETE:
697: case TS_STEP_PENDING:
698: h = ts->time_step;
699: t = (itime - ts->ptime)/h;
700: break;
701: case TS_STEP_COMPLETE:
702: h = ts->ptime - ts->ptime_prev;
703: t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
704: break;
705: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
706: }
707: PetscMalloc1(s,&b);
708: for (i=0; i<s; i++) b[i] = 0;
709: for (j=0,tt=t; j<p; j++,tt*=t) {
710: for (i=0; i<s; i++) {
711: b[i] += h * B[i*p+j] * tt;
712: }
713: }
715: VecCopy(rk->Y[0],X);
716: VecMAXPY(X,s,b,rk->YdotRHS);
718: PetscFree(b);
719: return(0);
720: }
722: /*------------------------------------------------------------*/
724: static PetscErrorCode TSRKTableauReset(TS ts)
725: {
726: TS_RK *rk = (TS_RK*)ts->data;
727: RKTableau tab = rk->tableau;
731: if (!tab) return(0);
732: PetscFree(rk->work);
733: VecDestroyVecs(tab->s,&rk->Y);
734: VecDestroyVecs(tab->s,&rk->YdotRHS);
735: VecDestroyVecs(tab->s*ts->numcost,&rk->VecDeltaLam);
736: VecDestroyVecs(tab->s*ts->numcost,&rk->VecDeltaMu);
737: return(0);
738: }
740: static PetscErrorCode TSReset_RK(TS ts)
741: {
742: TS_RK *rk = (TS_RK*)ts->data;
746: TSRKTableauReset(ts);
747: VecDestroy(&rk->VecCostIntegral0);
748: VecDestroyVecs(ts->numcost,&rk->VecSensiTemp);
749: return(0);
750: }
752: static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx)
753: {
755: return(0);
756: }
758: static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
759: {
761: return(0);
762: }
765: static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx)
766: {
768: return(0);
769: }
771: static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
772: {
775: return(0);
776: }
777: /*
778: static PetscErrorCode RKSetAdjCoe(RKTableau tab)
779: {
780: PetscReal *A,*b;
781: PetscInt s,i,j;
782: PetscErrorCode ierr;
785: s = tab->s;
786: PetscMalloc2(s*s,&A,s,&b);
788: for (i=0; i<s; i++)
789: for (j=0; j<s; j++) {
790: A[i*s+j] = (tab->b[s-1-i]==0) ? 0: -tab->A[s-1-i+(s-1-j)*s] * tab->b[s-1-j] / tab->b[s-1-i];
791: PetscPrintf(PETSC_COMM_WORLD,"Coefficients: A[%D][%D]=%.6f\n",i,j,A[i*s+j]);
792: }
794: for (i=0; i<s; i++) b[i] = (tab->b[s-1-i]==0)? 0: -tab->b[s-1-i];
796: PetscMemcpy(tab->A,A,s*s*sizeof(A[0]));
797: PetscMemcpy(tab->b,b,s*sizeof(b[0]));
798: PetscFree2(A,b);
799: return(0);
800: }
801: */
803: static PetscErrorCode TSRKTableauSetUp(TS ts)
804: {
805: TS_RK *rk = (TS_RK*)ts->data;
806: RKTableau tab = rk->tableau;
810: PetscMalloc1(tab->s,&rk->work);
811: VecDuplicateVecs(ts->vec_sol,tab->s,&rk->Y);
812: VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS);
813: return(0);
814: }
817: static PetscErrorCode TSSetUp_RK(TS ts)
818: {
819: TS_RK *rk = (TS_RK*)ts->data;
821: DM dm;
824: TSCheckImplicitTerm(ts);
825: TSRKTableauSetUp(ts);
826: if (!rk->VecCostIntegral0 && ts->vec_costintegral && ts->costintegralfwd) { /* back up cost integral */
827: VecDuplicate(ts->vec_costintegral,&rk->VecCostIntegral0);
828: }
829: TSGetDM(ts,&dm);
830: DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);
831: DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);
832: return(0);
833: }
836: /*------------------------------------------------------------*/
838: static PetscErrorCode TSSetFromOptions_RK(PetscOptionItems *PetscOptionsObject,TS ts)
839: {
840: TS_RK *rk = (TS_RK*)ts->data;
844: PetscOptionsHead(PetscOptionsObject,"RK ODE solver options");
845: {
846: RKTableauLink link;
847: PetscInt count,choice;
848: PetscBool flg;
849: const char **namelist;
850: for (link=RKTableauList,count=0; link; link=link->next,count++) ;
851: PetscMalloc1(count,&namelist);
852: for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
853: PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rk->tableau->name,&choice,&flg);
854: if (flg) {TSRKSetType(ts,namelist[choice]);}
855: PetscFree(namelist);
856: }
857: PetscOptionsTail();
858: return(0);
859: }
861: static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[])
862: {
864: PetscInt i;
865: size_t left,count;
866: char *p;
869: for (i=0,p=buf,left=len; i<n; i++) {
870: PetscSNPrintfCount(p,left,fmt,&count,x[i]);
871: if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer");
872: left -= count;
873: p += count;
874: *p++ = ' ';
875: }
876: p[i ? 0 : -1] = 0;
877: return(0);
878: }
880: static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer)
881: {
882: TS_RK *rk = (TS_RK*)ts->data;
883: PetscBool iascii;
887: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
888: if (iascii) {
889: RKTableau tab = rk->tableau;
890: TSRKType rktype;
891: char buf[512];
892: TSRKGetType(ts,&rktype);
893: PetscViewerASCIIPrintf(viewer," RK type %s\n",rktype);
894: PetscViewerASCIIPrintf(viewer," Order: %D\n",tab->order);
895: PetscViewerASCIIPrintf(viewer," FSAL property: %s\n",tab->FSAL ? "yes" : "no");
896: PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);
897: PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);
898: }
899: return(0);
900: }
902: static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer)
903: {
905: TSAdapt adapt;
908: TSGetAdapt(ts,&adapt);
909: TSAdaptLoad(adapt,viewer);
910: return(0);
911: }
913: /*@C
914: TSRKSetType - Set the type of RK scheme
916: Logically collective
918: Input Parameter:
919: + ts - timestepping context
920: - rktype - type of RK-scheme
922: Options Database:
923: . -ts_rk_type - <1fe,2a,3,3bs,4,5f,5dp,5bs>
925: Level: intermediate
927: .seealso: TSRKGetType(), TSRK, TSRKType, TSRK1FE, TSRK2A, TSRK3, TSRK3BS, TSRK4, TSRK5F, TSRK5DP, TSRK5BS
928: @*/
929: PetscErrorCode TSRKSetType(TS ts,TSRKType rktype)
930: {
936: PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));
937: return(0);
938: }
940: /*@C
941: TSRKGetType - Get the type of RK scheme
943: Logically collective
945: Input Parameter:
946: . ts - timestepping context
948: Output Parameter:
949: . rktype - type of RK-scheme
951: Level: intermediate
953: .seealso: TSRKGetType()
954: @*/
955: PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype)
956: {
961: PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));
962: return(0);
963: }
965: static PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype)
966: {
967: TS_RK *rk = (TS_RK*)ts->data;
970: *rktype = rk->tableau->name;
971: return(0);
972: }
973: static PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype)
974: {
975: TS_RK *rk = (TS_RK*)ts->data;
977: PetscBool match;
978: RKTableauLink link;
981: if (rk->tableau) {
982: PetscStrcmp(rk->tableau->name,rktype,&match);
983: if (match) return(0);
984: }
985: for (link = RKTableauList; link; link=link->next) {
986: PetscStrcmp(link->tab.name,rktype,&match);
987: if (match) {
988: if (ts->setupcalled) {TSRKTableauReset(ts);}
989: rk->tableau = &link->tab;
990: if (ts->setupcalled) {TSRKTableauSetUp(ts);}
991: ts->default_adapt_type = rk->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE;
992: return(0);
993: }
994: }
995: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype);
996: return(0);
997: }
999: static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y)
1000: {
1001: TS_RK *rk = (TS_RK*)ts->data;
1004: *ns = rk->tableau->s;
1005: if(Y) *Y = rk->Y;
1006: return(0);
1007: }
1009: static PetscErrorCode TSDestroy_RK(TS ts)
1010: {
1014: TSReset_RK(ts);
1015: if (ts->dm) {
1016: DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);
1017: DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);
1018: }
1019: PetscFree(ts->data);
1020: PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);
1021: PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);
1022: return(0);
1023: }
1025: /* ------------------------------------------------------------ */
1026: /*MC
1027: TSRK - ODE and DAE solver using Runge-Kutta schemes
1029: The user should provide the right hand side of the equation
1030: using TSSetRHSFunction().
1032: Notes:
1033: The default is TSRK3BS, it can be changed with TSRKSetType() or -ts_rk_type
1035: Level: beginner
1037: .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3,
1038: TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister()
1040: M*/
1041: PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts)
1042: {
1043: TS_RK *rk;
1047: TSRKInitializePackage();
1049: ts->ops->reset = TSReset_RK;
1050: ts->ops->destroy = TSDestroy_RK;
1051: ts->ops->view = TSView_RK;
1052: ts->ops->load = TSLoad_RK;
1053: ts->ops->setup = TSSetUp_RK;
1054: ts->ops->adjointsetup = TSAdjointSetUp_RK;
1055: ts->ops->step = TSStep_RK;
1056: ts->ops->interpolate = TSInterpolate_RK;
1057: ts->ops->evaluatestep = TSEvaluateStep_RK;
1058: ts->ops->rollback = TSRollBack_RK;
1059: ts->ops->setfromoptions = TSSetFromOptions_RK;
1060: ts->ops->getstages = TSGetStages_RK;
1061: ts->ops->adjointstep = TSAdjointStep_RK;
1063: ts->ops->adjointintegral = TSAdjointCostIntegral_RK;
1064: ts->ops->forwardintegral = TSForwardCostIntegral_RK;
1066: PetscNewLog(ts,&rk);
1067: ts->data = (void*)rk;
1069: PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);
1070: PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);
1072: TSRKSetType(ts,TSRKDefault);
1073: return(0);
1074: }