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 "src/ts/tsimpl.h" /*I "petscts.h" I*/ 11 #include "time.h" 12 13 typedef struct { 14 Vec y1,y2; /* work wectors for the two rk permuations */ 15 int 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 int p,s; /* variables to tell the size of the runge-kutta solver */ 24 clock_t start,end; /* variables to mesure cpu time */ 25 } TS_Rk; 26 27 EXTERN_C_BEGIN 28 #undef __FUNCT__ 29 #define __FUNCT__ "TSRKSetTolerance_RK" 30 PetscErrorCode TSRKSetTolerance_RK(TS ts,PetscReal aabs) 31 { 32 TS_Rk *rk = (TS_Rk*)ts->data; 33 34 PetscFunctionBegin; 35 rk->tolerance = aabs; 36 PetscFunctionReturn(0); 37 } 38 EXTERN_C_END 39 40 #undef __FUNCT__ 41 #define __FUNCT__ "TSRKSetTolerance" 42 /*@ 43 TSRKSetTolerance - Sets the total error the RK explicit time integrators 44 will allow over the given time interval. 45 46 Collective on TS 47 48 Input parameters: 49 + ts - the time-step context 50 - aabs - the absolute tolerance 51 52 Level: intermediate 53 54 .keywords: RK, tolerance 55 56 .seealso: TSPVodeSetTolerance() 57 58 @*/ 59 PetscErrorCode TSRKSetTolerance(TS ts,PetscReal aabs) 60 { 61 PetscErrorCode ierr,(*f)(TS,PetscReal); 62 63 PetscFunctionBegin; 64 ierr = PetscObjectQueryFunction((PetscObject)ts,"TSRKSetTolerance_C",(void (**)(void))&f);CHKERRQ(ierr); 65 if (f) { 66 ierr = (*f)(ts,aabs);CHKERRQ(ierr); 67 } 68 PetscFunctionReturn(0); 69 } 70 71 72 #undef __FUNCT__ 73 #define __FUNCT__ "TSSetUp_Rk" 74 static PetscErrorCode TSSetUp_Rk(TS ts) 75 { 76 TS_Rk *rk = (TS_Rk*)ts->data; 77 int ierr; 78 79 PetscFunctionBegin; 80 rk->nok = 0; 81 rk->nnok = 0; 82 rk->maxerror = rk->tolerance; 83 84 /* fixing maxerror: global vs local */ 85 rk->ferror = rk->maxerror / (ts->max_time - ts->ptime); 86 87 /* 34.0/45.0 gives double precision division */ 88 /* defining variables needed for Runge-Kutta computing */ 89 /* when changing below, please remember to change a, b1, b2 and c above! */ 90 /* Found in table on page 171: Dormand-Prince 5(4) */ 91 92 /* are these right? */ 93 rk->p=6; 94 rk->s=7; 95 96 rk->a[1][0]=1.0/5.0; 97 rk->a[2][0]=3.0/40.0; 98 rk->a[2][1]=9.0/40.0; 99 rk->a[3][0]=44.0/45.0; 100 rk->a[3][1]=-56.0/15.0; 101 rk->a[3][2]=32.0/9.0; 102 rk->a[4][0]=19372.0/6561.0; 103 rk->a[4][1]=-25360.0/2187.0; 104 rk->a[4][2]=64448.0/6561.0; 105 rk->a[4][3]=-212.0/729.0; 106 rk->a[5][0]=9017.0/3168.0; 107 rk->a[5][1]=-355.0/33.0; 108 rk->a[5][2]=46732.0/5247.0; 109 rk->a[5][3]=49.0/176.0; 110 rk->a[5][4]=-5103.0/18656.0; 111 rk->a[6][0]=35.0/384.0; 112 rk->a[6][1]=0.0; 113 rk->a[6][2]=500.0/1113.0; 114 rk->a[6][3]=125.0/192.0; 115 rk->a[6][4]=-2187.0/6784.0; 116 rk->a[6][5]=11.0/84.0; 117 118 119 rk->c[0]=0.0; 120 rk->c[1]=1.0/5.0; 121 rk->c[2]=3.0/10.0; 122 rk->c[3]=4.0/5.0; 123 rk->c[4]=8.0/9.0; 124 rk->c[5]=1.0; 125 rk->c[6]=1.0; 126 127 rk->b1[0]=35.0/384.0; 128 rk->b1[1]=0.0; 129 rk->b1[2]=500.0/1113.0; 130 rk->b1[3]=125.0/192.0; 131 rk->b1[4]=-2187.0/6784.0; 132 rk->b1[5]=11.0/84.0; 133 rk->b1[6]=0.0; 134 135 rk->b2[0]=5179.0/57600.0; 136 rk->b2[1]=0.0; 137 rk->b2[2]=7571.0/16695.0; 138 rk->b2[3]=393.0/640.0; 139 rk->b2[4]=-92097.0/339200.0; 140 rk->b2[5]=187.0/2100.0; 141 rk->b2[6]=1.0/40.0; 142 143 144 /* Found in table on page 170: Fehlberg 4(5) */ 145 /* 146 rk->p=5; 147 rk->s=6; 148 149 rk->a[1][0]=1.0/4.0; 150 rk->a[2][0]=3.0/32.0; 151 rk->a[2][1]=9.0/32.0; 152 rk->a[3][0]=1932.0/2197.0; 153 rk->a[3][1]=-7200.0/2197.0; 154 rk->a[3][2]=7296.0/2197.0; 155 rk->a[4][0]=439.0/216.0; 156 rk->a[4][1]=-8.0; 157 rk->a[4][2]=3680.0/513.0; 158 rk->a[4][3]=-845.0/4104.0; 159 rk->a[5][0]=-8.0/27.0; 160 rk->a[5][1]=2.0; 161 rk->a[5][2]=-3544.0/2565.0; 162 rk->a[5][3]=1859.0/4104.0; 163 rk->a[5][4]=-11.0/40.0; 164 165 rk->c[0]=0.0; 166 rk->c[1]=1.0/4.0; 167 rk->c[2]=3.0/8.0; 168 rk->c[3]=12.0/13.0; 169 rk->c[4]=1.0; 170 rk->c[5]=1.0/2.0; 171 172 rk->b1[0]=25.0/216.0; 173 rk->b1[1]=0.0; 174 rk->b1[2]=1408.0/2565.0; 175 rk->b1[3]=2197.0/4104.0; 176 rk->b1[4]=-1.0/5.0; 177 rk->b1[5]=0.0; 178 179 rk->b2[0]=16.0/135.0; 180 rk->b2[1]=0.0; 181 rk->b2[2]=6656.0/12825.0; 182 rk->b2[3]=28561.0/56430.0; 183 rk->b2[4]=-9.0/50.0; 184 rk->b2[5]=2.0/55.0; 185 */ 186 /* Found in table on page 169: Merson 4("5") */ 187 /* 188 rk->p=4; 189 rk->s=5; 190 rk->a[1][0] = 1.0/3.0; 191 rk->a[2][0] = 1.0/6.0; 192 rk->a[2][1] = 1.0/6.0; 193 rk->a[3][0] = 1.0/8.0; 194 rk->a[3][1] = 0.0; 195 rk->a[3][2] = 3.0/8.0; 196 rk->a[4][0] = 1.0/2.0; 197 rk->a[4][1] = 0.0; 198 rk->a[4][2] = -3.0/2.0; 199 rk->a[4][3] = 2.0; 200 201 rk->c[0] = 0.0; 202 rk->c[1] = 1.0/3.0; 203 rk->c[2] = 1.0/3.0; 204 rk->c[3] = 0.5; 205 rk->c[4] = 1.0; 206 207 rk->b1[0] = 1.0/2.0; 208 rk->b1[1] = 0.0; 209 rk->b1[2] = -3.0/2.0; 210 rk->b1[3] = 2.0; 211 rk->b1[4] = 0.0; 212 213 rk->b2[0] = 1.0/6.0; 214 rk->b2[1] = 0.0; 215 rk->b2[2] = 0.0; 216 rk->b2[3] = 2.0/3.0; 217 rk->b2[4] = 1.0/6.0; 218 */ 219 220 /* making b2 -> e=b1-b2 */ 221 /* 222 for(i=0;i<rk->s;i++){ 223 rk->b2[i] = (rk->b1[i]) - (rk->b2[i]); 224 } 225 */ 226 rk->b2[0]=71.0/57600.0; 227 rk->b2[1]=0.0; 228 rk->b2[2]=-71.0/16695.0; 229 rk->b2[3]=71.0/1920.0; 230 rk->b2[4]=-17253.0/339200.0; 231 rk->b2[5]=22.0/525.0; 232 rk->b2[6]=-1.0/40.0; 233 234 /* initializing vectors */ 235 ierr = VecDuplicate(ts->vec_sol,&rk->y1);CHKERRQ(ierr); 236 ierr = VecDuplicate(ts->vec_sol,&rk->y2);CHKERRQ(ierr); 237 ierr = VecDuplicate(rk->y1,&rk->tmp);CHKERRQ(ierr); 238 ierr = VecDuplicate(rk->y1,&rk->tmp_y);CHKERRQ(ierr); 239 ierr = VecDuplicateVecs(rk->y1,rk->s,&rk->k);CHKERRQ(ierr); 240 241 PetscFunctionReturn(0); 242 } 243 244 /*------------------------------------------------------------*/ 245 #undef __FUNCT__ 246 #define __FUNCT__ "TSRkqs" 247 PetscErrorCode TSRkqs(TS ts,PetscReal t,PetscReal h) 248 { 249 TS_Rk *rk = (TS_Rk*)ts->data; 250 PetscErrorCode ierr; 251 int j,l; 252 PetscReal tmp_t=t; 253 PetscScalar null=0.0,hh=h; 254 255 /* printf("h: %f, hh: %f",h,hh); */ 256 257 PetscFunctionBegin; 258 259 /* k[0]=0 */ 260 ierr = VecSet(&null,rk->k[0]);CHKERRQ(ierr); 261 262 /* k[0] = derivs(t,y1) */ 263 ierr = TSComputeRHSFunction(ts,t,rk->y1,rk->k[0]);CHKERRQ(ierr); 264 /* looping over runge-kutta variables */ 265 /* building the k - array of vectors */ 266 for(j = 1 ; j < rk->s ; j++){ 267 268 /* rk->tmp = 0 */ 269 ierr = VecSet(&null,rk->tmp);CHKERRQ(ierr); 270 271 for(l=0;l<j;l++){ 272 /* tmp += a(j,l)*k[l] */ 273 /* ierr = PetscPrintf(PETSC_COMM_WORLD,"a(%i,%i)=%f \n",j,l,rk->a[j][l]);CHKERRQ(ierr); */ 274 ierr = VecAXPY(&rk->a[j][l],rk->k[l],rk->tmp);CHKERRQ(ierr); 275 } 276 277 /* ierr = VecView(rk->tmp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 278 279 /* k[j] = derivs(t+c(j)*h,y1+h*tmp,k(j)) */ 280 /* I need the following helpers: 281 PetscScalar tmp_t=t+c(j)*h 282 Vec tmp_y=h*tmp+y1 283 */ 284 285 tmp_t = t + rk->c[j] * h; 286 287 /* tmp_y = h * tmp + y1 */ 288 ierr = VecWAXPY(&hh,rk->tmp,rk->y1,rk->tmp_y);CHKERRQ(ierr); 289 290 /* rk->k[j]=0 */ 291 ierr = VecSet(&null,rk->k[j]);CHKERRQ(ierr); 292 ierr = TSComputeRHSFunction(ts,tmp_t,rk->tmp_y,rk->k[j]);CHKERRQ(ierr); 293 } 294 295 /* tmp=0 and tmp_y=0 */ 296 ierr = VecSet(&null,rk->tmp);CHKERRQ(ierr); 297 ierr = VecSet(&null,rk->tmp_y);CHKERRQ(ierr); 298 299 for(j = 0 ; j < rk->s ; j++){ 300 /* tmp=b1[j]*k[j]+tmp */ 301 ierr = VecAXPY(&rk->b1[j],rk->k[j],rk->tmp);CHKERRQ(ierr); 302 /* tmp_y=b2[j]*k[j]+tmp_y */ 303 ierr = VecAXPY(&rk->b2[j],rk->k[j],rk->tmp_y);CHKERRQ(ierr); 304 } 305 306 /* y2 = hh * tmp_y */ 307 ierr = VecSet(&null,rk->y2);CHKERRQ(ierr); 308 ierr = VecAXPY(&hh,rk->tmp_y,rk->y2);CHKERRQ(ierr); 309 /* y1 = hh*tmp + y1 */ 310 ierr = VecAXPY(&hh,rk->tmp,rk->y1);CHKERRQ(ierr); 311 /* Finding difference between y1 and y2 */ 312 313 PetscFunctionReturn(0); 314 } 315 316 #undef __FUNCT__ 317 #define __FUNCT__ "TSStep_Rk" 318 static PetscErrorCode TSStep_Rk(TS ts,int *steps,PetscReal *ptime) 319 { 320 TS_Rk *rk = (TS_Rk*)ts->data; 321 int ierr; 322 PetscReal dt = 0.001; /* fixed first step guess */ 323 PetscReal norm=0.0,dt_fac=0.0,fac = 0.0/*,ttmp=0.0*/; 324 325 PetscFunctionBegin; 326 rk->start=clock(); 327 ierr=VecCopy(ts->vec_sol,rk->y1);CHKERRQ(ierr); 328 *steps = -ts->steps; 329 /* trying to save the vector */ 330 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 331 /* while loop to get from start to stop */ 332 while (ts->ptime < ts->max_time){ 333 /* calling rkqs */ 334 /* 335 -- input 336 ts - pointer to ts 337 ts->ptime - current time 338 dt - try this timestep 339 y1 - solution for this step 340 341 --output 342 y1 - suggested solution 343 y2 - check solution (runge - kutta second permutation) 344 */ 345 ierr = TSRkqs(ts,ts->ptime,dt);CHKERRQ(ierr); 346 /* checking for maxerror */ 347 /* comparing difference to maxerror */ 348 ierr = VecNorm(rk->y2,NORM_2,&norm);CHKERRQ(ierr); 349 /* modifying maxerror to satisfy this timestep */ 350 rk->maxerror = rk->ferror * dt; 351 /* ierr = PetscPrintf(PETSC_COMM_WORLD,"norm err: %f maxerror: %f dt: %f",norm,rk->maxerror,dt);CHKERRQ(ierr); */ 352 353 /* handling ok and not ok */ 354 if(norm < rk->maxerror){ 355 /* if ok: */ 356 ierr=VecCopy(rk->y1,ts->vec_sol);CHKERRQ(ierr); /* saves the suggested solution to current solution */ 357 ts->ptime += dt; /* storing the new current time */ 358 rk->nok++; 359 fac=5.0; 360 /* trying to save the vector */ 361 /* calling monitor */ 362 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 363 }else{ 364 /* if not OK */ 365 rk->nnok++; 366 fac=1.0; 367 ierr=VecCopy(ts->vec_sol,rk->y1);CHKERRQ(ierr); /* restores old solution */ 368 } 369 370 /*Computing next stepsize. See page 167 in Solving ODE 1 371 * 372 * h_new = h * min( facmax , max( facmin , fac * (tol/err)^(1/(p+1)) ) ) 373 * facmax set above 374 * facmin 375 */ 376 dt_fac = exp(log((rk->maxerror) / norm) / ((rk->p) + 1) ) * 0.9 ; 377 378 if(dt_fac > fac){ 379 /*ierr = PetscPrintf(PETSC_COMM_WORLD,"changing fac %f\n",fac);*/ 380 dt_fac = fac; 381 } 382 383 /* computing new dt */ 384 dt = dt * dt_fac; 385 386 if(ts->ptime+dt > ts->max_time){ 387 dt = ts->max_time - ts->ptime; 388 } 389 390 if(dt < 1e-14){ 391 ierr = PetscPrintf(PETSC_COMM_WORLD,"Very small steps: %f\n",dt);CHKERRQ(ierr); 392 dt = 1e-14; 393 } 394 395 /* trying to purify h */ 396 /* (did not give any visible result) */ 397 /* ttmp = ts->ptime + dt; 398 dt = ttmp - ts->ptime; */ 399 400 /* counting steps */ 401 ts->steps++; 402 } 403 404 ierr=VecCopy(rk->y1,ts->vec_sol);CHKERRQ(ierr); 405 *steps += ts->steps; 406 *ptime = ts->ptime; 407 rk->end=clock(); 408 PetscFunctionReturn(0); 409 } 410 411 /*------------------------------------------------------------*/ 412 #undef __FUNCT__ 413 #define __FUNCT__ "TSDestroy_Rk" 414 static PetscErrorCode TSDestroy_Rk(TS ts) 415 { 416 TS_Rk *rk = (TS_Rk*)ts->data; 417 PetscErrorCode ierr; 418 int i; 419 420 /* REMEMBER TO DESTROY ALL */ 421 422 PetscFunctionBegin; 423 if (rk->y1) {ierr = VecDestroy(rk->y1);CHKERRQ(ierr);} 424 if (rk->y2) {ierr = VecDestroy(rk->y2);CHKERRQ(ierr);} 425 if (rk->tmp) {ierr = VecDestroy(rk->tmp);CHKERRQ(ierr);} 426 if (rk->tmp_y) {ierr = VecDestroy(rk->tmp_y);CHKERRQ(ierr);} 427 for(i=0;i<rk->s;i++){ 428 if (rk->k[i]) {ierr = VecDestroy(rk->k[i]);CHKERRQ(ierr);} 429 } 430 ierr = PetscFree(rk);CHKERRQ(ierr); 431 PetscFunctionReturn(0); 432 } 433 /*------------------------------------------------------------*/ 434 435 #undef __FUNCT__ 436 #define __FUNCT__ "TSSetFromOptions_Rk" 437 static PetscErrorCode TSSetFromOptions_Rk(TS ts) 438 { 439 TS_Rk *rk = (TS_Rk*)ts->data; 440 PetscErrorCode ierr; 441 442 PetscFunctionBegin; 443 ierr = PetscOptionsHead("RK ODE solver options");CHKERRQ(ierr); 444 ierr = PetscOptionsReal("-ts_rk_tol","Tolerance for convergence","TSRKSetTolerance",rk->tolerance,&rk->tolerance,PETSC_NULL);CHKERRQ(ierr); 445 ierr = PetscOptionsTail();CHKERRQ(ierr); 446 PetscFunctionReturn(0); 447 } 448 449 #undef __FUNCT__ 450 #define __FUNCT__ "TSView_Rk" 451 static PetscErrorCode TSView_Rk(TS ts,PetscViewer viewer) 452 { 453 TS_Rk *rk = (TS_Rk*)ts->data; 454 PetscErrorCode ierr; 455 /*double elapsed;*/ 456 457 PetscFunctionBegin; 458 ierr = PetscPrintf(PETSC_COMM_WORLD," number of ok steps: %d\n",rk->nok);CHKERRQ(ierr); 459 ierr = PetscPrintf(PETSC_COMM_WORLD," number of rejected steps: %d\n",rk->nnok);CHKERRQ(ierr); 460 /* elapsed = ((double) (rk->end - rk->start)) / CLOCKS_PER_SEC; 461 462 ierr = PetscPrintf(PETSC_COMM_WORLD," CPU time used (in seconds): %f\n",elapsed);CHKERRQ(ierr); */ 463 PetscFunctionReturn(0); 464 } 465 466 /* ------------------------------------------------------------ */ 467 /*MC 468 TS_RK - ODE solver using the explicit Runge-Kutta methods 469 470 Options Database: 471 . -ts_rk_tol <tol> Tolerance for convergence 472 473 Contributed by: Asbjorn Hoiland Aarrestad, asbjorn@aarrestad.com, http://asbjorn.aarrestad.com/ 474 475 Level: beginner 476 477 .seealso: TSCreate(), TS, TSSetType(), TS_EULER, TSRKSetTolerance() 478 479 M*/ 480 481 EXTERN_C_BEGIN 482 #undef __FUNCT__ 483 #define __FUNCT__ "TSCreate_Rk" 484 PetscErrorCode TSCreate_Rk(TS ts) 485 { 486 TS_Rk *rk; 487 PetscErrorCode ierr; 488 489 PetscFunctionBegin; 490 ts->ops->setup = TSSetUp_Rk; 491 ts->ops->step = TSStep_Rk; 492 ts->ops->destroy = TSDestroy_Rk; 493 ts->ops->setfromoptions = TSSetFromOptions_Rk; 494 ts->ops->view = TSView_Rk; 495 496 ierr = PetscNew(TS_Rk,&rk);CHKERRQ(ierr); 497 PetscLogObjectMemory(ts,sizeof(TS_Rk)); 498 ts->data = (void*)rk; 499 500 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRKSetTolerance_C","TSRKSetTolerance_RK", 501 TSRKSetTolerance_RK);CHKERRQ(ierr); 502 503 PetscFunctionReturn(0); 504 } 505 EXTERN_C_END 506 507 508 509 510