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