Actual source code: rk.c

petsc-3.4.4 2014-03-13
  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;

 28: PetscErrorCode  TSRKSetTolerance_RK(TS ts,PetscReal aabs)
 29: {
 30:   TS_RK *rk = (TS_RK*)ts->data;

 33:   rk->tolerance = aabs;
 34:   return(0);
 35: }

 39: /*@
 40:    TSRKSetTolerance - Sets the total error the RK explicit time integrators
 41:                       will allow over the given time interval.

 43:    Logically Collective on TS

 45:    Input parameters:
 46: +    ts  - the time-step context
 47: -    aabs - the absolute tolerance

 49:    Level: intermediate

 51: .keywords: RK, tolerance

 53: .seealso: TSSundialsSetTolerance()

 55: @*/
 56: PetscErrorCode  TSRKSetTolerance(TS ts,PetscReal aabs)
 57: {

 62:   PetscTryMethod(ts,"TSRKSetTolerance_C",(TS,PetscReal),(ts,aabs));
 63:   return(0);
 64: }


 69: static PetscErrorCode TSSetUp_RK(TS ts)
 70: {
 71:   TS_RK          *rk = (TS_RK*)ts->data;

 75:   rk->nok      = 0;
 76:   rk->nnok     = 0;
 77:   rk->maxerror = rk->tolerance;

 79:   /* fixing maxerror: global vs local */
 80:   rk->ferror = rk->maxerror / (ts->max_time - ts->ptime);

 82:   /* 34.0/45.0 gives double precision division */
 83:   /* defining variables needed for Runge-Kutta computing */
 84:   /* when changing below, please remember to change a, b1, b2 and c above! */
 85:   /* Found in table on page 171: Dormand-Prince 5(4) */

 87:   /* are these right? */
 88:   rk->p=6;
 89:   rk->s=7;

 91:   rk->a[1][0]=1.0/5.0;
 92:   rk->a[2][0]=3.0/40.0;
 93:   rk->a[2][1]=9.0/40.0;
 94:   rk->a[3][0]=44.0/45.0;
 95:   rk->a[3][1]=-56.0/15.0;
 96:   rk->a[3][2]=32.0/9.0;
 97:   rk->a[4][0]=19372.0/6561.0;
 98:   rk->a[4][1]=-25360.0/2187.0;
 99:   rk->a[4][2]=64448.0/6561.0;
100:   rk->a[4][3]=-212.0/729.0;
101:   rk->a[5][0]=9017.0/3168.0;
102:   rk->a[5][1]=-355.0/33.0;
103:   rk->a[5][2]=46732.0/5247.0;
104:   rk->a[5][3]=49.0/176.0;
105:   rk->a[5][4]=-5103.0/18656.0;
106:   rk->a[6][0]=35.0/384.0;
107:   rk->a[6][1]=0.0;
108:   rk->a[6][2]=500.0/1113.0;
109:   rk->a[6][3]=125.0/192.0;
110:   rk->a[6][4]=-2187.0/6784.0;
111:   rk->a[6][5]=11.0/84.0;


114:   rk->c[0]=0.0;
115:   rk->c[1]=1.0/5.0;
116:   rk->c[2]=3.0/10.0;
117:   rk->c[3]=4.0/5.0;
118:   rk->c[4]=8.0/9.0;
119:   rk->c[5]=1.0;
120:   rk->c[6]=1.0;

122:   rk->b1[0]=35.0/384.0;
123:   rk->b1[1]=0.0;
124:   rk->b1[2]=500.0/1113.0;
125:   rk->b1[3]=125.0/192.0;
126:   rk->b1[4]=-2187.0/6784.0;
127:   rk->b1[5]=11.0/84.0;
128:   rk->b1[6]=0.0;

130:   rk->b2[0]=5179.0/57600.0;
131:   rk->b2[1]=0.0;
132:   rk->b2[2]=7571.0/16695.0;
133:   rk->b2[3]=393.0/640.0;
134:   rk->b2[4]=-92097.0/339200.0;
135:   rk->b2[5]=187.0/2100.0;
136:   rk->b2[6]=1.0/40.0;


139:   /* Found in table on page 170: Fehlberg 4(5) */
140:   /*
141:   rk->p=5;
142:   rk->s=6;

144:   rk->a[1][0]=1.0/4.0;
145:   rk->a[2][0]=3.0/32.0;
146:   rk->a[2][1]=9.0/32.0;
147:   rk->a[3][0]=1932.0/2197.0;
148:   rk->a[3][1]=-7200.0/2197.0;
149:   rk->a[3][2]=7296.0/2197.0;
150:   rk->a[4][0]=439.0/216.0;
151:   rk->a[4][1]=-8.0;
152:   rk->a[4][2]=3680.0/513.0;
153:   rk->a[4][3]=-845.0/4104.0;
154:   rk->a[5][0]=-8.0/27.0;
155:   rk->a[5][1]=2.0;
156:   rk->a[5][2]=-3544.0/2565.0;
157:   rk->a[5][3]=1859.0/4104.0;
158:   rk->a[5][4]=-11.0/40.0;

160:   rk->c[0]=0.0;
161:   rk->c[1]=1.0/4.0;
162:   rk->c[2]=3.0/8.0;
163:   rk->c[3]=12.0/13.0;
164:   rk->c[4]=1.0;
165:   rk->c[5]=1.0/2.0;

167:   rk->b1[0]=25.0/216.0;
168:   rk->b1[1]=0.0;
169:   rk->b1[2]=1408.0/2565.0;
170:   rk->b1[3]=2197.0/4104.0;
171:   rk->b1[4]=-1.0/5.0;
172:   rk->b1[5]=0.0;

174:   rk->b2[0]=16.0/135.0;
175:   rk->b2[1]=0.0;
176:   rk->b2[2]=6656.0/12825.0;
177:   rk->b2[3]=28561.0/56430.0;
178:   rk->b2[4]=-9.0/50.0;
179:   rk->b2[5]=2.0/55.0;
180:   */
181:   /* Found in table on page 169: Merson 4("5") */
182:   /*
183:   rk->p=4;
184:   rk->s=5;
185:   rk->a[1][0] = 1.0/3.0;
186:   rk->a[2][0] = 1.0/6.0;
187:   rk->a[2][1] = 1.0/6.0;
188:   rk->a[3][0] = 1.0/8.0;
189:   rk->a[3][1] = 0.0;
190:   rk->a[3][2] = 3.0/8.0;
191:   rk->a[4][0] = 1.0/2.0;
192:   rk->a[4][1] = 0.0;
193:   rk->a[4][2] = -3.0/2.0;
194:   rk->a[4][3] = 2.0;

196:   rk->c[0] = 0.0;
197:   rk->c[1] = 1.0/3.0;
198:   rk->c[2] = 1.0/3.0;
199:   rk->c[3] = 0.5;
200:   rk->c[4] = 1.0;

202:   rk->b1[0] = 1.0/2.0;
203:   rk->b1[1] = 0.0;
204:   rk->b1[2] = -3.0/2.0;
205:   rk->b1[3] = 2.0;
206:   rk->b1[4] = 0.0;

208:   rk->b2[0] = 1.0/6.0;
209:   rk->b2[1] = 0.0;
210:   rk->b2[2] = 0.0;
211:   rk->b2[3] = 2.0/3.0;
212:   rk->b2[4] = 1.0/6.0;
213:   */

215:   /* making b2 -> e=b1-b2 */
216:   /*
217:     for (i=0;i<rk->s;i++) {
218:      rk->b2[i] = (rk->b1[i]) - (rk->b2[i]);
219:   }
220:   */
221:   rk->b2[0]=71.0/57600.0;
222:   rk->b2[1]=0.0;
223:   rk->b2[2]=-71.0/16695.0;
224:   rk->b2[3]=71.0/1920.0;
225:   rk->b2[4]=-17253.0/339200.0;
226:   rk->b2[5]=22.0/525.0;
227:   rk->b2[6]=-1.0/40.0;

229:   /* initializing vectors */
230:   VecDuplicate(ts->vec_sol,&rk->y1);
231:   VecDuplicate(ts->vec_sol,&rk->y2);
232:   VecDuplicate(rk->y1,&rk->tmp);
233:   VecDuplicate(rk->y1,&rk->tmp_y);
234:   VecDuplicateVecs(rk->y1,rk->s,&rk->k);
235:   return(0);
236: }

