Actual source code: rk.c
petsc-3.5.4 2015-05-23
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: PetscScalar *work; /* Scalar work */
42: PetscReal stage_time;
43: TSStepStatus status;
44: } TS_RK;
46: /*MC
47: TSRK1 - First order forward Euler scheme.
49: This method has one stage.
51: Level: advanced
53: .seealso: TSRK
54: M*/
55: /*MC
56: TSRK2A - Second order RK scheme.
58: This method has two stages.
60: Level: advanced
62: .seealso: TSRK
63: M*/
64: /*MC
65: TSRK3 - Third order RK scheme.
67: This method has three stages.
69: Level: advanced
71: .seealso: TSRK
72: M*/
73: /*MC
74: TSRK3BS - Third order RK scheme of Bogacki-Shampine with 2nd order embedded method.
76: This method has four stages.
78: Level: advanced
80: .seealso: TSRK
81: M*/
82: /*MC
83: TSRK4 - Fourth order RK scheme.
85: This is the classical Runge-Kutta method with four stages.
87: Level: advanced
89: .seealso: TSRK
90: M*/
91: /*MC
92: TSRK5F - Fifth order Fehlberg RK scheme with a 4th order embedded method.
94: This method has six stages.
96: Level: advanced
98: .seealso: TSRK
99: M*/
100: /*MC
101: TSRK5DP - Fifth order Dormand-Prince RK scheme with the 4th order embedded method.
103: This method has seven stages.
105: Level: advanced
107: .seealso: TSRK
108: M*/
112: /*@C
113: TSRKRegisterAll - Registers all of the Runge-Kutta explicit methods in TSRK
115: Not Collective, but should be called by all processes which will need the schemes to be registered
117: Level: advanced
119: .keywords: TS, TSRK, register, all
121: .seealso: TSRKRegisterDestroy()
122: @*/
123: PetscErrorCode TSRKRegisterAll(void)
124: {
128: if (TSRKRegisterAllCalled) return(0);
129: TSRKRegisterAllCalled = PETSC_TRUE;
131: {
132: const PetscReal
133: A[1][1] = {{0.0}},
134: b[1] = {1.0};
135: TSRKRegister(TSRK1FE,1,1,&A[0][0],b,NULL,NULL,1,b);
136: }
137: {
138: const PetscReal
139: A[2][2] = {{0.0,0.0},
140: {1.0,0.0}},
141: b[2] = {0.5,0.5},
142: bembed[2] = {1.0,0};
143: TSRKRegister(TSRK2A,2,2,&A[0][0],b,NULL,bembed,2,b);
144: }
145: {
146: const PetscReal
147: A[3][3] = {{0,0,0},
148: {2.0/3.0,0,0},
149: {-1.0/3.0,1.0,0}},
150: b[3] = {0.25,0.5,0.25};
151: TSRKRegister(TSRK3,3,3,&A[0][0],b,NULL,NULL,3,b);
152: }
153: {
154: const PetscReal
155: A[4][4] = {{0,0,0,0},
156: {1.0/2.0,0,0,0},
157: {0,3.0/4.0,0,0},
158: {2.0/9.0,1.0/3.0,4.0/9.0,0}},
159: b[4] = {2.0/9.0,1.0/3.0,4.0/9.0,0},
160: bembed[4] = {7.0/24.0,1.0/4.0,1.0/3.0,1.0/8.0};
161: TSRKRegister(TSRK3BS,3,4,&A[0][0],b,NULL,bembed,3,b);
162: }
163: {
164: const PetscReal
165: A[4][4] = {{0,0,0,0},
166: {0.5,0,0,0},
167: {0,0.5,0,0},
168: {0,0,1.0,0}},
169: b[4] = {1.0/6.0,1.0/3.0,1.0/3.0,1.0/6.0};
170: TSRKRegister(TSRK4,4,4,&A[0][0],b,NULL,NULL,4,b);
171: }
172: {
173: const PetscReal
174: A[6][6] = {{0,0,0,0,0,0},
175: {0.25,0,0,0,0,0},
176: {3.0/32.0,9.0/32.0,0,0,0,0},
177: {1932.0/2197.0,-7200.0/2197.0,7296.0/2197.0,0,0,0},
178: {439.0/216.0,-8.0,3680.0/513.0,-845.0/4104.0,0,0},
179: {-8.0/27.0,2.0,-3544.0/2565.0,1859.0/4104.0,-11.0/40.0,0}},
180: b[6] = {16.0/135.0,0,6656.0/12825.0,28561.0/56430.0,-9.0/50.0,2.0/55.0},
181: bembed[6] = {25.0/216.0,0,1408.0/2565.0,2197.0/4104.0,-1.0/5.0,0};
182: TSRKRegister(TSRK5F,5,6,&A[0][0],b,NULL,bembed,5,b);
183: }
184: {
185: const PetscReal
186: A[7][7] = {{0,0,0,0,0,0,0},
187: {1.0/5.0,0,0,0,0,0,0},
188: {3.0/40.0,9.0/40.0,0,0,0,0,0},
189: {44.0/45.0,-56.0/15.0,32.0/9.0,0,0,0,0},
190: {19372.0/6561.0,-25360.0/2187.0,64448.0/6561.0,-212.0/729.0,0,0,0},
191: {9017.0/3168.0,-355.0/33.0,46732.0/5247.0,49.0/176.0,-5103.0/18656.0,0,0},
192: {35.0/384.0,0,500.0/1113.0,125.0/192.0,-2187.0/6784.0,11.0/84.0,0}},
193: b[7] = {35.0/384.0,0,500.0/1113.0,125.0/192.0,-2187.0/6784.0,11.0/84.0,0},
194: 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};
195: TSRKRegister(TSRK5DP,5,7,&A[0][0],b,NULL,bembed,5,b);
196: }
197: return(0);
198: }
202: /*@C
203: TSRKRegisterDestroy - Frees the list of schemes that were registered by TSRKRegister().
205: Not Collective
207: Level: advanced
209: .keywords: TSRK, register, destroy
210: .seealso: TSRKRegister(), TSRKRegisterAll()
211: @*/
212: PetscErrorCode TSRKRegisterDestroy(void)
213: {
215: RKTableauLink link;
218: while ((link = RKTableauList)) {
219: RKTableau t = &link->tab;
220: RKTableauList = link->next;
221: PetscFree3(t->A,t->b,t->c);
222: PetscFree (t->bembed);
223: PetscFree (t->binterp);
224: PetscFree (t->name);
225: PetscFree (link);
226: }
227: TSRKRegisterAllCalled = PETSC_FALSE;
228: return(0);
229: }
233: /*@C
234: TSRKInitializePackage - This function initializes everything in the TSRK package. It is called
235: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RK()
236: when using static libraries.
238: Level: developer
240: .keywords: TS, TSRK, initialize, package
241: .seealso: PetscInitialize()
242: @*/
243: PetscErrorCode TSRKInitializePackage(void)
244: {
248: if (TSRKPackageInitialized) return(0);
249: TSRKPackageInitialized = PETSC_TRUE;
250: TSRKRegisterAll();
251: PetscObjectComposedDataRegister(&explicit_stage_time_id);
252: PetscRegisterFinalize(TSRKFinalizePackage);
253: return(0);
254: }
258: /*@C
259: TSRKFinalizePackage - This function destroys everything in the TSRK package. It is
260: called from PetscFinalize().
262: Level: developer
264: .keywords: Petsc, destroy, package
265: .seealso: PetscFinalize()
266: @*/
267: PetscErrorCode TSRKFinalizePackage(void)
268: {
272: TSRKPackageInitialized = PETSC_FALSE;
273: TSRKRegisterDestroy();
274: return(0);
275: }
279: /*@C
280: TSRKRegister - register an RK scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation
282: Not Collective, but the same schemes should be registered on all processes on which they will be used
284: Input Parameters:
285: + name - identifier for method
286: . order - approximation order of method
287: . s - number of stages, this is the dimension of the matrices below
288: . A - stage coefficients (dimension s*s, row-major)
289: . b - step completion table (dimension s; NULL to use last row of A)
290: . c - abscissa (dimension s; NULL to use row sums of A)
291: . bembed - completion table for embedded method (dimension s; NULL if not available)
292: . pinterp - Order of the interpolation scheme, equal to the number of columns of binterp
293: - binterp - Coefficients of the interpolation formula (dimension s*pinterp; NULL to reuse binterpt)
295: Notes:
296: Several RK methods are provided, this function is only needed to create new methods.
298: Level: advanced
300: .keywords: TS, register
302: .seealso: TSRK
303: @*/
304: PetscErrorCode TSRKRegister(TSRKType name,PetscInt order,PetscInt s,
305: const PetscReal A[],const PetscReal b[],const PetscReal c[],
306: const PetscReal bembed[],
307: PetscInt pinterp,const PetscReal binterp[])
308: {
309: PetscErrorCode ierr;
310: RKTableauLink link;
311: RKTableau t;
312: PetscInt i,j;
315: PetscMalloc(sizeof(*link),&link);
316: PetscMemzero(link,sizeof(*link));
317: t = &link->tab;
318: PetscStrallocpy(name,&t->name);
319: t->order = order;
320: t->s = s;
321: PetscMalloc3(s*s,&t->A,s,&t->b,s,&t->c);
322: PetscMemcpy(t->A,A,s*s*sizeof(A[0]));
323: if (b) { PetscMemcpy(t->b,b,s*sizeof(b[0])); }
324: else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i];
325: if (c) { PetscMemcpy(t->c,c,s*sizeof(c[0])); }
326: else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j];
327: t->FSAL = PETSC_TRUE;
328: for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE;
329: if (bembed) {
330: PetscMalloc1(s,&t->bembed);
331: PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));
332: }
334: t->pinterp = pinterp;
335: PetscMalloc1(s*pinterp,&t->binterp);
336: PetscMemcpy(t->binterp,binterp,s*pinterp*sizeof(binterp[0]));
337: link->next = RKTableauList;
338: RKTableauList = link;
339: return(0);
340: }
344: /*
345: The step completion formula is
347: x1 = x0 + h b^T YdotRHS
349: This function can be called before or after ts->vec_sol has been updated.
350: Suppose we have a completion formula (b) and an embedded formula (be) of different order.
351: We can write
353: x1e = x0 + h be^T YdotRHS
354: = x1 - h b^T YdotRHS + h be^T YdotRHS
355: = x1 + h (be - b)^T YdotRHS
357: so we can evaluate the method with different order even after the step has been optimistically completed.
358: */
359: static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done)
360: {
361: TS_RK *rk = (TS_RK*)ts->data;
362: RKTableau tab = rk->tableau;
363: PetscScalar *w = rk->work;
364: PetscReal h;
365: PetscInt s = tab->s,j;
369: switch (rk->status) {
370: case TS_STEP_INCOMPLETE:
371: case TS_STEP_PENDING:
372: h = ts->time_step; break;
373: case TS_STEP_COMPLETE:
374: h = ts->time_step_prev; break;
375: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
376: }
377: if (order == tab->order) {
378: if (rk->status == TS_STEP_INCOMPLETE) {
379: VecCopy(ts->vec_sol,X);
380: for (j=0; j<s; j++) w[j] = h*tab->b[j];
381: VecMAXPY(X,s,w,rk->YdotRHS);
382: } else {VecCopy(ts->vec_sol,X);}
383: if (done) *done = PETSC_TRUE;
384: return(0);
385: } else if (order == tab->order-1) {
386: if (!tab->bembed) goto unavailable;
387: if (rk->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (be) */
388: VecCopy(ts->vec_sol,X);
389: for (j=0; j<s; j++) w[j] = h*tab->bembed[j];
390: VecMAXPY(X,s,w,rk->YdotRHS);
391: } else { /* Rollback and re-complete using (be-b) */
392: VecCopy(ts->vec_sol,X);
393: for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]);
394: VecMAXPY(X,s,w,rk->YdotRHS);
395: }
396: if (done) *done = PETSC_TRUE;
397: return(0);
398: }
399: unavailable:
400: if (done) *done = PETSC_FALSE;
401: else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"RK '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order);
402: return(0);
403: }
407: static PetscErrorCode TSStep_RK(TS ts)
408: {
409: TS_RK *rk = (TS_RK*)ts->data;
410: RKTableau tab = rk->tableau;
411: const PetscInt s = tab->s;
412: const PetscReal *A = tab->A,*b = tab->b,*c = tab->c;
413: PetscScalar *w = rk->work;
414: Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS;
415: TSAdapt adapt;
416: PetscInt i,j,reject,next_scheme;
417: PetscReal next_time_step;
418: PetscReal t;
419: PetscBool accept;
420: PetscErrorCode ierr;
424: next_time_step = ts->time_step;
425: t = ts->ptime;
426: accept = PETSC_TRUE;
427: rk->status = TS_STEP_INCOMPLETE;
430: for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) {
431: PetscReal h = ts->time_step;
432: TSPreStep(ts);
433: for (i=0; i<s; i++) {
434: rk->stage_time = t + h*c[i];
435: TSPreStage(ts,rk->stage_time);
436: VecCopy(ts->vec_sol,Y[i]);
437: for (j=0; j<i; j++) w[j] = h*A[i*s+j];
438: VecMAXPY(Y[i],i,w,YdotRHS);
439: TSPostStage(ts,rk->stage_time,i,Y);
440: TSGetAdapt(ts,&adapt);
441: TSAdaptCheckStage(adapt,ts,&accept);
442: if (!accept) goto reject_step;
443: TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);
444: }
445: TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);
446: rk->status = TS_STEP_PENDING;
448: /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */
449: TSGetAdapt(ts,&adapt);
450: TSAdaptCandidatesClear(adapt);
451: TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);
452: TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);
453: if (accept) {
454: /* ignore next_scheme for now */
455: ts->ptime += ts->time_step;
456: ts->time_step = next_time_step;
457: ts->steps++;
458: rk->status = TS_STEP_COMPLETE;
459: PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);
460: break;
461: } else { /* Roll back the current step */
462: for (j=0; j<s; j++) w[j] = -h*b[j];
463: VecMAXPY(ts->vec_sol,s,w,rk->YdotRHS);
464: ts->time_step = next_time_step;
465: rk->status = TS_STEP_INCOMPLETE;
466: }
467: reject_step: continue;
468: }
469: if (rk->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
470: return(0);
471: }
475: static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X)
476: {
477: TS_RK *rk = (TS_RK*)ts->data;
478: PetscInt s = rk->tableau->s,pinterp = rk->tableau->pinterp,i,j;
479: PetscReal h;
480: PetscReal tt,t;
481: PetscScalar *b;
482: const PetscReal *B = rk->tableau->binterp;
483: PetscErrorCode ierr;
486: if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name);
487: switch (rk->status) {
488: case TS_STEP_INCOMPLETE:
489: case TS_STEP_PENDING:
490: h = ts->time_step;
491: t = (itime - ts->ptime)/h;
492: break;
493: case TS_STEP_COMPLETE:
494: h = ts->time_step_prev;
495: t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
496: break;
497: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
498: }
499: PetscMalloc1(s,&b);
500: for (i=0; i<s; i++) b[i] = 0;
501: for (j=0,tt=t; j<pinterp; j++,tt*=t) {
502: for (i=0; i<s; i++) {
503: b[i] += h * B[i*pinterp+j] * tt;
504: }
505: }
506: VecCopy(rk->Y[0],X);
507: VecMAXPY(X,s,b,rk->YdotRHS);
508: PetscFree(b);
509: return(0);
510: }
512: /*------------------------------------------------------------*/
515: static PetscErrorCode TSReset_RK(TS ts)
516: {
517: TS_RK *rk = (TS_RK*)ts->data;
518: PetscInt s;
522: if (!rk->tableau) return(0);
523: s = rk->tableau->s;
524: VecDestroyVecs(s,&rk->Y);
525: VecDestroyVecs(s,&rk->YdotRHS);
526: PetscFree(rk->work);
527: return(0);
528: }
532: static PetscErrorCode TSDestroy_RK(TS ts)
533: {
537: TSReset_RK(ts);
538: PetscFree(ts->data);
539: PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);
540: PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);
541: return(0);
542: }
547: static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx)
548: {
550: return(0);
551: }
555: static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
556: {
558: return(0);
559: }
564: static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx)
565: {
567: return(0);
568: }
572: static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
573: {
576: return(0);
577: }
581: static PetscErrorCode TSSetUp_RK(TS ts)
582: {
583: TS_RK *rk = (TS_RK*)ts->data;
584: RKTableau tab;
585: PetscInt s;
587: DM dm;
590: if (!rk->tableau) {
591: TSRKSetType(ts,TSRKDefault);
592: }
593: tab = rk->tableau;
594: s = tab->s;
595: VecDuplicateVecs(ts->vec_sol,s,&rk->Y);
596: VecDuplicateVecs(ts->vec_sol,s,&rk->YdotRHS);
597: PetscMalloc1(s,&rk->work);
598: TSGetDM(ts,&dm);
599: if (dm) {
600: DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);
601: DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);
602: }
603: return(0);
604: }
605: /*------------------------------------------------------------*/
609: static PetscErrorCode TSSetFromOptions_RK(TS ts)
610: {
612: char rktype[256];
615: PetscOptionsHead("RK ODE solver options");
616: {
617: RKTableauLink link;
618: PetscInt count,choice;
619: PetscBool flg;
620: const char **namelist;
621: PetscStrncpy(rktype,TSRKDefault,sizeof(rktype));
622: for (link=RKTableauList,count=0; link; link=link->next,count++) ;
623: PetscMalloc1(count,&namelist);
624: for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
625: PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rktype,&choice,&flg);
626: TSRKSetType(ts,flg ? namelist[choice] : rktype);
627: PetscFree(namelist);
628: }
629: PetscOptionsTail();
630: return(0);
631: }
635: static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[])
636: {
638: PetscInt i;
639: size_t left,count;
640: char *p;
643: for (i=0,p=buf,left=len; i<n; i++) {
644: PetscSNPrintfCount(p,left,fmt,&count,x[i]);
645: if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer");
646: left -= count;
647: p += count;
648: *p++ = ' ';
649: }
650: p[i ? 0 : -1] = 0;
651: return(0);
652: }
656: static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer)
657: {
658: TS_RK *rk = (TS_RK*)ts->data;
659: RKTableau tab = rk->tableau;
660: PetscBool iascii;
662: TSAdapt adapt;
665: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
666: if (iascii) {
667: TSRKType rktype;
668: char buf[512];
669: TSRKGetType(ts,&rktype);
670: PetscViewerASCIIPrintf(viewer," RK %s\n",rktype);
671: PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);
672: PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);
673: PetscViewerASCIIPrintf(viewer,"FSAL: %s\n",tab->FSAL ? "yes" : "no");
674: }
675: TSGetAdapt(ts,&adapt);
676: TSAdaptView(adapt,viewer);
677: return(0);
678: }
682: static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer)
683: {
685: TSAdapt tsadapt;
688: TSGetAdapt(ts,&tsadapt);
689: TSAdaptLoad(tsadapt,viewer);
690: return(0);
691: }
695: /*@C
696: TSRKSetType - Set the type of RK scheme
698: Logically collective
700: Input Parameter:
701: + ts - timestepping context
702: - rktype - type of RK-scheme
704: Level: intermediate
706: .seealso: TSRKGetType(), TSRK, TSRK2, TSRK3, TSRKPRSSP2, TSRK4, TSRK5
707: @*/
708: PetscErrorCode TSRKSetType(TS ts,TSRKType rktype)
709: {
714: PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));
715: return(0);
716: }
720: /*@C
721: TSRKGetType - Get the type of RK scheme
723: Logically collective
725: Input Parameter:
726: . ts - timestepping context
728: Output Parameter:
729: . rktype - type of RK-scheme
731: Level: intermediate
733: .seealso: TSRKGetType()
734: @*/
735: PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype)
736: {
741: PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));
742: return(0);
743: }
747: PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype)
748: {
749: TS_RK *rk = (TS_RK*)ts->data;
753: if (!rk->tableau) {
754: TSRKSetType(ts,TSRKDefault);
755: }
756: *rktype = rk->tableau->name;
757: return(0);
758: }
761: PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype)
762: {
763: TS_RK *rk = (TS_RK*)ts->data;
765: PetscBool match;
766: RKTableauLink link;
769: if (rk->tableau) {
770: PetscStrcmp(rk->tableau->name,rktype,&match);
771: if (match) return(0);
772: }
773: for (link = RKTableauList; link; link=link->next) {
774: PetscStrcmp(link->tab.name,rktype,&match);
775: if (match) {
776: TSReset_RK(ts);
777: rk->tableau = &link->tab;
778: return(0);
779: }
780: }
781: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype);
782: return(0);
783: }
785: /* ------------------------------------------------------------ */
786: /*MC
787: TSRK - ODE and DAE solver using Runge-Kutta schemes
789: The user should provide the right hand side of the equation
790: using TSSetRHSFunction().
792: Notes:
793: The default is TSRK3, it can be changed with TSRKSetType() or -ts_rk_type
795: Level: beginner
797: .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3,
798: TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister()
800: M*/
803: PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts)
804: {
805: TS_RK *th;
809: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
810: TSRKInitializePackage();
811: #endif
813: ts->ops->reset = TSReset_RK;
814: ts->ops->destroy = TSDestroy_RK;
815: ts->ops->view = TSView_RK;
816: ts->ops->load = TSLoad_RK;
817: ts->ops->setup = TSSetUp_RK;
818: ts->ops->step = TSStep_RK;
819: ts->ops->interpolate = TSInterpolate_RK;
820: ts->ops->evaluatestep = TSEvaluateStep_RK;
821: ts->ops->setfromoptions = TSSetFromOptions_RK;
823: PetscNewLog(ts,&th);
824: ts->data = (void*)th;
826: PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);
827: PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);
828: return(0);
829: }