Actual source code: rk.c
petsc-3.3-p7 2013-05-11
1: /*
2: * Code for Timestepping with Runge Kutta
3: *
4: * Written by
5: * Asbjorn Hoiland Aarrestad
6: * asbjorn@aarrestad.com
7: * http://asbjorn.aarrestad.com/
8: *
9: */
10: #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
11: #include <time.h>
13: typedef struct {
14: Vec y1,y2; /* work wectors for the two rk permuations */
15: PetscInt nok,nnok; /* counters for ok and not ok steps */
16: PetscReal maxerror; /* variable to tell the maxerror allowed */
17: PetscReal ferror; /* variable to tell (global maxerror)/(total time) */
18: PetscReal tolerance; /* initial value set for maxerror by user */
19: Vec tmp,tmp_y,*k; /* two temp vectors and the k vectors for rk */
20: PetscScalar a[7][6]; /* rk scalars */
21: PetscScalar b1[7],b2[7]; /* rk scalars */
22: PetscReal c[7]; /* rk scalars */
23: PetscInt p,s; /* variables to tell the size of the runge-kutta solver */
24: } TS_RK;
26: EXTERN_C_BEGIN
29: PetscErrorCode TSRKSetTolerance_RK(TS ts,PetscReal aabs)
30: {
31: TS_RK *rk = (TS_RK*)ts->data;
34: rk->tolerance = aabs;
35: return(0);
36: }
37: EXTERN_C_END
41: /*@
42: TSRKSetTolerance - Sets the total error the RK explicit time integrators
43: will allow over the given time interval.
45: Logically Collective on TS
47: Input parameters:
48: + ts - the time-step context
49: - aabs - the absolute tolerance
51: Level: intermediate
53: .keywords: RK, tolerance
55: .seealso: TSSundialsSetTolerance()
57: @*/
58: PetscErrorCode TSRKSetTolerance(TS ts,PetscReal aabs)
59: {
64: PetscTryMethod(ts,"TSRKSetTolerance_C",(TS,PetscReal),(ts,aabs));
65: return(0);
66: }
71: static PetscErrorCode TSSetUp_RK(TS ts)
72: {
73: TS_RK *rk = (TS_RK*)ts->data;
77: rk->nok = 0;
78: rk->nnok = 0;
79: rk->maxerror = rk->tolerance;
81: /* fixing maxerror: global vs local */
82: rk->ferror = rk->maxerror / (ts->max_time - ts->ptime);
84: /* 34.0/45.0 gives double precision division */
85: /* defining variables needed for Runge-Kutta computing */
86: /* when changing below, please remember to change a, b1, b2 and c above! */
87: /* Found in table on page 171: Dormand-Prince 5(4) */
89: /* are these right? */
90: rk->p=6;
91: rk->s=7;
93: rk->a[1][0]=1.0/5.0;
94: rk->a[2][0]=3.0/40.0;
95: rk->a[2][1]=9.0/40.0;
96: rk->a[3][0]=44.0/45.0;
97: rk->a[3][1]=-56.0/15.0;
98: rk->a[3][2]=32.0/9.0;
99: rk->a[4][0]=19372.0/6561.0;
100: rk->a[4][1]=-25360.0/2187.0;
101: rk->a[4][2]=64448.0/6561.0;
102: rk->a[4][3]=-212.0/729.0;
103: rk->a[5][0]=9017.0/3168.0;
104: rk->a[5][1]=-355.0/33.0;
105: rk->a[5][2]=46732.0/5247.0;
106: rk->a[5][3]=49.0/176.0;
107: rk->a[5][4]=-5103.0/18656.0;
108: rk->a[6][0]=35.0/384.0;
109: rk->a[6][1]=0.0;
110: rk->a[6][2]=500.0/1113.0;
111: rk->a[6][3]=125.0/192.0;
112: rk->a[6][4]=-2187.0/6784.0;
113: rk->a[6][5]=11.0/84.0;
116: rk->c[0]=0.0;
117: rk->c[1]=1.0/5.0;
118: rk->c[2]=3.0/10.0;
119: rk->c[3]=4.0/5.0;
120: rk->c[4]=8.0/9.0;
121: rk->c[5]=1.0;
122: rk->c[6]=1.0;
124: rk->b1[0]=35.0/384.0;
125: rk->b1[1]=0.0;
126: rk->b1[2]=500.0/1113.0;
127: rk->b1[3]=125.0/192.0;
128: rk->b1[4]=-2187.0/6784.0;
129: rk->b1[5]=11.0/84.0;
130: rk->b1[6]=0.0;
132: rk->b2[0]=5179.0/57600.0;
133: rk->b2[1]=0.0;
134: rk->b2[2]=7571.0/16695.0;
135: rk->b2[3]=393.0/640.0;
136: rk->b2[4]=-92097.0/339200.0;
137: rk->b2[5]=187.0/2100.0;
138: rk->b2[6]=1.0/40.0;
141: /* Found in table on page 170: Fehlberg 4(5) */
142: /*
143: rk->p=5;
144: rk->s=6;
146: rk->a[1][0]=1.0/4.0;
147: rk->a[2][0]=3.0/32.0;
148: rk->a[2][1]=9.0/32.0;
149: rk->a[3][0]=1932.0/2197.0;
150: rk->a[3][1]=-7200.0/2197.0;
151: rk->a[3][2]=7296.0/2197.0;
152: rk->a[4][0]=439.0/216.0;
153: rk->a[4][1]=-8.0;
154: rk->a[4][2]=3680.0/513.0;
155: rk->a[4][3]=-845.0/4104.0;
156: rk->a[5][0]=-8.0/27.0;
157: rk->a[5][1]=2.0;
158: rk->a[5][2]=-3544.0/2565.0;
159: rk->a[5][3]=1859.0/4104.0;
160: rk->a[5][4]=-11.0/40.0;
162: rk->c[0]=0.0;
163: rk->c[1]=1.0/4.0;
164: rk->c[2]=3.0/8.0;
165: rk->c[3]=12.0/13.0;
166: rk->c[4]=1.0;
167: rk->c[5]=1.0/2.0;
169: rk->b1[0]=25.0/216.0;
170: rk->b1[1]=0.0;
171: rk->b1[2]=1408.0/2565.0;
172: rk->b1[3]=2197.0/4104.0;
173: rk->b1[4]=-1.0/5.0;
174: rk->b1[5]=0.0;
176: rk->b2[0]=16.0/135.0;
177: rk->b2[1]=0.0;
178: rk->b2[2]=6656.0/12825.0;
179: rk->b2[3]=28561.0/56430.0;
180: rk->b2[4]=-9.0/50.0;
181: rk->b2[5]=2.0/55.0;
182: */
183: /* Found in table on page 169: Merson 4("5") */
184: /*
185: rk->p=4;
186: rk->s=5;
187: rk->a[1][0] = 1.0/3.0;
188: rk->a[2][0] = 1.0/6.0;
189: rk->a[2][1] = 1.0/6.0;
190: rk->a[3][0] = 1.0/8.0;
191: rk->a[3][1] = 0.0;
192: rk->a[3][2] = 3.0/8.0;
193: rk->a[4][0] = 1.0/2.0;
194: rk->a[4][1] = 0.0;
195: rk->a[4][2] = -3.0/2.0;
196: rk->a[4][3] = 2.0;
198: rk->c[0] = 0.0;
199: rk->c[1] = 1.0/3.0;
200: rk->c[2] = 1.0/3.0;
201: rk->c[3] = 0.5;
202: rk->c[4] = 1.0;
204: rk->b1[0] = 1.0/2.0;
205: rk->b1[1] = 0.0;
206: rk->b1[2] = -3.0/2.0;
207: rk->b1[3] = 2.0;
208: rk->b1[4] = 0.0;
210: rk->b2[0] = 1.0/6.0;
211: rk->b2[1] = 0.0;
212: rk->b2[2] = 0.0;
213: rk->b2[3] = 2.0/3.0;
214: rk->b2[4] = 1.0/6.0;
215: */
217: /* making b2 -> e=b1-b2 */
218: /*
219: for(i=0;i<rk->s;i++){
220: rk->b2[i] = (rk->b1[i]) - (rk->b2[i]);
221: }
222: */
223: rk->b2[0]=71.0/57600.0;
224: rk->b2[1]=0.0;
225: rk->b2[2]=-71.0/16695.0;
226: rk->b2[3]=71.0/1920.0;
227: rk->b2[4]=-17253.0/339200.0;
228: rk->b2[5]=22.0/525.0;
229: rk->b2[6]=-1.0/40.0;
231: /* initializing vectors */
232: VecDuplicate(ts->vec_sol,&rk->y1);
233: VecDuplicate(ts->vec_sol,&rk->y2);
234: VecDuplicate(rk->y1,&rk->tmp);
235: VecDuplicate(rk->y1,&rk->tmp_y);
236: VecDuplicateVecs(rk->y1,rk->s,&rk->k);
238: return(0);
239: }
241: /*------------------------------------------------------------*/
244: PetscErrorCode TSRKqs(TS ts,PetscReal t,PetscReal h)
245: {
246: TS_RK *rk = (TS_RK*)ts->data;
248: PetscInt j,l;
249: PetscReal tmp_t=t;
250: PetscScalar hh=h;
253: /* k[0]=0 */
254: VecSet(rk->k[0],0.0);
256: /* k[0] = derivs(t,y1) */
257: TSComputeRHSFunction(ts,t,rk->y1,rk->k[0]);
258: /* looping over runge-kutta variables */
259: /* building the k - array of vectors */
260: for(j = 1 ; j < rk->s ; j++){
262: /* rk->tmp = 0 */
263: VecSet(rk->tmp,0.0);
265: for(l=0;l<j;l++){
266: /* tmp += a(j,l)*k[l] */
267: VecAXPY(rk->tmp,rk->a[j][l],rk->k[l]);
268: }
270: /* VecView(rk->tmp,PETSC_VIEWER_STDOUT_WORLD); */
272: /* k[j] = derivs(t+c(j)*h,y1+h*tmp,k(j)) */
273: /* I need the following helpers:
274: PetscScalar tmp_t=t+c(j)*h
275: Vec tmp_y=h*tmp+y1
276: */
278: tmp_t = t + rk->c[j] * h;
280: /* tmp_y = h * tmp + y1 */
281: VecWAXPY(rk->tmp_y,hh,rk->tmp,rk->y1);
283: /* rk->k[j]=0 */
284: VecSet(rk->k[j],0.0);
285: TSComputeRHSFunction(ts,tmp_t,rk->tmp_y,rk->k[j]);
286: }
288: /* tmp=0 and tmp_y=0 */
289: VecSet(rk->tmp,0.0);
290: VecSet(rk->tmp_y,0.0);
292: for(j = 0 ; j < rk->s ; j++){
293: /* tmp=b1[j]*k[j]+tmp */
294: VecAXPY(rk->tmp,rk->b1[j],rk->k[j]);
295: /* tmp_y=b2[j]*k[j]+tmp_y */
296: VecAXPY(rk->tmp_y,rk->b2[j],rk->k[j]);
297: }
299: /* y2 = hh * tmp_y */
300: VecSet(rk->y2,0.0);
301: VecAXPY(rk->y2,hh,rk->tmp_y);
302: /* y1 = hh*tmp + y1 */
303: VecAXPY(rk->y1,hh,rk->tmp);
304: /* Finding difference between y1 and y2 */
305: return(0);
306: }
310: static PetscErrorCode TSSolve_RK(TS ts)
311: {
312: TS_RK *rk = (TS_RK*)ts->data;
313: PetscReal norm=0.0,dt_fac=0.0,fac = 0.0/*,ttmp=0.0*/;
314: PetscInt i;
318: VecCopy(ts->vec_sol,rk->y1);
320: /* while loop to get from start to stop */
321: for (i = 0; i < ts->max_steps; i++) {
322: TSPreStep(ts); /* Note that this is called once per STEP, not once per STAGE. */
324: /* calling rkqs */
325: /*
326: -- input
327: ts - pointer to ts
328: ts->ptime - current time
329: ts->time_step - try this timestep
330: y1 - solution for this step
332: --output
333: y1 - suggested solution
334: y2 - check solution (runge - kutta second permutation)
335: */
336: TSRKqs(ts,ts->ptime,ts->time_step);
337: /* counting steps */
338: ts->steps++;
339: /* checking for maxerror */
340: /* comparing difference to maxerror */
341: VecNorm(rk->y2,NORM_2,&norm);
342: /* modifying maxerror to satisfy this timestep */
343: rk->maxerror = rk->ferror * ts->time_step;
344: /* PetscPrintf(PETSC_COMM_WORLD,"norm err: %f maxerror: %f dt: %f",norm,rk->maxerror,ts->time_step); */
346: /* handling ok and not ok */
347: if (norm < rk->maxerror){
348: /* if ok: */
349: ierr=VecCopy(rk->y1,ts->vec_sol); /* saves the suggested solution to current solution */
350: ts->ptime += ts->time_step; /* storing the new current time */
351: rk->nok++;
352: fac=5.0;
353: /* trying to save the vector */
354: TSPostStep(ts);
355: TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);
356: if (ts->ptime >= ts->max_time) break;
357: } else{
358: /* if not OK */
359: rk->nnok++;
360: fac=1.0;
361: ierr=VecCopy(ts->vec_sol,rk->y1); /* restores old solution */
362: }
364: /*Computing next stepsize. See page 167 in Solving ODE 1
365: *
366: * h_new = h * min( facmax , max( facmin , fac * (tol/err)^(1/(p+1)) ) )
367: * facmax set above
368: * facmin
369: */
370: dt_fac = exp(log((rk->maxerror) / norm) / ((rk->p) + 1) ) * 0.9 ;
372: if (dt_fac > fac){
373: /*PetscPrintf(PETSC_COMM_WORLD,"changing fac %f\n",fac);*/
374: dt_fac = fac;
375: }
377: /* computing new ts->time_step */
378: ts->time_step = ts->time_step * dt_fac;
380: if (ts->ptime+ts->time_step > ts->max_time){
381: ts->time_step = ts->max_time - ts->ptime;
382: }
384: if (ts->time_step < 1e-14){
385: PetscPrintf(PETSC_COMM_WORLD,"Very small steps: %f\n",ts->time_step);
386: ts->time_step = 1e-14;
387: }
389: /* trying to purify h */
390: /* (did not give any visible result) */
391: /* ttmp = ts->ptime + ts->time_step;
392: ts->time_step = ttmp - ts->ptime; */
394: }
396: ierr=VecCopy(rk->y1,ts->vec_sol);
397: return(0);
398: }
400: /*------------------------------------------------------------*/
403: static PetscErrorCode TSReset_RK(TS ts)
404: {
405: TS_RK *rk = (TS_RK*)ts->data;
409: VecDestroy(&rk->y1);
410: VecDestroy(&rk->y2);
411: VecDestroy(&rk->tmp);
412: VecDestroy(&rk->tmp_y);
413: if (rk->k) {VecDestroyVecs(rk->s,&rk->k);}
414: return(0);
415: }
419: static PetscErrorCode TSDestroy_RK(TS ts)
420: {
424: TSReset_RK(ts);
425: PetscFree(ts->data);
426: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRKSetTolerance_C","",PETSC_NULL);
427: return(0);
428: }
429: /*------------------------------------------------------------*/
433: static PetscErrorCode TSSetFromOptions_RK(TS ts)
434: {
435: TS_RK *rk = (TS_RK*)ts->data;
439: PetscOptionsHead("RK ODE solver options");
440: PetscOptionsReal("-ts_rk_tol","Tolerance for convergence","TSRKSetTolerance",rk->tolerance,&rk->tolerance,PETSC_NULL);
441: PetscOptionsTail();
442: return(0);
443: }
447: static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer)
448: {
449: TS_RK *rk = (TS_RK*)ts->data;
450: PetscBool iascii;
454: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
455: if (iascii) {
456: PetscViewerASCIIPrintf(viewer,"number of ok steps: %D\n",rk->nok);
457: PetscViewerASCIIPrintf(viewer,"number of rejected steps: %D\n",rk->nnok);
458: }
459: return(0);
460: }
462: /* ------------------------------------------------------------ */
463: /*MC
464: TSRK - ODE solver using the explicit Runge-Kutta methods
466: Options Database:
467: . -ts_rk_tol <tol> Tolerance for convergence
469: Contributed by: Asbjorn Hoiland Aarrestad, asbjorn@aarrestad.com, http://asbjorn.aarrestad.com/
471: Level: beginner
473: .seealso: TSCreate(), TS, TSSetType(), TSEULER, TSRKSetTolerance()
475: M*/
477: EXTERN_C_BEGIN
480: PetscErrorCode TSCreate_RK(TS ts)
481: {
482: TS_RK *rk;
486: ts->ops->setup = TSSetUp_RK;
487: ts->ops->solve = TSSolve_RK;
488: ts->ops->destroy = TSDestroy_RK;
489: ts->ops->setfromoptions = TSSetFromOptions_RK;
490: ts->ops->view = TSView_RK;
492: PetscNewLog(ts,TS_RK,&rk);
493: ts->data = (void*)rk;
495: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRKSetTolerance_C","TSRKSetTolerance_RK",TSRKSetTolerance_RK);
497: return(0);
498: }
499: EXTERN_C_END