238: /*------------------------------------------------------------*/
241: PetscErrorCode TSRKqs(TS ts,PetscReal t,PetscReal h)
242: {
243:   TS_RK          *rk = (TS_RK*)ts->data;
245:   PetscInt       j,l;
246:   PetscReal      tmp_t = t;
247:   PetscScalar    hh    = h;

250:   /* k[0]=0  */
251:   VecSet(rk->k[0],0.0);

253:   /* k[0] = derivs(t,y1) */
254:   TSComputeRHSFunction(ts,t,rk->y1,rk->k[0]);
255:   /* looping over runge-kutta variables */
256:   /* building the k - array of vectors */
257:   for (j = 1; j < rk->s; j++) {

259:     /* rk->tmp = 0 */
260:     VecSet(rk->tmp,0.0);

262:     for (l=0; l<j; l++) {
263:       /* tmp += a(j,l)*k[l] */
264:       VecAXPY(rk->tmp,rk->a[j][l],rk->k[l]);
265:     }

267:     /* VecView(rk->tmp,PETSC_VIEWER_STDOUT_WORLD); */

269:     /* k[j] = derivs(t+c(j)*h,y1+h*tmp,k(j)) */
270:     /* I need the following helpers:
271:        PetscScalar  tmp_t=t+c(j)*h
272:        Vec          tmp_y=h*tmp+y1
273:     */

275:     tmp_t = t + rk->c[j] * h;

277:     /* tmp_y = h * tmp + y1 */
278:     VecWAXPY(rk->tmp_y,hh,rk->tmp,rk->y1);

280:     /* rk->k[j]=0 */
281:     VecSet(rk->k[j],0.0);
282:     TSComputeRHSFunction(ts,tmp_t,rk->tmp_y,rk->k[j]);
283:   }

285:   /* tmp=0 and tmp_y=0 */
286:   VecSet(rk->tmp,0.0);
287:   VecSet(rk->tmp_y,0.0);

289:   for (j = 0; j < rk->s; j++) {
290:     /* tmp=b1[j]*k[j]+tmp  */
291:     VecAXPY(rk->tmp,rk->b1[j],rk->k[j]);
292:     /* tmp_y=b2[j]*k[j]+tmp_y */
293:     VecAXPY(rk->tmp_y,rk->b2[j],rk->k[j]);
294:   }

296:   /* y2 = hh * tmp_y */
297:   VecSet(rk->y2,0.0);
298:   VecAXPY(rk->y2,hh,rk->tmp_y);
299:   /* y1 = hh*tmp + y1 */
300:   VecAXPY(rk->y1,hh,rk->tmp);
301:   /* Finding difference between y1 and y2 */
302:   return(0);
303: }

