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