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 481 rk->status = TS_STEP_INCOMPLETE; 482 while (!ts->reason && rk->status != TS_STEP_COMPLETE) { 483 PetscReal t = ts->ptime; 484 PetscReal h = ts->time_step; 485 for (i=0; i<s; i++) { 486 rk->stage_time = t + h*c[i]; 487 ierr = TSPreStage(ts,rk->stage_time); CHKERRQ(ierr); 488 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 489 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 490 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 491 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 492 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 493 ierr = TSAdaptCheckStage(adapt,ts,rk->stage_time,Y[i],&stageok);CHKERRQ(ierr); 494 if (!stageok) goto reject_step; 495 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 496 } 497 498 rk->status = TS_STEP_INCOMPLETE; 499 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 500 rk->status = TS_STEP_PENDING; 501 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 502 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 503 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 504 ierr = TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr); 505 rk->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 506 if (!accept) { /* Roll back the current step */ 507 ierr = TSRollBack_RK(ts);CHKERRQ(ierr); 508 ts->time_step = next_time_step; 509 goto reject_step; 510 } 511 512 if (ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation*/ 513 rk->ptime = ts->ptime; 514 rk->time_step = ts->time_step; 515 } 516 517 ts->ptime += ts->time_step; 518 ts->time_step = next_time_step; 519 break; 520 521 reject_step: 522 ts->reject++; accept = PETSC_FALSE; 523 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 524 ts->reason = TS_DIVERGED_STEP_REJECTED; 525 ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr); 526 } 527 } 528 PetscFunctionReturn(0); 529 } 530 531 static PetscErrorCode TSAdjointSetUp_RK(TS ts) 532 { 533 TS_RK *rk = (TS_RK*)ts->data; 534 RKTableau tab = rk->tableau; 535 PetscInt s = tab->s; 536 PetscErrorCode ierr; 537 538 PetscFunctionBegin; 539 if (ts->adjointsetupcalled++) PetscFunctionReturn(0); 540 ierr = VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecDeltaLam);CHKERRQ(ierr); 541 if(ts->vecs_sensip) { 542 ierr = VecDuplicateVecs(ts->vecs_sensip[0],s*ts->numcost,&rk->VecDeltaMu);CHKERRQ(ierr); 543 } 544 ierr = VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&rk->VecSensiTemp);CHKERRQ(ierr); 545 PetscFunctionReturn(0); 546 } 547 548 static PetscErrorCode TSAdjointStep_RK(TS ts) 549 { 550 TS_RK *rk = (TS_RK*)ts->data; 551 RKTableau tab = rk->tableau; 552 const PetscInt s = tab->s; 553 const PetscReal *A = tab->A,*b = tab->b,*c = tab->c; 554 PetscScalar *w = rk->work; 555 Vec *Y = rk->Y,*VecDeltaLam = rk->VecDeltaLam,*VecDeltaMu = rk->VecDeltaMu,*VecSensiTemp = rk->VecSensiTemp; 556 PetscInt i,j,nadj; 557 PetscReal t = ts->ptime; 558 PetscErrorCode ierr; 559 PetscReal h = ts->time_step; 560 Mat J,Jp; 561 562 PetscFunctionBegin; 563 rk->status = TS_STEP_INCOMPLETE; 564 for (i=s-1; i>=0; i--) { 565 rk->stage_time = t + h*(1.0-c[i]); 566 for (nadj=0; nadj<ts->numcost; nadj++) { 567 ierr = VecCopy(ts->vecs_sensi[nadj],VecSensiTemp[nadj]);CHKERRQ(ierr); 568 ierr = VecScale(VecSensiTemp[nadj],-h*b[i]);CHKERRQ(ierr); 569 for (j=i+1; j<s; j++) { 570 ierr = VecAXPY(VecSensiTemp[nadj],-h*A[j*s+i],VecDeltaLam[nadj*s+j]);CHKERRQ(ierr); 571 } 572 } 573 /* Stage values of lambda */ 574 ierr = TSGetRHSJacobian(ts,&J,&Jp,NULL,NULL);CHKERRQ(ierr); 575 ierr = TSComputeRHSJacobian(ts,rk->stage_time,Y[i],J,Jp);CHKERRQ(ierr); 576 if (ts->vec_costintegral) { 577 ierr = TSAdjointComputeDRDYFunction(ts,rk->stage_time,Y[i],ts->vecs_drdy);CHKERRQ(ierr); 578 } 579 for (nadj=0; nadj<ts->numcost; nadj++) { 580 ierr = MatMultTranspose(J,VecSensiTemp[nadj],VecDeltaLam[nadj*s+i]);CHKERRQ(ierr); 581 if (ts->vec_costintegral) { 582 ierr = VecAXPY(VecDeltaLam[nadj*s+i],-h*b[i],ts->vecs_drdy[nadj]);CHKERRQ(ierr); 583 } 584 } 585 586 /* Stage values of mu */ 587 if(ts->vecs_sensip) { 588 ierr = TSAdjointComputeRHSJacobian(ts,rk->stage_time,Y[i],ts->Jacp);CHKERRQ(ierr); 589 if (ts->vec_costintegral) { 590 ierr = TSAdjointComputeDRDPFunction(ts,rk->stage_time,Y[i],ts->vecs_drdp);CHKERRQ(ierr); 591 } 592 593 for (nadj=0; nadj<ts->numcost; nadj++) { 594 ierr = MatMultTranspose(ts->Jacp,VecSensiTemp[nadj],VecDeltaMu[nadj*s+i]);CHKERRQ(ierr); 595 if (ts->vec_costintegral) { 596 ierr = VecAXPY(VecDeltaMu[nadj*s+i],-h*b[i],ts->vecs_drdp[nadj]);CHKERRQ(ierr); 597 } 598 } 599 } 600 } 601 602 for (j=0; j<s; j++) w[j] = 1.0; 603 for (nadj=0; nadj<ts->numcost; nadj++) { 604 ierr = VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecDeltaLam[nadj*s]);CHKERRQ(ierr); 605 if(ts->vecs_sensip) { 606 ierr = VecMAXPY(ts->vecs_sensip[nadj],s,w,&VecDeltaMu[nadj*s]);CHKERRQ(ierr); 607 } 608 } 609 rk->status = TS_STEP_COMPLETE; 610 PetscFunctionReturn(0); 611 } 612 613 static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X) 614 { 615 TS_RK *rk = (TS_RK*)ts->data; 616 PetscInt s = rk->tableau->s,pinterp = rk->tableau->pinterp,i,j; 617 PetscReal h; 618 PetscReal tt,t; 619 PetscScalar *b; 620 const PetscReal *B = rk->tableau->binterp; 621 PetscErrorCode ierr; 622 623 PetscFunctionBegin; 624 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 625 626 switch (rk->status) { 627 case TS_STEP_INCOMPLETE: 628 case TS_STEP_PENDING: 629 h = ts->time_step; 630 t = (itime - ts->ptime)/h; 631 break; 632 case TS_STEP_COMPLETE: 633 h = ts->ptime - ts->ptime_prev; 634 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 635 break; 636 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 637 } 638 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 639 for (i=0; i<s; i++) b[i] = 0; 640 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 641 for (i=0; i<s; i++) { 642 b[i] += h * B[i*pinterp+j] * tt; 643 } 644 } 645 646 ierr = VecCopy(rk->Y[0],X);CHKERRQ(ierr); 647 ierr = VecMAXPY(X,s,b,rk->YdotRHS);CHKERRQ(ierr); 648 649 ierr = PetscFree(b);CHKERRQ(ierr); 650 PetscFunctionReturn(0); 651 } 652 653 /*------------------------------------------------------------*/ 654 655 static PetscErrorCode TSRKTableauReset(TS ts) 656 { 657 TS_RK *rk = (TS_RK*)ts->data; 658 RKTableau tab = rk->tableau; 659 PetscErrorCode ierr; 660 661 PetscFunctionBegin; 662 if (!tab) PetscFunctionReturn(0); 663 ierr = PetscFree(rk->work);CHKERRQ(ierr); 664 ierr = VecDestroyVecs(tab->s,&rk->Y);CHKERRQ(ierr); 665 ierr = VecDestroyVecs(tab->s,&rk->YdotRHS);CHKERRQ(ierr); 666 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecDeltaLam);CHKERRQ(ierr); 667 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecDeltaMu);CHKERRQ(ierr); 668 PetscFunctionReturn(0); 669 } 670 671 static PetscErrorCode TSReset_RK(TS ts) 672 { 673 TS_RK *rk = (TS_RK*)ts->data; 674 PetscErrorCode ierr; 675 676 PetscFunctionBegin; 677 ierr = TSRKTableauReset(ts);CHKERRQ(ierr); 678 ierr = VecDestroy(&rk->VecCostIntegral0);CHKERRQ(ierr); 679 ierr = VecDestroyVecs(ts->numcost,&rk->VecSensiTemp);CHKERRQ(ierr); 680 PetscFunctionReturn(0); 681 } 682 683 static PetscErrorCode TSDestroy_RK(TS ts) 684 { 685 PetscErrorCode ierr; 686 687 PetscFunctionBegin; 688 ierr = TSReset_RK(ts);CHKERRQ(ierr); 689 ierr = PetscFree(ts->data);CHKERRQ(ierr); 690 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);CHKERRQ(ierr); 691 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);CHKERRQ(ierr); 692 PetscFunctionReturn(0); 693 } 694 695 696 static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx) 697 { 698 PetscFunctionBegin; 699 PetscFunctionReturn(0); 700 } 701 702 static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 703 { 704 PetscFunctionBegin; 705 PetscFunctionReturn(0); 706 } 707 708 709 static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx) 710 { 711 PetscFunctionBegin; 712 PetscFunctionReturn(0); 713 } 714 715 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 716 { 717 718 PetscFunctionBegin; 719 PetscFunctionReturn(0); 720 } 721 /* 722 static PetscErrorCode RKSetAdjCoe(RKTableau tab) 723 { 724 PetscReal *A,*b; 725 PetscInt s,i,j; 726 PetscErrorCode ierr; 727 728 PetscFunctionBegin; 729 s = tab->s; 730 ierr = PetscMalloc2(s*s,&A,s,&b);CHKERRQ(ierr); 731 732 for (i=0; i<s; i++) 733 for (j=0; j<s; j++) { 734 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]; 735 ierr = PetscPrintf(PETSC_COMM_WORLD,"Coefficients: A[%D][%D]=%.6f\n",i,j,A[i*s+j]);CHKERRQ(ierr); 736 } 737 738 for (i=0; i<s; i++) b[i] = (tab->b[s-1-i]==0)? 0: -tab->b[s-1-i]; 739 740 ierr = PetscMemcpy(tab->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 741 ierr = PetscMemcpy(tab->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 742 ierr = PetscFree2(A,b);CHKERRQ(ierr); 743 PetscFunctionReturn(0); 744 } 745 */ 746 747 static PetscErrorCode TSRKTableauSetUp(TS ts) 748 { 749 TS_RK *rk = (TS_RK*)ts->data; 750 RKTableau tab = rk->tableau; 751 PetscErrorCode ierr; 752 753 PetscFunctionBegin; 754 ierr = PetscMalloc1(tab->s,&rk->work);CHKERRQ(ierr); 755 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->Y);CHKERRQ(ierr); 756 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS);CHKERRQ(ierr); 757 PetscFunctionReturn(0); 758 } 759 760 761 static PetscErrorCode TSSetUp_RK(TS ts) 762 { 763 TS_RK *rk = (TS_RK*)ts->data; 764 PetscErrorCode ierr; 765 DM dm; 766 767 PetscFunctionBegin; 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