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 if (rk->MatFwdSensip0) { /* Unlike VecDestroy, MatDestroy does not check if the pointer is NULL */ 614 ierr = MatDestroy(&rk->MatFwdSensip0);CHKERRQ(ierr); 615 } 616 if (rk->MatsFwdStageSensip) { 617 for (i=0; i<tab->s; i++) { 618 ierr = MatDestroy(&rk->MatsFwdStageSensip[i]);CHKERRQ(ierr); 619 } 620 ierr = PetscFree(rk->MatsFwdStageSensip);CHKERRQ(ierr); 621 } 622 if (rk->MatsFwdSensipTemp) { 623 for (i=0; i<tab->s; i++) { 624 ierr = MatDestroy(&rk->MatsFwdSensipTemp[i]);CHKERRQ(ierr); 625 } 626 ierr = PetscFree(rk->MatsFwdSensipTemp);CHKERRQ(ierr); 627 } 628 ierr = VecDestroy(&rk->VecDeltaFwdSensipCol);CHKERRQ(ierr); 629 PetscFunctionReturn(0); 630 } 631 632 static PetscErrorCode TSStep_RK(TS ts) 633 { 634 TS_RK *rk = (TS_RK*)ts->data; 635 RKTableau tab = rk->tableau; 636 const PetscInt s = tab->s; 637 const PetscReal *A = tab->A,*c = tab->c; 638 PetscScalar *w = rk->work; 639 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 640 PetscBool FSAL = tab->FSAL; 641 TSAdapt adapt; 642 PetscInt i,j; 643 PetscInt rejections = 0; 644 PetscBool stageok,accept = PETSC_TRUE; 645 PetscReal next_time_step = ts->time_step; 646 PetscErrorCode ierr; 647 648 PetscFunctionBegin; 649 if (ts->steprollback || ts->steprestart) FSAL = PETSC_FALSE; 650 if (FSAL) { ierr = VecCopy(YdotRHS[s-1],YdotRHS[0]);CHKERRQ(ierr); } 651 652 rk->status = TS_STEP_INCOMPLETE; 653 while (!ts->reason && rk->status != TS_STEP_COMPLETE) { 654 PetscReal t = ts->ptime; 655 PetscReal h = ts->time_step; 656 for (i=0; i<s; i++) { 657 rk->stage_time = t + h*c[i]; 658 ierr = TSPreStage(ts,rk->stage_time); CHKERRQ(ierr); 659 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 660 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 661 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 662 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 663 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 664 ierr = TSAdaptCheckStage(adapt,ts,rk->stage_time,Y[i],&stageok);CHKERRQ(ierr); 665 if (!stageok) goto reject_step; 666 if (FSAL && !i) continue; 667 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 668 } 669 670 rk->status = TS_STEP_INCOMPLETE; 671 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 672 rk->status = TS_STEP_PENDING; 673 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 674 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 675 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);CHKERRQ(ierr); 676 ierr = TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr); 677 rk->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 678 if (!accept) { /* Roll back the current step */ 679 ierr = TSRollBack_RK(ts);CHKERRQ(ierr); 680 ts->time_step = next_time_step; 681 goto reject_step; 682 } 683 684 if (ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation */ 685 rk->ptime = ts->ptime; 686 rk->time_step = ts->time_step; 687 } 688 689 ts->ptime += ts->time_step; 690 ts->time_step = next_time_step; 691 break; 692 693 reject_step: 694 ts->reject++; accept = PETSC_FALSE; 695 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 696 ts->reason = TS_DIVERGED_STEP_REJECTED; 697 ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr); 698 } 699 } 700 PetscFunctionReturn(0); 701 } 702 703 static PetscErrorCode TSAdjointSetUp_RK(TS ts) 704 { 705 TS_RK *rk = (TS_RK*)ts->data; 706 RKTableau tab = rk->tableau; 707 PetscInt s = tab->s; 708 PetscErrorCode ierr; 709 710 PetscFunctionBegin; 711 if (ts->adjointsetupcalled++) PetscFunctionReturn(0); 712 ierr = VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecsDeltaLam);CHKERRQ(ierr); 713 ierr = VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&rk->VecsSensiTemp);CHKERRQ(ierr); 714 if(ts->vecs_sensip) { 715 ierr = VecDuplicate(ts->vecs_sensip[0],&rk->VecDeltaMu);CHKERRQ(ierr); 716 } 717 if (ts->vecs_sensi2) { 718 ierr = VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecsDeltaLam2);CHKERRQ(ierr); 719 ierr = VecDuplicateVecs(ts->vecs_sensi2[0],ts->numcost,&rk->VecsSensi2Temp);CHKERRQ(ierr); 720 } 721 if (ts->vecs_sensi2p) { 722 ierr = VecDuplicate(ts->vecs_sensi2p[0],&rk->VecDeltaMu2);CHKERRQ(ierr); 723 } 724 PetscFunctionReturn(0); 725 } 726 727 /* 728 Assumptions: 729 - TSStep_RK() always evaluates the step with b, not bembed. 730 */ 731 static PetscErrorCode TSAdjointStep_RK(TS ts) 732 { 733 TS_RK *rk = (TS_RK*)ts->data; 734 TS quadts = ts->quadraturets; 735 RKTableau tab = rk->tableau; 736 Mat J,Jquad; 737 const PetscInt s = tab->s; 738 const PetscReal *A = tab->A,*b = tab->b,*c = tab->c; 739 PetscScalar *w = rk->work,*xarr; 740 Vec *Y = rk->Y,*VecsDeltaLam = rk->VecsDeltaLam,VecDeltaMu = rk->VecDeltaMu,*VecsSensiTemp = rk->VecsSensiTemp; 741 Vec *VecsDeltaLam2 = rk->VecsDeltaLam2,VecDeltaMu2 = rk->VecDeltaMu2,*VecsSensi2Temp = rk->VecsSensi2Temp; 742 Vec VecDRDUTransCol = ts->vec_drdu_col,VecDRDPTransCol = ts->vec_drdp_col; 743 PetscInt i,j,nadj; 744 PetscReal t = ts->ptime; 745 PetscReal h = ts->time_step,stage_time; 746 PetscErrorCode ierr; 747 748 PetscFunctionBegin; 749 rk->status = TS_STEP_INCOMPLETE; 750 751 ierr = TSGetRHSJacobian(ts,&J,NULL,NULL,NULL);CHKERRQ(ierr); 752 if (quadts) { 753 ierr = TSGetRHSJacobian(quadts,&Jquad,NULL,NULL,NULL);CHKERRQ(ierr); 754 } 755 for (i=s-1; i>=0; i--) { 756 if (tab->FSAL && i == s-1) { 757 /* VecsDeltaLam[nadj*s+s-1] are initialized with zeros and the values never change.*/ 758 continue; 759 } 760 stage_time = t + h*(1.0-c[i]); 761 ierr = TSComputeRHSJacobian(ts,stage_time,Y[i],J,J);CHKERRQ(ierr); 762 if (quadts) { 763 ierr = TSComputeRHSJacobian(quadts,stage_time,Y[i],Jquad,Jquad);CHKERRQ(ierr); /* get r_u^T */ 764 } 765 if (ts->vecs_sensip) { 766 ierr = TSComputeRHSJacobianP(ts,stage_time,Y[i],ts->Jacprhs);CHKERRQ(ierr); /* get f_p */ 767 if (quadts) { 768 ierr = TSComputeRHSJacobianP(quadts,stage_time,Y[i],quadts->Jacprhs);CHKERRQ(ierr); /* get f_p for the quadrature */ 769 } 770 } 771 772 if (b[i]) { 773 for (j=i+1; j<s; j++) w[j-i-1] = A[j*s+i]/b[i]; /* coefficients for computing VecsSensiTemp */ 774 } else { 775 for (j=i+1; j<s; j++) w[j-i-1] = A[j*s+i]; /* coefficients for computing VecsSensiTemp */ 776 } 777 778 for (nadj=0; nadj<ts->numcost; nadj++) { 779 /* Stage values of lambda */ 780 if (b[i]) { 781 /* lambda_{n+1} + \sum_{j=i+1}^s a_{ji}/b[i]*lambda_{s,j} */ 782 ierr = VecCopy(ts->vecs_sensi[nadj],VecsSensiTemp[nadj]);CHKERRQ(ierr); /* VecDeltaLam is an vec array of size s by numcost */ 783 ierr = VecMAXPY(VecsSensiTemp[nadj],s-i-1,w,&VecsDeltaLam[nadj*s+i+1]);CHKERRQ(ierr); 784 ierr = MatMultTranspose(J,VecsSensiTemp[nadj],VecsDeltaLam[nadj*s+i]);CHKERRQ(ierr); /* VecsSensiTemp will be reused by 2nd-order adjoint */ 785 ierr = VecScale(VecsDeltaLam[nadj*s+i],-h*b[i]);CHKERRQ(ierr); 786 if (quadts) { 787 ierr = MatDenseGetColumn(Jquad,nadj,&xarr);CHKERRQ(ierr); 788 ierr = VecPlaceArray(VecDRDUTransCol,xarr);CHKERRQ(ierr); 789 ierr = VecAXPY(VecsDeltaLam[nadj*s+i],-h*b[i],VecDRDUTransCol);CHKERRQ(ierr); 790 ierr = VecResetArray(VecDRDUTransCol);CHKERRQ(ierr); 791 ierr = MatDenseRestoreColumn(Jquad,&xarr);CHKERRQ(ierr); 792 } 793 } else { 794 /* \sum_{j=i+1}^s a_{ji}*lambda_{s,j} */ 795 ierr = VecSet(VecsSensiTemp[nadj],0);CHKERRQ(ierr); 796 ierr = VecMAXPY(VecsSensiTemp[nadj],s-i-1,w,&VecsDeltaLam[nadj*s+i+1]);CHKERRQ(ierr); 797 ierr = MatMultTranspose(J,VecsSensiTemp[nadj],VecsDeltaLam[nadj*s+i]);CHKERRQ(ierr); 798 ierr = VecScale(VecsDeltaLam[nadj*s+i],-h);CHKERRQ(ierr); 799 } 800 801 /* Stage values of mu */ 802 if (ts->vecs_sensip) { 803 ierr = MatMultTranspose(ts->Jacprhs,VecsSensiTemp[nadj],VecDeltaMu);CHKERRQ(ierr); 804 if (b[i]) { 805 ierr = VecScale(VecDeltaMu,-h*b[i]);CHKERRQ(ierr); 806 if (quadts) { 807 ierr = MatDenseGetColumn(quadts->Jacprhs,nadj,&xarr);CHKERRQ(ierr); 808 ierr = VecPlaceArray(VecDRDPTransCol,xarr);CHKERRQ(ierr); 809 ierr = VecAXPY(VecDeltaMu,-h*b[i],VecDRDPTransCol);CHKERRQ(ierr); 810 ierr = VecResetArray(VecDRDPTransCol);CHKERRQ(ierr); 811 ierr = MatDenseRestoreColumn(quadts->Jacprhs,&xarr);CHKERRQ(ierr); 812 } 813 } else { 814 ierr = VecScale(VecDeltaMu,-h);CHKERRQ(ierr); 815 } 816 ierr = VecAXPY(ts->vecs_sensip[nadj],1.,VecDeltaMu);CHKERRQ(ierr); /* update sensip for each stage */ 817 } 818 } 819 820 if (ts->vecs_sensi2 && ts->forward_solve) { /* 2nd-order adjoint, TLM mode has to be turned on */ 821 /* Get w1 at t_{n+1} from TLM matrix */ 822 ierr = MatDenseGetColumn(rk->MatsFwdStageSensip[i],0,&xarr);CHKERRQ(ierr); 823 ierr = VecPlaceArray(ts->vec_sensip_col,xarr);CHKERRQ(ierr); 824 /* lambda_s^T F_UU w_1 */ 825 ierr = TSComputeRHSHessianProductFunction1(ts,stage_time,Y[i],VecsSensiTemp,ts->vec_sensip_col,ts->vecs_guu);CHKERRQ(ierr); 826 if (quadts) { 827 /* R_UU w_1 */ 828 ierr = TSComputeRHSHessianProductFunction1(quadts,stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_guu);CHKERRQ(ierr); 829 } 830 if (ts->vecs_sensip) { 831 /* lambda_s^T F_UP w_2 */ 832 ierr = TSComputeRHSHessianProductFunction2(ts,stage_time,Y[i],VecsSensiTemp,ts->vec_dir,ts->vecs_gup);CHKERRQ(ierr); 833 if (quadts) { 834 /* R_UP w_2 */ 835 ierr = TSComputeRHSHessianProductFunction2(quadts,stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_gup);CHKERRQ(ierr); 836 } 837 } 838 if (ts->vecs_sensi2p) { 839 /* lambda_s^T F_PU w_1 */ 840 ierr = TSComputeRHSHessianProductFunction3(ts,stage_time,Y[i],VecsSensiTemp,ts->vec_sensip_col,ts->vecs_gpu);CHKERRQ(ierr); 841 /* lambda_s^T F_PP w_2 */ 842 ierr = TSComputeRHSHessianProductFunction4(ts,stage_time,Y[i],VecsSensiTemp,ts->vec_dir,ts->vecs_gpp);CHKERRQ(ierr); 843 if (b[i] && quadts) { 844 /* R_PU w_1 */ 845 ierr = TSComputeRHSHessianProductFunction3(quadts,stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_gpu);CHKERRQ(ierr); 846 /* R_PP w_2 */ 847 ierr = TSComputeRHSHessianProductFunction4(quadts,stage_time,Y[i],NULL,ts->vec_dir,ts->vecs_gpp);CHKERRQ(ierr); 848 } 849 } 850 ierr = VecResetArray(ts->vec_sensip_col);CHKERRQ(ierr); 851 ierr = MatDenseRestoreColumn(rk->MatsFwdStageSensip[i],&xarr);CHKERRQ(ierr); 852 853 for (nadj=0; nadj<ts->numcost; nadj++) { 854 /* Stage values of lambda */ 855 if (b[i]) { 856 /* J_i^T*(Lambda_{n+1}+\sum_{j=i+1}^s a_{ji}/b_i*Lambda_{s,j} */ 857 ierr = VecCopy(ts->vecs_sensi2[nadj],VecsSensi2Temp[nadj]);CHKERRQ(ierr); 858 ierr = VecMAXPY(VecsSensi2Temp[nadj],s-i-1,w,&VecsDeltaLam2[nadj*s+i+1]);CHKERRQ(ierr); 859 ierr = MatMultTranspose(J,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj*s+i]);CHKERRQ(ierr); 860 ierr = VecScale(VecsDeltaLam2[nadj*s+i],-h*b[i]);CHKERRQ(ierr); 861 ierr = VecAXPY(VecsDeltaLam2[nadj*s+i],-h*b[i],ts->vecs_guu[nadj]);CHKERRQ(ierr); 862 if (ts->vecs_sensip) { 863 ierr = VecAXPY(VecsDeltaLam2[nadj*s+i],-h*b[i],ts->vecs_gup[nadj]);CHKERRQ(ierr); 864 } 865 } else { 866 /* \sum_{j=i+1}^s a_{ji}*Lambda_{s,j} */ 867 ierr = VecSet(VecsDeltaLam2[nadj*s+i],0);CHKERRQ(ierr); 868 ierr = VecMAXPY(VecsSensi2Temp[nadj],s-i-1,w,&VecsDeltaLam2[nadj*s+i+1]);CHKERRQ(ierr); 869 ierr = MatMultTranspose(J,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj*s+i]);CHKERRQ(ierr); 870 ierr = VecScale(VecsDeltaLam2[nadj*s+i],-h);CHKERRQ(ierr); 871 ierr = VecAXPY(VecsDeltaLam2[nadj*s+i],-h,ts->vecs_guu[nadj]);CHKERRQ(ierr); 872 if (ts->vecs_sensip) { 873 ierr = VecAXPY(VecsDeltaLam2[nadj*s+i],-h,ts->vecs_gup[nadj]);CHKERRQ(ierr); 874 } 875 } 876 if (ts->vecs_sensi2p) { /* 2nd-order adjoint for parameters */ 877 ierr = MatMultTranspose(ts->Jacprhs,VecsSensi2Temp[nadj],VecDeltaMu2);CHKERRQ(ierr); 878 if (b[i]) { 879 ierr = VecScale(VecDeltaMu2,-h*b[i]);CHKERRQ(ierr); 880 ierr = VecAXPY(VecDeltaMu2,-h*b[i],ts->vecs_gpu[nadj]);CHKERRQ(ierr); 881 ierr = VecAXPY(VecDeltaMu2,-h*b[i],ts->vecs_gpp[nadj]);CHKERRQ(ierr); 882 } else { 883 ierr = VecScale(VecDeltaMu2,-h);CHKERRQ(ierr); 884 ierr = VecAXPY(VecDeltaMu2,-h,ts->vecs_gpu[nadj]);CHKERRQ(ierr); 885 ierr = VecAXPY(VecDeltaMu2,-h,ts->vecs_gpp[nadj]);CHKERRQ(ierr); 886 } 887 ierr = VecAXPY(ts->vecs_sensi2p[nadj],1,VecDeltaMu2);CHKERRQ(ierr); /* update sensi2p for each stage */ 888 } 889 } 890 } 891 } 892 893 for (j=0; j<s; j++) w[j] = 1.0; 894 for (nadj=0; nadj<ts->numcost; nadj++) { /* no need to do this for mu's */ 895 ierr = VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecsDeltaLam[nadj*s]);CHKERRQ(ierr); 896 if (ts->vecs_sensi2) { 897 ierr = VecMAXPY(ts->vecs_sensi2[nadj],s,w,&VecsDeltaLam2[nadj*s]);CHKERRQ(ierr); 898 } 899 } 900 rk->status = TS_STEP_COMPLETE; 901 PetscFunctionReturn(0); 902 } 903 904 static PetscErrorCode TSAdjointReset_RK(TS ts) 905 { 906 TS_RK *rk = (TS_RK*)ts->data; 907 RKTableau tab = rk->tableau; 908 PetscErrorCode ierr; 909 910 PetscFunctionBegin; 911 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam);CHKERRQ(ierr); 912 ierr = VecDestroyVecs(ts->numcost,&rk->VecsSensiTemp);CHKERRQ(ierr); 913 ierr = VecDestroy(&rk->VecDeltaMu);CHKERRQ(ierr); 914 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam2);CHKERRQ(ierr); 915 ierr = VecDestroy(&rk->VecDeltaMu2);CHKERRQ(ierr); 916 ierr = VecDestroyVecs(ts->numcost,&rk->VecsSensi2Temp);CHKERRQ(ierr); 917 PetscFunctionReturn(0); 918 } 919 920 static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X) 921 { 922 TS_RK *rk = (TS_RK*)ts->data; 923 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 924 PetscReal h; 925 PetscReal tt,t; 926 PetscScalar *b; 927 const PetscReal *B = rk->tableau->binterp; 928 PetscErrorCode ierr; 929 930 PetscFunctionBegin; 931 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 932 933 switch (rk->status) { 934 case TS_STEP_INCOMPLETE: 935 case TS_STEP_PENDING: 936 h = ts->time_step; 937 t = (itime - ts->ptime)/h; 938 break; 939 case TS_STEP_COMPLETE: 940 h = ts->ptime - ts->ptime_prev; 941 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 942 break; 943 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 944 } 945 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 946 for (i=0; i<s; i++) b[i] = 0; 947 for (j=0,tt=t; j<p; j++,tt*=t) { 948 for (i=0; i<s; i++) { 949 b[i] += h * B[i*p+j] * tt; 950 } 951 } 952 ierr = VecCopy(rk->Y[0],X);CHKERRQ(ierr); 953 ierr = VecMAXPY(X,s,b,rk->YdotRHS);CHKERRQ(ierr); 954 ierr = PetscFree(b);CHKERRQ(ierr); 955 PetscFunctionReturn(0); 956 } 957 958 /*------------------------------------------------------------*/ 959 960 static PetscErrorCode TSRKTableauReset(TS ts) 961 { 962 TS_RK *rk = (TS_RK*)ts->data; 963 RKTableau tab = rk->tableau; 964 PetscErrorCode ierr; 965 966 PetscFunctionBegin; 967 if (!tab) PetscFunctionReturn(0); 968 ierr = PetscFree(rk->work);CHKERRQ(ierr); 969 ierr = VecDestroyVecs(tab->s,&rk->Y);CHKERRQ(ierr); 970 ierr = VecDestroyVecs(tab->s,&rk->YdotRHS);CHKERRQ(ierr); 971 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam);CHKERRQ(ierr); 972 ierr = VecDestroyVecs(ts->numcost,&rk->VecsSensiTemp);CHKERRQ(ierr); 973 ierr = VecDestroy(&rk->VecDeltaMu);CHKERRQ(ierr); 974 PetscFunctionReturn(0); 975 } 976 977 static PetscErrorCode TSReset_RK(TS ts) 978 { 979 PetscErrorCode ierr; 980 981 PetscFunctionBegin; 982 ierr = TSRKTableauReset(ts);CHKERRQ(ierr); 983 if (ts->use_splitrhsfunction) { 984 ierr = PetscTryMethod(ts,"TSReset_RK_MultirateSplit_C",(TS),(ts));CHKERRQ(ierr); 985 } else { 986 ierr = PetscTryMethod(ts,"TSReset_RK_MultirateNonsplit_C",(TS),(ts));CHKERRQ(ierr); 987 } 988 PetscFunctionReturn(0); 989 } 990 991 static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx) 992 { 993 PetscFunctionBegin; 994 PetscFunctionReturn(0); 995 } 996 997 static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 998 { 999 PetscFunctionBegin; 1000 PetscFunctionReturn(0); 1001 } 1002 1003 1004 static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx) 1005 { 1006 PetscFunctionBegin; 1007 PetscFunctionReturn(0); 1008 } 1009 1010 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1011 { 1012 1013 PetscFunctionBegin; 1014 PetscFunctionReturn(0); 1015 } 1016 /* 1017 static PetscErrorCode RKSetAdjCoe(RKTableau tab) 1018 { 1019 PetscReal *A,*b; 1020 PetscInt s,i,j; 1021 PetscErrorCode ierr; 1022 1023 PetscFunctionBegin; 1024 s = tab->s; 1025 ierr = PetscMalloc2(s*s,&A,s,&b);CHKERRQ(ierr); 1026 1027 for (i=0; i<s; i++) 1028 for (j=0; j<s; j++) { 1029 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]; 1030 ierr = PetscPrintf(PETSC_COMM_WORLD,"Coefficients: A[%D][%D]=%.6f\n",i,j,A[i*s+j]);CHKERRQ(ierr); 1031 } 1032 1033 for (i=0; i<s; i++) b[i] = (tab->b[s-1-i]==0)? 0: -tab->b[s-1-i]; 1034 1035 ierr = PetscMemcpy(tab->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 1036 ierr = PetscMemcpy(tab->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 1037 ierr = PetscFree2(A,b);CHKERRQ(ierr); 1038 PetscFunctionReturn(0); 1039 } 1040 */ 1041 1042 static PetscErrorCode TSRKTableauSetUp(TS ts) 1043 { 1044 TS_RK *rk = (TS_RK*)ts->data; 1045 RKTableau tab = rk->tableau; 1046 PetscErrorCode ierr; 1047 1048 PetscFunctionBegin; 1049 ierr = PetscMalloc1(tab->s,&rk->work);CHKERRQ(ierr); 1050 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->Y);CHKERRQ(ierr); 1051 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS);CHKERRQ(ierr); 1052 PetscFunctionReturn(0); 1053 } 1054 1055 static PetscErrorCode TSSetUp_RK(TS ts) 1056 { 1057 TS quadts = ts->quadraturets; 1058 PetscErrorCode ierr; 1059 DM dm; 1060 1061 PetscFunctionBegin; 1062 ierr = TSCheckImplicitTerm(ts);CHKERRQ(ierr); 1063 ierr = TSRKTableauSetUp(ts);CHKERRQ(ierr); 1064 if (quadts && ts->costintegralfwd) { 1065 Mat Jquad; 1066 ierr = TSGetRHSJacobian(quadts,&Jquad,NULL,NULL,NULL);CHKERRQ(ierr); 1067 } 1068 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1069 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 1070 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 1071 if (ts->use_splitrhsfunction) { 1072 ierr = PetscTryMethod(ts,"TSSetUp_RK_MultirateSplit_C",(TS),(ts));CHKERRQ(ierr); 1073 } else { 1074 ierr = PetscTryMethod(ts,"TSSetUp_RK_MultirateNonsplit_C",(TS),(ts));CHKERRQ(ierr); 1075 } 1076 PetscFunctionReturn(0); 1077 } 1078 1079 static PetscErrorCode TSSetFromOptions_RK(PetscOptionItems *PetscOptionsObject,TS ts) 1080 { 1081 TS_RK *rk = (TS_RK*)ts->data; 1082 PetscErrorCode ierr; 1083 1084 PetscFunctionBegin; 1085 ierr = PetscOptionsHead(PetscOptionsObject,"RK ODE solver options");CHKERRQ(ierr); 1086 { 1087 RKTableauLink link; 1088 PetscInt count,choice; 1089 PetscBool flg,use_multirate = PETSC_FALSE; 1090 const char **namelist; 1091 1092 for (link=RKTableauList,count=0; link; link=link->next,count++) ; 1093 ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 1094 for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1095 ierr = PetscOptionsBool("-ts_rk_multirate","Use interpolation-based multirate RK method","TSRKSetMultirate",rk->use_multirate,&use_multirate,&flg);CHKERRQ(ierr); 1096 if (flg) { 1097 ierr = TSRKSetMultirate(ts,use_multirate);CHKERRQ(ierr); 1098 } 1099 ierr = PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rk->tableau->name,&choice,&flg);CHKERRQ(ierr); 1100 if (flg) {ierr = TSRKSetType(ts,namelist[choice]);CHKERRQ(ierr);} 1101 ierr = PetscFree(namelist);CHKERRQ(ierr); 1102 } 1103 ierr = PetscOptionsTail();CHKERRQ(ierr); 1104 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Multirate methods options","");CHKERRQ(ierr); 1105 ierr = PetscOptionsInt("-ts_rk_dtratio","time step ratio between slow and fast","",rk->dtratio,&rk->dtratio,NULL);CHKERRQ(ierr); 1106 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1107 PetscFunctionReturn(0); 1108 } 1109 1110 static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer) 1111 { 1112 TS_RK *rk = (TS_RK*)ts->data; 1113 PetscBool iascii; 1114 PetscErrorCode ierr; 1115 1116 PetscFunctionBegin; 1117 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1118 if (iascii) { 1119 RKTableau tab = rk->tableau; 1120 TSRKType rktype; 1121 char buf[512]; 1122 ierr = TSRKGetType(ts,&rktype);CHKERRQ(ierr); 1123 ierr = PetscViewerASCIIPrintf(viewer," RK type %s\n",rktype);CHKERRQ(ierr); 1124 ierr = PetscViewerASCIIPrintf(viewer," Order: %D\n",tab->order);CHKERRQ(ierr); 1125 ierr = PetscViewerASCIIPrintf(viewer," FSAL property: %s\n",tab->FSAL ? "yes" : "no");CHKERRQ(ierr); 1126 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 1127 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 1128 } 1129 PetscFunctionReturn(0); 1130 } 1131 1132 static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer) 1133 { 1134 PetscErrorCode ierr; 1135 TSAdapt adapt; 1136 1137 PetscFunctionBegin; 1138 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1139 ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 1140 PetscFunctionReturn(0); 1141 } 1142 1143 /*@C 1144 TSRKSetType - Set the type of RK scheme 1145 1146 Logically collective 1147 1148 Input Parameter: 1149 + ts - timestepping context 1150 - rktype - type of RK-scheme 1151 1152 Options Database: 1153 . -ts_rk_type - <1fe,2a,3,3bs,4,5f,5dp,5bs> 1154 1155 Level: intermediate 1156 1157 .seealso: TSRKGetType(), TSRK, TSRKType, TSRK1FE, TSRK2A, TSRK3, TSRK3BS, TSRK4, TSRK5F, TSRK5DP, TSRK5BS 1158 @*/ 1159 PetscErrorCode TSRKSetType(TS ts,TSRKType rktype) 1160 { 1161 PetscErrorCode ierr; 1162 1163 PetscFunctionBegin; 1164 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1165 PetscValidCharPointer(rktype,2); 1166 ierr = PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));CHKERRQ(ierr); 1167 PetscFunctionReturn(0); 1168 } 1169 1170 /*@C 1171 TSRKGetType - Get the type of RK scheme 1172 1173 Logically collective 1174 1175 Input Parameter: 1176 . ts - timestepping context 1177 1178 Output Parameter: 1179 . rktype - type of RK-scheme 1180 1181 Level: intermediate 1182 1183 .seealso: TSRKGetType() 1184 @*/ 1185 PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype) 1186 { 1187 PetscErrorCode ierr; 1188 1189 PetscFunctionBegin; 1190 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1191 ierr = PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));CHKERRQ(ierr); 1192 PetscFunctionReturn(0); 1193 } 1194 1195 static PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype) 1196 { 1197 TS_RK *rk = (TS_RK*)ts->data; 1198 1199 PetscFunctionBegin; 1200 *rktype = rk->tableau->name; 1201 PetscFunctionReturn(0); 1202 } 1203 1204 static PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype) 1205 { 1206 TS_RK *rk = (TS_RK*)ts->data; 1207 PetscErrorCode ierr; 1208 PetscBool match; 1209 RKTableauLink link; 1210 1211 PetscFunctionBegin; 1212 if (rk->tableau) { 1213 ierr = PetscStrcmp(rk->tableau->name,rktype,&match);CHKERRQ(ierr); 1214 if (match) PetscFunctionReturn(0); 1215 } 1216 for (link = RKTableauList; link; link=link->next) { 1217 ierr = PetscStrcmp(link->tab.name,rktype,&match);CHKERRQ(ierr); 1218 if (match) { 1219 if (ts->setupcalled) {ierr = TSRKTableauReset(ts);CHKERRQ(ierr);} 1220 rk->tableau = &link->tab; 1221 if (ts->setupcalled) {ierr = TSRKTableauSetUp(ts);CHKERRQ(ierr);} 1222 ts->default_adapt_type = rk->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1223 PetscFunctionReturn(0); 1224 } 1225 } 1226 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype); 1227 PetscFunctionReturn(0); 1228 } 1229 1230 static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y) 1231 { 1232 TS_RK *rk = (TS_RK*)ts->data; 1233 1234 PetscFunctionBegin; 1235 if (ns) *ns = rk->tableau->s; 1236 if (Y) *Y = rk->Y; 1237 PetscFunctionReturn(0); 1238 } 1239 1240 static PetscErrorCode TSDestroy_RK(TS ts) 1241 { 1242 PetscErrorCode ierr; 1243 1244 PetscFunctionBegin; 1245 ierr = TSReset_RK(ts);CHKERRQ(ierr); 1246 if (ts->dm) { 1247 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 1248 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 1249 } 1250 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1251 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);CHKERRQ(ierr); 1252 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);CHKERRQ(ierr); 1253 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirate_C",NULL);CHKERRQ(ierr); 1254 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetMultirate_C",NULL);CHKERRQ(ierr); 1255 PetscFunctionReturn(0); 1256 } 1257 1258 /*@C 1259 TSRKSetMultirate - Use the interpolation-based multirate RK method 1260 1261 Logically collective 1262 1263 Input Parameter: 1264 + ts - timestepping context 1265 - 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 1266 1267 Options Database: 1268 . -ts_rk_multirate - <true,false> 1269 1270 Notes: 1271 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). 1272 1273 Level: intermediate 1274 1275 .seealso: TSRKGetMultirate() 1276 @*/ 1277 PetscErrorCode TSRKSetMultirate(TS ts,PetscBool use_multirate) 1278 { 1279 PetscErrorCode ierr; 1280 1281 PetscFunctionBegin; 1282 ierr = PetscTryMethod(ts,"TSRKSetMultirate_C",(TS,PetscBool),(ts,use_multirate));CHKERRQ(ierr); 1283 PetscFunctionReturn(0); 1284 } 1285 1286 /*@C 1287 TSRKGetMultirate - Gets whether to Use the interpolation-based multirate RK method 1288 1289 Not collective 1290 1291 Input Parameter: 1292 . ts - timestepping context 1293 1294 Output Parameter: 1295 . use_multirate - PETSC_TRUE if the multirate RK method is enabled, PETSC_FALSE otherwise 1296 1297 Level: intermediate 1298 1299 .seealso: TSRKSetMultirate() 1300 @*/ 1301 PetscErrorCode TSRKGetMultirate(TS ts,PetscBool *use_multirate) 1302 { 1303 PetscErrorCode ierr; 1304 1305 PetscFunctionBegin; 1306 ierr = PetscUseMethod(ts,"TSRKGetMultirate_C",(TS,PetscBool*),(ts,use_multirate));CHKERRQ(ierr); 1307 PetscFunctionReturn(0); 1308 } 1309 1310 /*MC 1311 TSRK - ODE and DAE solver using Runge-Kutta schemes 1312 1313 The user should provide the right hand side of the equation 1314 using TSSetRHSFunction(). 1315 1316 Notes: 1317 The default is TSRK3BS, it can be changed with TSRKSetType() or -ts_rk_type 1318 1319 Level: beginner 1320 1321 .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3, 1322 TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister(), TSRKSetMultirate(), TSRKGetMultirate() 1323 1324 M*/ 1325 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 1326 { 1327 TS_RK *rk; 1328 PetscErrorCode ierr; 1329 1330 PetscFunctionBegin; 1331 ierr = TSRKInitializePackage();CHKERRQ(ierr); 1332 1333 ts->ops->reset = TSReset_RK; 1334 ts->ops->destroy = TSDestroy_RK; 1335 ts->ops->view = TSView_RK; 1336 ts->ops->load = TSLoad_RK; 1337 ts->ops->setup = TSSetUp_RK; 1338 ts->ops->interpolate = TSInterpolate_RK; 1339 ts->ops->step = TSStep_RK; 1340 ts->ops->evaluatestep = TSEvaluateStep_RK; 1341 ts->ops->rollback = TSRollBack_RK; 1342 ts->ops->setfromoptions = TSSetFromOptions_RK; 1343 ts->ops->getstages = TSGetStages_RK; 1344 1345 ts->ops->adjointintegral = TSAdjointCostIntegral_RK; 1346 ts->ops->adjointsetup = TSAdjointSetUp_RK; 1347 ts->ops->adjointstep = TSAdjointStep_RK; 1348 ts->ops->adjointreset = TSAdjointReset_RK; 1349 1350 ts->ops->forwardintegral = TSForwardCostIntegral_RK; 1351 ts->ops->forwardsetup = TSForwardSetUp_RK; 1352 ts->ops->forwardreset = TSForwardReset_RK; 1353 ts->ops->forwardstep = TSForwardStep_RK; 1354 ts->ops->forwardgetstages= TSForwardGetStages_RK; 1355 1356 ierr = PetscNewLog(ts,&rk);CHKERRQ(ierr); 1357 ts->data = (void*)rk; 1358 1359 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);CHKERRQ(ierr); 1360 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);CHKERRQ(ierr); 1361 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirate_C",TSRKSetMultirate_RK);CHKERRQ(ierr); 1362 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetMultirate_C",TSRKGetMultirate_RK);CHKERRQ(ierr); 1363 1364 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 1365 rk->dtratio = 1; 1366 PetscFunctionReturn(0); 1367 } 1368