Actual source code: rk.c
petsc-3.6.1 2015-08-06
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;
16: static PetscInt explicit_stage_time_id;
18: typedef struct _RKTableau *RKTableau;
19: struct _RKTableau {
20: char *name;
21: PetscInt order; /* Classical approximation order of the method i */
22: PetscInt s; /* Number of stages */
23: PetscBool FSAL; /* flag to indicate if tableau is FSAL */
24: PetscInt pinterp; /* Interpolation order */
25: PetscReal *A,*b,*c; /* Tableau */
26: PetscReal *bembed; /* Embedded formula of order one less (order-1) */
27: PetscReal *binterp; /* Dense output formula */
28: PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */
29: };
30: typedef struct _RKTableauLink *RKTableauLink;
31: struct _RKTableauLink {
32: struct _RKTableau tab;
33: RKTableauLink next;
34: };
35: static RKTableauLink RKTableauList;
37: typedef struct {
38: RKTableau tableau;
39: Vec *Y; /* States computed during the step */
40: Vec *YdotRHS; /* Function evaluations for the non-stiff part */
41: Vec *VecDeltaLam; /* Increment of the adjoint sensitivity w.r.t IC at stage*/
42: Vec *VecDeltaMu; /* Increment of the adjoint sensitivity w.r.t P at stage*/
43: Vec *VecSensiTemp; /* Vector to be timed with Jacobian transpose*/
44: PetscScalar *work; /* Scalar work */
45: PetscReal stage_time;
46: TSStepStatus status;
47: } TS_RK;
49: /*MC
50: TSRK1 - First order forward Euler scheme.
52: This method has one stage.
54: Level: advanced
56: .seealso: TSRK
57: M*/
58: /*MC
59: TSRK2A - Second order RK scheme.
61: This method has two stages.
63: Level: advanced
65: .seealso: TSRK
66: M*/
67: /*MC
68: TSRK3 - Third order RK scheme.
70: This method has three stages.
72: Level: advanced
74: .seealso: TSRK
75: M*/
76: /*MC
77: TSRK3BS - Third order RK scheme of Bogacki-Shampine with 2nd order embedded method.
79: This method has four stages.
81: Level: advanced
83: .seealso: TSRK
84: M*/
85: /*MC
86: TSRK4 - Fourth order RK scheme.
88: This is the classical Runge-Kutta method with four stages.
90: Level: advanced
92: .seealso: TSRK
93: M*/
94: /*MC
95: TSRK5F - Fifth order Fehlberg RK scheme with a 4th order embedded method.
97: This method has six stages.
99: Level: advanced
101: .seealso: TSRK
102: M*/
103: /*MC
104: TSRK5DP - Fifth order Dormand-Prince RK scheme with the 4th order embedded method.
106: This method has seven stages.
108: Level: advanced
110: .seealso: TSRK
111: M*/
115: /*@C
116: TSRKRegisterAll - Registers all of the Runge-Kutta explicit methods in TSRK
118: Not Collective, but should be called by all processes which will need the schemes to be registered
120: Level: advanced
122: .keywords: TS, TSRK, register, all
124: .seealso: TSRKRegisterDestroy()
125: @*/
126: PetscErrorCode TSRKRegisterAll(void)
127: {
131: if (TSRKRegisterAllCalled) return(0);
132: TSRKRegisterAllCalled = PETSC_TRUE;
134: {
135: const PetscReal
136: A[1][1] = {{0.0}},
137: b[1] = {1.0};
138: TSRKRegister(TSRK1FE,1,1,&A[0][0],b,NULL,NULL,1,b);
139: }
140: {
141: const PetscReal
142: A[2][2] = {{0.0,0.0},
143: {1.0,0.0}},
144: b[2] = {0.5,0.5},
145: bembed[2] = {1.0,0};
146: TSRKRegister(TSRK2A,2,2,&A[0][0],b,NULL,bembed,2,b);
147: }
148: {
149: const PetscReal
150: A[3][3] = {{0,0,0},
151: {2.0/3.0,0,0},
152: {-1.0/3.0,1.0,0}},
153: b[3] = {0.25,0.5,0.25};
154: TSRKRegister(TSRK3,3,3,&A[0][0],b,NULL,NULL,3,b);
155: }
156: {
157: const PetscReal
158: A[4][4] = {{0,0,0,0},
159: {1.0/2.0,0,0,0},
160: {0,3.0/4.0,0,0},
161: {2.0/9.0,1.0/3.0,4.0/9.0,0}},
162: b[4] = {2.0/9.0,1.0/3.0,4.0/9.0,0},
163: bembed[4] = {7.0/24.0,1.0/4.0,1.0/3.0,1.0/8.0};
164: TSRKRegister(TSRK3BS,3,4,&A[0][0],b,NULL,bembed,3,b);
165: }
166: {
167: const PetscReal
168: A[4][4] = {{0,0,0,0},
169: {0.5,0,0,0},
170: {0,0.5,0,0},
171: {0,0,1.0,0}},
172: b[4] = {1.0/6.0,1.0/3.0,1.0/3.0,1.0/6.0};
173: TSRKRegister(TSRK4,4,4,&A[0][0],b,NULL,NULL,4,b);
174: }
175: {
176: const PetscReal
177: A[6][6] = {{0,0,0,0,0,0},
178: {0.25,0,0,0,0,0},
179: {3.0/32.0,9.0/32.0,0,0,0,0},
180: {1932.0/2197.0,-7200.0/2197.0,7296.0/2197.0,0,0,0},
181: {439.0/216.0,-8.0,3680.0/513.0,-845.0/4104.0,0,0},
182: {-8.0/27.0,2.0,-3544.0/2565.0,1859.0/4104.0,-11.0/40.0,0}},
183: b[6] = {16.0/135.0,0,6656.0/12825.0,28561.0/56430.0,-9.0/50.0,2.0/55.0},
184: bembed[6] = {25.0/216.0,0,1408.0/2565.0,2197.0/4104.0,-1.0/5.0,0};
185: TSRKRegister(TSRK5F,5,6,&A[0][0],b,NULL,bembed,5,b);
186: }
187: {
188: const PetscReal
189: A[7][7] = {{0,0,0,0,0,0,0},
190: {1.0/5.0,0,0,0,0,0,0},
191: {3.0/40.0,9.0/40.0,0,0,0,0,0},
192: {44.0/45.0,-56.0/15.0,32.0/9.0,0,0,0,0},
193: {19372.0/6561.0,-25360.0/2187.0,64448.0/6561.0,-212.0/729.0,0,0,0},
194: {9017.0/3168.0,-355.0/33.0,46732.0/5247.0,49.0/176.0,-5103.0/18656.0,0,0},
195: {35.0/384.0,0,500.0/1113.0,125.0/192.0,-2187.0/6784.0,11.0/84.0,0}},
196: b[7] = {35.0/384.0,0,500.0/1113.0,125.0/192.0,-2187.0/6784.0,11.0/84.0,0},
197: 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};
198: TSRKRegister(TSRK5DP,5,7,&A[0][0],b,NULL,bembed,5,b);
199: }
200: return(0);
201: }
205: /*@C
206: TSRKRegisterDestroy - Frees the list of schemes that were registered by TSRKRegister().
208: Not Collective
210: Level: advanced
212: .keywords: TSRK, register, destroy
213: .seealso: TSRKRegister(), TSRKRegisterAll()
214: @*/
215: PetscErrorCode TSRKRegisterDestroy(void)
216: {
218: RKTableauLink link;
221: while ((link = RKTableauList)) {
222: RKTableau t = &link->tab;
223: RKTableauList = link->next;
224: PetscFree3(t->A,t->b,t->c);
225: PetscFree (t->bembed);
226: PetscFree (t->binterp);
227: PetscFree (t->name);
228: PetscFree (link);
229: }
230: TSRKRegisterAllCalled = PETSC_FALSE;
231: return(0);
232: }
236: /*@C
237: TSRKInitializePackage - This function initializes everything in the TSRK package. It is called
238: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RK()
239: when using static libraries.
241: Level: developer
243: .keywords: TS, TSRK, initialize, package
244: .seealso: PetscInitialize()
245: @*/
246: PetscErrorCode TSRKInitializePackage(void)
247: {
251: if (TSRKPackageInitialized) return(0);
252: TSRKPackageInitialized = PETSC_TRUE;
253: TSRKRegisterAll();
254: PetscObjectComposedDataRegister(&explicit_stage_time_id);
255: PetscRegisterFinalize(TSRKFinalizePackage);
256: return(0);
257: }
261: /*@C
262: TSRKFinalizePackage - This function destroys everything in the TSRK package. It is
263: called from PetscFinalize().
265: Level: developer
267: .keywords: Petsc, destroy, package
268: .seealso: PetscFinalize()
269: @*/
270: PetscErrorCode TSRKFinalizePackage(void)
271: {
275: TSRKPackageInitialized = PETSC_FALSE;
276: TSRKRegisterDestroy();
277: return(0);
278: }
282: /*@C
283: TSRKRegister - register an RK scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation
285: Not Collective, but the same schemes should be registered on all processes on which they will be used
287: Input Parameters:
288: + name - identifier for method
289: . order - approximation order of method
290: . s - number of stages, this is the dimension of the matrices below
291: . A - stage coefficients (dimension s*s, row-major)
292: . b - step completion table (dimension s; NULL to use last row of A)
293: . c - abscissa (dimension s; NULL to use row sums of A)
294: . bembed - completion table for embedded method (dimension s; NULL if not available)
295: . pinterp - Order of the interpolation scheme, equal to the number of columns of binterp
296: - binterp - Coefficients of the interpolation formula (dimension s*pinterp; NULL to reuse binterpt)
298: Notes:
299: Several RK methods are provided, this function is only needed to create new methods.
301: Level: advanced
303: .keywords: TS, register
305: .seealso: TSRK
306: @*/
307: PetscErrorCode TSRKRegister(TSRKType name,PetscInt order,PetscInt s,
308: const PetscReal A[],const PetscReal b[],const PetscReal c[],
309: const PetscReal bembed[],
310: PetscInt pinterp,const PetscReal binterp[])
311: {
312: PetscErrorCode ierr;
313: RKTableauLink link;
314: RKTableau t;
315: PetscInt i,j;
318: PetscMalloc(sizeof(*link),&link);
319: PetscMemzero(link,sizeof(*link));
320: t = &link->tab;
321: PetscStrallocpy(name,&t->name);
322: t->order = order;
323: t->s = s;
324: PetscMalloc3(s*s,&t->A,s,&t->b,s,&t->c);
325: PetscMemcpy(t->A,A,s*s*sizeof(A[0]));
326: if (b) { PetscMemcpy(t->b,b,s*sizeof(b[0])); }
327: else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i];
328: if (c) { PetscMemcpy(t->c,c,s*sizeof(c[0])); }
329: else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j];
330: t->FSAL = PETSC_TRUE;
331: for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE;
332: if (bembed) {
333: PetscMalloc1(s,&t->bembed);
334: PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));
335: }
337: t->pinterp = pinterp;
338: PetscMalloc1(s*pinterp,&t->binterp);
339: PetscMemcpy(t->binterp,binterp,s*pinterp*sizeof(binterp[0]));
340: link->next = RKTableauList;
341: RKTableauList = link;
342: return(0);
343: }
347: /*
348: The step completion formula is
350: x1 = x0 + h b^T YdotRHS
352: This function can be called before or after ts->vec_sol has been updated.
353: Suppose we have a completion formula (b) and an embedded formula (be) of different order.
354: We can write
356: x1e = x0 + h be^T YdotRHS
357: = x1 - h b^T YdotRHS + h be^T YdotRHS
358: = x1 + h (be - b)^T YdotRHS
360: so we can evaluate the method with different order even after the step has been optimistically completed.
361: */
362: static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done)
363: {
364: TS_RK *rk = (TS_RK*)ts->data;
365: RKTableau tab = rk->tableau;
366: PetscScalar *w = rk->work;
367: PetscReal h;
368: PetscInt s = tab->s,j;
372: switch (rk->status) {
373: case TS_STEP_INCOMPLETE:
374: case TS_STEP_PENDING:
375: h = ts->time_step; break;
376: case TS_STEP_COMPLETE:
377: h = ts->time_step_prev; break;
378: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
379: }
380: if (order == tab->order) {
381: if (rk->status == TS_STEP_INCOMPLETE) {
382: VecCopy(ts->vec_sol,X);
383: for (j=0; j<s; j++) w[j] = h*tab->b[j];
384: VecMAXPY(X,s,w,rk->YdotRHS);
385: } else {VecCopy(ts->vec_sol,X);}
386: return(0);
387: } else if (order == tab->order-1) {
388: if (!tab->bembed) goto unavailable;
389: if (rk->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (be) */
390: VecCopy(ts->vec_sol,X);
391: for (j=0; j<s; j++) w[j] = h*tab->bembed[j];
392: VecMAXPY(X,s,w,rk->YdotRHS);
393: } else { /* Rollback and re-complete using (be-b) */
394: VecCopy(ts->vec_sol,X);
395: for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]);
396: VecMAXPY(X,s,w,rk->YdotRHS);
397: }
398: if (done) *done = PETSC_TRUE;
399: return(0);
400: }
401: unavailable:
402: if (done) *done = PETSC_FALSE;
403: else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"RK '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order);
404: return(0);
405: }
409: static PetscErrorCode TSStep_RK(TS ts)
410: {
411: TS_RK *rk = (TS_RK*)ts->data;
412: RKTableau tab = rk->tableau;
413: const PetscInt s = tab->s;
414: const PetscReal *A = tab->A,*b = tab->b,*c = tab->c;
415: PetscScalar *w = rk->work;
416: Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS;
417: TSAdapt adapt;
418: PetscInt i,j,reject,next_scheme;
419: PetscReal next_time_step;
420: PetscReal t;
421: PetscBool accept;
422: PetscErrorCode ierr;
426: next_time_step = ts->time_step;
427: t = ts->ptime;
428: accept = PETSC_TRUE;
429: rk->status = TS_STEP_INCOMPLETE;
432: for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) {
433: PetscReal h = ts->time_step;
434: TSPreStep(ts);
435: for (i=0; i<s; i++) {
436: rk->stage_time = t + h*c[i];
437: TSPreStage(ts,rk->stage_time);
438: VecCopy(ts->vec_sol,Y[i]);
439: for (j=0; j<i; j++) w[j] = h*A[i*s+j];
440: VecMAXPY(Y[i],i,w,YdotRHS);
441: TSPostStage(ts,rk->stage_time,i,Y);
442: TSGetAdapt(ts,&adapt);
443: TSAdaptCheckStage(adapt,ts,&accept);
444: if (!accept) goto reject_step;
445: TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);
446: }
447: TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);
448: rk->status = TS_STEP_PENDING;
450: /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */
451: TSGetAdapt(ts,&adapt);
452: TSAdaptCandidatesClear(adapt);
453: TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);
454: TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);
455: if (accept) {
456: if (ts->costintegralfwd) {
457: /* Evolve ts->vec_costintegral to compute integrals */
458: for (i=0; i<s; i++) {
459: TSAdjointComputeCostIntegrand(ts,t+h*c[i],Y[i],ts->vec_costintegrand);
460: VecAXPY(ts->vec_costintegral,h*b[i],ts->vec_costintegrand);
461: }
462: }
464: /* ignore next_scheme for now */
465: ts->ptime += ts->time_step;
466: ts->time_step = next_time_step;
467: ts->steps++;
468: rk->status = TS_STEP_COMPLETE;
469: PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);
470: break;
471: } else { /* Roll back the current step */
472: for (j=0; j<s; j++) w[j] = -h*b[j];
473: VecMAXPY(ts->vec_sol,s,w,rk->YdotRHS);
474: ts->time_step = next_time_step;
475: rk->status = TS_STEP_INCOMPLETE;
476: }
477: reject_step: continue;
478: }
479: if (rk->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
480: return(0);
481: }
485: static PetscErrorCode TSAdjointSetUp_RK(TS ts)
486: {
487: TS_RK *rk = (TS_RK*)ts->data;
488: RKTableau tab;
489: PetscInt s;
493: if (ts->adjointsetupcalled++) return(0);
494: tab = rk->tableau;
495: s = tab->s;
496: VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecDeltaLam);
497: if(ts->vecs_sensip) {
498: VecDuplicateVecs(ts->vecs_sensip[0],s*ts->numcost,&rk->VecDeltaMu);
499: }
500: VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&rk->VecSensiTemp);
501: return(0);
502: }
506: static PetscErrorCode TSAdjointStep_RK(TS ts)
507: {
508: TS_RK *rk = (TS_RK*)ts->data;
509: RKTableau tab = rk->tableau;
510: const PetscInt s = tab->s;
511: const PetscReal *A = tab->A,*b = tab->b,*c = tab->c;
512: PetscScalar *w = rk->work;
513: Vec *Y = rk->Y,*VecDeltaLam = rk->VecDeltaLam,*VecDeltaMu = rk->VecDeltaMu,*VecSensiTemp = rk->VecSensiTemp;
514: PetscInt i,j,nadj;
515: PetscReal t;
516: PetscErrorCode ierr;
517: PetscReal h = ts->time_step;
518: Mat J,Jp;
521: t = ts->ptime;
522: rk->status = TS_STEP_INCOMPLETE;
523: h = ts->time_step;
524: TSPreStep(ts);
525: for (i=s-1; i>=0; i--) {
526: rk->stage_time = t + h*(1.0-c[i]);
527: for (nadj=0; nadj<ts->numcost; nadj++) {
528: VecCopy(ts->vecs_sensi[nadj],VecSensiTemp[nadj]);
529: VecScale(VecSensiTemp[nadj],-h*b[i]);
530: for (j=i+1; j<s; j++) {
531: VecAXPY(VecSensiTemp[nadj],-h*A[j*s+i],VecDeltaLam[nadj*s+j]);
532: }
533: }
534: /* Stage values of lambda */
535: TSGetRHSJacobian(ts,&J,&Jp,NULL,NULL);
536: TSComputeRHSJacobian(ts,rk->stage_time,Y[i],J,Jp);
537: if (ts->vec_costintegral) {
538: TSAdjointComputeDRDYFunction(ts,rk->stage_time,Y[i],ts->vecs_drdy);
539: if (!ts->costintegralfwd) {
540: /* Evolve ts->vec_costintegral to compute integrals */
541: TSAdjointComputeCostIntegrand(ts,rk->stage_time,Y[i],ts->vec_costintegrand);
542: VecAXPY(ts->vec_costintegral,-h*b[i],ts->vec_costintegrand);
543: }
544: }
545: for (nadj=0; nadj<ts->numcost; nadj++) {
546: MatMultTranspose(J,VecSensiTemp[nadj],VecDeltaLam[nadj*s+i]);
547: if (ts->vec_costintegral) {
548: VecAXPY(VecDeltaLam[nadj*s+i],-h*b[i],ts->vecs_drdy[nadj]);
549: }
550: }
552: /* Stage values of mu */
553: if(ts->vecs_sensip) {
554: TSAdjointComputeRHSJacobian(ts,rk->stage_time,Y[i],ts->Jacp);
555: if (ts->vec_costintegral) {
556: TSAdjointComputeDRDPFunction(ts,rk->stage_time,Y[i],ts->vecs_drdp);
557: }
559: for (nadj=0; nadj<ts->numcost; nadj++) {
560: MatMultTranspose(ts->Jacp,VecSensiTemp[nadj],VecDeltaMu[nadj*s+i]);
561: if (ts->vec_costintegral) {
562: VecAXPY(VecDeltaMu[nadj*s+i],-h*b[i],ts->vecs_drdp[nadj]);
563: }
564: }
565: }
566: }
568: for (j=0; j<s; j++) w[j] = 1.0;
569: for (nadj=0; nadj<ts->numcost; nadj++) {
570: VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecDeltaLam[nadj*s]);
571: if(ts->vecs_sensip) {
572: VecMAXPY(ts->vecs_sensip[nadj],s,w,&VecDeltaMu[nadj*s]);
573: }
574: }
575: ts->ptime += ts->time_step;
576: ts->steps++;
577: rk->status = TS_STEP_COMPLETE;
578: return(0);
579: }
583: static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X)
584: {
585: TS_RK *rk = (TS_RK*)ts->data;
586: PetscInt s = rk->tableau->s,pinterp = rk->tableau->pinterp,i,j;
587: PetscReal h;
588: PetscReal tt,t;
589: PetscScalar *b;
590: const PetscReal *B = rk->tableau->binterp;
591: PetscErrorCode ierr;
594: if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name);
595: switch (rk->status) {
596: case TS_STEP_INCOMPLETE:
597: case TS_STEP_PENDING:
598: h = ts->time_step;
599: t = (itime - ts->ptime)/h;
600: break;
601: case TS_STEP_COMPLETE:
602: h = ts->time_step_prev;
603: t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
604: break;
605: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
606: }
607: PetscMalloc1(s,&b);
608: for (i=0; i<s; i++) b[i] = 0;
609: for (j=0,tt=t; j<pinterp; j++,tt*=t) {
610: for (i=0; i<s; i++) {
611: b[i] += h * B[i*pinterp+j] * tt;
612: }
613: }
614: VecCopy(rk->Y[0],X);
615: VecMAXPY(X,s,b,rk->YdotRHS);
616: PetscFree(b);
617: return(0);
618: }
620: /*------------------------------------------------------------*/
623: static PetscErrorCode TSReset_RK(TS ts)
624: {
625: TS_RK *rk = (TS_RK*)ts->data;
626: PetscInt s;
630: if (!rk->tableau) return(0);
631: s = rk->tableau->s;
632: VecDestroyVecs(s,&rk->Y);
633: VecDestroyVecs(s,&rk->YdotRHS);
634: VecDestroyVecs(s*ts->numcost,&rk->VecDeltaLam);
635: VecDestroyVecs(s*ts->numcost,&rk->VecDeltaMu);
636: VecDestroyVecs(ts->numcost,&rk->VecSensiTemp);
637: PetscFree(rk->work);
638: return(0);
639: }
643: static PetscErrorCode TSDestroy_RK(TS ts)
644: {
648: TSReset_RK(ts);
649: PetscFree(ts->data);
650: PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);
651: PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);
652: return(0);
653: }
658: static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx)
659: {
661: return(0);
662: }
666: static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
667: {
669: return(0);
670: }
675: static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx)
676: {
678: return(0);
679: }
683: static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
684: {
687: return(0);
688: }
689: /*
692: static PetscErrorCode RKSetAdjCoe(RKTableau tab)
693: {
694: PetscReal *A,*b;
695: PetscInt s,i,j;
696: PetscErrorCode ierr;
699: s = tab->s;
700: PetscMalloc2(s*s,&A,s,&b);
702: for (i=0; i<s; i++)
703: for (j=0; j<s; j++) {
704: 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];
705: PetscPrintf(PETSC_COMM_WORLD,"Coefficients: A[%D][%D]=%.6f\n",i,j,A[i*s+j]);
706: }
708: for (i=0; i<s; i++) b[i] = (tab->b[s-1-i]==0)? 0: -tab->b[s-1-i];
710: PetscMemcpy(tab->A,A,s*s*sizeof(A[0]));
711: PetscMemcpy(tab->b,b,s*sizeof(b[0]));
712: PetscFree2(A,b);
713: return(0);
714: }
715: */
718: static PetscErrorCode TSSetUp_RK(TS ts)
719: {
720: TS_RK *rk = (TS_RK*)ts->data;
721: RKTableau tab;
722: PetscInt s;
724: DM dm;
727: if (!rk->tableau) {
728: TSRKSetType(ts,TSRKDefault);
729: }
730: tab = rk->tableau;
731: s = tab->s;
732: VecDuplicateVecs(ts->vec_sol,s,&rk->Y);
733: VecDuplicateVecs(ts->vec_sol,s,&rk->YdotRHS);
734: PetscMalloc1(s,&rk->work);
735: TSGetDM(ts,&dm);
736: if (dm) {
737: DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);
738: DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);
739: }
740: return(0);
741: }
744: /*------------------------------------------------------------*/
748: static PetscErrorCode TSSetFromOptions_RK(PetscOptions *PetscOptionsObject,TS ts)
749: {
751: char rktype[256];
754: PetscOptionsHead(PetscOptionsObject,"RK ODE solver options");
755: {
756: RKTableauLink link;
757: PetscInt count,choice;
758: PetscBool flg;
759: const char **namelist;
760: PetscStrncpy(rktype,TSRKDefault,sizeof(rktype));
761: for (link=RKTableauList,count=0; link; link=link->next,count++) ;
762: PetscMalloc1(count,&namelist);
763: for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
764: PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rktype,&choice,&flg);
765: TSRKSetType(ts,flg ? namelist[choice] : rktype);
766: PetscFree(namelist);
767: }
768: PetscOptionsTail();
769: return(0);
770: }
774: static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[])
775: {
777: PetscInt i;
778: size_t left,count;
779: char *p;
782: for (i=0,p=buf,left=len; i<n; i++) {
783: PetscSNPrintfCount(p,left,fmt,&count,x[i]);
784: if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer");
785: left -= count;
786: p += count;
787: *p++ = ' ';
788: }
789: p[i ? 0 : -1] = 0;
790: return(0);
791: }
795: static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer)
796: {
797: TS_RK *rk = (TS_RK*)ts->data;
798: RKTableau tab = rk->tableau;
799: PetscBool iascii;
801: TSAdapt adapt;
804: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
805: if (iascii) {
806: TSRKType rktype;
807: char buf[512];
808: TSRKGetType(ts,&rktype);
809: PetscViewerASCIIPrintf(viewer," RK %s\n",rktype);
810: PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);
811: PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);
812: PetscViewerASCIIPrintf(viewer,"FSAL: %s\n",tab->FSAL ? "yes" : "no");
813: }
814: TSGetAdapt(ts,&adapt);
815: TSAdaptView(adapt,viewer);
816: return(0);
817: }
821: static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer)
822: {
824: TSAdapt tsadapt;
827: TSGetAdapt(ts,&tsadapt);
828: TSAdaptLoad(tsadapt,viewer);
829: return(0);
830: }
834: /*@C
835: TSRKSetType - Set the type of RK scheme
837: Logically collective
839: Input Parameter:
840: + ts - timestepping context
841: - rktype - type of RK-scheme
843: Level: intermediate
845: .seealso: TSRKGetType(), TSRK, TSRK2, TSRK3, TSRKPRSSP2, TSRK4, TSRK5
846: @*/
847: PetscErrorCode TSRKSetType(TS ts,TSRKType rktype)
848: {
853: PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));
854: return(0);
855: }
859: /*@C
860: TSRKGetType - Get the type of RK scheme
862: Logically collective
864: Input Parameter:
865: . ts - timestepping context
867: Output Parameter:
868: . rktype - type of RK-scheme
870: Level: intermediate
872: .seealso: TSRKGetType()
873: @*/
874: PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype)
875: {
880: PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));
881: return(0);
882: }
886: PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype)
887: {
888: TS_RK *rk = (TS_RK*)ts->data;
892: if (!rk->tableau) {
893: TSRKSetType(ts,TSRKDefault);
894: }
895: *rktype = rk->tableau->name;
896: return(0);
897: }
900: PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype)
901: {
902: TS_RK *rk = (TS_RK*)ts->data;
904: PetscBool match;
905: RKTableauLink link;
908: if (rk->tableau) {
909: PetscStrcmp(rk->tableau->name,rktype,&match);
910: if (match) return(0);
911: }
912: for (link = RKTableauList; link; link=link->next) {
913: PetscStrcmp(link->tab.name,rktype,&match);
914: if (match) {
915: TSReset_RK(ts);
916: rk->tableau = &link->tab;
917: return(0);
918: }
919: }
920: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype);
921: return(0);
922: }
926: static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y)
927: {
928: TS_RK *rk = (TS_RK*)ts->data;
931: *ns = rk->tableau->s;
932: if(Y) *Y = rk->Y;
933: return(0);
934: }
937: /* ------------------------------------------------------------ */
938: /*MC
939: TSRK - ODE and DAE solver using Runge-Kutta schemes
941: The user should provide the right hand side of the equation
942: using TSSetRHSFunction().
944: Notes:
945: The default is TSRK3, it can be changed with TSRKSetType() or -ts_rk_type
947: Level: beginner
949: .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3,
950: TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister()
952: M*/
955: PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts)
956: {
957: TS_RK *th;
961: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
962: TSRKInitializePackage();
963: #endif
965: ts->ops->reset = TSReset_RK;
966: ts->ops->destroy = TSDestroy_RK;
967: ts->ops->view = TSView_RK;
968: ts->ops->load = TSLoad_RK;
969: ts->ops->setup = TSSetUp_RK;
970: ts->ops->adjointsetup = TSAdjointSetUp_RK;
971: ts->ops->step = TSStep_RK;
972: ts->ops->interpolate = TSInterpolate_RK;
973: ts->ops->evaluatestep = TSEvaluateStep_RK;
974: ts->ops->setfromoptions = TSSetFromOptions_RK;
975: ts->ops->getstages = TSGetStages_RK;
976: ts->ops->adjointstep = TSAdjointStep_RK;
978: PetscNewLog(ts,&th);
979: ts->data = (void*)th;
981: PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);
982: PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);
983: return(0);
984: }