1 /* 2 Code for time stepping with the Runge-Kutta method 3 4 Notes: 5 The general system is written as 6 7 Udot = F(t,U) 8 9 */ 10 #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/ 11 #include <petscdm.h> 12 13 static TSRKType TSRKDefault = TSRK3BS; 14 static PetscBool TSRKRegisterAllCalled; 15 static PetscBool TSRKPackageInitialized; 16 static PetscInt explicit_stage_time_id; 17 18 typedef struct _RKTableau *RKTableau; 19 struct _RKTableau { 20 char *name; 21 PetscInt order; /* Classical approximation order of the method i */ 22 PetscInt s; /* Number of stages */ 23 PetscBool FSAL; /* flag to indicate if tableau is FSAL */ 24 PetscInt pinterp; /* Interpolation order */ 25 PetscReal *A,*b,*c; /* Tableau */ 26 PetscReal *bembed; /* Embedded formula of order one less (order-1) */ 27 PetscReal *binterp; /* Dense output formula */ 28 PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 29 }; 30 typedef struct _RKTableauLink *RKTableauLink; 31 struct _RKTableauLink { 32 struct _RKTableau tab; 33 RKTableauLink next; 34 }; 35 static RKTableauLink RKTableauList; 36 37 typedef struct { 38 RKTableau tableau; 39 Vec *Y; /* States computed during the step */ 40 Vec *YdotRHS; /* Function evaluations for the non-stiff part */ 41 PetscScalar *work; /* Scalar work */ 42 PetscReal stage_time; 43 TSStepStatus status; 44 } TS_RK; 45 46 /*MC 47 TSRK1 - First order forward Euler scheme. 48 49 This method has one stage. 50 51 Level: advanced 52 53 .seealso: TSRK 54 M*/ 55 /*MC 56 TSRK2A - Second order RK scheme. 57 58 This method has two stages. 59 60 Level: advanced 61 62 .seealso: TSRK 63 M*/ 64 /*MC 65 TSRK3 - Third order RK scheme. 66 67 This method has three stages. 68 69 Level: advanced 70 71 .seealso: TSRK 72 M*/ 73 /*MC 74 TSRK3BS - Third order RK scheme of Bogacki-Shampine with 2nd order embedded method. 75 76 This method has four stages. 77 78 Level: advanced 79 80 .seealso: TSRK 81 M*/ 82 /*MC 83 TSRK4 - Fourth order RK scheme. 84 85 This is the classical Runge-Kutta method with four stages. 86 87 Level: advanced 88 89 .seealso: TSRK 90 M*/ 91 /*MC 92 TSRK5F - Fifth order Fehlberg RK scheme with a 4th order embedded method. 93 94 This method has six stages. 95 96 Level: advanced 97 98 .seealso: TSRK 99 M*/ 100 /*MC 101 TSRK5DP - Fifth order Dormand-Prince RK scheme with the 4th order embedded method. 102 103 This method has seven stages. 104 105 Level: advanced 106 107 .seealso: TSRK 108 M*/ 109 110 #undef __FUNCT__ 111 #define __FUNCT__ "TSRKRegisterAll" 112 /*@C 113 TSRKRegisterAll - Registers all of the Runge-Kutta explicit methods in TSRK 114 115 Not Collective, but should be called by all processes which will need the schemes to be registered 116 117 Level: advanced 118 119 .keywords: TS, TSRK, register, all 120 121 .seealso: TSRKRegisterDestroy() 122 @*/ 123 PetscErrorCode TSRKRegisterAll(void) 124 { 125 PetscErrorCode ierr; 126 127 PetscFunctionBegin; 128 if (TSRKRegisterAllCalled) PetscFunctionReturn(0); 129 TSRKRegisterAllCalled = PETSC_TRUE; 130 131 { 132 const PetscReal 133 A[1][1] = {{0.0}}, 134 b[1] = {1.0}; 135 ierr = TSRKRegister(TSRK1FE,1,1,&A[0][0],b,NULL,NULL,1,b);CHKERRQ(ierr); 136 } 137 { 138 const PetscReal 139 A[2][2] = {{0.0,0.0}, 140 {1.0,0.0}}, 141 b[2] = {0.5,0.5}, 142 bembed[2] = {1.0,0}; 143 ierr = TSRKRegister(TSRK2A,2,2,&A[0][0],b,NULL,bembed,2,b);CHKERRQ(ierr); 144 } 145 { 146 const PetscReal 147 A[3][3] = {{0,0,0}, 148 {2.0/3.0,0,0}, 149 {-1.0/3.0,1.0,0}}, 150 b[3] = {0.25,0.5,0.25}; 151 ierr = TSRKRegister(TSRK3,3,3,&A[0][0],b,NULL,NULL,3,b);CHKERRQ(ierr); 152 } 153 { 154 const PetscReal 155 A[4][4] = {{0,0,0,0}, 156 {1.0/2.0,0,0,0}, 157 {0,3.0/4.0,0,0}, 158 {2.0/9.0,1.0/3.0,4.0/9.0,0}}, 159 b[4] = {2.0/9.0,1.0/3.0,4.0/9.0,0}, 160 bembed[4] = {7.0/24.0,1.0/4.0,1.0/3.0,1.0/8.0}; 161 ierr = TSRKRegister(TSRK3BS,3,4,&A[0][0],b,NULL,bembed,3,b);CHKERRQ(ierr); 162 } 163 { 164 const PetscReal 165 A[4][4] = {{0,0,0,0}, 166 {0.5,0,0,0}, 167 {0,0.5,0,0}, 168 {0,0,1.0,0}}, 169 b[4] = {1.0/6.0,1.0/3.0,1.0/3.0,1.0/6.0}; 170 ierr = TSRKRegister(TSRK4,4,4,&A[0][0],b,NULL,NULL,4,b);CHKERRQ(ierr); 171 } 172 { 173 const PetscReal 174 A[6][6] = {{0,0,0,0,0,0}, 175 {0.25,0,0,0,0,0}, 176 {3.0/32.0,9.0/32.0,0,0,0,0}, 177 {1932.0/2197.0,-7200.0/2197.0,7296.0/2197.0,0,0,0}, 178 {439.0/216.0,-8.0,3680.0/513.0,-845.0/4104.0,0,0}, 179 {-8.0/27.0,2.0,-3544.0/2565.0,1859.0/4104.0,-11.0/40.0,0}}, 180 b[6] = {16.0/135.0,0,6656.0/12825.0,28561.0/56430.0,-9.0/50.0,2.0/55.0}, 181 bembed[6] = {25.0/216.0,0,1408.0/2565.0,2197.0/4104.0,-1.0/5.0,0}; 182 ierr = TSRKRegister(TSRK5F,5,6,&A[0][0],b,NULL,bembed,5,b);CHKERRQ(ierr); 183 } 184 { 185 const PetscReal 186 A[7][7] = {{0,0,0,0,0,0,0}, 187 {1.0/5.0,0,0,0,0,0,0}, 188 {3.0/40.0,9.0/40.0,0,0,0,0,0}, 189 {44.0/45.0,-56.0/15.0,32.0/9.0,0,0,0,0}, 190 {19372.0/6561.0,-25360.0/2187.0,64448.0/6561.0,-212.0/729.0,0,0,0}, 191 {9017.0/3168.0,-355.0/33.0,46732.0/5247.0,49.0/176.0,-5103.0/18656.0,0,0}, 192 {35.0/384.0,0,500.0/1113.0,125.0/192.0,-2187.0/6784.0,11.0/84.0,0}}, 193 b[7] = {35.0/384.0,0,500.0/1113.0,125.0/192.0,-2187.0/6784.0,11.0/84.0,0}, 194 bembed[7] = {5179.0/57600.0,0,7571.0/16695.0,393.0/640.0,-92097.0/339200.0,187.0/2100.0,1.0/40.0}; 195 ierr = TSRKRegister(TSRK5DP,5,7,&A[0][0],b,NULL,bembed,5,b);CHKERRQ(ierr); 196 } 197 PetscFunctionReturn(0); 198 } 199 200 #undef __FUNCT__ 201 #define __FUNCT__ "TSRKRegisterDestroy" 202 /*@C 203 TSRKRegisterDestroy - Frees the list of schemes that were registered by TSRKRegister(). 204 205 Not Collective 206 207 Level: advanced 208 209 .keywords: TSRK, register, destroy 210 .seealso: TSRKRegister(), TSRKRegisterAll() 211 @*/ 212 PetscErrorCode TSRKRegisterDestroy(void) 213 { 214 PetscErrorCode ierr; 215 RKTableauLink link; 216 217 PetscFunctionBegin; 218 while ((link = RKTableauList)) { 219 RKTableau t = &link->tab; 220 RKTableauList = link->next; 221 ierr = PetscFree3(t->A,t->b,t->c); CHKERRQ(ierr); 222 ierr = PetscFree (t->bembed); CHKERRQ(ierr); 223 ierr = PetscFree (t->binterp); CHKERRQ(ierr); 224 ierr = PetscFree (t->name); CHKERRQ(ierr); 225 ierr = PetscFree (link); CHKERRQ(ierr); 226 } 227 TSRKRegisterAllCalled = PETSC_FALSE; 228 PetscFunctionReturn(0); 229 } 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "TSRKInitializePackage" 233 /*@C 234 TSRKInitializePackage - This function initializes everything in the TSRK package. It is called 235 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RK() 236 when using static libraries. 237 238 Level: developer 239 240 .keywords: TS, TSRK, initialize, package 241 .seealso: PetscInitialize() 242 @*/ 243 PetscErrorCode TSRKInitializePackage(void) 244 { 245 PetscErrorCode ierr; 246 247 PetscFunctionBegin; 248 if (TSRKPackageInitialized) PetscFunctionReturn(0); 249 TSRKPackageInitialized = PETSC_TRUE; 250 ierr = TSRKRegisterAll();CHKERRQ(ierr); 251 ierr = PetscObjectComposedDataRegister(&explicit_stage_time_id);CHKERRQ(ierr); 252 ierr = PetscRegisterFinalize(TSRKFinalizePackage);CHKERRQ(ierr); 253 PetscFunctionReturn(0); 254 } 255 256 #undef __FUNCT__ 257 #define __FUNCT__ "TSRKFinalizePackage" 258 /*@C 259 TSRKFinalizePackage - This function destroys everything in the TSRK package. It is 260 called from PetscFinalize(). 261 262 Level: developer 263 264 .keywords: Petsc, destroy, package 265 .seealso: PetscFinalize() 266 @*/ 267 PetscErrorCode TSRKFinalizePackage(void) 268 { 269 PetscErrorCode ierr; 270 271 PetscFunctionBegin; 272 TSRKPackageInitialized = PETSC_FALSE; 273 ierr = TSRKRegisterDestroy();CHKERRQ(ierr); 274 PetscFunctionReturn(0); 275 } 276 277 #undef __FUNCT__ 278 #define __FUNCT__ "TSRKRegister" 279 /*@C 280 TSRKRegister - register an RK scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 281 282 Not Collective, but the same schemes should be registered on all processes on which they will be used 283 284 Input Parameters: 285 + name - identifier for method 286 . order - approximation order of method 287 . s - number of stages, this is the dimension of the matrices below 288 . A - stage coefficients (dimension s*s, row-major) 289 . b - step completion table (dimension s; NULL to use last row of A) 290 . c - abscissa (dimension s; NULL to use row sums of A) 291 . bembed - completion table for embedded method (dimension s; NULL if not available) 292 . pinterp - Order of the interpolation scheme, equal to the number of columns of binterp 293 - binterp - Coefficients of the interpolation formula (dimension s*pinterp; NULL to reuse binterpt) 294 295 Notes: 296 Several RK methods are provided, this function is only needed to create new methods. 297 298 Level: advanced 299 300 .keywords: TS, register 301 302 .seealso: TSRK 303 @*/ 304 PetscErrorCode TSRKRegister(TSRKType name,PetscInt order,PetscInt s, 305 const PetscReal A[],const PetscReal b[],const PetscReal c[], 306 const PetscReal bembed[], 307 PetscInt pinterp,const PetscReal binterp[]) 308 { 309 PetscErrorCode ierr; 310 RKTableauLink link; 311 RKTableau t; 312 PetscInt i,j; 313 314 PetscFunctionBegin; 315 ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 316 ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr); 317 t = &link->tab; 318 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 319 t->order = order; 320 t->s = s; 321 ierr = PetscMalloc3(s*s,&t->A,s,&t->b,s,&t->c);CHKERRQ(ierr); 322 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 323 if (b) { ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); } 324 else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i]; 325 if (c) { ierr = PetscMemcpy(t->c,c,s*sizeof(c[0]));CHKERRQ(ierr); } 326 else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 327 t->FSAL = PETSC_TRUE; 328 for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE; 329 if (bembed) { 330 ierr = PetscMalloc1(s,&t->bembed);CHKERRQ(ierr); 331 ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr); 332 } 333 334 t->pinterp = pinterp; 335 ierr = PetscMalloc1(s*pinterp,&t->binterp);CHKERRQ(ierr); 336 ierr = PetscMemcpy(t->binterp,binterp,s*pinterp*sizeof(binterp[0]));CHKERRQ(ierr); 337 link->next = RKTableauList; 338 RKTableauList = link; 339 PetscFunctionReturn(0); 340 } 341 342 #undef __FUNCT__ 343 #define __FUNCT__ "TSEvaluateStep_RK" 344 /* 345 The step completion formula is 346 347 x1 = x0 + h b^T YdotRHS 348 349 This function can be called before or after ts->vec_sol has been updated. 350 Suppose we have a completion formula (b) and an embedded formula (be) of different order. 351 We can write 352 353 x1e = x0 + h be^T YdotRHS 354 = x1 - h b^T YdotRHS + h be^T YdotRHS 355 = x1 + h (be - b)^T YdotRHS 356 357 so we can evaluate the method with different order even after the step has been optimistically completed. 358 */ 359 static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done) 360 { 361 TS_RK *rk = (TS_RK*)ts->data; 362 RKTableau tab = rk->tableau; 363 PetscScalar *w = rk->work; 364 PetscReal h; 365 PetscInt s = tab->s,j; 366 PetscErrorCode ierr; 367 368 PetscFunctionBegin; 369 switch (rk->status) { 370 case TS_STEP_INCOMPLETE: 371 case TS_STEP_PENDING: 372 h = ts->time_step; break; 373 case TS_STEP_COMPLETE: 374 h = ts->time_step_prev; break; 375 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 376 } 377 if (order == tab->order) { 378 if (rk->status == TS_STEP_INCOMPLETE) { 379 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 380 for (j=0; j<s; j++) w[j] = h*tab->b[j]; 381 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 382 } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 383 if (done) *done = PETSC_TRUE; 384 PetscFunctionReturn(0); 385 } else if (order == tab->order-1) { 386 if (!tab->bembed) goto unavailable; 387 if (rk->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (be) */ 388 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 389 for (j=0; j<s; j++) w[j] = h*tab->bembed[j]; 390 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 391 } else { /* Rollback and re-complete using (be-b) */ 392 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 393 for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]); 394 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 395 } 396 if (done) *done = PETSC_TRUE; 397 PetscFunctionReturn(0); 398 } 399 unavailable: 400 if (done) *done = PETSC_FALSE; 401 else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"RK '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order); 402 PetscFunctionReturn(0); 403 } 404 405 #undef __FUNCT__ 406 #define __FUNCT__ "TSStep_RK" 407 static PetscErrorCode TSStep_RK(TS ts) 408 { 409 TS_RK *rk = (TS_RK*)ts->data; 410 RKTableau tab = rk->tableau; 411 const PetscInt s = tab->s; 412 const PetscReal *A = tab->A,*b = tab->b,*c = tab->c; 413 PetscScalar *w = rk->work; 414 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 415 TSAdapt adapt; 416 PetscInt i,j,reject,next_scheme; 417 PetscReal next_time_step; 418 PetscReal t; 419 PetscBool accept; 420 PetscErrorCode ierr; 421 422 PetscFunctionBegin; 423 424 next_time_step = ts->time_step; 425 t = ts->ptime; 426 accept = PETSC_TRUE; 427 rk->status = TS_STEP_INCOMPLETE; 428 429 430 for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 431 PetscReal h = ts->time_step; 432 ierr = TSPreStep(ts);CHKERRQ(ierr); 433 for (i=0; i<s; i++) { 434 rk->stage_time = t + h*c[i]; 435 ierr = TSPreStage(ts,rk->stage_time); CHKERRQ(ierr); 436 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 437 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 438 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 439 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 440 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 441 ierr = TSAdaptCheckStage(adapt,ts,&accept);CHKERRQ(ierr); 442 if (!accept) goto reject_step; 443 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 444 } 445 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 446 rk->status = TS_STEP_PENDING; 447 448 /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 449 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 450 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 451 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 452 ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 453 if (accept) { 454 /* ignore next_scheme for now */ 455 ts->ptime += ts->time_step; 456 ts->time_step = next_time_step; 457 ts->steps++; 458 rk->status = TS_STEP_COMPLETE; 459 ierr = PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);CHKERRQ(ierr); 460 break; 461 } else { /* Roll back the current step */ 462 for (j=0; j<s; j++) w[j] = -h*b[j]; 463 ierr = VecMAXPY(ts->vec_sol,s,w,rk->YdotRHS);CHKERRQ(ierr); 464 ts->time_step = next_time_step; 465 rk->status = TS_STEP_INCOMPLETE; 466 } 467 reject_step: continue; 468 } 469 if (rk->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 470 PetscFunctionReturn(0); 471 } 472 473 #undef __FUNCT__ 474 #define __FUNCT__ "TSStepAdj_RK" 475 static PetscErrorCode TSStepAdj_RK(TS tsadj) 476 { 477 PetscErrorCode ierr; 478 PetscFunctionBegin; 479 480 ierr = TSStep_RK(tsadj);CHKERRQ(ierr); 481 482 PetscFunctionReturn(0); 483 } 484 485 #undef __FUNCT__ 486 #define __FUNCT__ "TSInterpolate_RK" 487 static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X) 488 { 489 TS_RK *rk = (TS_RK*)ts->data; 490 PetscInt s = rk->tableau->s,pinterp = rk->tableau->pinterp,i,j; 491 PetscReal h; 492 PetscReal tt,t; 493 PetscScalar *b; 494 const PetscReal *B = rk->tableau->binterp; 495 PetscErrorCode ierr; 496 497 PetscFunctionBegin; 498 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 499 switch (rk->status) { 500 case TS_STEP_INCOMPLETE: 501 case TS_STEP_PENDING: 502 h = ts->time_step; 503 t = (itime - ts->ptime)/h; 504 break; 505 case TS_STEP_COMPLETE: 506 h = ts->time_step_prev; 507 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 508 break; 509 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 510 } 511 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 512 for (i=0; i<s; i++) b[i] = 0; 513 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 514 for (i=0; i<s; i++) { 515 b[i] += h * B[i*pinterp+j] * tt; 516 } 517 } 518 ierr = VecCopy(rk->Y[0],X);CHKERRQ(ierr); 519 ierr = VecMAXPY(X,s,b,rk->YdotRHS);CHKERRQ(ierr); 520 ierr = PetscFree(b);CHKERRQ(ierr); 521 PetscFunctionReturn(0); 522 } 523 524 /*------------------------------------------------------------*/ 525 #undef __FUNCT__ 526 #define __FUNCT__ "TSReset_RK" 527 static PetscErrorCode TSReset_RK(TS ts) 528 { 529 TS_RK *rk = (TS_RK*)ts->data; 530 PetscInt s; 531 PetscErrorCode ierr; 532 533 PetscFunctionBegin; 534 if (!rk->tableau) PetscFunctionReturn(0); 535 s = rk->tableau->s; 536 ierr = VecDestroyVecs(s,&rk->Y);CHKERRQ(ierr); 537 ierr = VecDestroyVecs(s,&rk->YdotRHS);CHKERRQ(ierr); 538 ierr = PetscFree(rk->work);CHKERRQ(ierr); 539 PetscFunctionReturn(0); 540 } 541 542 #undef __FUNCT__ 543 #define __FUNCT__ "TSDestroy_RK" 544 static PetscErrorCode TSDestroy_RK(TS ts) 545 { 546 PetscErrorCode ierr; 547 548 PetscFunctionBegin; 549 ierr = TSReset_RK(ts);CHKERRQ(ierr); 550 ierr = PetscFree(ts->data);CHKERRQ(ierr); 551 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);CHKERRQ(ierr); 552 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);CHKERRQ(ierr); 553 PetscFunctionReturn(0); 554 } 555 556 557 #undef __FUNCT__ 558 #define __FUNCT__ "DMCoarsenHook_TSRK" 559 static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx) 560 { 561 PetscFunctionBegin; 562 PetscFunctionReturn(0); 563 } 564 565 #undef __FUNCT__ 566 #define __FUNCT__ "DMRestrictHook_TSRK" 567 static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 568 { 569 PetscFunctionBegin; 570 PetscFunctionReturn(0); 571 } 572 573 574 #undef __FUNCT__ 575 #define __FUNCT__ "DMSubDomainHook_TSRK" 576 static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx) 577 { 578 PetscFunctionBegin; 579 PetscFunctionReturn(0); 580 } 581 582 #undef __FUNCT__ 583 #define __FUNCT__ "DMSubDomainRestrictHook_TSRK" 584 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 585 { 586 587 PetscFunctionBegin; 588 PetscFunctionReturn(0); 589 } 590 591 #undef __FUNCT__ 592 #define __FUNCT__ "RKSetAdjCoe" 593 static PetscErrorCode RKSetAdjCoe(RKTableau tab) 594 { 595 PetscReal *A,*b; 596 PetscInt s,i,j; 597 PetscErrorCode ierr; 598 599 PetscFunctionBegin; 600 s = tab->s; 601 ierr = PetscMalloc2(s*s,&A,s,&b);CHKERRQ(ierr); 602 603 for (i=0; i<s; i++) 604 for (j=0; j<s; j++) { 605 A[i*s+j] = (tab->b[s-1-i]==0) ? 0: -tab->A[s-1-i+(s-1-j)*s] * tab->b[s-1-j] / tab->b[s-1-i]; 606 ierr = PetscPrintf(PETSC_COMM_WORLD,"Coefficients: A[%D][%D]=%.6f\n",i,j,A[i*s+j]);CHKERRQ(ierr); 607 } 608 609 for (i=0; i<s; i++) b[i] = (tab->b[s-1-i]==0)? 0: -tab->b[s-1-i]; 610 611 ierr = PetscMemcpy(tab->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 612 ierr = PetscMemcpy(tab->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 613 ierr = PetscFree2(A,b);CHKERRQ(ierr); 614 PetscFunctionReturn(0); 615 } 616 617 #undef __FUNCT__ 618 #define __FUNCT__ "TSSetUp_RK" 619 static PetscErrorCode TSSetUp_RK(TS ts) 620 { 621 TS_RK *rk = (TS_RK*)ts->data; 622 RKTableau tab; 623 PetscInt s; 624 PetscErrorCode ierr; 625 DM dm; 626 627 PetscFunctionBegin; 628 if (!rk->tableau) { 629 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 630 } 631 tab = rk->tableau; 632 s = tab->s; 633 if (ts->reverse_mode) { 634 ierr = RKSetAdjCoe(tab);CHKERRQ(ierr); 635 } 636 ierr = VecDuplicateVecs(ts->vec_sol,s,&rk->Y);CHKERRQ(ierr); 637 ierr = VecDuplicateVecs(ts->vec_sol,s,&rk->YdotRHS);CHKERRQ(ierr); 638 ierr = PetscMalloc1(s,&rk->work);CHKERRQ(ierr); 639 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 640 if (dm) { 641 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 642 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 643 } 644 PetscFunctionReturn(0); 645 } 646 647 /*------------------------------------------------------------*/ 648 649 #undef __FUNCT__ 650 #define __FUNCT__ "TSSetFromOptions_RK" 651 static PetscErrorCode TSSetFromOptions_RK(PetscOptions *PetscOptionsObject,TS ts) 652 { 653 PetscErrorCode ierr; 654 char rktype[256]; 655 656 PetscFunctionBegin; 657 ierr = PetscOptionsHead(PetscOptionsObject,"RK ODE solver options");CHKERRQ(ierr); 658 { 659 RKTableauLink link; 660 PetscInt count,choice; 661 PetscBool flg; 662 const char **namelist; 663 ierr = PetscStrncpy(rktype,TSRKDefault,sizeof(rktype));CHKERRQ(ierr); 664 for (link=RKTableauList,count=0; link; link=link->next,count++) ; 665 ierr = PetscMalloc1(count,&namelist);CHKERRQ(ierr); 666 for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 667 ierr = PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rktype,&choice,&flg);CHKERRQ(ierr); 668 ierr = TSRKSetType(ts,flg ? namelist[choice] : rktype);CHKERRQ(ierr); 669 ierr = PetscFree(namelist);CHKERRQ(ierr); 670 } 671 ierr = PetscOptionsTail();CHKERRQ(ierr); 672 PetscFunctionReturn(0); 673 } 674 675 #undef __FUNCT__ 676 #define __FUNCT__ "PetscFormatRealArray" 677 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 678 { 679 PetscErrorCode ierr; 680 PetscInt i; 681 size_t left,count; 682 char *p; 683 684 PetscFunctionBegin; 685 for (i=0,p=buf,left=len; i<n; i++) { 686 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 687 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 688 left -= count; 689 p += count; 690 *p++ = ' '; 691 } 692 p[i ? 0 : -1] = 0; 693 PetscFunctionReturn(0); 694 } 695 696 #undef __FUNCT__ 697 #define __FUNCT__ "TSView_RK" 698 static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer) 699 { 700 TS_RK *rk = (TS_RK*)ts->data; 701 RKTableau tab = rk->tableau; 702 PetscBool iascii; 703 PetscErrorCode ierr; 704 TSAdapt adapt; 705 706 PetscFunctionBegin; 707 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 708 if (iascii) { 709 TSRKType rktype; 710 char buf[512]; 711 ierr = TSRKGetType(ts,&rktype);CHKERRQ(ierr); 712 ierr = PetscViewerASCIIPrintf(viewer," RK %s\n",rktype);CHKERRQ(ierr); 713 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 714 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 715 ierr = PetscViewerASCIIPrintf(viewer,"FSAL: %s\n",tab->FSAL ? "yes" : "no");CHKERRQ(ierr); 716 } 717 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 718 ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 719 PetscFunctionReturn(0); 720 } 721 722 #undef __FUNCT__ 723 #define __FUNCT__ "TSLoad_RK" 724 static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer) 725 { 726 PetscErrorCode ierr; 727 TSAdapt tsadapt; 728 729 PetscFunctionBegin; 730 ierr = TSGetAdapt(ts,&tsadapt);CHKERRQ(ierr); 731 ierr = TSAdaptLoad(tsadapt,viewer);CHKERRQ(ierr); 732 PetscFunctionReturn(0); 733 } 734 735 #undef __FUNCT__ 736 #define __FUNCT__ "TSRKSetType" 737 /*@C 738 TSRKSetType - Set the type of RK scheme 739 740 Logically collective 741 742 Input Parameter: 743 + ts - timestepping context 744 - rktype - type of RK-scheme 745 746 Level: intermediate 747 748 .seealso: TSRKGetType(), TSRK, TSRK2, TSRK3, TSRKPRSSP2, TSRK4, TSRK5 749 @*/ 750 PetscErrorCode TSRKSetType(TS ts,TSRKType rktype) 751 { 752 PetscErrorCode ierr; 753 754 PetscFunctionBegin; 755 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 756 ierr = PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));CHKERRQ(ierr); 757 PetscFunctionReturn(0); 758 } 759 760 #undef __FUNCT__ 761 #define __FUNCT__ "TSRKGetType" 762 /*@C 763 TSRKGetType - Get the type of RK scheme 764 765 Logically collective 766 767 Input Parameter: 768 . ts - timestepping context 769 770 Output Parameter: 771 . rktype - type of RK-scheme 772 773 Level: intermediate 774 775 .seealso: TSRKGetType() 776 @*/ 777 PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype) 778 { 779 PetscErrorCode ierr; 780 781 PetscFunctionBegin; 782 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 783 ierr = PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));CHKERRQ(ierr); 784 PetscFunctionReturn(0); 785 } 786 787 #undef __FUNCT__ 788 #define __FUNCT__ "TSRKGetType_RK" 789 PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype) 790 { 791 TS_RK *rk = (TS_RK*)ts->data; 792 PetscErrorCode ierr; 793 794 PetscFunctionBegin; 795 if (!rk->tableau) { 796 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 797 } 798 *rktype = rk->tableau->name; 799 PetscFunctionReturn(0); 800 } 801 #undef __FUNCT__ 802 #define __FUNCT__ "TSRKSetType_RK" 803 PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype) 804 { 805 TS_RK *rk = (TS_RK*)ts->data; 806 PetscErrorCode ierr; 807 PetscBool match; 808 RKTableauLink link; 809 810 PetscFunctionBegin; 811 if (rk->tableau) { 812 ierr = PetscStrcmp(rk->tableau->name,rktype,&match);CHKERRQ(ierr); 813 if (match) PetscFunctionReturn(0); 814 } 815 for (link = RKTableauList; link; link=link->next) { 816 ierr = PetscStrcmp(link->tab.name,rktype,&match);CHKERRQ(ierr); 817 if (match) { 818 ierr = TSReset_RK(ts);CHKERRQ(ierr); 819 rk->tableau = &link->tab; 820 PetscFunctionReturn(0); 821 } 822 } 823 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype); 824 PetscFunctionReturn(0); 825 } 826 827 #undef __FUNCT__ 828 #define __FUNCT__ "TSGetStages_RK" 829 static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y) 830 { 831 TS_RK *rk = (TS_RK*)ts->data; 832 833 PetscFunctionBegin; 834 *ns = rk->tableau->s; 835 if(Y) *Y = rk->Y; 836 PetscFunctionReturn(0); 837 } 838 839 840 /* ------------------------------------------------------------ */ 841 /*MC 842 TSRK - ODE and DAE solver using Runge-Kutta schemes 843 844 The user should provide the right hand side of the equation 845 using TSSetRHSFunction(). 846 847 Notes: 848 The default is TSRK3, it can be changed with TSRKSetType() or -ts_rk_type 849 850 Level: beginner 851 852 .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3, 853 TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister() 854 855 M*/ 856 #undef __FUNCT__ 857 #define __FUNCT__ "TSCreate_RK" 858 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 859 { 860 TS_RK *th; 861 PetscErrorCode ierr; 862 863 PetscFunctionBegin; 864 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 865 ierr = TSRKInitializePackage();CHKERRQ(ierr); 866 #endif 867 868 ts->ops->reset = TSReset_RK; 869 ts->ops->destroy = TSDestroy_RK; 870 ts->ops->view = TSView_RK; 871 ts->ops->load = TSLoad_RK; 872 ts->ops->setup = TSSetUp_RK; 873 ts->ops->step = TSStep_RK; 874 ts->ops->interpolate = TSInterpolate_RK; 875 ts->ops->evaluatestep = TSEvaluateStep_RK; 876 ts->ops->setfromoptions = TSSetFromOptions_RK; 877 ts->ops->getstages = TSGetStages_RK; 878 ts->ops->stepadj = TSStepAdj_RK; 879 880 ierr = PetscNewLog(ts,&th);CHKERRQ(ierr); 881 ts->data = (void*)th; 882 883 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);CHKERRQ(ierr); 884 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);CHKERRQ(ierr); 885 PetscFunctionReturn(0); 886 } 887