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> 12 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; 25 26 #undef __FUNCT__ 27 #define __FUNCT__ "TSRKSetTolerance_RK" 28 PetscErrorCode TSRKSetTolerance_RK(TS ts,PetscReal aabs) 29 { 30 TS_RK *rk = (TS_RK*)ts->data; 31 32 PetscFunctionBegin; 33 rk->tolerance = aabs; 34 PetscFunctionReturn(0); 35 } 36 37 #undef __FUNCT__ 38 #define __FUNCT__ "TSRKSetTolerance" 39 /*@ 40 TSRKSetTolerance - Sets the total error the RK explicit time integrators 41 will allow over the given time interval. 42 43 Logically Collective on TS 44 45 Input parameters: 46 + ts - the time-step context 47 - aabs - the absolute tolerance 48 49 Level: intermediate 50 51 .keywords: RK, tolerance 52 53 .seealso: TSSundialsSetTolerance() 54 55 @*/ 56 PetscErrorCode TSRKSetTolerance(TS ts,PetscReal aabs) 57 { 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 PetscValidLogicalCollectiveReal(ts,aabs,2); 62 ierr = PetscTryMethod(ts,"TSRKSetTolerance_C",(TS,PetscReal),(ts,aabs));CHKERRQ(ierr); 63 PetscFunctionReturn(0); 64 } 65 66 67 #undef __FUNCT__ 68 #define __FUNCT__ "TSSetUp_RK" 69 static PetscErrorCode TSSetUp_RK(TS ts) 70 { 71 TS_RK *rk = (TS_RK*)ts->data; 72 PetscErrorCode ierr; 73 74 PetscFunctionBegin; 75 rk->nok = 0; 76 rk->nnok = 0; 77 rk->maxerror = rk->tolerance; 78 79 /* fixing maxerror: global vs local */ 80 rk->ferror = rk->maxerror / (ts->max_time - ts->ptime); 81 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) */ 86 87 /* are these right? */ 88 rk->p=6; 89 rk->s=7; 90 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; 112 113 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; 121 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; 129 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; 137 138 139 /* Found in table on page 170: Fehlberg 4(5) */ 140 /* 141 rk->p=5; 142 rk->s=6; 143 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; 159 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; 166 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; 173 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; 195 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; 201 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; 207 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 */ 214 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; 228 229 /* initializing vectors */ 230 ierr = VecDuplicate(ts->vec_sol,&rk->y1);CHKERRQ(ierr); 231 ierr = VecDuplicate(ts->vec_sol,&rk->y2);CHKERRQ(ierr); 232 ierr = VecDuplicate(rk->y1,&rk->tmp);CHKERRQ(ierr); 233 ierr = VecDuplicate(rk->y1,&rk->tmp_y);CHKERRQ(ierr); 234 ierr = VecDuplicateVecs(rk->y1,rk->s,&rk->k);CHKERRQ(ierr); 235 PetscFunctionReturn(0); 236 } 237 238 /*------------------------------------------------------------*/ 239 #undef __FUNCT__ 240 #define __FUNCT__ "TSRKqs" 241 PetscErrorCode TSRKqs(TS ts,PetscReal t,PetscReal h) 242 { 243 TS_RK *rk = (TS_RK*)ts->data; 244 PetscErrorCode ierr; 245 PetscInt j,l; 246 PetscReal tmp_t = t; 247 PetscScalar hh = h; 248 249 PetscFunctionBegin; 250 /* k[0]=0 */ 251 ierr = VecSet(rk->k[0],0.0);CHKERRQ(ierr); 252 253 /* k[0] = derivs(t,y1) */ 254 ierr = TSComputeRHSFunction(ts,t,rk->y1,rk->k[0]);CHKERRQ(ierr); 255 /* looping over runge-kutta variables */ 256 /* building the k - array of vectors */ 257 for (j = 1; j < rk->s; j++) { 258 259 /* rk->tmp = 0 */ 260 ierr = VecSet(rk->tmp,0.0);CHKERRQ(ierr); 261 262 for (l=0; l<j; l++) { 263 /* tmp += a(j,l)*k[l] */ 264 ierr = VecAXPY(rk->tmp,rk->a[j][l],rk->k[l]);CHKERRQ(ierr); 265 } 266 267 /* ierr = VecView(rk->tmp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 268 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 */ 274 275 tmp_t = t + rk->c[j] * h; 276 277 /* tmp_y = h * tmp + y1 */ 278 ierr = VecWAXPY(rk->tmp_y,hh,rk->tmp,rk->y1);CHKERRQ(ierr); 279 280 /* rk->k[j]=0 */ 281 ierr = VecSet(rk->k[j],0.0);CHKERRQ(ierr); 282 ierr = TSComputeRHSFunction(ts,tmp_t,rk->tmp_y,rk->k[j]);CHKERRQ(ierr); 283 } 284 285 /* tmp=0 and tmp_y=0 */ 286 ierr = VecSet(rk->tmp,0.0);CHKERRQ(ierr); 287 ierr = VecSet(rk->tmp_y,0.0);CHKERRQ(ierr); 288 289 for (j = 0; j < rk->s; j++) { 290 /* tmp=b1[j]*k[j]+tmp */ 291 ierr = VecAXPY(rk->tmp,rk->b1[j],rk->k[j]);CHKERRQ(ierr); 292 /* tmp_y=b2[j]*k[j]+tmp_y */ 293 ierr = VecAXPY(rk->tmp_y,rk->b2[j],rk->k[j]);CHKERRQ(ierr); 294 } 295 296 /* y2 = hh * tmp_y */ 297 ierr = VecSet(rk->y2,0.0);CHKERRQ(ierr); 298 ierr = VecAXPY(rk->y2,hh,rk->tmp_y);CHKERRQ(ierr); 299 /* y1 = hh*tmp + y1 */ 300 ierr = VecAXPY(rk->y1,hh,rk->tmp);CHKERRQ(ierr); 301 /* Finding difference between y1 and y2 */ 302 PetscFunctionReturn(0); 303 } 304 305 #undef __FUNCT__ 306 #define __FUNCT__ "TSSolve_RK" 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; 312 PetscErrorCode ierr; 313 314 PetscFunctionBegin; 315 ierr = VecCopy(ts->vec_sol,rk->y1);CHKERRQ(ierr); 316 317 /* while loop to get from start to stop */ 318 for (i = 0; i < ts->max_steps; i++) { 319 ierr = TSPreStep(ts);CHKERRQ(ierr); /* Note that this is called once per STEP, not once per STAGE. */ 320 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 328 329 --output 330 y1 - suggested solution 331 y2 - check solution (runge - kutta second permutation) 332 */ 333 ierr = TSRKqs(ts,ts->ptime,ts->time_step);CHKERRQ(ierr); 334 /* counting steps */ 335 ts->steps++; 336 /* checking for maxerror */ 337 /* comparing difference to maxerror */ 338 ierr = VecNorm(rk->y2,NORM_2,&norm);CHKERRQ(ierr); 339 /* modifying maxerror to satisfy this timestep */ 340 rk->maxerror = rk->ferror * ts->time_step; 341 /* ierr = PetscPrintf(PETSC_COMM_WORLD,"norm err: %f maxerror: %f dt: %f",norm,rk->maxerror,ts->time_step);CHKERRQ(ierr); */ 342 343 /* handling ok and not ok */ 344 if (norm < rk->maxerror) { 345 /* if ok: */ 346 ierr = VecCopy(rk->y1,ts->vec_sol);CHKERRQ(ierr); /* 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 ierr = TSPostStep(ts);CHKERRQ(ierr); 352 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 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);CHKERRQ(ierr); /* restores old solution */ 359 } 360 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; 368 369 if (dt_fac > fac) dt_fac = fac; 370 371 372 /* computing new ts->time_step */ 373 ts->time_step = ts->time_step * dt_fac; 374 375 if (ts->ptime+ts->time_step > ts->max_time) ts->time_step = ts->max_time - ts->ptime; 376 377 if (ts->time_step < 1e-14) { 378 ierr = PetscPrintf(PETSC_COMM_WORLD,"Very small steps: %f\n",ts->time_step);CHKERRQ(ierr); 379 ts->time_step = 1e-14; 380 } 381 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; */ 386 387 } 388 389 ierr=VecCopy(rk->y1,ts->vec_sol);CHKERRQ(ierr); 390 PetscFunctionReturn(0); 391 } 392 393 /*------------------------------------------------------------*/ 394 #undef __FUNCT__ 395 #define __FUNCT__ "TSReset_RK" 396 static PetscErrorCode TSReset_RK(TS ts) 397 { 398 TS_RK *rk = (TS_RK*)ts->data; 399 PetscErrorCode ierr; 400 401 PetscFunctionBegin; 402 ierr = VecDestroy(&rk->y1);CHKERRQ(ierr); 403 ierr = VecDestroy(&rk->y2);CHKERRQ(ierr); 404 ierr = VecDestroy(&rk->tmp);CHKERRQ(ierr); 405 ierr = VecDestroy(&rk->tmp_y);CHKERRQ(ierr); 406 if (rk->k) {ierr = VecDestroyVecs(rk->s,&rk->k);CHKERRQ(ierr);} 407 PetscFunctionReturn(0); 408 } 409 410 #undef __FUNCT__ 411 #define __FUNCT__ "TSDestroy_RK" 412 static PetscErrorCode TSDestroy_RK(TS ts) 413 { 414 PetscErrorCode ierr; 415 416 PetscFunctionBegin; 417 ierr = TSReset_RK(ts);CHKERRQ(ierr); 418 ierr = PetscFree(ts->data);CHKERRQ(ierr); 419 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetTolerance_C",NULL);CHKERRQ(ierr); 420 PetscFunctionReturn(0); 421 } 422 /*------------------------------------------------------------*/ 423 424 #undef __FUNCT__ 425 #define __FUNCT__ "TSSetFromOptions_RK" 426 static PetscErrorCode TSSetFromOptions_RK(TS ts) 427 { 428 TS_RK *rk = (TS_RK*)ts->data; 429 PetscErrorCode ierr; 430 431 PetscFunctionBegin; 432 ierr = PetscOptionsHead("RK ODE solver options");CHKERRQ(ierr); 433 ierr = PetscOptionsReal("-ts_rk_tol","Tolerance for convergence","TSRKSetTolerance",rk->tolerance,&rk->tolerance,NULL);CHKERRQ(ierr); 434 ierr = PetscOptionsTail();CHKERRQ(ierr); 435 PetscFunctionReturn(0); 436 } 437 438 #undef __FUNCT__ 439 #define __FUNCT__ "TSView_RK" 440 static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer) 441 { 442 TS_RK *rk = (TS_RK*)ts->data; 443 PetscBool iascii; 444 PetscErrorCode ierr; 445 446 PetscFunctionBegin; 447 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 448 if (iascii) { 449 ierr = PetscViewerASCIIPrintf(viewer,"number of ok steps: %D\n",rk->nok);CHKERRQ(ierr); 450 ierr = PetscViewerASCIIPrintf(viewer,"number of rejected steps: %D\n",rk->nnok);CHKERRQ(ierr); 451 } 452 PetscFunctionReturn(0); 453 } 454 455 /* ------------------------------------------------------------ */ 456 /*MC 457 TSRK - ODE solver using the explicit Runge-Kutta methods 458 459 Options Database: 460 . -ts_rk_tol <tol> Tolerance for convergence 461 462 Contributed by: Asbjorn Hoiland Aarrestad, asbjorn@aarrestad.com, http://asbjorn.aarrestad.com/ 463 464 Level: beginner 465 466 .seealso: TSCreate(), TS, TSSetType(), TSEULER, TSRKSetTolerance() 467 468 M*/ 469 470 #undef __FUNCT__ 471 #define __FUNCT__ "TSCreate_RK" 472 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 473 { 474 TS_RK *rk; 475 PetscErrorCode ierr; 476 477 PetscFunctionBegin; 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; 483 484 ierr = PetscNewLog(ts,TS_RK,&rk);CHKERRQ(ierr); 485 ts->data = (void*)rk; 486 487 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetTolerance_C",TSRKSetTolerance_RK);CHKERRQ(ierr); 488 PetscFunctionReturn(0); 489 } 490