Actual source code: rk.c
petsc-3.7.5 2017-01-01
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> /*I "petscts.h" I*/
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: PetscBool FSAL; /* flag to indicate if tableau is FSAL */
23: PetscInt pinterp; /* Interpolation order */
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: TSRK1 - First order forward Euler scheme.
54: This method has one stage.
56: Level: advanced
58: .seealso: TSRK
59: M*/
60: /*MC
61: TSRK2A - Second order RK scheme.
63: This method has two stages.
65: Level: advanced
67: .seealso: TSRK
68: M*/
69: /*MC
70: TSRK3 - Third order RK scheme.
72: This method has three stages.
74: Level: advanced
76: .seealso: TSRK
77: M*/
78: /*MC
79: TSRK3BS - Third order RK scheme of Bogacki-Shampine with 2nd order embedded method.
81: This method has four stages.
83: Level: advanced
85: .seealso: TSRK
86: M*/
87: /*MC
88: TSRK4 - Fourth order RK scheme.
90: This is the classical Runge-Kutta method with four stages.
92: Level: advanced
94: .seealso: TSRK
95: M*/
96: /*MC
97: TSRK5F - Fifth order Fehlberg RK scheme with a 4th order embedded method.
99: This method has six stages.
101: Level: advanced
103: .seealso: TSRK
104: M*/
105: /*MC
106: TSRK5DP - Fifth order Dormand-Prince RK scheme with the 4th order embedded method.
108: This method has seven stages.
110: Level: advanced
112: .seealso: TSRK
113: M*/
117: /*@C
118: TSRKRegisterAll - Registers all of the Runge-Kutta explicit methods in TSRK
120: Not Collective, but should be called by all processes which will need the schemes to be registered
122: Level: advanced
124: .keywords: TS, TSRK, register, all
126: .seealso: TSRKRegisterDestroy()
127: @*/
128: PetscErrorCode TSRKRegisterAll(void)
129: {
133: if (TSRKRegisterAllCalled) return(0);
134: TSRKRegisterAllCalled = PETSC_TRUE;
136: {
137: const PetscReal
138: A[1][1] = {{0.0}},
139: b[1] = {1.0};
140: TSRKRegister(TSRK1FE,1,1,&A[0][0],b,NULL,NULL,1,b);
141: }
142: {
143: const PetscReal
144: A[2][2] = {{0.0,0.0},
145: {1.0,0.0}},
146: b[2] = {0.5,0.5},
147: bembed[2] = {1.0,0};
148: TSRKRegister(TSRK2A,2,2,&A[0][0],b,NULL,bembed,1,b);
149: }
150: {
151: const PetscReal
152: A[3][3] = {{0,0,0},
153: {2.0/3.0,0,0},
154: {-1.0/3.0,1.0,0}},
155: b[3] = {0.25,0.5,0.25};
156: TSRKRegister(TSRK3,3,3,&A[0][0],b,NULL,NULL,1,b);
157: }
158: {
159: const PetscReal
160: A[4][4] = {{0,0,0,0},
161: {1.0/2.0,0,0,0},
162: {0,3.0/4.0,0,0},
163: {2.0/9.0,1.0/3.0,4.0/9.0,0}},
164: b[4] = {2.0/9.0,1.0/3.0,4.0/9.0,0},
165: bembed[4] = {7.0/24.0,1.0/4.0,1.0/3.0,1.0/8.0};
166: TSRKRegister(TSRK3BS,3,4,&A[0][0],b,NULL,bembed,1,b);
167: }
168: {
169: const PetscReal
170: A[4][4] = {{0,0,0,0},
171: {0.5,0,0,0},
172: {0,0.5,0,0},
173: {0,0,1.0,0}},
174: b[4] = {1.0/6.0,1.0/3.0,1.0/3.0,1.0/6.0};
175: TSRKRegister(TSRK4,4,4,&A[0][0],b,NULL,NULL,1,b);
176: }
177: {
178: const PetscReal
179: A[6][6] = {{0,0,0,0,0,0},
180: {0.25,0,0,0,0,0},
181: {3.0/32.0,9.0/32.0,0,0,0,0},
182: {1932.0/2197.0,-7200.0/2197.0,7296.0/2197.0,0,0,0},
183: {439.0/216.0,-8.0,3680.0/513.0,-845.0/4104.0,0,0},
184: {-8.0/27.0,2.0,-3544.0/2565.0,1859.0/4104.0,-11.0/40.0,0}},
185: b[6] = {16.0/135.0,0,6656.0/12825.0,28561.0/56430.0,-9.0/50.0,2.0/55.0},
186: bembed[6] = {25.0/216.0,0,1408.0/2565.0,2197.0/4104.0,-1.0/5.0,0};
187: TSRKRegister(TSRK5F,5,6,&A[0][0],b,NULL,bembed,1,b);
188: }
189: {
190: const PetscReal
191: A[7][7] = {{0,0,0,0,0,0,0},
192: {1.0/5.0,0,0,0,0,0,0},
193: {3.0/40.0,9.0/40.0,0,0,0,0,0},
194: {44.0/45.0,-56.0/15.0,32.0/9.0,0,0,0,0},
195: {19372.0/6561.0,-25360.0/2187.0,64448.0/6561.0,-212.0/729.0,0,0,0},
196: {9017.0/3168.0,-355.0/33.0,46732.0/5247.0,49.0/176.0,-5103.0/18656.0,0,0},
197: {35.0/384.0,0,500.0/1113.0,125.0/192.0,-2187.0/6784.0,11.0/84.0,0}},
198: b[7] = {35.0/384.0,0,500.0/1113.0,125.0/192.0,-2187.0/6784.0,11.0/84.0,0},
199: bembed[7] = {5179.0/57600.0,0,7571.0/16695.0,393.0/640.0,-92097.0/339200.0,187.0/2100.0,1.0/40.0};
200: TSRKRegister(TSRK5DP,5,7,&A[0][0],b,NULL,bembed,1,b);
201: }
202: return(0);
203: }
207: /*@C
208: TSRKRegisterDestroy - Frees the list of schemes that were registered by TSRKRegister().
210: Not Collective
212: Level: advanced
214: .keywords: TSRK, register, destroy
215: .seealso: TSRKRegister(), TSRKRegisterAll()
216: @*/
217: PetscErrorCode TSRKRegisterDestroy(void)
218: {
220: RKTableauLink link;
223: while ((link = RKTableauList)) {
224: RKTableau t = &link->tab;
225: RKTableauList = link->next;
226: PetscFree3(t->A,t->b,t->c);
227: PetscFree (t->bembed);
228: PetscFree (t->binterp);
229: PetscFree (t->name);
230: PetscFree (link);
231: }
232: TSRKRegisterAllCalled = PETSC_FALSE;
233: return(0);
234: }
238: /*@C
239: TSRKInitializePackage - This function initializes everything in the TSRK package. It is called
240: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RK()
241: when using static libraries.
243: Level: developer
245: .keywords: TS, TSRK, initialize, package
246: .seealso: PetscInitialize()
247: @*/
248: PetscErrorCode TSRKInitializePackage(void)
249: {
253: if (TSRKPackageInitialized) return(0);
254: TSRKPackageInitialized = PETSC_TRUE;
255: TSRKRegisterAll();
256: PetscRegisterFinalize(TSRKFinalizePackage);
257: return(0);
258: }
262: /*@C
263: TSRKFinalizePackage - This function destroys everything in the TSRK package. It is
264: called from PetscFinalize().
266: Level: developer
268: .keywords: Petsc, destroy, package
269: .seealso: PetscFinalize()
270: @*/
271: PetscErrorCode TSRKFinalizePackage(void)
272: {
276: TSRKPackageInitialized = PETSC_FALSE;
277: TSRKRegisterDestroy();
278: return(0);
279: }
283: /*@C
284: TSRKRegister - register an RK scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation
286: Not Collective, but the same schemes should be registered on all processes on which they will be used
288: Input Parameters:
289: + name - identifier for method
290: . order - approximation order of method
291: . s - number of stages, this is the dimension of the matrices below
292: . A - stage coefficients (dimension s*s, row-major)
293: . b - step completion table (dimension s; NULL to use last row of A)
294: . c - abscissa (dimension s; NULL to use row sums of A)
295: . bembed - completion table for embedded method (dimension s; NULL if not available)
296: . pinterp - Order of the interpolation scheme, equal to the number of columns of binterp
297: - binterp - Coefficients of the interpolation formula (dimension s*pinterp; NULL to reuse binterpt)
299: Notes:
300: Several RK methods are provided, this function is only needed to create new methods.
302: Level: advanced
304: .keywords: TS, register
306: .seealso: TSRK
307: @*/
308: PetscErrorCode TSRKRegister(TSRKType name,PetscInt order,PetscInt s,
309: const PetscReal A[],const PetscReal b[],const PetscReal c[],
310: const PetscReal bembed[],
311: PetscInt pinterp,const PetscReal binterp[])
312: {
313: PetscErrorCode ierr;
314: RKTableauLink link;
315: RKTableau t;
316: PetscInt i,j;
319: PetscMalloc(sizeof(*link),&link);
320: PetscMemzero(link,sizeof(*link));
321: t = &link->tab;
322: PetscStrallocpy(name,&t->name);
323: t->order = order;
324: t->s = s;
325: PetscMalloc3(s*s,&t->A,s,&t->b,s,&t->c);
326: PetscMemcpy(t->A,A,s*s*sizeof(A[0]));
327: if (b) { PetscMemcpy(t->b,b,s*sizeof(b[0])); }
328: else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i];
329: if (c) { PetscMemcpy(t->c,c,s*sizeof(c[0])); }
330: else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j];
331: t->FSAL = PETSC_TRUE;
332: for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE;
333: if (bembed) {
334: PetscMalloc1(s,&t->bembed);
335: PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));
336: }
338: t->pinterp = pinterp;
339: PetscMalloc1(s*pinterp,&t->binterp);
340: PetscMemcpy(t->binterp,binterp,s*pinterp*sizeof(binterp[0]));
341: link->next = RKTableauList;
342: RKTableauList = link;
343: return(0);
344: }
348: /*
349: The step completion formula is
351: x1 = x0 + h b^T YdotRHS
353: This function can be called before or after ts->vec_sol has been updated.
354: Suppose we have a completion formula (b) and an embedded formula (be) of different order.
355: We can write
357: x1e = x0 + h be^T YdotRHS
358: = x1 - h b^T YdotRHS + h be^T YdotRHS
359: = x1 + h (be - b)^T YdotRHS
361: so we can evaluate the method with different order even after the step has been optimistically completed.
362: */
363: static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done)
364: {
365: TS_RK *rk = (TS_RK*)ts->data;
366: RKTableau tab = rk->tableau;
367: PetscScalar *w = rk->work;
368: PetscReal h;
369: PetscInt s = tab->s,j;
373: switch (rk->status) {
374: case TS_STEP_INCOMPLETE:
375: case TS_STEP_PENDING:
376: h = ts->time_step; break;
377: case TS_STEP_COMPLETE:
378: h = ts->ptime - ts->ptime_prev; break;
379: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
380: }
381: if (order == tab->order) {
382: if (rk->status == TS_STEP_INCOMPLETE) {
383: VecCopy(ts->vec_sol,X);
384: for (j=0; j<s; j++) w[j] = h*tab->b[j];
385: VecMAXPY(X,s,w,rk->YdotRHS);
386: } else {VecCopy(ts->vec_sol,X);}
387: return(0);
388: } else if (order == tab->order-1) {
389: if (!tab->bembed) goto unavailable;
390: if (rk->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (be) */
391: VecCopy(ts->vec_sol,X);
392: for (j=0; j<s; j++) w[j] = h*tab->bembed[j];
393: VecMAXPY(X,s,w,rk->YdotRHS);
394: } else { /* Rollback and re-complete using (be-b) */
395: VecCopy(ts->vec_sol,X);
396: for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]);
397: VecMAXPY(X,s,w,rk->YdotRHS);
398: if (ts->vec_costintegral && ts->costintegralfwd) {
399: VecCopy(rk->VecCostIntegral0,ts->vec_costintegral);
400: }
401: }
402: if (done) *done = PETSC_TRUE;
403: return(0);
404: }
405: unavailable:
406: if (done) *done = PETSC_FALSE;
407: 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);
408: return(0);
409: }
413: static PetscErrorCode TSForwardCostIntegral_RK(TS ts)
414: {
415: TS_RK *rk = (TS_RK*)ts->data;
416: RKTableau tab = rk->tableau;
417: const PetscInt s = tab->s;
418: const PetscReal *b = tab->b,*c = tab->c;
419: Vec *Y = rk->Y;
420: PetscInt i;
421: PetscErrorCode ierr;
424: /* backup cost integral */
425: VecCopy(ts->vec_costintegral,rk->VecCostIntegral0);
426: for (i=s-1; i>=0; i--) {
427: /* Evolve ts->vec_costintegral to compute integrals */
428: TSAdjointComputeCostIntegrand(ts,rk->ptime+rk->time_step*(1.0-c[i]),Y[i],ts->vec_costintegrand);
429: VecAXPY(ts->vec_costintegral,rk->time_step*b[i],ts->vec_costintegrand);
430: }
431: return(0);
432: }
436: static PetscErrorCode TSAdjointCostIntegral_RK(TS ts)
437: {
438: TS_RK *rk = (TS_RK*)ts->data;
439: RKTableau tab = rk->tableau;
440: const PetscInt s = tab->s;
441: const PetscReal *b = tab->b,*c = tab->c;
442: Vec *Y = rk->Y;
443: PetscInt i;
444: PetscErrorCode ierr;
447: for (i=s-1; i>=0; i--) {
448: /* Evolve ts->vec_costintegral to compute integrals */
449: TSAdjointComputeCostIntegrand(ts,ts->ptime-ts->time_step*(1.0-c[i]),Y[i],ts->vec_costintegrand);
450: VecAXPY(ts->vec_costintegral,-ts->time_step*b[i],ts->vec_costintegrand);
451: }
452: return(0);
453: }
457: static PetscErrorCode TSRollBack_RK(TS ts)
458: {
459: TS_RK *rk = (TS_RK*)ts->data;
460: RKTableau tab = rk->tableau;
461: const PetscInt s = tab->s;
462: const PetscReal *b = tab->b;
463: PetscScalar *w = rk->work;
464: Vec *YdotRHS = rk->YdotRHS;
465: PetscInt j;
466: PetscReal h;
467: PetscErrorCode ierr;
470: switch (rk->status) {
471: case TS_STEP_INCOMPLETE:
472: case TS_STEP_PENDING:
473: h = ts->time_step; break;
474: case TS_STEP_COMPLETE:
475: h = ts->ptime - ts->ptime_prev; break;
476: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
477: }
478: for (j=0; j<s; j++) w[j] = -h*b[j];
479: VecMAXPY(ts->vec_sol,s,w,YdotRHS);
480: return(0);
481: }
485: static PetscErrorCode TSStep_RK(TS ts)
486: {
487: TS_RK *rk = (TS_RK*)ts->data;
488: RKTableau tab = rk->tableau;
489: const PetscInt s = tab->s;
490: const PetscReal *A = tab->A,*c = tab->c;
491: PetscScalar *w = rk->work;
492: Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS;
493: TSAdapt adapt;
494: PetscInt i,j;
495: PetscInt rejections = 0;
496: PetscBool stageok,accept = PETSC_TRUE;
497: PetscReal next_time_step = ts->time_step;
498: PetscErrorCode ierr;
502: rk->status = TS_STEP_INCOMPLETE;
503: while (!ts->reason && rk->status != TS_STEP_COMPLETE) {
504: PetscReal t = ts->ptime;
505: PetscReal h = ts->time_step;
506: for (i=0; i<s; i++) {
507: rk->stage_time = t + h*c[i];
508: TSPreStage(ts,rk->stage_time);
509: VecCopy(ts->vec_sol,Y[i]);
510: for (j=0; j<i; j++) w[j] = h*A[i*s+j];
511: VecMAXPY(Y[i],i,w,YdotRHS);
512: TSPostStage(ts,rk->stage_time,i,Y);
513: TSGetAdapt(ts,&adapt);
514: TSAdaptCheckStage(adapt,ts,rk->stage_time,Y[i],&stageok);
515: if (!stageok) goto reject_step;
516: TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);
517: }
519: rk->status = TS_STEP_INCOMPLETE;
520: TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);
521: rk->status = TS_STEP_PENDING;
522: TSGetAdapt(ts,&adapt);
523: TSAdaptCandidatesClear(adapt);
524: TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);
525: TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
526: rk->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
527: if (!accept) { /* Roll back the current step */
528: TSRollBack_RK(ts);
529: ts->time_step = next_time_step;
530: goto reject_step;
531: }
533: if (ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation*/
534: rk->ptime = ts->ptime;
535: rk->time_step = ts->time_step;
536: }
538: ts->ptime += ts->time_step;
539: ts->time_step = next_time_step;
540: break;
542: reject_step:
543: ts->reject++; accept = PETSC_FALSE;
544: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
545: ts->reason = TS_DIVERGED_STEP_REJECTED;
546: PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
547: }
548: }
549: return(0);
550: }
554: static PetscErrorCode TSAdjointSetUp_RK(TS ts)
555: {
556: TS_RK *rk = (TS_RK*)ts->data;
557: RKTableau tab = rk->tableau;
558: PetscInt s = tab->s;
562: if (ts->adjointsetupcalled++) return(0);
563: VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecDeltaLam);
564: if(ts->vecs_sensip) {
565: VecDuplicateVecs(ts->vecs_sensip[0],s*ts->numcost,&rk->VecDeltaMu);
566: }
567: VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&rk->VecSensiTemp);
568: return(0);
569: }
573: static PetscErrorCode TSAdjointStep_RK(TS ts)
574: {
575: TS_RK *rk = (TS_RK*)ts->data;
576: RKTableau tab = rk->tableau;
577: const PetscInt s = tab->s;
578: const PetscReal *A = tab->A,*b = tab->b,*c = tab->c;
579: PetscScalar *w = rk->work;
580: Vec *Y = rk->Y,*VecDeltaLam = rk->VecDeltaLam,*VecDeltaMu = rk->VecDeltaMu,*VecSensiTemp = rk->VecSensiTemp;
581: PetscInt i,j,nadj;
582: PetscReal t = ts->ptime;
583: PetscErrorCode ierr;
584: PetscReal h = ts->time_step;
585: Mat J,Jp;
588: rk->status = TS_STEP_INCOMPLETE;
589: for (i=s-1; i>=0; i--) {
590: rk->stage_time = t + h*(1.0-c[i]);
591: for (nadj=0; nadj<ts->numcost; nadj++) {
592: VecCopy(ts->vecs_sensi[nadj],VecSensiTemp[nadj]);
593: VecScale(VecSensiTemp[nadj],-h*b[i]);
594: for (j=i+1; j<s; j++) {
595: VecAXPY(VecSensiTemp[nadj],-h*A[j*s+i],VecDeltaLam[nadj*s+j]);
596: }
597: }
598: /* Stage values of lambda */
599: TSGetRHSJacobian(ts,&J,&Jp,NULL,NULL);
600: TSComputeRHSJacobian(ts,rk->stage_time,Y[i],J,Jp);
601: if (ts->vec_costintegral) {
602: TSAdjointComputeDRDYFunction(ts,rk->stage_time,Y[i],ts->vecs_drdy);
603: }
604: for (nadj=0; nadj<ts->numcost; nadj++) {
605: MatMultTranspose(J,VecSensiTemp[nadj],VecDeltaLam[nadj*s+i]);
606: if (ts->vec_costintegral) {
607: VecAXPY(VecDeltaLam[nadj*s+i],-h*b[i],ts->vecs_drdy[nadj]);
608: }
609: }
611: /* Stage values of mu */
612: if(ts->vecs_sensip) {
613: TSAdjointComputeRHSJacobian(ts,rk->stage_time,Y[i],ts->Jacp);
614: if (ts->vec_costintegral) {
615: TSAdjointComputeDRDPFunction(ts,rk->stage_time,Y[i],ts->vecs_drdp);
616: }
618: for (nadj=0; nadj<ts->numcost; nadj++) {
619: MatMultTranspose(ts->Jacp,VecSensiTemp[nadj],VecDeltaMu[nadj*s+i]);
620: if (ts->vec_costintegral) {
621: VecAXPY(VecDeltaMu[nadj*s+i],-h*b[i],ts->vecs_drdp[nadj]);
622: }
623: }
624: }
625: }
627: for (j=0; j<s; j++) w[j] = 1.0;
628: for (nadj=0; nadj<ts->numcost; nadj++) {
629: VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecDeltaLam[nadj*s]);
630: if(ts->vecs_sensip) {
631: VecMAXPY(ts->vecs_sensip[nadj],s,w,&VecDeltaMu[nadj*s]);
632: }
633: }
634: rk->status = TS_STEP_COMPLETE;
635: return(0);
636: }
640: static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X)
641: {
642: TS_RK *rk = (TS_RK*)ts->data;
643: PetscInt s = rk->tableau->s,pinterp = rk->tableau->pinterp,i,j;
644: PetscReal h;
645: PetscReal tt,t;
646: PetscScalar *b;
647: const PetscReal *B = rk->tableau->binterp;
648: PetscErrorCode ierr;
651: if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name);
653: switch (rk->status) {
654: case TS_STEP_INCOMPLETE:
655: case TS_STEP_PENDING:
656: h = ts->time_step;
657: t = (itime - ts->ptime)/h;
658: break;
659: case TS_STEP_COMPLETE:
660: h = ts->ptime - ts->ptime_prev;
661: t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
662: break;
663: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
664: }
665: PetscMalloc1(s,&b);
666: for (i=0; i<s; i++) b[i] = 0;
667: for (j=0,tt=t; j<pinterp; j++,tt*=t) {
668: for (i=0; i<s; i++) {
669: b[i] += h * B[i*pinterp+j] * tt;
670: }
671: }
673: VecCopy(rk->Y[0],X);
674: VecMAXPY(X,s,b,rk->YdotRHS);
676: PetscFree(b);
677: return(0);
678: }
680: /*------------------------------------------------------------*/
684: static PetscErrorCode TSRKTableauReset(TS ts)
685: {
686: TS_RK *rk = (TS_RK*)ts->data;
687: RKTableau tab = rk->tableau;
691: if (!tab) return(0);
692: PetscFree(rk->work);
693: VecDestroyVecs(tab->s,&rk->Y);
694: VecDestroyVecs(tab->s,&rk->YdotRHS);
695: VecDestroyVecs(tab->s*ts->numcost,&rk->VecDeltaLam);
696: VecDestroyVecs(tab->s*ts->numcost,&rk->VecDeltaMu);
697: return(0);
698: }
702: static PetscErrorCode TSReset_RK(TS ts)
703: {
704: TS_RK *rk = (TS_RK*)ts->data;
708: TSRKTableauReset(ts);
709: VecDestroy(&rk->VecCostIntegral0);
710: VecDestroyVecs(ts->numcost,&rk->VecSensiTemp);
711: return(0);
712: }
716: static PetscErrorCode TSDestroy_RK(TS ts)
717: {
721: TSReset_RK(ts);
722: PetscFree(ts->data);
723: PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);
724: PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);
725: return(0);
726: }
731: static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx)
732: {
734: return(0);
735: }
739: static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
740: {
742: return(0);
743: }
748: static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx)
749: {
751: return(0);
752: }
756: static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
757: {
760: return(0);
761: }
762: /*
765: static PetscErrorCode RKSetAdjCoe(RKTableau tab)
766: {
767: PetscReal *A,*b;
768: PetscInt s,i,j;
769: PetscErrorCode ierr;
772: s = tab->s;
773: PetscMalloc2(s*s,&A,s,&b);
775: for (i=0; i<s; i++)
776: for (j=0; j<s; j++) {
777: 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];
778: PetscPrintf(PETSC_COMM_WORLD,"Coefficients: A[%D][%D]=%.6f\n",i,j,A[i*s+j]);
779: }
781: for (i=0; i<s; i++) b[i] = (tab->b[s-1-i]==0)? 0: -tab->b[s-1-i];
783: PetscMemcpy(tab->A,A,s*s*sizeof(A[0]));
784: PetscMemcpy(tab->b,b,s*sizeof(b[0]));
785: PetscFree2(A,b);
786: return(0);
787: }
788: */
792: static PetscErrorCode TSRKTableauSetUp(TS ts)
793: {
794: TS_RK *rk = (TS_RK*)ts->data;
795: RKTableau tab = rk->tableau;
799: PetscMalloc1(tab->s,&rk->work);
800: VecDuplicateVecs(ts->vec_sol,tab->s,&rk->Y);
801: VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS);
802: return(0);
803: }
808: static PetscErrorCode TSSetUp_RK(TS ts)
809: {
810: TS_RK *rk = (TS_RK*)ts->data;
812: DM dm;
815: TSRKTableauSetUp(ts);
816: if (!rk->VecCostIntegral0 && ts->vec_costintegral && ts->costintegralfwd) { /* back up cost integral */
817: VecDuplicate(ts->vec_costintegral,&rk->VecCostIntegral0);
818: }
819: TSGetDM(ts,&dm);
820: if (dm) {
821: DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);
822: DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);
823: }
824: return(0);
825: }
828: /*------------------------------------------------------------*/
832: static PetscErrorCode TSSetFromOptions_RK(PetscOptionItems *PetscOptionsObject,TS ts)
833: {
834: TS_RK *rk = (TS_RK*)ts->data;
838: PetscOptionsHead(PetscOptionsObject,"RK ODE solver options");
839: {
840: RKTableauLink link;
841: PetscInt count,choice;
842: PetscBool flg;
843: const char **namelist;
844: for (link=RKTableauList,count=0; link; link=link->next,count++) ;
845: PetscMalloc1(count,&namelist);
846: for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
847: PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rk->tableau->name,&choice,&flg);
848: if (flg) {TSRKSetType(ts,namelist[choice]);}
849: PetscFree(namelist);
850: }
851: PetscOptionsTail();
852: return(0);
853: }
857: static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[])
858: {
860: PetscInt i;
861: size_t left,count;
862: char *p;
865: for (i=0,p=buf,left=len; i<n; i++) {
866: PetscSNPrintfCount(p,left,fmt,&count,x[i]);
867: if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer");
868: left -= count;
869: p += count;
870: *p++ = ' ';
871: }
872: p[i ? 0 : -1] = 0;
873: return(0);
874: }
878: static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer)
879: {
880: TS_RK *rk = (TS_RK*)ts->data;
881: PetscBool iascii;
885: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
886: if (iascii) {
887: RKTableau tab = rk->tableau;
888: TSRKType rktype;
889: char buf[512];
890: TSRKGetType(ts,&rktype);
891: PetscViewerASCIIPrintf(viewer," RK %s\n",rktype);
892: PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);
893: PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);
894: PetscViewerASCIIPrintf(viewer,"FSAL: %s\n",tab->FSAL ? "yes" : "no");
895: }
896: if (ts->adapt) {TSAdaptView(ts->adapt,viewer);}
897: return(0);
898: }
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: }
915: /*@C
916: TSRKSetType - Set the type of RK scheme
918: Logically collective
920: Input Parameter:
921: + ts - timestepping context
922: - rktype - type of RK-scheme
924: Level: intermediate
926: .seealso: TSRKGetType(), TSRK, TSRK2, TSRK3, TSRKPRSSP2, TSRK4, TSRK5
927: @*/
928: PetscErrorCode TSRKSetType(TS ts,TSRKType rktype)
929: {
934: PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));
935: return(0);
936: }
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: }
967: static PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype)
968: {
969: TS_RK *rk = (TS_RK*)ts->data;
972: *rktype = rk->tableau->name;
973: return(0);
974: }
977: static PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype)
978: {
979: TS_RK *rk = (TS_RK*)ts->data;
981: PetscBool match;
982: RKTableauLink link;
985: if (rk->tableau) {
986: PetscStrcmp(rk->tableau->name,rktype,&match);
987: if (match) return(0);
988: }
989: for (link = RKTableauList; link; link=link->next) {
990: PetscStrcmp(link->tab.name,rktype,&match);
991: if (match) {
992: if (ts->setupcalled) {TSRKTableauReset(ts);}
993: rk->tableau = &link->tab;
994: if (ts->setupcalled) {TSRKTableauSetUp(ts);}
995: return(0);
996: }
997: }
998: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype);
999: return(0);
1000: }
1004: static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y)
1005: {
1006: TS_RK *rk = (TS_RK*)ts->data;
1009: *ns = rk->tableau->s;
1010: if(Y) *Y = rk->Y;
1011: return(0);
1012: }
1015: /* ------------------------------------------------------------ */
1016: /*MC
1017: TSRK - ODE and DAE solver using Runge-Kutta schemes
1019: The user should provide the right hand side of the equation
1020: using TSSetRHSFunction().
1022: Notes:
1023: The default is TSRK3, it can be changed with TSRKSetType() or -ts_rk_type
1025: Level: beginner
1027: .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3,
1028: TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister()
1030: M*/
1033: PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts)
1034: {
1035: TS_RK *rk;
1039: TSRKInitializePackage();
1041: ts->ops->reset = TSReset_RK;
1042: ts->ops->destroy = TSDestroy_RK;
1043: ts->ops->view = TSView_RK;
1044: ts->ops->load = TSLoad_RK;
1045: ts->ops->setup = TSSetUp_RK;
1046: ts->ops->adjointsetup = TSAdjointSetUp_RK;
1047: ts->ops->step = TSStep_RK;
1048: ts->ops->interpolate = TSInterpolate_RK;
1049: ts->ops->evaluatestep = TSEvaluateStep_RK;
1050: ts->ops->rollback = TSRollBack_RK;
1051: ts->ops->setfromoptions = TSSetFromOptions_RK;
1052: ts->ops->getstages = TSGetStages_RK;
1053: ts->ops->adjointstep = TSAdjointStep_RK;
1055: ts->ops->adjointintegral = TSAdjointCostIntegral_RK;
1056: ts->ops->forwardintegral = TSForwardCostIntegral_RK;
1058: PetscNewLog(ts,&rk);
1059: ts->data = (void*)rk;
1061: PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);
1062: PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);
1064: TSRKSetType(ts,TSRKDefault);
1065: return(0);
1066: }