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