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 17 typedef struct _RKTableau *RKTableau; 18 struct _RKTableau { 19 char *name; 20 PetscInt order; /* Classical approximation order of the method i */ 21 PetscInt s; /* Number of stages */ 22 PetscInt p; /* Interpolation order */ 23 PetscBool FSAL; /* flag to indicate if tableau is FSAL */ 24 PetscReal *A,*b,*c; /* Tableau */ 25 PetscReal *bembed; /* Embedded formula of order one less (order-1) */ 26 PetscReal *binterp; /* Dense output formula */ 27 PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 28 }; 29 typedef struct _RKTableauLink *RKTableauLink; 30 struct _RKTableauLink { 31 struct _RKTableau tab; 32 RKTableauLink next; 33 }; 34 static RKTableauLink RKTableauList; 35 36 typedef struct { 37 RKTableau tableau; 38 Vec *Y; /* States computed during the step */ 39 Vec *YdotRHS; /* Function evaluations for the non-stiff part */ 40 Vec *VecDeltaLam; /* Increment of the adjoint sensitivity w.r.t IC at stage */ 41 Vec *VecDeltaMu; /* Increment of the adjoint sensitivity w.r.t P at stage */ 42 Vec *VecSensiTemp; /* Vector to be timed with Jacobian transpose */ 43 Vec VecCostIntegral0; /* backup for roll-backs due to events */ 44 PetscScalar *work; /* Scalar work */ 45 PetscReal stage_time; 46 TSStepStatus status; 47 PetscReal ptime; 48 PetscReal time_step; 49 } TS_RK; 50 51 /*MC 52 TSRK1 - First order forward Euler scheme. 53 54 This method has one stage. 55 56 Level: advanced 57 58 .seealso: TSRK 59 M*/ 60 /*MC 61 TSRK2A - Second order RK scheme. 62 63 This method has two stages. 64 65 Level: advanced 66 67 .seealso: TSRK 68 M*/ 69 /*MC 70 TSRK3 - Third order RK scheme. 71 72 This method has three stages. 73 74 Level: advanced 75 76 .seealso: TSRK 77 M*/ 78 /*MC 79 TSRK3BS - Third order RK scheme of Bogacki-Shampine with 2nd order embedded method. 80 81 This method has four stages with the First Same As Last (FSAL) property. 82 83 Level: advanced 84 85 References: https://doi.org/10.1016/0893-9659(89)90079-7 86 87 .seealso: TSRK 88 M*/ 89 /*MC 90 TSRK4 - Fourth order RK scheme. 91 92 This is the classical Runge-Kutta method with four stages. 93 94 Level: advanced 95 96 .seealso: TSRK 97 M*/ 98 /*MC 99 TSRK5F - Fifth order Fehlberg RK scheme with a 4th order embedded method. 100 101 This method has six stages. 102 103 Level: advanced 104 105 .seealso: TSRK 106 M*/ 107 /*MC 108 TSRK5DP - Fifth order Dormand-Prince RK scheme with the 4th order embedded method. 109 110 This method has seven stages with the First Same As Last (FSAL) property. 111 112 Level: advanced 113 114 References: https://doi.org/10.1016/0771-050X(80)90013-3 115 116 .seealso: TSRK 117 M*/ 118 /*MC 119 TSRK5BS - Fifth order Bogacki-Shampine RK scheme with 4th order embedded method. 120 121 This method has eight stages with the First Same As Last (FSAL) property. 122 123 Level: advanced 124 125 References: https://doi.org/10.1016/0898-1221(96)00141-1 126 127 .seealso: TSRK 128 M*/ 129 130 /*@C 131 TSRKRegisterAll - Registers all of the Runge-Kutta explicit methods in TSRK 132 133 Not Collective, but should be called by all processes which will need the schemes to be registered 134 135 Level: advanced 136 137 .keywords: TS, TSRK, register, all 138 139 .seealso: TSRKRegisterDestroy() 140 @*/ 141 PetscErrorCode TSRKRegisterAll(void) 142 { 143 PetscErrorCode ierr; 144 145 PetscFunctionBegin; 146 if (TSRKRegisterAllCalled) PetscFunctionReturn(0); 147 TSRKRegisterAllCalled = PETSC_TRUE; 148 149 #define RC PetscRealConstant 150 { 151 const PetscReal 152 A[1][1] = {{0}}, 153 b[1] = {RC(1.0)}; 154 ierr = TSRKRegister(TSRK1FE,1,1,&A[0][0],b,NULL,NULL,0,NULL);CHKERRQ(ierr); 155 } 156 { 157 const PetscReal 158 A[2][2] = {{0,0}, 159 {RC(1.0),0}}, 160 b[2] = {RC(0.5),RC(0.5)}, 161 bembed[2] = {RC(1.0),0}; 162 ierr = TSRKRegister(TSRK2A,2,2,&A[0][0],b,NULL,bembed,0,NULL);CHKERRQ(ierr); 163 } 164 { 165 const PetscReal 166 A[3][3] = {{0,0,0}, 167 {RC(2.0)/RC(3.0),0,0}, 168 {RC(-1.0)/RC(3.0),RC(1.0),0}}, 169 b[3] = {RC(0.25),RC(0.5),RC(0.25)}; 170 ierr = TSRKRegister(TSRK3,3,3,&A[0][0],b,NULL,NULL,0,NULL);CHKERRQ(ierr); 171 } 172 { 173 const PetscReal 174 A[4][4] = {{0,0,0,0}, 175 {RC(1.0)/RC(2.0),0,0,0}, 176 {0,RC(3.0)/RC(4.0),0,0}, 177 {RC(2.0)/RC(9.0),RC(1.0)/RC(3.0),RC(4.0)/RC(9.0),0}}, 178 b[4] = {RC(2.0)/RC(9.0),RC(1.0)/RC(3.0),RC(4.0)/RC(9.0),0}, 179 bembed[4] = {RC(7.0)/RC(24.0),RC(1.0)/RC(4.0),RC(1.0)/RC(3.0),RC(1.0)/RC(8.0)}; 180 ierr = TSRKRegister(TSRK3BS,3,4,&A[0][0],b,NULL,bembed,0,NULL);CHKERRQ(ierr); 181 } 182 { 183 const PetscReal 184 A[4][4] = {{0,0,0,0}, 185 {RC(0.5),0,0,0}, 186 {0,RC(0.5),0,0}, 187 {0,0,RC(1.0),0}}, 188 b[4] = {RC(1.0)/RC(6.0),RC(1.0)/RC(3.0),RC(1.0)/RC(3.0),RC(1.0)/RC(6.0)}; 189 ierr = TSRKRegister(TSRK4,4,4,&A[0][0],b,NULL,NULL,0,NULL);CHKERRQ(ierr); 190 } 191 { 192 const PetscReal 193 A[6][6] = {{0,0,0,0,0,0}, 194 {RC(0.25),0,0,0,0,0}, 195 {RC(3.0)/RC(32.0),RC(9.0)/RC(32.0),0,0,0,0}, 196 {RC(1932.0)/RC(2197.0),RC(-7200.0)/RC(2197.0),RC(7296.0)/RC(2197.0),0,0,0}, 197 {RC(439.0)/RC(216.0),RC(-8.0),RC(3680.0)/RC(513.0),RC(-845.0)/RC(4104.0),0,0}, 198 {RC(-8.0)/RC(27.0),RC(2.0),RC(-3544.0)/RC(2565.0),RC(1859.0)/RC(4104.0),RC(-11.0)/RC(40.0),0}}, 199 b[6] = {RC(16.0)/RC(135.0),0,RC(6656.0)/RC(12825.0),RC(28561.0)/RC(56430.0),RC(-9.0)/RC(50.0),RC(2.0)/RC(55.0)}, 200 bembed[6] = {RC(25.0)/RC(216.0),0,RC(1408.0)/RC(2565.0),RC(2197.0)/RC(4104.0),RC(-1.0)/RC(5.0),0}; 201 ierr = TSRKRegister(TSRK5F,5,6,&A[0][0],b,NULL,bembed,0,NULL);CHKERRQ(ierr); 202 } 203 { 204 const PetscReal 205 A[7][7] = {{0,0,0,0,0,0,0}, 206 {RC(1.0)/RC(5.0),0,0,0,0,0,0}, 207 {RC(3.0)/RC(40.0),RC(9.0)/RC(40.0),0,0,0,0,0}, 208 {RC(44.0)/RC(45.0),RC(-56.0)/RC(15.0),RC(32.0)/RC(9.0),0,0,0,0}, 209 {RC(19372.0)/RC(6561.0),RC(-25360.0)/RC(2187.0),RC(64448.0)/RC(6561.0),RC(-212.0)/RC(729.0),0,0,0}, 210 {RC(9017.0)/RC(3168.0),RC(-355.0)/RC(33.0),RC(46732.0)/RC(5247.0),RC(49.0)/RC(176.0),RC(-5103.0)/RC(18656.0),0,0}, 211 {RC(35.0)/RC(384.0),0,RC(500.0)/RC(1113.0),RC(125.0)/RC(192.0),RC(-2187.0)/RC(6784.0),RC(11.0)/RC(84.0),0}}, 212 b[7] = {RC(35.0)/RC(384.0),0,RC(500.0)/RC(1113.0),RC(125.0)/RC(192.0),RC(-2187.0)/RC(6784.0),RC(11.0)/RC(84.0),0}, 213 bembed[7] = {RC(5179.0)/RC(57600.0),0,RC(7571.0)/RC(16695.0),RC(393.0)/RC(640.0),RC(-92097.0)/RC(339200.0),RC(187.0)/RC(2100.0),RC(1.0)/RC(40.0)}; 214 ierr = TSRKRegister(TSRK5DP,5,7,&A[0][0],b,NULL,bembed,0,NULL);CHKERRQ(ierr); 215 } 216 { 217 const PetscReal 218 A[8][8] = {{0,0,0,0,0,0,0,0}, 219 {RC(1.0)/RC(6.0),0,0,0,0,0,0,0}, 220 {RC(2.0)/RC(27.0),RC(4.0)/RC(27.0),0,0,0,0,0,0}, 221 {RC(183.0)/RC(1372.0),RC(-162.0)/RC(343.0),RC(1053.0)/RC(1372.0),0,0,0,0,0}, 222 {RC(68.0)/RC(297.0),RC(-4.0)/RC(11.0),RC(42.0)/RC(143.0),RC(1960.0)/RC(3861.0),0,0,0,0}, 223 {RC(597.0)/RC(22528.0),RC(81.0)/RC(352.0),RC(63099.0)/RC(585728.0),RC(58653.0)/RC(366080.0),RC(4617.0)/RC(20480.0),0,0,0}, 224 {RC(174197.0)/RC(959244.0),RC(-30942.0)/RC(79937.0),RC(8152137.0)/RC(19744439.0),RC(666106.0)/RC(1039181.0),RC(-29421.0)/RC(29068.0),RC(482048.0)/RC(414219.0),0,0}, 225 {RC(587.0)/RC(8064.0),0,RC(4440339.0)/RC(15491840.0),RC(24353.0)/RC(124800.0),RC(387.0)/RC(44800.0),RC(2152.0)/RC(5985.0),RC(7267.0)/RC(94080.0),0}}, 226 b[8] = {RC(587.0)/RC(8064.0),0,RC(4440339.0)/RC(15491840.0),RC(24353.0)/RC(124800.0),RC(387.0)/RC(44800.0),RC(2152.0)/RC(5985.0),RC(7267.0)/RC(94080.0),0}, 227 bembed[8] = {RC(2479.0)/RC(34992.0),0,RC(123.0)/RC(416.0),RC(612941.0)/RC(3411720.0),RC(43.0)/RC(1440.0),RC(2272.0)/RC(6561.0),RC(79937.0)/RC(1113912.0),RC(3293.0)/RC(556956.0)}; 228 ierr = TSRKRegister(TSRK5BS,5,8,&A[0][0],b,NULL,bembed,0,NULL);CHKERRQ(ierr); 229 } 230 #undef RC 231 PetscFunctionReturn(0); 232 } 233 234 /*@C 235 TSRKRegisterDestroy - Frees the list of schemes that were registered by TSRKRegister(). 236 237 Not Collective 238 239 Level: advanced 240 241 .keywords: TSRK, register, destroy 242 .seealso: TSRKRegister(), TSRKRegisterAll() 243 @*/ 244 PetscErrorCode TSRKRegisterDestroy(void) 245 { 246 PetscErrorCode ierr; 247 RKTableauLink link; 248 249 PetscFunctionBegin; 250 while ((link = RKTableauList)) { 251 RKTableau t = &link->tab; 252 RKTableauList = link->next; 253 ierr = PetscFree3(t->A,t->b,t->c); CHKERRQ(ierr); 254 ierr = PetscFree (t->bembed); CHKERRQ(ierr); 255 ierr = PetscFree (t->binterp); CHKERRQ(ierr); 256 ierr = PetscFree (t->name); CHKERRQ(ierr); 257 ierr = PetscFree (link); CHKERRQ(ierr); 258 } 259 TSRKRegisterAllCalled = PETSC_FALSE; 260 PetscFunctionReturn(0); 261 } 262 263 /*@C 264 TSRKInitializePackage - This function initializes everything in the TSRK package. It is called 265 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RK() 266 when using static libraries. 267 268 Level: developer 269 270 .keywords: TS, TSRK, initialize, package 271 .seealso: PetscInitialize() 272 @*/ 273 PetscErrorCode TSRKInitializePackage(void) 274 { 275 PetscErrorCode ierr; 276 277 PetscFunctionBegin; 278 if (TSRKPackageInitialized) PetscFunctionReturn(0); 279 TSRKPackageInitialized = PETSC_TRUE; 280 ierr = TSRKRegisterAll();CHKERRQ(ierr); 281 ierr = PetscRegisterFinalize(TSRKFinalizePackage);CHKERRQ(ierr); 282 PetscFunctionReturn(0); 283 } 284 285 /*@C 286 TSRKFinalizePackage - This function destroys everything in the TSRK package. It is 287 called from PetscFinalize(). 288 289 Level: developer 290 291 .keywords: Petsc, destroy, package 292 .seealso: PetscFinalize() 293 @*/ 294 PetscErrorCode TSRKFinalizePackage(void) 295 { 296 PetscErrorCode ierr; 297 298 PetscFunctionBegin; 299 TSRKPackageInitialized = PETSC_FALSE; 300 ierr = TSRKRegisterDestroy();CHKERRQ(ierr); 301 PetscFunctionReturn(0); 302 } 303 304 /*@C 305 TSRKRegister - register an RK scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 306 307 Not Collective, but the same schemes should be registered on all processes on which they will be used 308 309 Input Parameters: 310 + name - identifier for method 311 . order - approximation order of method 312 . s - number of stages, this is the dimension of the matrices below 313 . A - stage coefficients (dimension s*s, row-major) 314 . b - step completion table (dimension s; NULL to use last row of A) 315 . c - abscissa (dimension s; NULL to use row sums of A) 316 . bembed - completion table for embedded method (dimension s; NULL if not available) 317 . p - Order of the interpolation scheme, equal to the number of columns of binterp 318 - binterp - Coefficients of the interpolation formula (dimension s*p; NULL to reuse b with p=1) 319 320 Notes: 321 Several RK methods are provided, this function is only needed to create new methods. 322 323 Level: advanced 324 325 .keywords: TS, register 326 327 .seealso: TSRK 328 @*/ 329 PetscErrorCode TSRKRegister(TSRKType name,PetscInt order,PetscInt s, 330 const PetscReal A[],const PetscReal b[],const PetscReal c[], 331 const PetscReal bembed[],PetscInt p,const PetscReal binterp[]) 332 { 333 PetscErrorCode ierr; 334 RKTableauLink link; 335 RKTableau t; 336 PetscInt i,j; 337 338 PetscFunctionBegin; 339 PetscValidCharPointer(name,1); 340 PetscValidRealPointer(A,4); 341 if (b) PetscValidRealPointer(b,5); 342 if (c) PetscValidRealPointer(c,6); 343 if (bembed) PetscValidRealPointer(bembed,7); 344 if (binterp || p > 1) PetscValidRealPointer(binterp,9); 345 346 ierr = PetscNew(&link);CHKERRQ(ierr); 347 t = &link->tab; 348 349 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 350 t->order = order; 351 t->s = s; 352 ierr = PetscMalloc3(s*s,&t->A,s,&t->b,s,&t->c);CHKERRQ(ierr); 353 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 354 if (b) { ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); } 355 else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i]; 356 if (c) { ierr = PetscMemcpy(t->c,c,s*sizeof(c[0]));CHKERRQ(ierr); } 357 else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 358 t->FSAL = PETSC_TRUE; 359 for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE; 360 361 if (bembed) { 362 ierr = PetscMalloc1(s,&t->bembed);CHKERRQ(ierr); 363 ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr); 364 } 365 366 if (!binterp) { p = 1; binterp = t->b; } 367 t->p = p; 368 ierr = PetscMalloc1(s*p,&t->binterp);CHKERRQ(ierr); 369 ierr = PetscMemcpy(t->binterp,binterp,s*p*sizeof(binterp[0]));CHKERRQ(ierr); 370 371 link->next = RKTableauList; 372 RKTableauList = link; 373 PetscFunctionReturn(0); 374 } 375 376 /* 377 The step completion formula is 378 379 x1 = x0 + h b^T YdotRHS 380 381 This function can be called before or after ts->vec_sol has been updated. 382 Suppose we have a completion formula (b) and an embedded formula (be) of different order. 383 We can write 384 385 x1e = x0 + h be^T YdotRHS 386 = x1 - h b^T YdotRHS + h be^T YdotRHS 387 = x1 + h (be - b)^T YdotRHS 388 389 so we can evaluate the method with different order even after the step has been optimistically completed. 390 */ 391 static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done) 392 { 393 TS_RK *rk = (TS_RK*)ts->data; 394 RKTableau tab = rk->tableau; 395 PetscScalar *w = rk->work; 396 PetscReal h; 397 PetscInt s = tab->s,j; 398 PetscErrorCode ierr; 399 400 PetscFunctionBegin; 401 switch (rk->status) { 402 case TS_STEP_INCOMPLETE: 403 case TS_STEP_PENDING: 404 h = ts->time_step; break; 405 case TS_STEP_COMPLETE: 406 h = ts->ptime - ts->ptime_prev; break; 407 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 408 } 409 if (order == tab->order) { 410 if (rk->status == TS_STEP_INCOMPLETE) { 411 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 412 for (j=0; j<s; j++) w[j] = h*tab->b[j]; 413 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 414 } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 415 PetscFunctionReturn(0); 416 } else if (order == tab->order-1) { 417 if (!tab->bembed) goto unavailable; 418 if (rk->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (be) */ 419 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 420 for (j=0; j<s; j++) w[j] = h*tab->bembed[j]; 421 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 422 } else { /* Rollback and re-complete using (be-b) */ 423 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 424 for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]); 425 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 426 if (ts->vec_costintegral && ts->costintegralfwd) { 427 ierr = VecCopy(rk->VecCostIntegral0,ts->vec_costintegral);CHKERRQ(ierr); 428 } 429 } 430 if (done) *done = PETSC_TRUE; 431 PetscFunctionReturn(0); 432 } 433 unavailable: 434 if (done) *done = PETSC_FALSE; 435 else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"RK '%s' of order %D cannot evaluate step at order %D. Consider using -ts_adapt_type none or a different method that has an embedded estimate.",tab->name,tab->order,order); 436 PetscFunctionReturn(0); 437 } 438 439 static PetscErrorCode TSForwardCostIntegral_RK(TS ts) 440 { 441 TS_RK *rk = (TS_RK*)ts->data; 442 RKTableau tab = rk->tableau; 443 const PetscInt s = tab->s; 444 const PetscReal *b = tab->b,*c = tab->c; 445 Vec *Y = rk->Y; 446 PetscInt i; 447 PetscErrorCode ierr; 448 449 PetscFunctionBegin; 450 /* backup cost integral */ 451 ierr = VecCopy(ts->vec_costintegral,rk->VecCostIntegral0);CHKERRQ(ierr); 452 for (i=s-1; i>=0; i--) { 453 /* Evolve ts->vec_costintegral to compute integrals */ 454 ierr = TSAdjointComputeCostIntegrand(ts,rk->ptime+rk->time_step*(1.0-c[i]),Y[i],ts->vec_costintegrand);CHKERRQ(ierr); 455 ierr = VecAXPY(ts->vec_costintegral,rk->time_step*b[i],ts->vec_costintegrand);CHKERRQ(ierr); 456 } 457 PetscFunctionReturn(0); 458 } 459 460 static PetscErrorCode TSAdjointCostIntegral_RK(TS ts) 461 { 462 TS_RK *rk = (TS_RK*)ts->data; 463 RKTableau tab = rk->tableau; 464 const PetscInt s = tab->s; 465 const PetscReal *b = tab->b,*c = tab->c; 466 Vec *Y = rk->Y; 467 PetscInt i; 468 PetscErrorCode ierr; 469 470 PetscFunctionBegin; 471 for (i=s-1; i>=0; i--) { 472 /* Evolve ts->vec_costintegral to compute integrals */ 473 ierr = TSAdjointComputeCostIntegrand(ts,ts->ptime-ts->time_step*(1.0-c[i]),Y[i],ts->vec_costintegrand);CHKERRQ(ierr); 474 ierr = VecAXPY(ts->vec_costintegral,-ts->time_step*b[i],ts->vec_costintegrand);CHKERRQ(ierr); 475 } 476 PetscFunctionReturn(0); 477 } 478 479 static PetscErrorCode TSRollBack_RK(TS ts) 480 { 481 TS_RK *rk = (TS_RK*)ts->data; 482 RKTableau tab = rk->tableau; 483 const PetscInt s = tab->s; 484 const PetscReal *b = tab->b; 485 PetscScalar *w = rk->work; 486 Vec *YdotRHS = rk->YdotRHS; 487 PetscInt j; 488 PetscReal h; 489 PetscErrorCode ierr; 490 491 PetscFunctionBegin; 492 switch (rk->status) { 493 case TS_STEP_INCOMPLETE: 494 case TS_STEP_PENDING: 495 h = ts->time_step; break; 496 case TS_STEP_COMPLETE: 497 h = ts->ptime - ts->ptime_prev; break; 498 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 499 } 500 for (j=0; j<s; j++) w[j] = -h*b[j]; 501 ierr = VecMAXPY(ts->vec_sol,s,w,YdotRHS);CHKERRQ(ierr); 502 PetscFunctionReturn(0); 503 } 504 505 static PetscErrorCode TSStep_RK(TS ts) 506 { 507 TS_RK *rk = (TS_RK*)ts->data; 508 RKTableau tab = rk->tableau; 509 const PetscInt s = tab->s; 510 const PetscReal *A = tab->A,*c = tab->c; 511 PetscScalar *w = rk->work; 512 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 513 PetscBool FSAL = tab->FSAL; 514 TSAdapt adapt; 515 PetscInt i,j; 516 PetscInt rejections = 0; 517 PetscBool stageok,accept = PETSC_TRUE; 518 PetscReal next_time_step = ts->time_step; 519 PetscErrorCode ierr; 520 521 PetscFunctionBegin; 522 if (ts->steprollback || ts->steprestart) FSAL = PETSC_FALSE; 523 if (FSAL) { ierr = VecCopy(YdotRHS[s-1],YdotRHS[0]);CHKERRQ(ierr); } 524 525 rk->status = TS_STEP_INCOMPLETE; 526 while (!ts->reason && rk->status != TS_STEP_COMPLETE) { 527 PetscReal t = ts->ptime; 528 PetscReal h = ts->time_step; 529 for (i=0; i<s; i++) { 530 rk->stage_time = t + h*c[i]; 531 ierr = TSPreStage(ts,rk->stage_time); CHKERRQ(ierr); 532 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 533 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 534 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 535 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 536 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 537 ierr = TSAdaptCheckStage(adapt,ts,rk->stage_time,Y[i],&stageok);CHKERRQ(ierr); 538 if (!stageok) goto reject_step; 539 if (FSAL && !i) continue; 540 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 541 } 542 543 rk->status = TS_STEP_INCOMPLETE; 544 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 545 rk->status = TS_STEP_PENDING; 546 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 547 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 548 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);CHKERRQ(ierr); 549 ierr = TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr); 550 rk->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 551 if (!accept) { /* Roll back the current step */ 552 ierr = TSRollBack_RK(ts);CHKERRQ(ierr); 553 ts->time_step = next_time_step; 554 goto reject_step; 555 } 556 557 if (ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation*/ 558 rk->ptime = ts->ptime; 559 rk->time_step = ts->time_step; 560 } 561 562 ts->ptime += ts->time_step; 563 ts->time_step = next_time_step; 564 break; 565 566 reject_step: 567 ts->reject++; accept = PETSC_FALSE; 568 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 569 ts->reason = TS_DIVERGED_STEP_REJECTED; 570 ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr); 571 } 572 } 573 PetscFunctionReturn(0); 574 } 575 576 static PetscErrorCode TSAdjointSetUp_RK(TS ts) 577 { 578 TS_RK *rk = (TS_RK*)ts->data; 579 RKTableau tab = rk->tableau; 580 PetscInt s = tab->s; 581 PetscErrorCode ierr; 582 583 PetscFunctionBegin; 584 if (ts->adjointsetupcalled++) PetscFunctionReturn(0); 585 ierr = VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecDeltaLam);CHKERRQ(ierr); 586 if(ts->vecs_sensip) { 587 ierr = VecDuplicateVecs(ts->vecs_sensip[0],s*ts->numcost,&rk->VecDeltaMu);CHKERRQ(ierr); 588 } 589 ierr = VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&rk->VecSensiTemp);CHKERRQ(ierr); 590 PetscFunctionReturn(0); 591 } 592 593 static PetscErrorCode TSAdjointStep_RK(TS ts) 594 { 595 TS_RK *rk = (TS_RK*)ts->data; 596 RKTableau tab = rk->tableau; 597 const PetscInt s = tab->s; 598 const PetscReal *A = tab->A,*b = tab->b,*c = tab->c; 599 PetscScalar *w = rk->work; 600 Vec *Y = rk->Y,*VecDeltaLam = rk->VecDeltaLam,*VecDeltaMu = rk->VecDeltaMu,*VecSensiTemp = rk->VecSensiTemp; 601 PetscInt i,j,nadj; 602 PetscReal t = ts->ptime; 603 PetscErrorCode ierr; 604 PetscReal h = ts->time_step; 605 Mat J,Jp; 606 607 PetscFunctionBegin; 608 rk->status = TS_STEP_INCOMPLETE; 609 for (i=s-1; i>=0; i--) { 610 rk->stage_time = t + h*(1.0-c[i]); 611 for (nadj=0; nadj<ts->numcost; nadj++) { 612 ierr = VecCopy(ts->vecs_sensi[nadj],VecSensiTemp[nadj]);CHKERRQ(ierr); 613 ierr = VecScale(VecSensiTemp[nadj],-h*b[i]);CHKERRQ(ierr); 614 for (j=i+1; j<s; j++) { 615 ierr = VecAXPY(VecSensiTemp[nadj],-h*A[j*s+i],VecDeltaLam[nadj*s+j]);CHKERRQ(ierr); 616 } 617 } 618 /* Stage values of lambda */ 619 ierr = TSGetRHSJacobian(ts,&J,&Jp,NULL,NULL);CHKERRQ(ierr); 620 ierr = TSComputeRHSJacobian(ts,rk->stage_time,Y[i],J,Jp);CHKERRQ(ierr); 621 if (ts->vec_costintegral) { 622 ierr = TSAdjointComputeDRDYFunction(ts,rk->stage_time,Y[i],ts->vecs_drdy);CHKERRQ(ierr); 623 } 624 for (nadj=0; nadj<ts->numcost; nadj++) { 625 ierr = MatMultTranspose(J,VecSensiTemp[nadj],VecDeltaLam[nadj*s+i]);CHKERRQ(ierr); 626 if (ts->vec_costintegral) { 627 ierr = VecAXPY(VecDeltaLam[nadj*s+i],-h*b[i],ts->vecs_drdy[nadj]);CHKERRQ(ierr); 628 } 629 } 630 631 /* Stage values of mu */ 632 if(ts->vecs_sensip) { 633 ierr = TSAdjointComputeRHSJacobian(ts,rk->stage_time,Y[i],ts->Jacp);CHKERRQ(ierr); 634 if (ts->vec_costintegral) { 635 ierr = TSAdjointComputeDRDPFunction(ts,rk->stage_time,Y[i],ts->vecs_drdp);CHKERRQ(ierr); 636 } 637 638 for (nadj=0; nadj<ts->numcost; nadj++) { 639 ierr = MatMultTranspose(ts->Jacp,VecSensiTemp[nadj],VecDeltaMu[nadj*s+i]);CHKERRQ(ierr); 640 if (ts->vec_costintegral) { 641 ierr = VecAXPY(VecDeltaMu[nadj*s+i],-h*b[i],ts->vecs_drdp[nadj]);CHKERRQ(ierr); 642 } 643 } 644 } 645 } 646 647 for (j=0; j<s; j++) w[j] = 1.0; 648 for (nadj=0; nadj<ts->numcost; nadj++) { 649 ierr = VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecDeltaLam[nadj*s]);CHKERRQ(ierr); 650 if(ts->vecs_sensip) { 651 ierr = VecMAXPY(ts->vecs_sensip[nadj],s,w,&VecDeltaMu[nadj*s]);CHKERRQ(ierr); 652 } 653 } 654 rk->status = TS_STEP_COMPLETE; 655 PetscFunctionReturn(0); 656 } 657 658 static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X) 659 { 660 TS_RK *rk = (TS_RK*)ts->data; 661 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 662 PetscReal h; 663 PetscReal tt,t; 664 PetscScalar *b; 665 const PetscReal *B = rk->tableau->binterp; 666 PetscErrorCode ierr; 667 668 PetscFunctionBegin; 669 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 670 671 switch (rk->status) { 672 case TS_STEP_INCOMPLETE: 673 case TS_STEP_PENDING: 674 h = ts->time_step; 675 t = (itime - ts->ptime)/h; 676 break; 677 case TS_STEP_COMPLETE: 678 h = ts->ptime - ts->ptime_prev; 679 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 680 break; 681 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 682 } 683 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 684 for (i=0; i<s; i++) b[i] = 0; 685 for (j=0,tt=t; j<p; j++,tt*=t) { 686 for (i=0; i<s; i++) { 687 b[i] += h * B[i*p+j] * tt; 688 } 689 } 690 691 ierr = VecCopy(rk->Y[0],X);CHKERRQ(ierr); 692 ierr = VecMAXPY(X,s,b,rk->YdotRHS);CHKERRQ(ierr); 693 694 ierr = PetscFree(b);CHKERRQ(ierr); 695 PetscFunctionReturn(0); 696 } 697 698 /*------------------------------------------------------------*/ 699 700 static PetscErrorCode TSRKTableauReset(TS ts) 701 { 702 TS_RK *rk = (TS_RK*)ts->data; 703 RKTableau tab = rk->tableau; 704 PetscErrorCode ierr; 705 706 PetscFunctionBegin; 707 if (!tab) PetscFunctionReturn(0); 708 ierr = PetscFree(rk->work);CHKERRQ(ierr); 709 ierr = VecDestroyVecs(tab->s,&rk->Y);CHKERRQ(ierr); 710 ierr = VecDestroyVecs(tab->s,&rk->YdotRHS);CHKERRQ(ierr); 711 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecDeltaLam);CHKERRQ(ierr); 712 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecDeltaMu);CHKERRQ(ierr); 713 PetscFunctionReturn(0); 714 } 715 716 static PetscErrorCode TSReset_RK(TS ts) 717 { 718 TS_RK *rk = (TS_RK*)ts->data; 719 PetscErrorCode ierr; 720 721 PetscFunctionBegin; 722 ierr = TSRKTableauReset(ts);CHKERRQ(ierr); 723 ierr = VecDestroy(&rk->VecCostIntegral0);CHKERRQ(ierr); 724 ierr = VecDestroyVecs(ts->numcost,&rk->VecSensiTemp);CHKERRQ(ierr); 725 PetscFunctionReturn(0); 726 } 727 728 static PetscErrorCode TSDestroy_RK(TS ts) 729 { 730 PetscErrorCode ierr; 731 732 PetscFunctionBegin; 733 ierr = TSReset_RK(ts);CHKERRQ(ierr); 734 ierr = PetscFree(ts->data);CHKERRQ(ierr); 735 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);CHKERRQ(ierr); 736 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);CHKERRQ(ierr); 737 PetscFunctionReturn(0); 738 } 739 740 741 static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx) 742 { 743 PetscFunctionBegin; 744 PetscFunctionReturn(0); 745 } 746 747 static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 748 { 749 PetscFunctionBegin; 750 PetscFunctionReturn(0); 751 } 752 753 754 static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx) 755 { 756 PetscFunctionBegin; 757 PetscFunctionReturn(0); 758 } 759 760 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 761 { 762 763 PetscFunctionBegin; 764 PetscFunctionReturn(0); 765 } 766 /* 767 static PetscErrorCode RKSetAdjCoe(RKTableau tab) 768 { 769 PetscReal *A,*b; 770 PetscInt s,i,j; 771 PetscErrorCode ierr; 772 773 PetscFunctionBegin; 774 s = tab->s; 775 ierr = PetscMalloc2(s*s,&A,s,&b);CHKERRQ(ierr); 776 777 for (i=0; i<s; i++) 778 for (j=0; j<s; j++) { 779 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]; 780 ierr = PetscPrintf(PETSC_COMM_WORLD,"Coefficients: A[%D][%D]=%.6f\n",i,j,A[i*s+j]);CHKERRQ(ierr); 781 } 782 783 for (i=0; i<s; i++) b[i] = (tab->b[s-1-i]==0)? 0: -tab->b[s-1-i]; 784 785 ierr = PetscMemcpy(tab->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 786 ierr = PetscMemcpy(tab->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 787 ierr = PetscFree2(A,b);CHKERRQ(ierr); 788 PetscFunctionReturn(0); 789 } 790 */ 791 792 static PetscErrorCode TSRKTableauSetUp(TS ts) 793 { 794 TS_RK *rk = (TS_RK*)ts->data; 795 RKTableau tab = rk->tableau; 796 PetscErrorCode ierr; 797 798 PetscFunctionBegin; 799 ierr = PetscMalloc1(tab->s,&rk->work);CHKERRQ(ierr); 800 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->Y);CHKERRQ(ierr); 801 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS);CHKERRQ(ierr); 802 PetscFunctionReturn(0); 803 } 804 805 806 static PetscErrorCode TSSetUp_RK(TS ts) 807 { 808 TS_RK *rk = (TS_RK*)ts->data; 809 PetscErrorCode ierr; 810 DM dm; 811 812 PetscFunctionBegin; 813 ierr = TSCheckImplicitTerm(ts);CHKERRQ(ierr); 814 ierr = TSRKTableauSetUp(ts);CHKERRQ(ierr); 815 if (!rk->VecCostIntegral0 && ts->vec_costintegral && ts->costintegralfwd) { /* back up cost integral */ 816 ierr = VecDuplicate(ts->vec_costintegral,&rk->VecCostIntegral0);CHKERRQ(ierr); 817 } 818 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 819 if (dm) { 820 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 821 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 822 } 823 PetscFunctionReturn(0); 824 } 825 826 827 /*------------------------------------------------------------*/ 828 829 static PetscErrorCode TSSetFromOptions_RK(PetscOptionItems *PetscOptionsObject,TS ts) 830 { 831 TS_RK *rk = (TS_RK*)ts->data; 832 PetscErrorCode ierr; 833 834 PetscFunctionBegin; 835 ierr = PetscOptionsHead(PetscOptionsObject,"RK ODE solver options");CHKERRQ(ierr); 836 { 837 RKTableauLink link; 838 PetscInt count,choice; 839 PetscBool flg; 840 const char **namelist; 841 for (link=RKTableauList,count=0; link; link=link->next,count++) ; 842 ierr = PetscMalloc1(count,&namelist);CHKERRQ(ierr); 843 for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 844 ierr = PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rk->tableau->name,&choice,&flg);CHKERRQ(ierr); 845 if (flg) {ierr = TSRKSetType(ts,namelist[choice]);CHKERRQ(ierr);} 846 ierr = PetscFree(namelist);CHKERRQ(ierr); 847 } 848 ierr = PetscOptionsTail();CHKERRQ(ierr); 849 PetscFunctionReturn(0); 850 } 851 852 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 853 { 854 PetscErrorCode ierr; 855 PetscInt i; 856 size_t left,count; 857 char *p; 858 859 PetscFunctionBegin; 860 for (i=0,p=buf,left=len; i<n; i++) { 861 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 862 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 863 left -= count; 864 p += count; 865 *p++ = ' '; 866 } 867 p[i ? 0 : -1] = 0; 868 PetscFunctionReturn(0); 869 } 870 871 static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer) 872 { 873 TS_RK *rk = (TS_RK*)ts->data; 874 PetscBool iascii; 875 PetscErrorCode ierr; 876 877 PetscFunctionBegin; 878 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 879 if (iascii) { 880 RKTableau tab = rk->tableau; 881 TSRKType rktype; 882 char buf[512]; 883 ierr = TSRKGetType(ts,&rktype);CHKERRQ(ierr); 884 ierr = PetscViewerASCIIPrintf(viewer," RK type %s\n",rktype);CHKERRQ(ierr); 885 ierr = PetscViewerASCIIPrintf(viewer," Order: %D\n",tab->order);CHKERRQ(ierr); 886 ierr = PetscViewerASCIIPrintf(viewer," FSAL property: %s\n",tab->FSAL ? "yes" : "no");CHKERRQ(ierr); 887 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 888 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 889 } 890 PetscFunctionReturn(0); 891 } 892 893 static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer) 894 { 895 PetscErrorCode ierr; 896 TSAdapt adapt; 897 898 PetscFunctionBegin; 899 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 900 ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 901 PetscFunctionReturn(0); 902 } 903 904 /*@C 905 TSRKSetType - Set the type of RK scheme 906 907 Logically collective 908 909 Input Parameter: 910 + ts - timestepping context 911 - rktype - type of RK-scheme 912 913 Level: intermediate 914 915 .seealso: TSRKGetType(), TSRK, TSRK2, TSRK3, TSRKPRSSP2, TSRK4, TSRK5 916 @*/ 917 PetscErrorCode TSRKSetType(TS ts,TSRKType rktype) 918 { 919 PetscErrorCode ierr; 920 921 PetscFunctionBegin; 922 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 923 PetscValidCharPointer(rktype,2); 924 ierr = PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));CHKERRQ(ierr); 925 PetscFunctionReturn(0); 926 } 927 928 /*@C 929 TSRKGetType - Get the type of RK scheme 930 931 Logically collective 932 933 Input Parameter: 934 . ts - timestepping context 935 936 Output Parameter: 937 . rktype - type of RK-scheme 938 939 Level: intermediate 940 941 .seealso: TSRKGetType() 942 @*/ 943 PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype) 944 { 945 PetscErrorCode ierr; 946 947 PetscFunctionBegin; 948 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 949 ierr = PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));CHKERRQ(ierr); 950 PetscFunctionReturn(0); 951 } 952 953 static PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype) 954 { 955 TS_RK *rk = (TS_RK*)ts->data; 956 957 PetscFunctionBegin; 958 *rktype = rk->tableau->name; 959 PetscFunctionReturn(0); 960 } 961 static PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype) 962 { 963 TS_RK *rk = (TS_RK*)ts->data; 964 PetscErrorCode ierr; 965 PetscBool match; 966 RKTableauLink link; 967 968 PetscFunctionBegin; 969 if (rk->tableau) { 970 ierr = PetscStrcmp(rk->tableau->name,rktype,&match);CHKERRQ(ierr); 971 if (match) PetscFunctionReturn(0); 972 } 973 for (link = RKTableauList; link; link=link->next) { 974 ierr = PetscStrcmp(link->tab.name,rktype,&match);CHKERRQ(ierr); 975 if (match) { 976 if (ts->setupcalled) {ierr = TSRKTableauReset(ts);CHKERRQ(ierr);} 977 rk->tableau = &link->tab; 978 if (ts->setupcalled) {ierr = TSRKTableauSetUp(ts);CHKERRQ(ierr);} 979 ts->default_adapt_type = rk->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 980 PetscFunctionReturn(0); 981 } 982 } 983 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype); 984 PetscFunctionReturn(0); 985 } 986 987 static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y) 988 { 989 TS_RK *rk = (TS_RK*)ts->data; 990 991 PetscFunctionBegin; 992 *ns = rk->tableau->s; 993 if(Y) *Y = rk->Y; 994 PetscFunctionReturn(0); 995 } 996 997 998 /* ------------------------------------------------------------ */ 999 /*MC 1000 TSRK - ODE and DAE solver using Runge-Kutta schemes 1001 1002 The user should provide the right hand side of the equation 1003 using TSSetRHSFunction(). 1004 1005 Notes: 1006 The default is TSRK3BS, it can be changed with TSRKSetType() or -ts_rk_type 1007 1008 Level: beginner 1009 1010 .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3, 1011 TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister() 1012 1013 M*/ 1014 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 1015 { 1016 TS_RK *rk; 1017 PetscErrorCode ierr; 1018 1019 PetscFunctionBegin; 1020 ierr = TSRKInitializePackage();CHKERRQ(ierr); 1021 1022 ts->ops->reset = TSReset_RK; 1023 ts->ops->destroy = TSDestroy_RK; 1024 ts->ops->view = TSView_RK; 1025 ts->ops->load = TSLoad_RK; 1026 ts->ops->setup = TSSetUp_RK; 1027 ts->ops->adjointsetup = TSAdjointSetUp_RK; 1028 ts->ops->step = TSStep_RK; 1029 ts->ops->interpolate = TSInterpolate_RK; 1030 ts->ops->evaluatestep = TSEvaluateStep_RK; 1031 ts->ops->rollback = TSRollBack_RK; 1032 ts->ops->setfromoptions = TSSetFromOptions_RK; 1033 ts->ops->getstages = TSGetStages_RK; 1034 ts->ops->adjointstep = TSAdjointStep_RK; 1035 1036 ts->ops->adjointintegral = TSAdjointCostIntegral_RK; 1037 ts->ops->forwardintegral = TSForwardCostIntegral_RK; 1038 1039 ierr = PetscNewLog(ts,&rk);CHKERRQ(ierr); 1040 ts->data = (void*)rk; 1041 1042 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);CHKERRQ(ierr); 1043 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);CHKERRQ(ierr); 1044 1045 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 1046 PetscFunctionReturn(0); 1047 } 1048