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 11 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 12 #include <petscdm.h> 13 #include <../src/ts/impls/explicit/rk/rk.h> 14 #include <../src/ts/impls/explicit/rk/mrk.h> 15 16 static TSRKType TSRKDefault = TSRK3BS; 17 static PetscBool TSRKRegisterAllCalled; 18 static PetscBool TSRKPackageInitialized; 19 20 static RKTableauLink RKTableauList; 21 22 /*MC 23 TSRK1FE - First order forward Euler scheme. 24 25 This method has one stage. 26 27 Options database: 28 . -ts_rk_type 1fe 29 30 Level: advanced 31 32 .seealso: TSRK, TSRKType, TSRKSetType() 33 M*/ 34 /*MC 35 TSRK2A - Second order RK scheme. 36 37 This method has two stages. 38 39 Options database: 40 . -ts_rk_type 2a 41 42 Level: advanced 43 44 .seealso: TSRK, TSRKType, TSRKSetType() 45 M*/ 46 /*MC 47 TSRK3 - Third order RK scheme. 48 49 This method has three stages. 50 51 Options database: 52 . -ts_rk_type 3 53 54 Level: advanced 55 56 .seealso: TSRK, TSRKType, TSRKSetType() 57 M*/ 58 /*MC 59 TSRK3BS - Third order RK scheme of Bogacki-Shampine with 2nd order embedded method. 60 61 This method has four stages with the First Same As Last (FSAL) property. 62 63 Options database: 64 . -ts_rk_type 3bs 65 66 Level: advanced 67 68 References: https://doi.org/10.1016/0893-9659(89)90079-7 69 70 .seealso: TSRK, TSRKType, TSRKSetType() 71 M*/ 72 /*MC 73 TSRK4 - Fourth order RK scheme. 74 75 This is the classical Runge-Kutta method with four stages. 76 77 Options database: 78 . -ts_rk_type 4 79 80 Level: advanced 81 82 .seealso: TSRK, TSRKType, TSRKSetType() 83 M*/ 84 /*MC 85 TSRK5F - Fifth order Fehlberg RK scheme with a 4th order embedded method. 86 87 This method has six stages. 88 89 Options database: 90 . -ts_rk_type 5f 91 92 Level: advanced 93 94 .seealso: TSRK, TSRKType, TSRKSetType() 95 M*/ 96 /*MC 97 TSRK5DP - Fifth order Dormand-Prince RK scheme with the 4th order embedded method. 98 99 This method has seven stages with the First Same As Last (FSAL) property. 100 101 Options database: 102 . -ts_rk_type 5dp 103 104 Level: advanced 105 106 References: https://doi.org/10.1016/0771-050X(80)90013-3 107 108 .seealso: TSRK, TSRKType, TSRKSetType() 109 M*/ 110 /*MC 111 TSRK5BS - Fifth order Bogacki-Shampine RK scheme with 4th order embedded method. 112 113 This method has eight stages with the First Same As Last (FSAL) property. 114 115 Options database: 116 . -ts_rk_type 5bs 117 118 Level: advanced 119 120 References: https://doi.org/10.1016/0898-1221(96)00141-1 121 122 .seealso: TSRK, TSRKType, TSRKSetType() 123 M*/ 124 125 /*@C 126 TSRKRegisterAll - Registers all of the Runge-Kutta explicit methods in TSRK 127 128 Not Collective, but should be called by all processes which will need the schemes to be registered 129 130 Level: advanced 131 132 .keywords: TS, TSRK, register, all 133 134 .seealso: TSRKRegisterDestroy() 135 @*/ 136 PetscErrorCode TSRKRegisterAll(void) 137 { 138 PetscErrorCode ierr; 139 140 PetscFunctionBegin; 141 if (TSRKRegisterAllCalled) PetscFunctionReturn(0); 142 TSRKRegisterAllCalled = PETSC_TRUE; 143 144 #define RC PetscRealConstant 145 { 146 const PetscReal 147 A[1][1] = {{0}}, 148 b[1] = {RC(1.0)}; 149 ierr = TSRKRegister(TSRK1FE,1,1,&A[0][0],b,NULL,NULL,0,NULL);CHKERRQ(ierr); 150 } 151 { 152 const PetscReal 153 A[2][2] = {{0,0}, 154 {RC(1.0),0}}, 155 b[2] = {RC(0.5),RC(0.5)}, 156 bembed[2] = {RC(1.0),0}; 157 ierr = TSRKRegister(TSRK2A,2,2,&A[0][0],b,NULL,bembed,0,NULL);CHKERRQ(ierr); 158 } 159 { 160 const PetscReal 161 A[3][3] = {{0,0,0}, 162 {RC(2.0)/RC(3.0),0,0}, 163 {RC(-1.0)/RC(3.0),RC(1.0),0}}, 164 b[3] = {RC(0.25),RC(0.5),RC(0.25)}; 165 ierr = TSRKRegister(TSRK3,3,3,&A[0][0],b,NULL,NULL,0,NULL);CHKERRQ(ierr); 166 } 167 { 168 const PetscReal 169 A[4][4] = {{0,0,0,0}, 170 {RC(1.0)/RC(2.0),0,0,0}, 171 {0,RC(3.0)/RC(4.0),0,0}, 172 {RC(2.0)/RC(9.0),RC(1.0)/RC(3.0),RC(4.0)/RC(9.0),0}}, 173 b[4] = {RC(2.0)/RC(9.0),RC(1.0)/RC(3.0),RC(4.0)/RC(9.0),0}, 174 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)}; 175 ierr = TSRKRegister(TSRK3BS,3,4,&A[0][0],b,NULL,bembed,0,NULL);CHKERRQ(ierr); 176 } 177 { 178 const PetscReal 179 A[4][4] = {{0,0,0,0}, 180 {RC(0.5),0,0,0}, 181 {0,RC(0.5),0,0}, 182 {0,0,RC(1.0),0}}, 183 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)}; 184 ierr = TSRKRegister(TSRK4,4,4,&A[0][0],b,NULL,NULL,0,NULL);CHKERRQ(ierr); 185 } 186 { 187 const PetscReal 188 A[6][6] = {{0,0,0,0,0,0}, 189 {RC(0.25),0,0,0,0,0}, 190 {RC(3.0)/RC(32.0),RC(9.0)/RC(32.0),0,0,0,0}, 191 {RC(1932.0)/RC(2197.0),RC(-7200.0)/RC(2197.0),RC(7296.0)/RC(2197.0),0,0,0}, 192 {RC(439.0)/RC(216.0),RC(-8.0),RC(3680.0)/RC(513.0),RC(-845.0)/RC(4104.0),0,0}, 193 {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}}, 194 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)}, 195 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}; 196 ierr = TSRKRegister(TSRK5F,5,6,&A[0][0],b,NULL,bembed,0,NULL);CHKERRQ(ierr); 197 } 198 { 199 const PetscReal 200 A[7][7] = {{0,0,0,0,0,0,0}, 201 {RC(1.0)/RC(5.0),0,0,0,0,0,0}, 202 {RC(3.0)/RC(40.0),RC(9.0)/RC(40.0),0,0,0,0,0}, 203 {RC(44.0)/RC(45.0),RC(-56.0)/RC(15.0),RC(32.0)/RC(9.0),0,0,0,0}, 204 {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}, 205 {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}, 206 {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}}, 207 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}, 208 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)}, 209 binterp[7][5] = {{RC(1.0),RC(-4034104133.0)/RC(1410260304.0),RC(105330401.0)/RC(33982176.0),RC(-13107642775.0)/RC(11282082432.0),RC(6542295.0)/RC(470086768.0)}, 210 {0,0,0,0,0}, 211 {0,RC(132343189600.0)/RC(32700410799.0),RC(-833316000.0)/RC(131326951.0),RC(91412856700.0)/RC(32700410799.0),RC(-523383600.0)/RC(10900136933.0)}, 212 {0,RC(-115792950.0)/RC(29380423.0),RC(185270875.0)/RC(16991088.0),RC(-12653452475.0)/RC(1880347072.0),RC(98134425.0)/RC(235043384.0)}, 213 {0,RC(70805911779.0)/RC(24914598704.0),RC(-4531260609.0)/RC(600351776.0),RC(988140236175.0)/RC(199316789632.0),RC(-14307999165.0)/RC(24914598704.0)}, 214 {0,RC(-331320693.0)/RC(205662961.0),RC(31361737.0)/RC(7433601.0),RC(-2426908385.0)/RC(822651844.0),RC(97305120.0)/RC(205662961.0)}, 215 {0,RC(44764047.0)/RC(29380423.0),RC(-1532549.0)/RC(353981.0),RC(90730570.0)/RC(29380423.0),RC(-8293050.0)/RC(29380423.0)}}; 216 217 ierr = TSRKRegister(TSRK5DP,5,7,&A[0][0],b,NULL,bembed,5,binterp[0]);CHKERRQ(ierr); 218 } 219 { 220 const PetscReal 221 A[8][8] = {{0,0,0,0,0,0,0,0}, 222 {RC(1.0)/RC(6.0),0,0,0,0,0,0,0}, 223 {RC(2.0)/RC(27.0),RC(4.0)/RC(27.0),0,0,0,0,0,0}, 224 {RC(183.0)/RC(1372.0),RC(-162.0)/RC(343.0),RC(1053.0)/RC(1372.0),0,0,0,0,0}, 225 {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}, 226 {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}, 227 {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}, 228 {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}}, 229 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}, 230 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)}; 231 ierr = TSRKRegister(TSRK5BS,5,8,&A[0][0],b,NULL,bembed,0,NULL);CHKERRQ(ierr); 232 } 233 #undef RC 234 PetscFunctionReturn(0); 235 } 236 237 /*@C 238 TSRKRegisterDestroy - Frees the list of schemes that were registered by TSRKRegister(). 239 240 Not Collective 241 242 Level: advanced 243 244 .keywords: TSRK, register, destroy 245 .seealso: TSRKRegister(), TSRKRegisterAll() 246 @*/ 247 PetscErrorCode TSRKRegisterDestroy(void) 248 { 249 PetscErrorCode ierr; 250 RKTableauLink link; 251 252 PetscFunctionBegin; 253 while ((link = RKTableauList)) { 254 RKTableau t = &link->tab; 255 RKTableauList = link->next; 256 ierr = PetscFree3(t->A,t->b,t->c); CHKERRQ(ierr); 257 ierr = PetscFree (t->bembed); CHKERRQ(ierr); 258 ierr = PetscFree (t->binterp); CHKERRQ(ierr); 259 ierr = PetscFree (t->name); CHKERRQ(ierr); 260 ierr = PetscFree (link); CHKERRQ(ierr); 261 } 262 TSRKRegisterAllCalled = PETSC_FALSE; 263 PetscFunctionReturn(0); 264 } 265 266 /*@C 267 TSRKInitializePackage - This function initializes everything in the TSRK package. It is called 268 from TSInitializePackage(). 269 270 Level: developer 271 272 .keywords: TS, TSRK, initialize, package 273 .seealso: PetscInitialize() 274 @*/ 275 PetscErrorCode TSRKInitializePackage(void) 276 { 277 PetscErrorCode ierr; 278 279 PetscFunctionBegin; 280 if (TSRKPackageInitialized) PetscFunctionReturn(0); 281 TSRKPackageInitialized = PETSC_TRUE; 282 ierr = TSRKRegisterAll();CHKERRQ(ierr); 283 ierr = PetscRegisterFinalize(TSRKFinalizePackage);CHKERRQ(ierr); 284 PetscFunctionReturn(0); 285 } 286 287 /*@C 288 TSRKFinalizePackage - This function destroys everything in the TSRK package. It is 289 called from PetscFinalize(). 290 291 Level: developer 292 293 .keywords: Petsc, destroy, package 294 .seealso: PetscFinalize() 295 @*/ 296 PetscErrorCode TSRKFinalizePackage(void) 297 { 298 PetscErrorCode ierr; 299 300 PetscFunctionBegin; 301 TSRKPackageInitialized = PETSC_FALSE; 302 ierr = TSRKRegisterDestroy();CHKERRQ(ierr); 303 PetscFunctionReturn(0); 304 } 305 306 /*@C 307 TSRKRegister - register an RK scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 308 309 Not Collective, but the same schemes should be registered on all processes on which they will be used 310 311 Input Parameters: 312 + name - identifier for method 313 . order - approximation order of method 314 . s - number of stages, this is the dimension of the matrices below 315 . A - stage coefficients (dimension s*s, row-major) 316 . b - step completion table (dimension s; NULL to use last row of A) 317 . c - abscissa (dimension s; NULL to use row sums of A) 318 . bembed - completion table for embedded method (dimension s; NULL if not available) 319 . p - Order of the interpolation scheme, equal to the number of columns of binterp 320 - binterp - Coefficients of the interpolation formula (dimension s*p; NULL to reuse b with p=1) 321 322 Notes: 323 Several RK methods are provided, this function is only needed to create new methods. 324 325 Level: advanced 326 327 .keywords: TS, register 328 329 .seealso: TSRK 330 @*/ 331 PetscErrorCode TSRKRegister(TSRKType name,PetscInt order,PetscInt s, 332 const PetscReal A[],const PetscReal b[],const PetscReal c[], 333 const PetscReal bembed[],PetscInt p,const PetscReal binterp[]) 334 { 335 PetscErrorCode ierr; 336 RKTableauLink link; 337 RKTableau t; 338 PetscInt i,j; 339 340 PetscFunctionBegin; 341 PetscValidCharPointer(name,1); 342 PetscValidRealPointer(A,4); 343 if (b) PetscValidRealPointer(b,5); 344 if (c) PetscValidRealPointer(c,6); 345 if (bembed) PetscValidRealPointer(bembed,7); 346 if (binterp || p > 1) PetscValidRealPointer(binterp,9); 347 348 ierr = TSRKInitializePackage();CHKERRQ(ierr); 349 ierr = PetscNew(&link);CHKERRQ(ierr); 350 t = &link->tab; 351 352 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 353 t->order = order; 354 t->s = s; 355 ierr = PetscMalloc3(s*s,&t->A,s,&t->b,s,&t->c);CHKERRQ(ierr); 356 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 357 if (b) { ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); } 358 else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i]; 359 if (c) { ierr = PetscMemcpy(t->c,c,s*sizeof(c[0]));CHKERRQ(ierr); } 360 else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 361 t->FSAL = PETSC_TRUE; 362 for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE; 363 364 if (bembed) { 365 ierr = PetscMalloc1(s,&t->bembed);CHKERRQ(ierr); 366 ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr); 367 } 368 369 if (!binterp) { p = 1; binterp = t->b; } 370 t->p = p; 371 ierr = PetscMalloc1(s*p,&t->binterp);CHKERRQ(ierr); 372 ierr = PetscMemcpy(t->binterp,binterp,s*p*sizeof(binterp[0]));CHKERRQ(ierr); 373 374 link->next = RKTableauList; 375 RKTableauList = link; 376 PetscFunctionReturn(0); 377 } 378 379 /* 380 This is for single-step RK method 381 The step completion formula is 382 383 x1 = x0 + h b^T YdotRHS 384 385 This function can be called before or after ts->vec_sol has been updated. 386 Suppose we have a completion formula (b) and an embedded formula (be) of different order. 387 We can write 388 389 x1e = x0 + h be^T YdotRHS 390 = x1 - h b^T YdotRHS + h be^T YdotRHS 391 = x1 + h (be - b)^T YdotRHS 392 393 so we can evaluate the method with different order even after the step has been optimistically completed. 394 */ 395 static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done) 396 { 397 TS_RK *rk = (TS_RK*)ts->data; 398 RKTableau tab = rk->tableau; 399 PetscScalar *w = rk->work; 400 PetscReal h; 401 PetscInt s = tab->s,j; 402 PetscErrorCode ierr; 403 404 PetscFunctionBegin; 405 switch (rk->status) { 406 case TS_STEP_INCOMPLETE: 407 case TS_STEP_PENDING: 408 h = ts->time_step; break; 409 case TS_STEP_COMPLETE: 410 h = ts->ptime - ts->ptime_prev; break; 411 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 412 } 413 if (order == tab->order) { 414 if (rk->status == TS_STEP_INCOMPLETE) { 415 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 416 for (j=0; j<s; j++) w[j] = h*tab->b[j]/rk->dtratio; 417 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 418 } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 419 PetscFunctionReturn(0); 420 } else if (order == tab->order-1) { 421 if (!tab->bembed) goto unavailable; 422 if (rk->status == TS_STEP_INCOMPLETE) { /*Complete with the embedded method (be)*/ 423 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 424 for (j=0; j<s; j++) w[j] = h*tab->bembed[j]; 425 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 426 } else { /*Rollback and re-complete using (be-b) */ 427 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 428 for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]); 429 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 430 } 431 if (done) *done = PETSC_TRUE; 432 PetscFunctionReturn(0); 433 } 434 unavailable: 435 if (done) *done = PETSC_FALSE; 436 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); 437 PetscFunctionReturn(0); 438 } 439 440 static PetscErrorCode TSForwardCostIntegral_RK(TS ts) 441 { 442 TS_RK *rk = (TS_RK*)ts->data; 443 TS quadts = ts->quadraturets; 444 RKTableau tab = rk->tableau; 445 const PetscInt s = tab->s; 446 const PetscReal *b = tab->b,*c = tab->c; 447 Vec *Y = rk->Y; 448 PetscInt i; 449 PetscErrorCode ierr; 450 451 PetscFunctionBegin; 452 /* No need to backup quadts->vec_sol since it can be reverted in TSRollBack_RK */ 453 for (i=s-1; i>=0; i--) { 454 /* Evolve quadrature TS solution to compute integrals */ 455 ierr = TSComputeRHSFunction(quadts,rk->ptime+rk->time_step*c[i],Y[i],ts->vec_costintegrand);CHKERRQ(ierr); 456 ierr = VecAXPY(quadts->vec_sol,rk->time_step*b[i],ts->vec_costintegrand);CHKERRQ(ierr); 457 } 458 PetscFunctionReturn(0); 459 } 460 461 static PetscErrorCode TSAdjointCostIntegral_RK(TS ts) 462 { 463 TS_RK *rk = (TS_RK*)ts->data; 464 RKTableau tab = rk->tableau; 465 TS quadts = ts->quadraturets; 466 const PetscInt s = tab->s; 467 const PetscReal *b = tab->b,*c = tab->c; 468 Vec *Y = rk->Y; 469 PetscInt i; 470 PetscErrorCode ierr; 471 472 PetscFunctionBegin; 473 for (i=s-1; i>=0; i--) { 474 /* Evolve quadrature TS solution to compute integrals */ 475 ierr = TSComputeRHSFunction(quadts,ts->ptime+ts->time_step*(1.0-c[i]),Y[i],ts->vec_costintegrand);CHKERRQ(ierr); 476 ierr = VecAXPY(quadts->vec_sol,-ts->time_step*b[i],ts->vec_costintegrand);CHKERRQ(ierr); 477 } 478 PetscFunctionReturn(0); 479 } 480 481 static PetscErrorCode TSRollBack_RK(TS ts) 482 { 483 TS_RK *rk = (TS_RK*)ts->data; 484 TS quadts = ts->quadraturets; 485 RKTableau tab = rk->tableau; 486 const PetscInt s = tab->s; 487 const PetscReal *b = tab->b,*c = tab->c; 488 PetscScalar *w = rk->work; 489 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 490 PetscInt j; 491 PetscReal h; 492 PetscErrorCode ierr; 493 494 PetscFunctionBegin; 495 switch (rk->status) { 496 case TS_STEP_INCOMPLETE: 497 case TS_STEP_PENDING: 498 h = ts->time_step; break; 499 case TS_STEP_COMPLETE: 500 h = ts->ptime - ts->ptime_prev; break; 501 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 502 } 503 for (j=0; j<s; j++) w[j] = -h*b[j]; 504 ierr = VecMAXPY(ts->vec_sol,s,w,YdotRHS);CHKERRQ(ierr); 505 if (quadts && ts->costintegralfwd) { 506 for (j=0; j<s; j++) { 507 /* Revert the quadrature TS solution */ 508 ierr = TSComputeRHSFunction(quadts,rk->ptime+h*c[j],Y[j],ts->vec_costintegrand);CHKERRQ(ierr); 509 ierr = VecAXPY(quadts->vec_sol,-h*b[j],ts->vec_costintegrand);CHKERRQ(ierr); 510 } 511 } 512 PetscFunctionReturn(0); 513 } 514 515 static PetscErrorCode TSForwardStep_RK(TS ts) 516 { 517 TS_RK *rk = (TS_RK*)ts->data; 518 RKTableau tab = rk->tableau; 519 Mat J,*MatsFwdSensipTemp = rk->MatsFwdSensipTemp; 520 const PetscInt s = tab->s; 521 const PetscReal *A = tab->A,*c = tab->c,*b = tab->b; 522 Vec *Y = rk->Y; 523 PetscInt i,j; 524 PetscReal stage_time,h = ts->time_step; 525 PetscBool zero; 526 PetscErrorCode ierr; 527 528 PetscFunctionBegin; 529 ierr = MatCopy(ts->mat_sensip,rk->MatFwdSensip0,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 530 ierr = TSGetRHSJacobian(ts,&J,NULL,NULL,NULL);CHKERRQ(ierr); 531 532 for (i=0; i<s; i++) { 533 stage_time = ts->ptime + h*c[i]; 534 zero = PETSC_FALSE; 535 if (b[i] == 0 && i == s-1) zero = PETSC_TRUE; 536 /* TLM Stage values */ 537 if(!i) { 538 ierr = MatCopy(ts->mat_sensip,rk->MatsFwdStageSensip[i],SAME_NONZERO_PATTERN);CHKERRQ(ierr); 539 } else if (!zero) { 540 ierr = MatZeroEntries(rk->MatsFwdStageSensip[i]);CHKERRQ(ierr); 541 for (j=0; j<i; j++) { 542 ierr = MatAXPY(rk->MatsFwdStageSensip[i],h*A[i*s+j],MatsFwdSensipTemp[j],SAME_NONZERO_PATTERN);CHKERRQ(ierr); 543 } 544 ierr = MatAXPY(rk->MatsFwdStageSensip[i],1.,ts->mat_sensip,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 545 } else { 546 ierr = MatZeroEntries(rk->MatsFwdStageSensip[i]);CHKERRQ(ierr); 547 } 548 549 ierr = TSComputeRHSJacobian(ts,stage_time,Y[i],J,J);CHKERRQ(ierr); 550 ierr = MatMatMult(J,rk->MatsFwdStageSensip[i],MAT_REUSE_MATRIX,PETSC_DEFAULT,&MatsFwdSensipTemp[i]);CHKERRQ(ierr); 551 if (ts->Jacprhs) { 552 ierr = TSComputeRHSJacobianP(ts,stage_time,Y[i],ts->Jacprhs);CHKERRQ(ierr); /* get f_p */ 553 if (ts->vecs_sensi2p) { /* TLM used for 2nd-order adjoint */ 554 PetscScalar *xarr; 555 ierr = MatDenseGetColumn(MatsFwdSensipTemp[i],0,&xarr);CHKERRQ(ierr); 556 ierr = VecPlaceArray(rk->VecDeltaFwdSensipCol,xarr);CHKERRQ(ierr); 557 ierr = MatMultAdd(ts->Jacprhs,ts->vec_dir,rk->VecDeltaFwdSensipCol,rk->VecDeltaFwdSensipCol);CHKERRQ(ierr); 558 ierr = VecResetArray(rk->VecDeltaFwdSensipCol);CHKERRQ(ierr);CHKERRQ(ierr); 559 ierr = MatDenseRestoreColumn(MatsFwdSensipTemp[i],&xarr);CHKERRQ(ierr); 560 } else { 561 ierr = MatAXPY(MatsFwdSensipTemp[i],1.,ts->Jacprhs,SUBSET_NONZERO_PATTERN);CHKERRQ(ierr); 562 } 563 } 564 } 565 566 for (i=0; i<s; i++) { 567 ierr = MatAXPY(ts->mat_sensip,h*b[i],rk->MatsFwdSensipTemp[i],SAME_NONZERO_PATTERN);CHKERRQ(ierr); 568 } 569 rk->status = TS_STEP_COMPLETE; 570 PetscFunctionReturn(0); 571 } 572 573 static PetscErrorCode TSForwardGetStages_RK(TS ts,PetscInt *ns,Mat **stagesensip) 574 { 575 TS_RK *rk = (TS_RK*)ts->data; 576 RKTableau tab = rk->tableau; 577 578 PetscFunctionBegin; 579 if (ns) *ns = tab->s; 580 if (stagesensip) *stagesensip = rk->MatsFwdStageSensip; 581 PetscFunctionReturn(0); 582 } 583 584 static PetscErrorCode TSForwardSetUp_RK(TS ts) 585 { 586 TS_RK *rk = (TS_RK*)ts->data; 587 RKTableau tab = rk->tableau; 588 PetscInt i; 589 PetscErrorCode ierr; 590 591 PetscFunctionBegin; 592 /* backup sensitivity results for roll-backs */ 593 ierr = MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&rk->MatFwdSensip0);CHKERRQ(ierr); 594 595 ierr = PetscMalloc1(tab->s,&rk->MatsFwdStageSensip);CHKERRQ(ierr); 596 ierr = PetscMalloc1(tab->s,&rk->MatsFwdSensipTemp);CHKERRQ(ierr); 597 for(i=0; i<tab->s; i++) { 598 ierr = MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&rk->MatsFwdStageSensip[i]);CHKERRQ(ierr); 599 ierr = MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&rk->MatsFwdSensipTemp[i]);CHKERRQ(ierr); 600 } 601 ierr = VecDuplicate(ts->vec_sol,&rk->VecDeltaFwdSensipCol);CHKERRQ(ierr); 602 PetscFunctionReturn(0); 603 } 604 605 static PetscErrorCode TSForwardReset_RK(TS ts) 606 { 607 TS_RK *rk = (TS_RK*)ts->data; 608 RKTableau tab = rk->tableau; 609 PetscInt i; 610 PetscErrorCode ierr; 611 612 PetscFunctionBegin; 613 ierr = MatDestroy(&rk->MatFwdSensip0);CHKERRQ(ierr); 614 if (rk->MatsFwdStageSensip) { 615 for (i=0; i<tab->s; i++) { 616 ierr = MatDestroy(&rk->MatsFwdStageSensip[i]);CHKERRQ(ierr); 617 } 618 ierr = PetscFree(rk->MatsFwdStageSensip);CHKERRQ(ierr); 619 } 620 if (rk->MatsFwdSensipTemp) { 621 for (i=0; i<tab->s; i++) { 622 ierr = MatDestroy(&rk->MatsFwdSensipTemp[i]);CHKERRQ(ierr); 623 } 624 ierr = PetscFree(rk->MatsFwdSensipTemp);CHKERRQ(ierr); 625 } 626 ierr = VecDestroy(&rk->VecDeltaFwdSensipCol);CHKERRQ(ierr); 627 PetscFunctionReturn(0); 628 } 629 630 static PetscErrorCode TSStep_RK(TS ts) 631 { 632 TS_RK *rk = (TS_RK*)ts->data; 633 RKTableau tab = rk->tableau; 634 const PetscInt s = tab->s; 635 const PetscReal *A = tab->A,*c = tab->c; 636 PetscScalar *w = rk->work; 637 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 638 PetscBool FSAL = tab->FSAL; 639 TSAdapt adapt; 640 PetscInt i,j; 641 PetscInt rejections = 0; 642 PetscBool stageok,accept = PETSC_TRUE; 643 PetscReal next_time_step = ts->time_step; 644 PetscErrorCode ierr; 645 646 PetscFunctionBegin; 647 if (ts->steprollback || ts->steprestart) FSAL = PETSC_FALSE; 648 if (FSAL) { ierr = VecCopy(YdotRHS[s-1],YdotRHS[0]);CHKERRQ(ierr); } 649 650 rk->status = TS_STEP_INCOMPLETE; 651 while (!ts->reason && rk->status != TS_STEP_COMPLETE) { 652 PetscReal t = ts->ptime; 653 PetscReal h = ts->time_step; 654 for (i=0; i<s; i++) { 655 rk->stage_time = t + h*c[i]; 656 ierr = TSPreStage(ts,rk->stage_time); CHKERRQ(ierr); 657 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 658 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 659 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 660 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 661 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 662 ierr = TSAdaptCheckStage(adapt,ts,rk->stage_time,Y[i],&stageok);CHKERRQ(ierr); 663 if (!stageok) goto reject_step; 664 if (FSAL && !i) continue; 665 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 666 } 667 668 rk->status = TS_STEP_INCOMPLETE; 669 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 670 rk->status = TS_STEP_PENDING; 671 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 672 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 673 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);CHKERRQ(ierr); 674 ierr = TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr); 675 rk->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 676 if (!accept) { /* Roll back the current step */ 677 ierr = TSRollBack_RK(ts);CHKERRQ(ierr); 678 ts->time_step = next_time_step; 679 goto reject_step; 680 } 681 682 if (ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation */ 683 rk->ptime = ts->ptime; 684 rk->time_step = ts->time_step; 685 } 686 687 ts->ptime += ts->time_step; 688 ts->time_step = next_time_step; 689 break; 690 691 reject_step: 692 ts->reject++; accept = PETSC_FALSE; 693 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 694 ts->reason = TS_DIVERGED_STEP_REJECTED; 695 ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr); 696 } 697 } 698 PetscFunctionReturn(0); 699 } 700 701 static PetscErrorCode TSAdjointSetUp_RK(TS ts) 702 { 703 TS_RK *rk = (TS_RK*)ts->data; 704 RKTableau tab = rk->tableau; 705 PetscInt s = tab->s; 706 PetscErrorCode ierr; 707 708 PetscFunctionBegin; 709 if (ts->adjointsetupcalled++) PetscFunctionReturn(0); 710 ierr = VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecsDeltaLam);CHKERRQ(ierr); 711 ierr = VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&rk->VecsSensiTemp);CHKERRQ(ierr); 712 if(ts->vecs_sensip) { 713 ierr = VecDuplicate(ts->vecs_sensip[0],&rk->VecDeltaMu);CHKERRQ(ierr); 714 } 715 if (ts->vecs_sensi2) { 716 ierr = VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecsDeltaLam2);CHKERRQ(ierr); 717 ierr = VecDuplicateVecs(ts->vecs_sensi2[0],ts->numcost,&rk->VecsSensi2Temp);CHKERRQ(ierr); 718 } 719 if (ts->vecs_sensi2p) { 720 ierr = VecDuplicate(ts->vecs_sensi2p[0],&rk->VecDeltaMu2);CHKERRQ(ierr); 721 } 722 PetscFunctionReturn(0); 723 } 724 725 /* 726 Assumptions: 727 - TSStep_RK() always evaluates the step with b, not bembed. 728 */ 729 static PetscErrorCode TSAdjointStep_RK(TS ts) 730 { 731 TS_RK *rk = (TS_RK*)ts->data; 732 TS quadts = ts->quadraturets; 733 RKTableau tab = rk->tableau; 734 Mat J,Jquad; 735 const PetscInt s = tab->s; 736 const PetscReal *A = tab->A,*b = tab->b,*c = tab->c; 737 PetscScalar *w = rk->work,*xarr; 738 Vec *Y = rk->Y,*VecsDeltaLam = rk->VecsDeltaLam,VecDeltaMu = rk->VecDeltaMu,*VecsSensiTemp = rk->VecsSensiTemp; 739 Vec *VecsDeltaLam2 = rk->VecsDeltaLam2,VecDeltaMu2 = rk->VecDeltaMu2,*VecsSensi2Temp = rk->VecsSensi2Temp; 740 Vec VecDRDUTransCol = ts->vec_drdu_col,VecDRDPTransCol = ts->vec_drdp_col; 741 PetscInt i,j,nadj; 742 PetscReal t = ts->ptime; 743 PetscReal h = ts->time_step,stage_time; 744 PetscErrorCode ierr; 745 746 PetscFunctionBegin; 747 rk->status = TS_STEP_INCOMPLETE; 748 749 ierr = TSGetRHSJacobian(ts,&J,NULL,NULL,NULL);CHKERRQ(ierr); 750 if (quadts) { 751 ierr = TSGetRHSJacobian(quadts,&Jquad,NULL,NULL,NULL);CHKERRQ(ierr); 752 } 753 for (i=s-1; i>=0; i--) { 754 if (tab->FSAL && i == s-1) { 755 /* VecsDeltaLam[nadj*s+s-1] are initialized with zeros and the values never change.*/ 756 continue; 757 } 758 stage_time = t + h*(1.0-c[i]); 759 ierr = TSComputeRHSJacobian(ts,stage_time,Y[i],J,J);CHKERRQ(ierr); 760 if (quadts) { 761 ierr = TSComputeRHSJacobian(quadts,stage_time,Y[i],Jquad,Jquad);CHKERRQ(ierr); /* get r_u^T */ 762 } 763 if (ts->vecs_sensip) { 764 ierr = TSComputeRHSJacobianP(ts,stage_time,Y[i],ts->Jacprhs);CHKERRQ(ierr); /* get f_p */ 765 if (quadts) { 766 ierr = TSComputeRHSJacobianP(quadts,stage_time,Y[i],quadts->Jacprhs);CHKERRQ(ierr); /* get f_p for the quadrature */ 767 } 768 } 769 770 if (b[i]) { 771 for (j=i+1; j<s; j++) w[j-i-1] = A[j*s+i]/b[i]; /* coefficients for computing VecsSensiTemp */ 772 } else { 773 for (j=i+1; j<s; j++) w[j-i-1] = A[j*s+i]; /* coefficients for computing VecsSensiTemp */ 774 } 775 776 for (nadj=0; nadj<ts->numcost; nadj++) { 777 /* Stage values of lambda */ 778 if (b[i]) { 779 /* lambda_{n+1} + \sum_{j=i+1}^s a_{ji}/b[i]*lambda_{s,j} */ 780 ierr = VecCopy(ts->vecs_sensi[nadj],VecsSensiTemp[nadj]);CHKERRQ(ierr); /* VecDeltaLam is an vec array of size s by numcost */ 781 ierr = VecMAXPY(VecsSensiTemp[nadj],s-i-1,w,&VecsDeltaLam[nadj*s+i+1]);CHKERRQ(ierr); 782 ierr = MatMultTranspose(J,VecsSensiTemp[nadj],VecsDeltaLam[nadj*s+i]);CHKERRQ(ierr); /* VecsSensiTemp will be reused by 2nd-order adjoint */ 783 ierr = VecScale(VecsDeltaLam[nadj*s+i],-h*b[i]);CHKERRQ(ierr); 784 if (quadts) { 785 ierr = MatDenseGetColumn(Jquad,nadj,&xarr);CHKERRQ(ierr); 786 ierr = VecPlaceArray(VecDRDUTransCol,xarr);CHKERRQ(ierr); 787 ierr = VecAXPY(VecsDeltaLam[nadj*s+i],-h*b[i],VecDRDUTransCol);CHKERRQ(ierr); 788 ierr = VecResetArray(VecDRDUTransCol);CHKERRQ(ierr); 789 ierr = MatDenseRestoreColumn(Jquad,&xarr);CHKERRQ(ierr); 790 } 791 } else { 792 /* \sum_{j=i+1}^s a_{ji}*lambda_{s,j} */ 793 ierr = VecSet(VecsSensiTemp[nadj],0);CHKERRQ(ierr); 794 ierr = VecMAXPY(VecsSensiTemp[nadj],s-i-1,w,&VecsDeltaLam[nadj*s+i+1]);CHKERRQ(ierr); 795 ierr = MatMultTranspose(J,VecsSensiTemp[nadj],VecsDeltaLam[nadj*s+i]);CHKERRQ(ierr); 796 ierr = VecScale(VecsDeltaLam[nadj*s+i],-h);CHKERRQ(ierr); 797 } 798 799 /* Stage values of mu */ 800 if (ts->vecs_sensip) { 801 ierr = MatMultTranspose(ts->Jacprhs,VecsSensiTemp[nadj],VecDeltaMu);CHKERRQ(ierr); 802 if (b[i]) { 803 ierr = VecScale(VecDeltaMu,-h*b[i]);CHKERRQ(ierr); 804 if (quadts) { 805 ierr = MatDenseGetColumn(quadts->Jacprhs,nadj,&xarr);CHKERRQ(ierr); 806 ierr = VecPlaceArray(VecDRDPTransCol,xarr);CHKERRQ(ierr); 807 ierr = VecAXPY(VecDeltaMu,-h*b[i],VecDRDPTransCol);CHKERRQ(ierr); 808 ierr = VecResetArray(VecDRDPTransCol);CHKERRQ(ierr); 809 ierr = MatDenseRestoreColumn(quadts->Jacprhs,&xarr);CHKERRQ(ierr); 810 } 811 } else { 812 ierr = VecScale(VecDeltaMu,-h);CHKERRQ(ierr); 813 } 814 ierr = VecAXPY(ts->vecs_sensip[nadj],1.,VecDeltaMu);CHKERRQ(ierr); /* update sensip for each stage */ 815 } 816 } 817 818 if (ts->vecs_sensi2 && ts->forward_solve) { /* 2nd-order adjoint, TLM mode has to be turned on */ 819 /* Get w1 at t_{n+1} from TLM matrix */ 820 ierr = MatDenseGetColumn(rk->MatsFwdStageSensip[i],0,&xarr);CHKERRQ(ierr); 821 ierr = VecPlaceArray(ts->vec_sensip_col,xarr);CHKERRQ(ierr); 822 /* lambda_s^T F_UU w_1 */ 823 ierr = TSComputeRHSHessianProductFunctionUU(ts,stage_time,Y[i],VecsSensiTemp,ts->vec_sensip_col,ts->vecs_guu);CHKERRQ(ierr); 824 if (quadts) { 825 /* R_UU w_1 */ 826 ierr = TSComputeRHSHessianProductFunctionUU(quadts,stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_guu);CHKERRQ(ierr); 827 } 828 if (ts->vecs_sensip) { 829 /* lambda_s^T F_UP w_2 */ 830 ierr = TSComputeRHSHessianProductFunctionUP(ts,stage_time,Y[i],VecsSensiTemp,ts->vec_dir,ts->vecs_gup);CHKERRQ(ierr); 831 if (quadts) { 832 /* R_UP w_2 */ 833 ierr = TSComputeRHSHessianProductFunctionUP(quadts,stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_gup);CHKERRQ(ierr); 834 } 835 } 836 if (ts->vecs_sensi2p) { 837 /* lambda_s^T F_PU w_1 */ 838 ierr = TSComputeRHSHessianProductFunctionPU(ts,stage_time,Y[i],VecsSensiTemp,ts->vec_sensip_col,ts->vecs_gpu);CHKERRQ(ierr); 839 /* lambda_s^T F_PP w_2 */ 840 ierr = TSComputeRHSHessianProductFunctionPP(ts,stage_time,Y[i],VecsSensiTemp,ts->vec_dir,ts->vecs_gpp);CHKERRQ(ierr); 841 if (b[i] && quadts) { 842 /* R_PU w_1 */ 843 ierr = TSComputeRHSHessianProductFunctionPU(quadts,stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_gpu);CHKERRQ(ierr); 844 /* R_PP w_2 */ 845 ierr = TSComputeRHSHessianProductFunctionPP(quadts,stage_time,Y[i],NULL,ts->vec_dir,ts->vecs_gpp);CHKERRQ(ierr); 846 } 847 } 848 ierr = VecResetArray(ts->vec_sensip_col);CHKERRQ(ierr); 849 ierr = MatDenseRestoreColumn(rk->MatsFwdStageSensip[i],&xarr);CHKERRQ(ierr); 850 851 for (nadj=0; nadj<ts->numcost; nadj++) { 852 /* Stage values of lambda */ 853 if (b[i]) { 854 /* J_i^T*(Lambda_{n+1}+\sum_{j=i+1}^s a_{ji}/b_i*Lambda_{s,j} */ 855 ierr = VecCopy(ts->vecs_sensi2[nadj],VecsSensi2Temp[nadj]);CHKERRQ(ierr); 856 ierr = VecMAXPY(VecsSensi2Temp[nadj],s-i-1,w,&VecsDeltaLam2[nadj*s+i+1]);CHKERRQ(ierr); 857 ierr = MatMultTranspose(J,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj*s+i]);CHKERRQ(ierr); 858 ierr = VecScale(VecsDeltaLam2[nadj*s+i],-h*b[i]);CHKERRQ(ierr); 859 ierr = VecAXPY(VecsDeltaLam2[nadj*s+i],-h*b[i],ts->vecs_guu[nadj]);CHKERRQ(ierr); 860 if (ts->vecs_sensip) { 861 ierr = VecAXPY(VecsDeltaLam2[nadj*s+i],-h*b[i],ts->vecs_gup[nadj]);CHKERRQ(ierr); 862 } 863 } else { 864 /* \sum_{j=i+1}^s a_{ji}*Lambda_{s,j} */ 865 ierr = VecSet(VecsDeltaLam2[nadj*s+i],0);CHKERRQ(ierr); 866 ierr = VecMAXPY(VecsSensi2Temp[nadj],s-i-1,w,&VecsDeltaLam2[nadj*s+i+1]);CHKERRQ(ierr); 867 ierr = MatMultTranspose(J,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj*s+i]);CHKERRQ(ierr); 868 ierr = VecScale(VecsDeltaLam2[nadj*s+i],-h);CHKERRQ(ierr); 869 ierr = VecAXPY(VecsDeltaLam2[nadj*s+i],-h,ts->vecs_guu[nadj]);CHKERRQ(ierr); 870 if (ts->vecs_sensip) { 871 ierr = VecAXPY(VecsDeltaLam2[nadj*s+i],-h,ts->vecs_gup[nadj]);CHKERRQ(ierr); 872 } 873 } 874 if (ts->vecs_sensi2p) { /* 2nd-order adjoint for parameters */ 875 ierr = MatMultTranspose(ts->Jacprhs,VecsSensi2Temp[nadj],VecDeltaMu2);CHKERRQ(ierr); 876 if (b[i]) { 877 ierr = VecScale(VecDeltaMu2,-h*b[i]);CHKERRQ(ierr); 878 ierr = VecAXPY(VecDeltaMu2,-h*b[i],ts->vecs_gpu[nadj]);CHKERRQ(ierr); 879 ierr = VecAXPY(VecDeltaMu2,-h*b[i],ts->vecs_gpp[nadj]);CHKERRQ(ierr); 880 } else { 881 ierr = VecScale(VecDeltaMu2,-h);CHKERRQ(ierr); 882 ierr = VecAXPY(VecDeltaMu2,-h,ts->vecs_gpu[nadj]);CHKERRQ(ierr); 883 ierr = VecAXPY(VecDeltaMu2,-h,ts->vecs_gpp[nadj]);CHKERRQ(ierr); 884 } 885 ierr = VecAXPY(ts->vecs_sensi2p[nadj],1,VecDeltaMu2);CHKERRQ(ierr); /* update sensi2p for each stage */ 886 } 887 } 888 } 889 } 890 891 for (j=0; j<s; j++) w[j] = 1.0; 892 for (nadj=0; nadj<ts->numcost; nadj++) { /* no need to do this for mu's */ 893 ierr = VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecsDeltaLam[nadj*s]);CHKERRQ(ierr); 894 if (ts->vecs_sensi2) { 895 ierr = VecMAXPY(ts->vecs_sensi2[nadj],s,w,&VecsDeltaLam2[nadj*s]);CHKERRQ(ierr); 896 } 897 } 898 rk->status = TS_STEP_COMPLETE; 899 PetscFunctionReturn(0); 900 } 901 902 static PetscErrorCode TSAdjointReset_RK(TS ts) 903 { 904 TS_RK *rk = (TS_RK*)ts->data; 905 RKTableau tab = rk->tableau; 906 PetscErrorCode ierr; 907 908 PetscFunctionBegin; 909 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam);CHKERRQ(ierr); 910 ierr = VecDestroyVecs(ts->numcost,&rk->VecsSensiTemp);CHKERRQ(ierr); 911 ierr = VecDestroy(&rk->VecDeltaMu);CHKERRQ(ierr); 912 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam2);CHKERRQ(ierr); 913 ierr = VecDestroy(&rk->VecDeltaMu2);CHKERRQ(ierr); 914 ierr = VecDestroyVecs(ts->numcost,&rk->VecsSensi2Temp);CHKERRQ(ierr); 915 PetscFunctionReturn(0); 916 } 917 918 static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X) 919 { 920 TS_RK *rk = (TS_RK*)ts->data; 921 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 922 PetscReal h; 923 PetscReal tt,t; 924 PetscScalar *b; 925 const PetscReal *B = rk->tableau->binterp; 926 PetscErrorCode ierr; 927 928 PetscFunctionBegin; 929 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 930 931 switch (rk->status) { 932 case TS_STEP_INCOMPLETE: 933 case TS_STEP_PENDING: 934 h = ts->time_step; 935 t = (itime - ts->ptime)/h; 936 break; 937 case TS_STEP_COMPLETE: 938 h = ts->ptime - ts->ptime_prev; 939 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 940 break; 941 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 942 } 943 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 944 for (i=0; i<s; i++) b[i] = 0; 945 for (j=0,tt=t; j<p; j++,tt*=t) { 946 for (i=0; i<s; i++) { 947 b[i] += h * B[i*p+j] * tt; 948 } 949 } 950 ierr = VecCopy(rk->Y[0],X);CHKERRQ(ierr); 951 ierr = VecMAXPY(X,s,b,rk->YdotRHS);CHKERRQ(ierr); 952 ierr = PetscFree(b);CHKERRQ(ierr); 953 PetscFunctionReturn(0); 954 } 955 956 /*------------------------------------------------------------*/ 957 958 static PetscErrorCode TSRKTableauReset(TS ts) 959 { 960 TS_RK *rk = (TS_RK*)ts->data; 961 RKTableau tab = rk->tableau; 962 PetscErrorCode ierr; 963 964 PetscFunctionBegin; 965 if (!tab) PetscFunctionReturn(0); 966 ierr = PetscFree(rk->work);CHKERRQ(ierr); 967 ierr = VecDestroyVecs(tab->s,&rk->Y);CHKERRQ(ierr); 968 ierr = VecDestroyVecs(tab->s,&rk->YdotRHS);CHKERRQ(ierr); 969 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam);CHKERRQ(ierr); 970 ierr = VecDestroyVecs(ts->numcost,&rk->VecsSensiTemp);CHKERRQ(ierr); 971 ierr = VecDestroy(&rk->VecDeltaMu);CHKERRQ(ierr); 972 PetscFunctionReturn(0); 973 } 974 975 static PetscErrorCode TSReset_RK(TS ts) 976 { 977 PetscErrorCode ierr; 978 979 PetscFunctionBegin; 980 ierr = TSRKTableauReset(ts);CHKERRQ(ierr); 981 if (ts->use_splitrhsfunction) { 982 ierr = PetscTryMethod(ts,"TSReset_RK_MultirateSplit_C",(TS),(ts));CHKERRQ(ierr); 983 } else { 984 ierr = PetscTryMethod(ts,"TSReset_RK_MultirateNonsplit_C",(TS),(ts));CHKERRQ(ierr); 985 } 986 PetscFunctionReturn(0); 987 } 988 989 static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx) 990 { 991 PetscFunctionBegin; 992 PetscFunctionReturn(0); 993 } 994 995 static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 996 { 997 PetscFunctionBegin; 998 PetscFunctionReturn(0); 999 } 1000 1001 1002 static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx) 1003 { 1004 PetscFunctionBegin; 1005 PetscFunctionReturn(0); 1006 } 1007 1008 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1009 { 1010 1011 PetscFunctionBegin; 1012 PetscFunctionReturn(0); 1013 } 1014 /* 1015 static PetscErrorCode RKSetAdjCoe(RKTableau tab) 1016 { 1017 PetscReal *A,*b; 1018 PetscInt s,i,j; 1019 PetscErrorCode ierr; 1020 1021 PetscFunctionBegin; 1022 s = tab->s; 1023 ierr = PetscMalloc2(s*s,&A,s,&b);CHKERRQ(ierr); 1024 1025 for (i=0; i<s; i++) 1026 for (j=0; j<s; j++) { 1027 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]; 1028 ierr = PetscPrintf(PETSC_COMM_WORLD,"Coefficients: A[%D][%D]=%.6f\n",i,j,A[i*s+j]);CHKERRQ(ierr); 1029 } 1030 1031 for (i=0; i<s; i++) b[i] = (tab->b[s-1-i]==0)? 0: -tab->b[s-1-i]; 1032 1033 ierr = PetscMemcpy(tab->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 1034 ierr = PetscMemcpy(tab->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 1035 ierr = PetscFree2(A,b);CHKERRQ(ierr); 1036 PetscFunctionReturn(0); 1037 } 1038 */ 1039 1040 static PetscErrorCode TSRKTableauSetUp(TS ts) 1041 { 1042 TS_RK *rk = (TS_RK*)ts->data; 1043 RKTableau tab = rk->tableau; 1044 PetscErrorCode ierr; 1045 1046 PetscFunctionBegin; 1047 ierr = PetscMalloc1(tab->s,&rk->work);CHKERRQ(ierr); 1048 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->Y);CHKERRQ(ierr); 1049 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS);CHKERRQ(ierr); 1050 PetscFunctionReturn(0); 1051 } 1052 1053 static PetscErrorCode TSSetUp_RK(TS ts) 1054 { 1055 TS quadts = ts->quadraturets; 1056 PetscErrorCode ierr; 1057 DM dm; 1058 1059 PetscFunctionBegin; 1060 ierr = TSCheckImplicitTerm(ts);CHKERRQ(ierr); 1061 ierr = TSRKTableauSetUp(ts);CHKERRQ(ierr); 1062 if (quadts && ts->costintegralfwd) { 1063 Mat Jquad; 1064 ierr = TSGetRHSJacobian(quadts,&Jquad,NULL,NULL,NULL);CHKERRQ(ierr); 1065 } 1066 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1067 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 1068 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 1069 if (ts->use_splitrhsfunction) { 1070 ierr = PetscTryMethod(ts,"TSSetUp_RK_MultirateSplit_C",(TS),(ts));CHKERRQ(ierr); 1071 } else { 1072 ierr = PetscTryMethod(ts,"TSSetUp_RK_MultirateNonsplit_C",(TS),(ts));CHKERRQ(ierr); 1073 } 1074 PetscFunctionReturn(0); 1075 } 1076 1077 static PetscErrorCode TSSetFromOptions_RK(PetscOptionItems *PetscOptionsObject,TS ts) 1078 { 1079 TS_RK *rk = (TS_RK*)ts->data; 1080 PetscErrorCode ierr; 1081 1082 PetscFunctionBegin; 1083 ierr = PetscOptionsHead(PetscOptionsObject,"RK ODE solver options");CHKERRQ(ierr); 1084 { 1085 RKTableauLink link; 1086 PetscInt count,choice; 1087 PetscBool flg,use_multirate = PETSC_FALSE; 1088 const char **namelist; 1089 1090 for (link=RKTableauList,count=0; link; link=link->next,count++) ; 1091 ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 1092 for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1093 ierr = PetscOptionsBool("-ts_rk_multirate","Use interpolation-based multirate RK method","TSRKSetMultirate",rk->use_multirate,&use_multirate,&flg);CHKERRQ(ierr); 1094 if (flg) { 1095 ierr = TSRKSetMultirate(ts,use_multirate);CHKERRQ(ierr); 1096 } 1097 ierr = PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rk->tableau->name,&choice,&flg);CHKERRQ(ierr); 1098 if (flg) {ierr = TSRKSetType(ts,namelist[choice]);CHKERRQ(ierr);} 1099 ierr = PetscFree(namelist);CHKERRQ(ierr); 1100 } 1101 ierr = PetscOptionsTail();CHKERRQ(ierr); 1102 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Multirate methods options","");CHKERRQ(ierr); 1103 ierr = PetscOptionsInt("-ts_rk_dtratio","time step ratio between slow and fast","",rk->dtratio,&rk->dtratio,NULL);CHKERRQ(ierr); 1104 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1105 PetscFunctionReturn(0); 1106 } 1107 1108 static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer) 1109 { 1110 TS_RK *rk = (TS_RK*)ts->data; 1111 PetscBool iascii; 1112 PetscErrorCode ierr; 1113 1114 PetscFunctionBegin; 1115 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1116 if (iascii) { 1117 RKTableau tab = rk->tableau; 1118 TSRKType rktype; 1119 char buf[512]; 1120 ierr = TSRKGetType(ts,&rktype);CHKERRQ(ierr); 1121 ierr = PetscViewerASCIIPrintf(viewer," RK type %s\n",rktype);CHKERRQ(ierr); 1122 ierr = PetscViewerASCIIPrintf(viewer," Order: %D\n",tab->order);CHKERRQ(ierr); 1123 ierr = PetscViewerASCIIPrintf(viewer," FSAL property: %s\n",tab->FSAL ? "yes" : "no");CHKERRQ(ierr); 1124 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 1125 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 1126 } 1127 PetscFunctionReturn(0); 1128 } 1129 1130 static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer) 1131 { 1132 PetscErrorCode ierr; 1133 TSAdapt adapt; 1134 1135 PetscFunctionBegin; 1136 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1137 ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 1138 PetscFunctionReturn(0); 1139 } 1140 1141 /*@C 1142 TSRKSetType - Set the type of RK scheme 1143 1144 Logically collective 1145 1146 Input Parameter: 1147 + ts - timestepping context 1148 - rktype - type of RK-scheme 1149 1150 Options Database: 1151 . -ts_rk_type - <1fe,2a,3,3bs,4,5f,5dp,5bs> 1152 1153 Level: intermediate 1154 1155 .seealso: TSRKGetType(), TSRK, TSRKType, TSRK1FE, TSRK2A, TSRK3, TSRK3BS, TSRK4, TSRK5F, TSRK5DP, TSRK5BS 1156 @*/ 1157 PetscErrorCode TSRKSetType(TS ts,TSRKType rktype) 1158 { 1159 PetscErrorCode ierr; 1160 1161 PetscFunctionBegin; 1162 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1163 PetscValidCharPointer(rktype,2); 1164 ierr = PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));CHKERRQ(ierr); 1165 PetscFunctionReturn(0); 1166 } 1167 1168 /*@C 1169 TSRKGetType - Get the type of RK scheme 1170 1171 Logically collective 1172 1173 Input Parameter: 1174 . ts - timestepping context 1175 1176 Output Parameter: 1177 . rktype - type of RK-scheme 1178 1179 Level: intermediate 1180 1181 .seealso: TSRKGetType() 1182 @*/ 1183 PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype) 1184 { 1185 PetscErrorCode ierr; 1186 1187 PetscFunctionBegin; 1188 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1189 ierr = PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));CHKERRQ(ierr); 1190 PetscFunctionReturn(0); 1191 } 1192 1193 static PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype) 1194 { 1195 TS_RK *rk = (TS_RK*)ts->data; 1196 1197 PetscFunctionBegin; 1198 *rktype = rk->tableau->name; 1199 PetscFunctionReturn(0); 1200 } 1201 1202 static PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype) 1203 { 1204 TS_RK *rk = (TS_RK*)ts->data; 1205 PetscErrorCode ierr; 1206 PetscBool match; 1207 RKTableauLink link; 1208 1209 PetscFunctionBegin; 1210 if (rk->tableau) { 1211 ierr = PetscStrcmp(rk->tableau->name,rktype,&match);CHKERRQ(ierr); 1212 if (match) PetscFunctionReturn(0); 1213 } 1214 for (link = RKTableauList; link; link=link->next) { 1215 ierr = PetscStrcmp(link->tab.name,rktype,&match);CHKERRQ(ierr); 1216 if (match) { 1217 if (ts->setupcalled) {ierr = TSRKTableauReset(ts);CHKERRQ(ierr);} 1218 rk->tableau = &link->tab; 1219 if (ts->setupcalled) {ierr = TSRKTableauSetUp(ts);CHKERRQ(ierr);} 1220 ts->default_adapt_type = rk->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1221 PetscFunctionReturn(0); 1222 } 1223 } 1224 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype); 1225 PetscFunctionReturn(0); 1226 } 1227 1228 static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y) 1229 { 1230 TS_RK *rk = (TS_RK*)ts->data; 1231 1232 PetscFunctionBegin; 1233 if (ns) *ns = rk->tableau->s; 1234 if (Y) *Y = rk->Y; 1235 PetscFunctionReturn(0); 1236 } 1237 1238 static PetscErrorCode TSDestroy_RK(TS ts) 1239 { 1240 PetscErrorCode ierr; 1241 1242 PetscFunctionBegin; 1243 ierr = TSReset_RK(ts);CHKERRQ(ierr); 1244 if (ts->dm) { 1245 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 1246 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 1247 } 1248 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1249 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);CHKERRQ(ierr); 1250 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);CHKERRQ(ierr); 1251 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirate_C",NULL);CHKERRQ(ierr); 1252 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetMultirate_C",NULL);CHKERRQ(ierr); 1253 PetscFunctionReturn(0); 1254 } 1255 1256 /*@C 1257 TSRKSetMultirate - Use the interpolation-based multirate RK method 1258 1259 Logically collective 1260 1261 Input Parameter: 1262 + ts - timestepping context 1263 - use_multirate - PETSC_TRUE enables the multirate RK method, sets the basic method to be RK2A and sets the ratio between slow stepsize and fast stepsize to be 2 1264 1265 Options Database: 1266 . -ts_rk_multirate - <true,false> 1267 1268 Notes: 1269 The multirate method requires interpolation. The default interpolation works for 1st- and 2nd- order RK, but not for high-order RKs except TSRK5DP which comes with the interpolation coeffcients (binterp). 1270 1271 Level: intermediate 1272 1273 .seealso: TSRKGetMultirate() 1274 @*/ 1275 PetscErrorCode TSRKSetMultirate(TS ts,PetscBool use_multirate) 1276 { 1277 PetscErrorCode ierr; 1278 1279 PetscFunctionBegin; 1280 ierr = PetscTryMethod(ts,"TSRKSetMultirate_C",(TS,PetscBool),(ts,use_multirate));CHKERRQ(ierr); 1281 PetscFunctionReturn(0); 1282 } 1283 1284 /*@C 1285 TSRKGetMultirate - Gets whether to Use the interpolation-based multirate RK method 1286 1287 Not collective 1288 1289 Input Parameter: 1290 . ts - timestepping context 1291 1292 Output Parameter: 1293 . use_multirate - PETSC_TRUE if the multirate RK method is enabled, PETSC_FALSE otherwise 1294 1295 Level: intermediate 1296 1297 .seealso: TSRKSetMultirate() 1298 @*/ 1299 PetscErrorCode TSRKGetMultirate(TS ts,PetscBool *use_multirate) 1300 { 1301 PetscErrorCode ierr; 1302 1303 PetscFunctionBegin; 1304 ierr = PetscUseMethod(ts,"TSRKGetMultirate_C",(TS,PetscBool*),(ts,use_multirate));CHKERRQ(ierr); 1305 PetscFunctionReturn(0); 1306 } 1307 1308 /*MC 1309 TSRK - ODE and DAE solver using Runge-Kutta schemes 1310 1311 The user should provide the right hand side of the equation 1312 using TSSetRHSFunction(). 1313 1314 Notes: 1315 The default is TSRK3BS, it can be changed with TSRKSetType() or -ts_rk_type 1316 1317 Level: beginner 1318 1319 .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3, 1320 TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister(), TSRKSetMultirate(), TSRKGetMultirate() 1321 1322 M*/ 1323 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 1324 { 1325 TS_RK *rk; 1326 PetscErrorCode ierr; 1327 1328 PetscFunctionBegin; 1329 ierr = TSRKInitializePackage();CHKERRQ(ierr); 1330 1331 ts->ops->reset = TSReset_RK; 1332 ts->ops->destroy = TSDestroy_RK; 1333 ts->ops->view = TSView_RK; 1334 ts->ops->load = TSLoad_RK; 1335 ts->ops->setup = TSSetUp_RK; 1336 ts->ops->interpolate = TSInterpolate_RK; 1337 ts->ops->step = TSStep_RK; 1338 ts->ops->evaluatestep = TSEvaluateStep_RK; 1339 ts->ops->rollback = TSRollBack_RK; 1340 ts->ops->setfromoptions = TSSetFromOptions_RK; 1341 ts->ops->getstages = TSGetStages_RK; 1342 1343 ts->ops->adjointintegral = TSAdjointCostIntegral_RK; 1344 ts->ops->adjointsetup = TSAdjointSetUp_RK; 1345 ts->ops->adjointstep = TSAdjointStep_RK; 1346 ts->ops->adjointreset = TSAdjointReset_RK; 1347 1348 ts->ops->forwardintegral = TSForwardCostIntegral_RK; 1349 ts->ops->forwardsetup = TSForwardSetUp_RK; 1350 ts->ops->forwardreset = TSForwardReset_RK; 1351 ts->ops->forwardstep = TSForwardStep_RK; 1352 ts->ops->forwardgetstages= TSForwardGetStages_RK; 1353 1354 ierr = PetscNewLog(ts,&rk);CHKERRQ(ierr); 1355 ts->data = (void*)rk; 1356 1357 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);CHKERRQ(ierr); 1358 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);CHKERRQ(ierr); 1359 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirate_C",TSRKSetMultirate_RK);CHKERRQ(ierr); 1360 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetMultirate_C",TSRKGetMultirate_RK);CHKERRQ(ierr); 1361 1362 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 1363 rk->dtratio = 1; 1364 PetscFunctionReturn(0); 1365 } 1366