307: static PetscErrorCode TSSolve_RK(TS ts)
308: {
309:   TS_RK          *rk = (TS_RK*)ts->data;
310:   PetscReal      norm=0.0,dt_fac=0.0,fac = 0.0 /*,ttmp=0.0*/;
311:   PetscInt       i;

315:   VecCopy(ts->vec_sol,rk->y1);

317:   /* while loop to get from start to stop */
318:   for (i = 0; i < ts->max_steps; i++) {
319:     TSPreStep(ts); /* Note that this is called once per STEP, not once per STAGE. */

321:     /* calling rkqs */
322:     /*
323:       -- input
324:       ts        - pointer to ts
325:       ts->ptime - current time
326:       ts->time_step        - try this timestep
327:       y1        - solution for this step

329:       --output
330:       y1        - suggested solution
331:       y2        - check solution (runge - kutta second permutation)
332:     */
333:     TSRKqs(ts,ts->ptime,ts->time_step);
334:     /* counting steps */
335:     ts->steps++;
336:     /* checking for maxerror */
337:     /* comparing difference to maxerror */
338:     VecNorm(rk->y2,NORM_2,&norm);
339:     /* modifying maxerror to satisfy this timestep */
340:     rk->maxerror = rk->ferror * ts->time_step;
341:     /* PetscPrintf(PETSC_COMM_WORLD,"norm err: %f maxerror: %f dt: %f",norm,rk->maxerror,ts->time_step); */

343:     /* handling ok and not ok */
344:     if (norm < rk->maxerror) {
345:       /* if ok: */
346:       VecCopy(rk->y1,ts->vec_sol); /* saves the suggested solution to current solution */
347:       ts->ptime += ts->time_step;   /* storing the new current time */
348:       rk->nok++;
349:       fac=5.0;
350:       /* trying to save the vector */
351:       TSPostStep(ts);
352:       TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);
353:       if (ts->ptime >= ts->max_time) break;
354:     } else {
355:       /* if not OK */
356:       rk->nnok++;
357:       fac =1.0;
358:       ierr=VecCopy(ts->vec_sol,rk->y1);    /* restores old solution */
359:     }

361:     /*Computing next stepsize. See page 167 in Solving ODE 1
362:      *
363:      * h_new = h * min(facmax , max(facmin , fac * (tol/err)^(1/(p+1))))
364:      * facmax set above
365:      * facmin
366:      */
367:     dt_fac = exp(log((rk->maxerror) / norm) / ((rk->p) + 1)) * 0.9;

369:     if (dt_fac > fac) dt_fac = fac;
370: 

372:     /* computing new ts->time_step */
373:     ts->time_step = ts->time_step * dt_fac;

375:     if (ts->ptime+ts->time_step > ts->max_time) ts->time_step = ts->max_time - ts->ptime;

377:     if (ts->time_step < 1e-14) {
378:       PetscPrintf(PETSC_COMM_WORLD,"Very small steps: %f\n",ts->time_step);
379:       ts->time_step = 1e-14;
380:     }

382:     /* trying to purify h */
383:     /* (did not give any visible result) */
384:     /* ttmp = ts->ptime + ts->time_step;
385:        ts->time_step = ttmp - ts->ptime; */

387:   }

