1 /* 2 Code for timestepping with Runge-Kutta method 3 4 Notes: 5 The general system is written as 6 7 F(t,U,Udot) = G(t,U) 8 9 where F represents the stiff part of the physics and G represents the non-stiff part. 10 11 */ 12 #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/ 13 #include <petscdm.h> 14 15 static TSRKType TSRKDefault = TSRK3; 16 static PetscBool TSRKRegisterAllCalled; 17 static PetscBool TSRKPackageInitialized; 18 static PetscInt explicit_stage_time_id; 19 20 typedef struct _RKTableau *RKTableau; 21 struct _RKTableau { 22 char *name; 23 PetscInt order; /* Classical approximation order of the method i */ 24 PetscInt s; /* Number of stages */ 25 PetscBool FSAL; /* flag to indicate if tableau is FSAL */ 26 PetscInt pinterp; /* Interpolation order */ 27 PetscReal *A,*b,*c; /* Tableau */ 28 PetscReal *bembed; /* Embedded formula of order one less (order-1) */ 29 PetscReal *binterp; /* Dense output formula */ 30 PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 31 }; 32 typedef struct _RKTableauLink *RKTableauLink; 33 struct _RKTableauLink { 34 struct _RKTableau tab; 35 RKTableauLink next; 36 }; 37 static RKTableauLink RKTableauList; 38 39 typedef struct { 40 RKTableau tableau; 41 Vec *Y; /* States computed during the step */ 42 Vec *YdotRHS; /* Function evaluations for the non-stiff part */ 43 PetscScalar *work; /* Scalar work */ 44 PetscReal stage_time; 45 TSStepStatus status; 46 } TS_RK; 47 48 /*MC 49 TSRK1 - First order forward Euler scheme. 50 51 This method has one stage. 52 53 Level: advanced 54 55 .seealso: TSRK 56 M*/ 57 /*MC 58 TSRK2A - Second order RK scheme. 59 60 This method has two stages. 61 62 Level: advanced 63 64 .seealso: TSRK 65 M*/ 66 /*MC 67 TSRK3 - Third order RK scheme. 68 69 This method has three stages. 70 71 Level: advanced 72 73 .seealso: TSRK 74 M*/ 75 /*MC 76 TSRK3BS - Third order RK scheme of Bogacki-Shampine with 2nd order embedded method. 77 78 This method has four stages. 79 80 Level: advanced 81 82 .seealso: TSRK 83 M*/ 84 /*MC 85 TSRK4 - Fourth order RK scheme. 86 87 This method has four stages. 88 89 Level: advanced 90 91 .seealso: TSRK 92 M*/ 93 /*MC 94 TSRK5F - Fifth order Fehlberg RK scheme with 4th order embedded method. 95 96 This method has six stages. 97 98 Level: advanced 99 100 .seealso: TSRK 101 M*/ 102 /*MC 103 TSRK5DP - Fifth order Dormand-Prince RK scheme with 4th order embedded method. 104 105 This method has seven stages. 106 107 Level: advanced 108 109 .seealso: TSRK 110 M*/ 111 112 #undef __FUNCT__ 113 #define __FUNCT__ "TSRKRegisterAll" 114 /*@C 115 TSRKRegisterAll - Registers all of the Runge-Kutta explicit methods in TSRK 116 117 Not Collective, but should be called by all processes which will need the schemes to be registered 118 119 Level: advanced 120 121 .keywords: TS, TSRK, register, all 122 123 .seealso: TSRKRegisterDestroy() 124 @*/ 125 PetscErrorCode TSRKRegisterAll(void) 126 { 127 PetscErrorCode ierr; 128 129 PetscFunctionBegin; 130 if (TSRKRegisterAllCalled) PetscFunctionReturn(0); 131 TSRKRegisterAllCalled = PETSC_TRUE; 132 133 { 134 const PetscReal 135 A[1][1] = {{0.0}}, 136 b[1] = {1.0}; 137 ierr = TSRKRegister(TSRK1FE,1,1,&A[0][0],b,NULL,NULL,1,b);CHKERRQ(ierr); 138 } 139 { 140 const PetscReal 141 A[2][2] = {{0.0,0.0}, 142 {1.0,0.0}}, 143 b[2] = {0.5,0.5}, 144 bembed[2] = {1.0,0}; 145 ierr = TSRKRegister(TSRK2A,2,2,&A[0][0],b,NULL,bembed,2,b);CHKERRQ(ierr); 146 } 147 { 148 const PetscReal 149 A[3][3] = {{0,0,0}, 150 {2.0/3.0,0,0}, 151 {-1.0/3.0,1.0,0}}, 152 b[3] = {0.25,0.5,0.25}; 153 ierr = TSRKRegister(TSRK3,3,3,&A[0][0],b,NULL,NULL,3,b);CHKERRQ(ierr); 154 } 155 { 156 const PetscReal 157 A[4][4] = {{0,0,0,0}, 158 {1.0/2.0,0,0,0}, 159 {0,3.0/4.0,0,0}, 160 {2.0/9.0,1.0/3.0,4.0/9.0,0}}, 161 b[4] = {2.0/9.0,1.0/3.0,4.0/9.0,0}, 162 bembed[4] = {7.0/24.0,1.0/4.0,1.0/3.0,1.0/8.0}; 163 ierr = TSRKRegister(TSRK3BS,3,4,&A[0][0],b,NULL,bembed,3,b);CHKERRQ(ierr); 164 } 165 { 166 const PetscReal 167 A[4][4] = {{0,0,0,0}, 168 {0.5,0,0,0}, 169 {0,0.5,0,0}, 170 {0,0,1.0,0}}, 171 b[4] = {1.0/6.0,1.0/3.0,1.0/3.0,1.0/6.0}; 172 ierr = TSRKRegister(TSRK4,4,4,&A[0][0],b,NULL,NULL,4,b);CHKERRQ(ierr); 173 } 174 { 175 const PetscReal 176 A[6][6] = {{0,0,0,0,0,0}, 177 {0.25,0,0,0,0,0}, 178 {3.0/32.0,9.0/32.0,0,0,0,0}, 179 {1932.0/2197.0,-7200.0/2197.0,7296.0/2197.0,0,0,0}, 180 {439.0/216.0,-8.0,3680.0/513.0,-845.0/4104.0,0,0}, 181 {-8.0/27.0,2.0,-3544.0/2565.0,1859.0/4104.0,-11.0/40.0,0}}, 182 b[6] = {16.0/135.0,0,6656.0/12825.0,28561.0/56430.0,-9.0/50.0,2.0/55.0}, 183 bembed[6] = {25.0/216.0,0,1408.0/2565.0,2197.0/4104.0,-1.0/5.0,0}; 184 ierr = TSRKRegister(TSRK5F,5,6,&A[0][0],b,NULL,bembed,5,b);CHKERRQ(ierr); 185 } 186 { 187 const PetscReal 188 A[7][7] = {{0,0,0,0,0,0,0}, 189 {1.0/5.0,0,0,0,0,0,0}, 190 {3.0/40.0,9.0/40.0,0,0,0,0,0}, 191 {44.0/45.0,-56.0/15.0,32.0/9.0,0,0,0,0}, 192 {19372.0/6561.0,-25360.0/2187.0,64448.0/6561.0,-212.0/729.0,0,0,0}, 193 {9017.0/3168.0,-355.0/33.0,46732.0/5247.0,49.0/176.0,-5103.0/18656.0,0,0}, 194 {35.0/384.0,0,500.0/1113.0,125.0/192.0,-2187.0/6784.0,11.0/84.0,0}}, 195 b[7] = {35.0/384.0,0,500.0/1113.0,125.0/192.0,-2187.0/6784.0,11.0/84.0,0}, 196 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}; 197 ierr = TSRKRegister(TSRK5DP,5,7,&A[0][0],b,NULL,bembed,5,b);CHKERRQ(ierr); 198 } 199 PetscFunctionReturn(0); 200 } 201 202 #undef __FUNCT__ 203 #define __FUNCT__ "TSRKRegisterDestroy" 204 /*@C 205 TSRKRegisterDestroy - Frees the list of schemes that were registered by TSRKRegister(). 206 207 Not Collective 208 209 Level: advanced 210 211 .keywords: TSRK, register, destroy 212 .seealso: TSRKRegister(), TSRKRegisterAll() 213 @*/ 214 PetscErrorCode TSRKRegisterDestroy(void) 215 { 216 PetscErrorCode ierr; 217 RKTableauLink link; 218 219 PetscFunctionBegin; 220 while ((link = RKTableauList)) { 221 RKTableau t = &link->tab; 222 RKTableauList = link->next; 223 ierr = PetscFree3(t->A,t->b,t->c); CHKERRQ(ierr); 224 ierr = PetscFree (t->bembed); CHKERRQ(ierr); 225 ierr = PetscFree (t->binterp); CHKERRQ(ierr); 226 ierr = PetscFree (t->name); CHKERRQ(ierr); 227 ierr = PetscFree (link); CHKERRQ(ierr); 228 } 229 TSRKRegisterAllCalled = PETSC_FALSE; 230 PetscFunctionReturn(0); 231 } 232 233 #undef __FUNCT__ 234 #define __FUNCT__ "TSRKInitializePackage" 235 /*@C 236 TSRKInitializePackage - This function initializes everything in the TSRK package. It is called 237 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RK() 238 when using static libraries. 239 240 Level: developer 241 242 .keywords: TS, TSRK, initialize, package 243 .seealso: PetscInitialize() 244 @*/ 245 PetscErrorCode TSRKInitializePackage(void) 246 { 247 PetscErrorCode ierr; 248 249 PetscFunctionBegin; 250 if (TSRKPackageInitialized) PetscFunctionReturn(0); 251 TSRKPackageInitialized = PETSC_TRUE; 252 ierr = TSRKRegisterAll();CHKERRQ(ierr); 253 ierr = PetscObjectComposedDataRegister(&explicit_stage_time_id);CHKERRQ(ierr); 254 ierr = PetscRegisterFinalize(TSRKFinalizePackage);CHKERRQ(ierr); 255 PetscFunctionReturn(0); 256 } 257 258 #undef __FUNCT__ 259 #define __FUNCT__ "TSRKFinalizePackage" 260 /*@C 261 TSRKFinalizePackage - This function destroys everything in the TSRK package. It is 262 called from PetscFinalize(). 263 264 Level: developer 265 266 .keywords: Petsc, destroy, package 267 .seealso: PetscFinalize() 268 @*/ 269 PetscErrorCode TSRKFinalizePackage(void) 270 { 271 PetscErrorCode ierr; 272 273 PetscFunctionBegin; 274 TSRKPackageInitialized = PETSC_FALSE; 275 ierr = TSRKRegisterDestroy();CHKERRQ(ierr); 276 PetscFunctionReturn(0); 277 } 278 279 #undef __FUNCT__ 280 #define __FUNCT__ "TSRKRegister" 281 /*@C 282 TSRKRegister - register an RK scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 283 284 Not Collective, but the same schemes should be registered on all processes on which they will be used 285 286 Input Parameters: 287 + name - identifier for method 288 . order - approximation order of method 289 . s - number of stages, this is the dimension of the matrices below 290 . A - stage coefficients (dimension s*s, row-major) 291 . b - step completion table (dimension s; NULL to use last row of A) 292 . c - abscissa (dimension s; NULL to use row sums of A) 293 . bembed - completion table for embedded method (dimension s; NULL if not available) 294 . pinterp - Order of the interpolation scheme, equal to the number of columns of binterp 295 - binterp - Coefficients of the interpolation formula (dimension s*pinterp; NULL to reuse binterpt) 296 297 Notes: 298 Several RK methods are provided, this function is only needed to create new methods. 299 300 Level: advanced 301 302 .keywords: TS, register 303 304 .seealso: TSRK 305 @*/ 306 PetscErrorCode TSRKRegister(TSRKType name,PetscInt order,PetscInt s, 307 const PetscReal A[],const PetscReal b[],const PetscReal c[], 308 const PetscReal bembed[], 309 PetscInt pinterp,const PetscReal binterp[]) 310 { 311 PetscErrorCode ierr; 312 RKTableauLink link; 313 RKTableau t; 314 PetscInt i,j; 315 316 PetscFunctionBegin; 317 ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 318 ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr); 319 t = &link->tab; 320 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 321 t->order = order; 322 t->s = s; 323 ierr = PetscMalloc3(s*s,PetscReal,&t->A,s,PetscReal,&t->b,s,PetscReal,&t->c);CHKERRQ(ierr); 324 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 325 if (b) { ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); } 326 else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i]; 327 if (c) { ierr = PetscMemcpy(t->c,c,s*sizeof(c[0]));CHKERRQ(ierr); } 328 else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 329 t->FSAL = PETSC_TRUE; 330 for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE; 331 if (bembed) { 332 ierr = PetscMalloc(s*sizeof(PetscReal),&t->bembed);CHKERRQ(ierr); 333 ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr); 334 } 335 336 t->pinterp = pinterp; 337 ierr = PetscMalloc(s*pinterp*sizeof(PetscReal),&t->binterp);CHKERRQ(ierr); 338 ierr = PetscMemcpy(t->binterp,binterp,s*pinterp*sizeof(binterp[0]));CHKERRQ(ierr); 339 link->next = RKTableauList; 340 RKTableauList = link; 341 PetscFunctionReturn(0); 342 } 343 344 #undef __FUNCT__ 345 #define __FUNCT__ "TSEvaluateStep_RK" 346 /* 347 The step completion formula is 348 349 x1 = x0 + h b^T YdotRHS 350 351 This function can be called before or after ts->vec_sol has been updated. 352 Suppose we have a completion formula (b) and an embedded formula (be) of different order. 353 We can write 354 355 x1e = x0 + h be^T YdotRHS 356 = x1 - h b^T YdotRHS + h be^T YdotRHS 357 = x1 + h (be - b)^T YdotRHS 358 359 so we can evaluate the method with different order even after the step has been optimistically completed. 360 */ 361 static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done) 362 { 363 TS_RK *rk = (TS_RK*)ts->data; 364 RKTableau tab = rk->tableau; 365 PetscScalar *w = rk->work; 366 PetscReal h; 367 PetscInt s = tab->s,j; 368 PetscErrorCode ierr; 369 370 PetscFunctionBegin; 371 switch (rk->status) { 372 case TS_STEP_INCOMPLETE: 373 case TS_STEP_PENDING: 374 h = ts->time_step; break; 375 case TS_STEP_COMPLETE: 376 h = ts->time_step_prev; break; 377 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 378 } 379 if (order == tab->order) { 380 if (rk->status == TS_STEP_INCOMPLETE) { 381 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 382 for (j=0; j<s; j++) w[j] = h*tab->b[j]; 383 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 384 } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 385 if (done) *done = PETSC_TRUE; 386 PetscFunctionReturn(0); 387 } else if (order == tab->order-1) { 388 if (!tab->bembed) goto unavailable; 389 if (rk->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (be) */ 390 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 391 for (j=0; j<s; j++) w[j] = h*tab->bembed[j]; 392 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 393 } else { /* Rollback and re-complete using (be-b) */ 394 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 395 for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]); 396 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 397 } 398 if (done) *done = PETSC_TRUE; 399 PetscFunctionReturn(0); 400 } 401 unavailable: 402 if (done) *done = PETSC_FALSE; 403 else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"RK '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order); 404 PetscFunctionReturn(0); 405 } 406 407 #undef __FUNCT__ 408 #define __FUNCT__ "TSStep_RK" 409 static PetscErrorCode TSStep_RK(TS ts) 410 { 411 TS_RK *rk = (TS_RK*)ts->data; 412 RKTableau tab = rk->tableau; 413 const PetscInt s = tab->s; 414 const PetscReal *A = tab->A,*b = tab->b,*c = tab->c; 415 PetscScalar *w = rk->work; 416 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 417 TSAdapt adapt; 418 PetscInt i,j,reject,next_scheme; 419 PetscReal next_time_step; 420 PetscReal t; 421 PetscBool accept; 422 PetscErrorCode ierr; 423 424 PetscFunctionBegin; 425 426 next_time_step = ts->time_step; 427 t = ts->ptime; 428 accept = PETSC_TRUE; 429 rk->status = TS_STEP_INCOMPLETE; 430 431 432 for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 433 PetscReal h = ts->time_step; 434 ierr = TSPreStep(ts);CHKERRQ(ierr); 435 for (i=0; i<s; i++) { 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 = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 440 ierr = TSAdaptCheckStage(adapt,ts,&accept);CHKERRQ(ierr); 441 if (!accept) goto reject_step; 442 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 443 } 444 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 445 rk->status = TS_STEP_PENDING; 446 447 /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 448 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 449 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 450 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 451 ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 452 if (accept) { 453 /* ignore next_scheme for now */ 454 ts->ptime += ts->time_step; 455 ts->time_step = next_time_step; 456 ts->steps++; 457 rk->status = TS_STEP_COMPLETE; 458 ierr = PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);CHKERRQ(ierr); 459 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 = PetscMalloc(s*sizeof(PetscScalar),&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 = PetscMalloc(s*sizeof(rk->work[0]),&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 = PetscMalloc(count*sizeof(char*),&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,TS_RK,&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