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