389:   ierr=VecCopy(rk->y1,ts->vec_sol);
390:   return(0);
391: }

393: /*------------------------------------------------------------*/
396: static PetscErrorCode TSReset_RK(TS ts)
397: {
398:   TS_RK          *rk = (TS_RK*)ts->data;

402:   VecDestroy(&rk->y1);
403:   VecDestroy(&rk->y2);
404:   VecDestroy(&rk->tmp);
405:   VecDestroy(&rk->tmp_y);
406:   if (rk->k) {VecDestroyVecs(rk->s,&rk->k);}
407:   return(0);
408: }

412: static PetscErrorCode TSDestroy_RK(TS ts)
413: {

417:   TSReset_RK(ts);
418:   PetscFree(ts->data);
419:   PetscObjectComposeFunction((PetscObject)ts,"TSRKSetTolerance_C",NULL);
420:   return(0);
421: }
422: /*------------------------------------------------------------*/

426: static PetscErrorCode TSSetFromOptions_RK(TS ts)
427: {
428:   TS_RK          *rk = (TS_RK*)ts->data;

432:   PetscOptionsHead("RK ODE solver options");
433:   PetscOptionsReal("-ts_rk_tol","Tolerance for convergence","TSRKSetTolerance",rk->tolerance,&rk->tolerance,NULL);
434:   PetscOptionsTail();
435:   return(0);
436: }

440: static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer)
441: {
442:   TS_RK          *rk = (TS_RK*)ts->data;
443:   PetscBool      iascii;

447:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
448:   if (iascii) {
449:     PetscViewerASCIIPrintf(viewer,"number of ok steps: %D\n",rk->nok);
450:     PetscViewerASCIIPrintf(viewer,"number of rejected steps: %D\n",rk->nnok);
451:   }
452:   return(0);
453: }

455: /* ------------------------------------------------------------ */
456: /*MC
457:       TSRK - ODE solver using the explicit Runge-Kutta methods

459:    Options Database:
460: .  -ts_rk_tol <tol> Tolerance for convergence

462:   Contributed by: Asbjorn Hoiland Aarrestad, asbjorn@aarrestad.com, http://asbjorn.aarrestad.com/

464:   Level: beginner

466: .seealso:  TSCreate(), TS, TSSetType(), TSEULER, TSRKSetTolerance()

468: M*/

472: PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts)
473: {
474:   TS_RK          *rk;

478:   ts->ops->setup          = TSSetUp_RK;
479:   ts->ops->solve          = TSSolve_RK;
480:   ts->ops->destroy        = TSDestroy_RK;
481:   ts->ops->setfromoptions = TSSetFromOptions_RK;
482:   ts->ops->view           = TSView_RK;

484:   PetscNewLog(ts,TS_RK,&rk);
485:   ts->data = (void*)rk;

487:   PetscObjectComposeFunction((PetscObject)ts,"TSRKSetTolerance_C",TSRKSetTolerance_RK);
488:   return(0);
489: }