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__ "TSInterpolate_RK" 475 static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X) 476 { 477 TS_RK *rk = (TS_RK*)ts->data; 478 PetscInt s = rk->tableau->s,pinterp = rk->tableau->pinterp,i,j; 479 PetscReal h; 480 PetscReal tt,t; 481 PetscScalar *b; 482 const PetscReal *B = rk->tableau->binterp; 483 PetscErrorCode ierr; 484 485 PetscFunctionBegin; 486 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 487 switch (rk->status) { 488 case TS_STEP_INCOMPLETE: 489 case TS_STEP_PENDING: 490 h = ts->time_step; 491 t = (itime - ts->ptime)/h; 492 break; 493 case TS_STEP_COMPLETE: 494 h = ts->time_step_prev; 495 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 496 break; 497 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 498 } 499 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 500 for (i=0; i<s; i++) b[i] = 0; 501 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 502 for (i=0; i<s; i++) { 503 b[i] += h * B[i*pinterp+j] * tt; 504 } 505 } 506 ierr = VecCopy(rk->Y[0],X);CHKERRQ(ierr); 507 ierr = VecMAXPY(X,s,b,rk->YdotRHS);CHKERRQ(ierr); 508 ierr = PetscFree(b);CHKERRQ(ierr); 509 PetscFunctionReturn(0); 510 } 511 512 /*------------------------------------------------------------*/ 513 #undef __FUNCT__ 514 #define __FUNCT__ "TSReset_RK" 515 static PetscErrorCode TSReset_RK(TS ts) 516 { 517 TS_RK *rk = (TS_RK*)ts->data; 518 PetscInt s; 519 PetscErrorCode ierr; 520 521 PetscFunctionBegin; 522 if (!rk->tableau) PetscFunctionReturn(0); 523 s = rk->tableau->s; 524 ierr = VecDestroyVecs(s,&rk->Y);CHKERRQ(ierr); 525 ierr = VecDestroyVecs(s,&rk->YdotRHS);CHKERRQ(ierr); 526 ierr = PetscFree(rk->work);CHKERRQ(ierr); 527 PetscFunctionReturn(0); 528 } 529 530 #undef __FUNCT__ 531 #define __FUNCT__ "TSDestroy_RK" 532 static PetscErrorCode TSDestroy_RK(TS ts) 533 { 534 PetscErrorCode ierr; 535 536 PetscFunctionBegin; 537 ierr = TSReset_RK(ts);CHKERRQ(ierr); 538 ierr = PetscFree(ts->data);CHKERRQ(ierr); 539 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);CHKERRQ(ierr); 540 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);CHKERRQ(ierr); 541 PetscFunctionReturn(0); 542 } 543 544 545 #undef __FUNCT__ 546 #define __FUNCT__ "DMCoarsenHook_TSRK" 547 static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx) 548 { 549 PetscFunctionBegin; 550 PetscFunctionReturn(0); 551 } 552 553 #undef __FUNCT__ 554 #define __FUNCT__ "DMRestrictHook_TSRK" 555 static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 556 { 557 PetscFunctionBegin; 558 PetscFunctionReturn(0); 559 } 560 561 562 #undef __FUNCT__ 563 #define __FUNCT__ "DMSubDomainHook_TSRK" 564 static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx) 565 { 566 PetscFunctionBegin; 567 PetscFunctionReturn(0); 568 } 569 570 #undef __FUNCT__ 571 #define __FUNCT__ "DMSubDomainRestrictHook_TSRK" 572 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 573 { 574 575 PetscFunctionBegin; 576 PetscFunctionReturn(0); 577 } 578 579 #undef __FUNCT__ 580 #define __FUNCT__ "TSSetUp_RK" 581 static PetscErrorCode TSSetUp_RK(TS ts) 582 { 583 TS_RK *rk = (TS_RK*)ts->data; 584 RKTableau tab; 585 PetscInt s; 586 PetscErrorCode ierr; 587 DM dm; 588 589 PetscFunctionBegin; 590 if (!rk->tableau) { 591 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 592 } 593 tab = rk->tableau; 594 s = tab->s; 595 ierr = VecDuplicateVecs(ts->vec_sol,s,&rk->Y);CHKERRQ(ierr); 596 ierr = VecDuplicateVecs(ts->vec_sol,s,&rk->YdotRHS);CHKERRQ(ierr); 597 ierr = PetscMalloc1(s,&rk->work);CHKERRQ(ierr); 598 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 599 if (dm) { 600 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 601 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 602 } 603 PetscFunctionReturn(0); 604 } 605 /*------------------------------------------------------------*/ 606 607 #undef __FUNCT__ 608 #define __FUNCT__ "TSSetFromOptions_RK" 609 static PetscErrorCode TSSetFromOptions_RK(TS ts) 610 { 611 PetscErrorCode ierr; 612 char rktype[256]; 613 614 PetscFunctionBegin; 615 ierr = PetscOptionsHead("RK ODE solver options");CHKERRQ(ierr); 616 { 617 RKTableauLink link; 618 PetscInt count,choice; 619 PetscBool flg; 620 const char **namelist; 621 ierr = PetscStrncpy(rktype,TSRKDefault,sizeof(rktype));CHKERRQ(ierr); 622 for (link=RKTableauList,count=0; link; link=link->next,count++) ; 623 ierr = PetscMalloc1(count,&namelist);CHKERRQ(ierr); 624 for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 625 ierr = PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rktype,&choice,&flg);CHKERRQ(ierr); 626 ierr = TSRKSetType(ts,flg ? namelist[choice] : rktype);CHKERRQ(ierr); 627 ierr = PetscFree(namelist);CHKERRQ(ierr); 628 } 629 ierr = PetscOptionsTail();CHKERRQ(ierr); 630 PetscFunctionReturn(0); 631 } 632 633 #undef __FUNCT__ 634 #define __FUNCT__ "PetscFormatRealArray" 635 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 636 { 637 PetscErrorCode ierr; 638 PetscInt i; 639 size_t left,count; 640 char *p; 641 642 PetscFunctionBegin; 643 for (i=0,p=buf,left=len; i<n; i++) { 644 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 645 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 646 left -= count; 647 p += count; 648 *p++ = ' '; 649 } 650 p[i ? 0 : -1] = 0; 651 PetscFunctionReturn(0); 652 } 653 654 #undef __FUNCT__ 655 #define __FUNCT__ "TSView_RK" 656 static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer) 657 { 658 TS_RK *rk = (TS_RK*)ts->data; 659 RKTableau tab = rk->tableau; 660 PetscBool iascii; 661 PetscErrorCode ierr; 662 TSAdapt adapt; 663 664 PetscFunctionBegin; 665 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 666 if (iascii) { 667 TSRKType rktype; 668 char buf[512]; 669 ierr = TSRKGetType(ts,&rktype);CHKERRQ(ierr); 670 ierr = PetscViewerASCIIPrintf(viewer," RK %s\n",rktype);CHKERRQ(ierr); 671 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 672 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 673 ierr = PetscViewerASCIIPrintf(viewer,"FSAL: %s\n",tab->FSAL ? "yes" : "no");CHKERRQ(ierr); 674 } 675 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 676 ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 677 PetscFunctionReturn(0); 678 } 679 680 #undef __FUNCT__ 681 #define __FUNCT__ "TSLoad_RK" 682 static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer) 683 { 684 PetscErrorCode ierr; 685 TSAdapt tsadapt; 686 687 PetscFunctionBegin; 688 ierr = TSGetAdapt(ts,&tsadapt);CHKERRQ(ierr); 689 ierr = TSAdaptLoad(tsadapt,viewer);CHKERRQ(ierr); 690 PetscFunctionReturn(0); 691 } 692 693 #undef __FUNCT__ 694 #define __FUNCT__ "TSRKSetType" 695 /*@C 696 TSRKSetType - Set the type of RK scheme 697 698 Logically collective 699 700 Input Parameter: 701 + ts - timestepping context 702 - rktype - type of RK-scheme 703 704 Level: intermediate 705 706 .seealso: TSRKGetType(), TSRK, TSRK2, TSRK3, TSRKPRSSP2, TSRK4, TSRK5 707 @*/ 708 PetscErrorCode TSRKSetType(TS ts,TSRKType rktype) 709 { 710 PetscErrorCode ierr; 711 712 PetscFunctionBegin; 713 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 714 ierr = PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));CHKERRQ(ierr); 715 PetscFunctionReturn(0); 716 } 717 718 #undef __FUNCT__ 719 #define __FUNCT__ "TSRKGetType" 720 /*@C 721 TSRKGetType - Get the type of RK scheme 722 723 Logically collective 724 725 Input Parameter: 726 . ts - timestepping context 727 728 Output Parameter: 729 . rktype - type of RK-scheme 730 731 Level: intermediate 732 733 .seealso: TSRKGetType() 734 @*/ 735 PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype) 736 { 737 PetscErrorCode ierr; 738 739 PetscFunctionBegin; 740 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 741 ierr = PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));CHKERRQ(ierr); 742 PetscFunctionReturn(0); 743 } 744 745 #undef __FUNCT__ 746 #define __FUNCT__ "TSRKGetType_RK" 747 PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype) 748 { 749 TS_RK *rk = (TS_RK*)ts->data; 750 PetscErrorCode ierr; 751 752 PetscFunctionBegin; 753 if (!rk->tableau) { 754 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 755 } 756 *rktype = rk->tableau->name; 757 PetscFunctionReturn(0); 758 } 759 #undef __FUNCT__ 760 #define __FUNCT__ "TSRKSetType_RK" 761 PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype) 762 { 763 TS_RK *rk = (TS_RK*)ts->data; 764 PetscErrorCode ierr; 765 PetscBool match; 766 RKTableauLink link; 767 768 PetscFunctionBegin; 769 if (rk->tableau) { 770 ierr = PetscStrcmp(rk->tableau->name,rktype,&match);CHKERRQ(ierr); 771 if (match) PetscFunctionReturn(0); 772 } 773 for (link = RKTableauList; link; link=link->next) { 774 ierr = PetscStrcmp(link->tab.name,rktype,&match);CHKERRQ(ierr); 775 if (match) { 776 ierr = TSReset_RK(ts);CHKERRQ(ierr); 777 rk->tableau = &link->tab; 778 PetscFunctionReturn(0); 779 } 780 } 781 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype); 782 PetscFunctionReturn(0); 783 } 784 785 /* ------------------------------------------------------------ */ 786 /*MC 787 TSRK - ODE and DAE solver using Runge-Kutta schemes 788 789 The user should provide the right hand side of the equation 790 using TSSetRHSFunction(). 791 792 Notes: 793 The default is TSRK3, it can be changed with TSRKSetType() or -ts_rk_type 794 795 Level: beginner 796 797 .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3, 798 TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister() 799 800 M*/ 801 #undef __FUNCT__ 802 #define __FUNCT__ "TSCreate_RK" 803 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 804 { 805 TS_RK *th; 806 PetscErrorCode ierr; 807 808 PetscFunctionBegin; 809 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 810 ierr = TSRKInitializePackage();CHKERRQ(ierr); 811 #endif 812 813 ts->ops->reset = TSReset_RK; 814 ts->ops->destroy = TSDestroy_RK; 815 ts->ops->view = TSView_RK; 816 ts->ops->load = TSLoad_RK; 817 ts->ops->setup = TSSetUp_RK; 818 ts->ops->step = TSStep_RK; 819 ts->ops->interpolate = TSInterpolate_RK; 820 ts->ops->evaluatestep = TSEvaluateStep_RK; 821 ts->ops->setfromoptions = TSSetFromOptions_RK; 822 823 ierr = PetscNewLog(ts,&th);CHKERRQ(ierr); 824 ts->data = (void*)th; 825 826 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);CHKERRQ(ierr); 827 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);CHKERRQ(ierr); 828 PetscFunctionReturn(0); 829 } 830