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 PetscInt p; /* Interpolation order */ 23 PetscBool FSAL; /* flag to indicate if tableau is FSAL */ 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,0,NULL);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,0,NULL);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,0,NULL);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,0,NULL);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,0,NULL);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,0,NULL);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,0,NULL);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 . p - Order of the interpolation scheme, equal to the number of columns of binterp 287 - binterp - Coefficients of the interpolation formula (dimension s*p; NULL to reuse b with p=1) 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[],PetscInt p,const PetscReal binterp[]) 301 { 302 PetscErrorCode ierr; 303 RKTableauLink link; 304 RKTableau t; 305 PetscInt i,j; 306 307 PetscFunctionBegin; 308 PetscValidCharPointer(name,1); 309 PetscValidRealPointer(A,4); 310 if (b) PetscValidRealPointer(b,5); 311 if (c) PetscValidRealPointer(c,6); 312 if (bembed) PetscValidRealPointer(bembed,7); 313 if (binterp || p > 1) PetscValidRealPointer(binterp,9); 314 315 ierr = PetscNew(&link);CHKERRQ(ierr); 316 t = &link->tab; 317 318 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 319 t->order = order; 320 t->s = s; 321 ierr = PetscMalloc3(s*s,&t->A,s,&t->b,s,&t->c);CHKERRQ(ierr); 322 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 323 if (b) { ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); } 324 else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i]; 325 if (c) { ierr = PetscMemcpy(t->c,c,s*sizeof(c[0]));CHKERRQ(ierr); } 326 else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 327 t->FSAL = PETSC_TRUE; 328 for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE; 329 330 if (bembed) { 331 ierr = PetscMalloc1(s,&t->bembed);CHKERRQ(ierr); 332 ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr); 333 } 334 335 if (!binterp) { p = 1; binterp = t->b; } 336 t->p = p; 337 ierr = PetscMalloc1(s*p,&t->binterp);CHKERRQ(ierr); 338 ierr = PetscMemcpy(t->binterp,binterp,s*p*sizeof(binterp[0]));CHKERRQ(ierr); 339 340 link->next = RKTableauList; 341 RKTableauList = link; 342 PetscFunctionReturn(0); 343 } 344 345 /* 346 The step completion formula is 347 348 x1 = x0 + h b^T YdotRHS 349 350 This function can be called before or after ts->vec_sol has been updated. 351 Suppose we have a completion formula (b) and an embedded formula (be) of different order. 352 We can write 353 354 x1e = x0 + h be^T YdotRHS 355 = x1 - h b^T YdotRHS + h be^T YdotRHS 356 = x1 + h (be - b)^T YdotRHS 357 358 so we can evaluate the method with different order even after the step has been optimistically completed. 359 */ 360 static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done) 361 { 362 TS_RK *rk = (TS_RK*)ts->data; 363 RKTableau tab = rk->tableau; 364 PetscScalar *w = rk->work; 365 PetscReal h; 366 PetscInt s = tab->s,j; 367 PetscErrorCode ierr; 368 369 PetscFunctionBegin; 370 switch (rk->status) { 371 case TS_STEP_INCOMPLETE: 372 case TS_STEP_PENDING: 373 h = ts->time_step; break; 374 case TS_STEP_COMPLETE: 375 h = ts->ptime - ts->ptime_prev; break; 376 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 377 } 378 if (order == tab->order) { 379 if (rk->status == TS_STEP_INCOMPLETE) { 380 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 381 for (j=0; j<s; j++) w[j] = h*tab->b[j]; 382 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 383 } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 384 PetscFunctionReturn(0); 385 } else if (order == tab->order-1) { 386 if (!tab->bembed) goto unavailable; 387 if (rk->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (be) */ 388 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 389 for (j=0; j<s; j++) w[j] = h*tab->bembed[j]; 390 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 391 } else { /* Rollback and re-complete using (be-b) */ 392 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 393 for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]); 394 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 395 if (ts->vec_costintegral && ts->costintegralfwd) { 396 ierr = VecCopy(rk->VecCostIntegral0,ts->vec_costintegral);CHKERRQ(ierr); 397 } 398 } 399 if (done) *done = PETSC_TRUE; 400 PetscFunctionReturn(0); 401 } 402 unavailable: 403 if (done) *done = PETSC_FALSE; 404 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); 405 PetscFunctionReturn(0); 406 } 407 408 static PetscErrorCode TSForwardCostIntegral_RK(TS ts) 409 { 410 TS_RK *rk = (TS_RK*)ts->data; 411 RKTableau tab = rk->tableau; 412 const PetscInt s = tab->s; 413 const PetscReal *b = tab->b,*c = tab->c; 414 Vec *Y = rk->Y; 415 PetscInt i; 416 PetscErrorCode ierr; 417 418 PetscFunctionBegin; 419 /* backup cost integral */ 420 ierr = VecCopy(ts->vec_costintegral,rk->VecCostIntegral0);CHKERRQ(ierr); 421 for (i=s-1; i>=0; i--) { 422 /* Evolve ts->vec_costintegral to compute integrals */ 423 ierr = TSAdjointComputeCostIntegrand(ts,rk->ptime+rk->time_step*(1.0-c[i]),Y[i],ts->vec_costintegrand);CHKERRQ(ierr); 424 ierr = VecAXPY(ts->vec_costintegral,rk->time_step*b[i],ts->vec_costintegrand);CHKERRQ(ierr); 425 } 426 PetscFunctionReturn(0); 427 } 428 429 static PetscErrorCode TSAdjointCostIntegral_RK(TS ts) 430 { 431 TS_RK *rk = (TS_RK*)ts->data; 432 RKTableau tab = rk->tableau; 433 const PetscInt s = tab->s; 434 const PetscReal *b = tab->b,*c = tab->c; 435 Vec *Y = rk->Y; 436 PetscInt i; 437 PetscErrorCode ierr; 438 439 PetscFunctionBegin; 440 for (i=s-1; i>=0; i--) { 441 /* Evolve ts->vec_costintegral to compute integrals */ 442 ierr = TSAdjointComputeCostIntegrand(ts,ts->ptime-ts->time_step*(1.0-c[i]),Y[i],ts->vec_costintegrand);CHKERRQ(ierr); 443 ierr = VecAXPY(ts->vec_costintegral,-ts->time_step*b[i],ts->vec_costintegrand);CHKERRQ(ierr); 444 } 445 PetscFunctionReturn(0); 446 } 447 448 static PetscErrorCode TSRollBack_RK(TS ts) 449 { 450 TS_RK *rk = (TS_RK*)ts->data; 451 RKTableau tab = rk->tableau; 452 const PetscInt s = tab->s; 453 const PetscReal *b = tab->b; 454 PetscScalar *w = rk->work; 455 Vec *YdotRHS = rk->YdotRHS; 456 PetscInt j; 457 PetscReal h; 458 PetscErrorCode ierr; 459 460 PetscFunctionBegin; 461 switch (rk->status) { 462 case TS_STEP_INCOMPLETE: 463 case TS_STEP_PENDING: 464 h = ts->time_step; break; 465 case TS_STEP_COMPLETE: 466 h = ts->ptime - ts->ptime_prev; break; 467 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 468 } 469 for (j=0; j<s; j++) w[j] = -h*b[j]; 470 ierr = VecMAXPY(ts->vec_sol,s,w,YdotRHS);CHKERRQ(ierr); 471 PetscFunctionReturn(0); 472 } 473 474 static PetscErrorCode TSStep_RK(TS ts) 475 { 476 TS_RK *rk = (TS_RK*)ts->data; 477 RKTableau tab = rk->tableau; 478 const PetscInt s = tab->s; 479 const PetscReal *A = tab->A,*c = tab->c; 480 PetscScalar *w = rk->work; 481 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 482 TSAdapt adapt; 483 PetscInt i,j; 484 PetscInt rejections = 0; 485 PetscBool stageok,accept = PETSC_TRUE; 486 PetscReal next_time_step = ts->time_step; 487 PetscErrorCode ierr; 488 489 PetscFunctionBegin; 490 rk->status = TS_STEP_INCOMPLETE; 491 while (!ts->reason && rk->status != TS_STEP_COMPLETE) { 492 PetscReal t = ts->ptime; 493 PetscReal h = ts->time_step; 494 for (i=0; i<s; i++) { 495 rk->stage_time = t + h*c[i]; 496 ierr = TSPreStage(ts,rk->stage_time); CHKERRQ(ierr); 497 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 498 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 499 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 500 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 501 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 502 ierr = TSAdaptCheckStage(adapt,ts,rk->stage_time,Y[i],&stageok);CHKERRQ(ierr); 503 if (!stageok) goto reject_step; 504 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 505 } 506 507 rk->status = TS_STEP_INCOMPLETE; 508 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 509 rk->status = TS_STEP_PENDING; 510 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 511 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 512 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 513 ierr = TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr); 514 rk->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 515 if (!accept) { /* Roll back the current step */ 516 ierr = TSRollBack_RK(ts);CHKERRQ(ierr); 517 ts->time_step = next_time_step; 518 goto reject_step; 519 } 520 521 if (ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation*/ 522 rk->ptime = ts->ptime; 523 rk->time_step = ts->time_step; 524 } 525 526 ts->ptime += ts->time_step; 527 ts->time_step = next_time_step; 528 break; 529 530 reject_step: 531 ts->reject++; accept = PETSC_FALSE; 532 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 533 ts->reason = TS_DIVERGED_STEP_REJECTED; 534 ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr); 535 } 536 } 537 PetscFunctionReturn(0); 538 } 539 540 static PetscErrorCode TSAdjointSetUp_RK(TS ts) 541 { 542 TS_RK *rk = (TS_RK*)ts->data; 543 RKTableau tab = rk->tableau; 544 PetscInt s = tab->s; 545 PetscErrorCode ierr; 546 547 PetscFunctionBegin; 548 if (ts->adjointsetupcalled++) PetscFunctionReturn(0); 549 ierr = VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecDeltaLam);CHKERRQ(ierr); 550 if(ts->vecs_sensip) { 551 ierr = VecDuplicateVecs(ts->vecs_sensip[0],s*ts->numcost,&rk->VecDeltaMu);CHKERRQ(ierr); 552 } 553 ierr = VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&rk->VecSensiTemp);CHKERRQ(ierr); 554 PetscFunctionReturn(0); 555 } 556 557 static PetscErrorCode TSAdjointStep_RK(TS ts) 558 { 559 TS_RK *rk = (TS_RK*)ts->data; 560 RKTableau tab = rk->tableau; 561 const PetscInt s = tab->s; 562 const PetscReal *A = tab->A,*b = tab->b,*c = tab->c; 563 PetscScalar *w = rk->work; 564 Vec *Y = rk->Y,*VecDeltaLam = rk->VecDeltaLam,*VecDeltaMu = rk->VecDeltaMu,*VecSensiTemp = rk->VecSensiTemp; 565 PetscInt i,j,nadj; 566 PetscReal t = ts->ptime; 567 PetscErrorCode ierr; 568 PetscReal h = ts->time_step; 569 Mat J,Jp; 570 571 PetscFunctionBegin; 572 rk->status = TS_STEP_INCOMPLETE; 573 for (i=s-1; i>=0; i--) { 574 rk->stage_time = t + h*(1.0-c[i]); 575 for (nadj=0; nadj<ts->numcost; nadj++) { 576 ierr = VecCopy(ts->vecs_sensi[nadj],VecSensiTemp[nadj]);CHKERRQ(ierr); 577 ierr = VecScale(VecSensiTemp[nadj],-h*b[i]);CHKERRQ(ierr); 578 for (j=i+1; j<s; j++) { 579 ierr = VecAXPY(VecSensiTemp[nadj],-h*A[j*s+i],VecDeltaLam[nadj*s+j]);CHKERRQ(ierr); 580 } 581 } 582 /* Stage values of lambda */ 583 ierr = TSGetRHSJacobian(ts,&J,&Jp,NULL,NULL);CHKERRQ(ierr); 584 ierr = TSComputeRHSJacobian(ts,rk->stage_time,Y[i],J,Jp);CHKERRQ(ierr); 585 if (ts->vec_costintegral) { 586 ierr = TSAdjointComputeDRDYFunction(ts,rk->stage_time,Y[i],ts->vecs_drdy);CHKERRQ(ierr); 587 } 588 for (nadj=0; nadj<ts->numcost; nadj++) { 589 ierr = MatMultTranspose(J,VecSensiTemp[nadj],VecDeltaLam[nadj*s+i]);CHKERRQ(ierr); 590 if (ts->vec_costintegral) { 591 ierr = VecAXPY(VecDeltaLam[nadj*s+i],-h*b[i],ts->vecs_drdy[nadj]);CHKERRQ(ierr); 592 } 593 } 594 595 /* Stage values of mu */ 596 if(ts->vecs_sensip) { 597 ierr = TSAdjointComputeRHSJacobian(ts,rk->stage_time,Y[i],ts->Jacp);CHKERRQ(ierr); 598 if (ts->vec_costintegral) { 599 ierr = TSAdjointComputeDRDPFunction(ts,rk->stage_time,Y[i],ts->vecs_drdp);CHKERRQ(ierr); 600 } 601 602 for (nadj=0; nadj<ts->numcost; nadj++) { 603 ierr = MatMultTranspose(ts->Jacp,VecSensiTemp[nadj],VecDeltaMu[nadj*s+i]);CHKERRQ(ierr); 604 if (ts->vec_costintegral) { 605 ierr = VecAXPY(VecDeltaMu[nadj*s+i],-h*b[i],ts->vecs_drdp[nadj]);CHKERRQ(ierr); 606 } 607 } 608 } 609 } 610 611 for (j=0; j<s; j++) w[j] = 1.0; 612 for (nadj=0; nadj<ts->numcost; nadj++) { 613 ierr = VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecDeltaLam[nadj*s]);CHKERRQ(ierr); 614 if(ts->vecs_sensip) { 615 ierr = VecMAXPY(ts->vecs_sensip[nadj],s,w,&VecDeltaMu[nadj*s]);CHKERRQ(ierr); 616 } 617 } 618 rk->status = TS_STEP_COMPLETE; 619 PetscFunctionReturn(0); 620 } 621 622 static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X) 623 { 624 TS_RK *rk = (TS_RK*)ts->data; 625 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 626 PetscReal h; 627 PetscReal tt,t; 628 PetscScalar *b; 629 const PetscReal *B = rk->tableau->binterp; 630 PetscErrorCode ierr; 631 632 PetscFunctionBegin; 633 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 634 635 switch (rk->status) { 636 case TS_STEP_INCOMPLETE: 637 case TS_STEP_PENDING: 638 h = ts->time_step; 639 t = (itime - ts->ptime)/h; 640 break; 641 case TS_STEP_COMPLETE: 642 h = ts->ptime - ts->ptime_prev; 643 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 644 break; 645 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 646 } 647 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 648 for (i=0; i<s; i++) b[i] = 0; 649 for (j=0,tt=t; j<p; j++,tt*=t) { 650 for (i=0; i<s; i++) { 651 b[i] += h * B[i*p+j] * tt; 652 } 653 } 654 655 ierr = VecCopy(rk->Y[0],X);CHKERRQ(ierr); 656 ierr = VecMAXPY(X,s,b,rk->YdotRHS);CHKERRQ(ierr); 657 658 ierr = PetscFree(b);CHKERRQ(ierr); 659 PetscFunctionReturn(0); 660 } 661 662 /*------------------------------------------------------------*/ 663 664 static PetscErrorCode TSRKTableauReset(TS ts) 665 { 666 TS_RK *rk = (TS_RK*)ts->data; 667 RKTableau tab = rk->tableau; 668 PetscErrorCode ierr; 669 670 PetscFunctionBegin; 671 if (!tab) PetscFunctionReturn(0); 672 ierr = PetscFree(rk->work);CHKERRQ(ierr); 673 ierr = VecDestroyVecs(tab->s,&rk->Y);CHKERRQ(ierr); 674 ierr = VecDestroyVecs(tab->s,&rk->YdotRHS);CHKERRQ(ierr); 675 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecDeltaLam);CHKERRQ(ierr); 676 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecDeltaMu);CHKERRQ(ierr); 677 PetscFunctionReturn(0); 678 } 679 680 static PetscErrorCode TSReset_RK(TS ts) 681 { 682 TS_RK *rk = (TS_RK*)ts->data; 683 PetscErrorCode ierr; 684 685 PetscFunctionBegin; 686 ierr = TSRKTableauReset(ts);CHKERRQ(ierr); 687 ierr = VecDestroy(&rk->VecCostIntegral0);CHKERRQ(ierr); 688 ierr = VecDestroyVecs(ts->numcost,&rk->VecSensiTemp);CHKERRQ(ierr); 689 PetscFunctionReturn(0); 690 } 691 692 static PetscErrorCode TSDestroy_RK(TS ts) 693 { 694 PetscErrorCode ierr; 695 696 PetscFunctionBegin; 697 ierr = TSReset_RK(ts);CHKERRQ(ierr); 698 ierr = PetscFree(ts->data);CHKERRQ(ierr); 699 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);CHKERRQ(ierr); 700 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);CHKERRQ(ierr); 701 PetscFunctionReturn(0); 702 } 703 704 705 static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx) 706 { 707 PetscFunctionBegin; 708 PetscFunctionReturn(0); 709 } 710 711 static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 712 { 713 PetscFunctionBegin; 714 PetscFunctionReturn(0); 715 } 716 717 718 static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx) 719 { 720 PetscFunctionBegin; 721 PetscFunctionReturn(0); 722 } 723 724 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 725 { 726 727 PetscFunctionBegin; 728 PetscFunctionReturn(0); 729 } 730 /* 731 static PetscErrorCode RKSetAdjCoe(RKTableau tab) 732 { 733 PetscReal *A,*b; 734 PetscInt s,i,j; 735 PetscErrorCode ierr; 736 737 PetscFunctionBegin; 738 s = tab->s; 739 ierr = PetscMalloc2(s*s,&A,s,&b);CHKERRQ(ierr); 740 741 for (i=0; i<s; i++) 742 for (j=0; j<s; j++) { 743 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]; 744 ierr = PetscPrintf(PETSC_COMM_WORLD,"Coefficients: A[%D][%D]=%.6f\n",i,j,A[i*s+j]);CHKERRQ(ierr); 745 } 746 747 for (i=0; i<s; i++) b[i] = (tab->b[s-1-i]==0)? 0: -tab->b[s-1-i]; 748 749 ierr = PetscMemcpy(tab->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 750 ierr = PetscMemcpy(tab->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 751 ierr = PetscFree2(A,b);CHKERRQ(ierr); 752 PetscFunctionReturn(0); 753 } 754 */ 755 756 static PetscErrorCode TSRKTableauSetUp(TS ts) 757 { 758 TS_RK *rk = (TS_RK*)ts->data; 759 RKTableau tab = rk->tableau; 760 PetscErrorCode ierr; 761 762 PetscFunctionBegin; 763 ierr = PetscMalloc1(tab->s,&rk->work);CHKERRQ(ierr); 764 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->Y);CHKERRQ(ierr); 765 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS);CHKERRQ(ierr); 766 PetscFunctionReturn(0); 767 } 768 769 770 static PetscErrorCode TSSetUp_RK(TS ts) 771 { 772 TS_RK *rk = (TS_RK*)ts->data; 773 PetscErrorCode ierr; 774 DM dm; 775 776 PetscFunctionBegin; 777 ierr = TSCheckImplicitTerm(ts);CHKERRQ(ierr); 778 ierr = TSRKTableauSetUp(ts);CHKERRQ(ierr); 779 if (!rk->VecCostIntegral0 && ts->vec_costintegral && ts->costintegralfwd) { /* back up cost integral */ 780 ierr = VecDuplicate(ts->vec_costintegral,&rk->VecCostIntegral0);CHKERRQ(ierr); 781 } 782 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 783 if (dm) { 784 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 785 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 786 } 787 PetscFunctionReturn(0); 788 } 789 790 791 /*------------------------------------------------------------*/ 792 793 static PetscErrorCode TSSetFromOptions_RK(PetscOptionItems *PetscOptionsObject,TS ts) 794 { 795 TS_RK *rk = (TS_RK*)ts->data; 796 PetscErrorCode ierr; 797 798 PetscFunctionBegin; 799 ierr = PetscOptionsHead(PetscOptionsObject,"RK ODE solver options");CHKERRQ(ierr); 800 { 801 RKTableauLink link; 802 PetscInt count,choice; 803 PetscBool flg; 804 const char **namelist; 805 for (link=RKTableauList,count=0; link; link=link->next,count++) ; 806 ierr = PetscMalloc1(count,&namelist);CHKERRQ(ierr); 807 for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 808 ierr = PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rk->tableau->name,&choice,&flg);CHKERRQ(ierr); 809 if (flg) {ierr = TSRKSetType(ts,namelist[choice]);CHKERRQ(ierr);} 810 ierr = PetscFree(namelist);CHKERRQ(ierr); 811 } 812 ierr = PetscOptionsTail();CHKERRQ(ierr); 813 PetscFunctionReturn(0); 814 } 815 816 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 817 { 818 PetscErrorCode ierr; 819 PetscInt i; 820 size_t left,count; 821 char *p; 822 823 PetscFunctionBegin; 824 for (i=0,p=buf,left=len; i<n; i++) { 825 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 826 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 827 left -= count; 828 p += count; 829 *p++ = ' '; 830 } 831 p[i ? 0 : -1] = 0; 832 PetscFunctionReturn(0); 833 } 834 835 static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer) 836 { 837 TS_RK *rk = (TS_RK*)ts->data; 838 PetscBool iascii; 839 PetscErrorCode ierr; 840 841 PetscFunctionBegin; 842 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 843 if (iascii) { 844 RKTableau tab = rk->tableau; 845 TSRKType rktype; 846 char buf[512]; 847 ierr = TSRKGetType(ts,&rktype);CHKERRQ(ierr); 848 ierr = PetscViewerASCIIPrintf(viewer," RK: %s\n",rktype);CHKERRQ(ierr); 849 ierr = PetscViewerASCIIPrintf(viewer," Order: %D\n",tab->order);CHKERRQ(ierr); 850 ierr = PetscViewerASCIIPrintf(viewer," FSAL property: %s\n",tab->FSAL ? "yes" : "no");CHKERRQ(ierr); 851 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 852 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 853 } 854 if (ts->adapt) {ierr = TSAdaptView(ts->adapt,viewer);CHKERRQ(ierr);} 855 PetscFunctionReturn(0); 856 } 857 858 static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer) 859 { 860 PetscErrorCode ierr; 861 TSAdapt adapt; 862 863 PetscFunctionBegin; 864 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 865 ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 866 PetscFunctionReturn(0); 867 } 868 869 /*@C 870 TSRKSetType - Set the type of RK scheme 871 872 Logically collective 873 874 Input Parameter: 875 + ts - timestepping context 876 - rktype - type of RK-scheme 877 878 Level: intermediate 879 880 .seealso: TSRKGetType(), TSRK, TSRK2, TSRK3, TSRKPRSSP2, TSRK4, TSRK5 881 @*/ 882 PetscErrorCode TSRKSetType(TS ts,TSRKType rktype) 883 { 884 PetscErrorCode ierr; 885 886 PetscFunctionBegin; 887 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 888 ierr = PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));CHKERRQ(ierr); 889 PetscFunctionReturn(0); 890 } 891 892 /*@C 893 TSRKGetType - Get the type of RK scheme 894 895 Logically collective 896 897 Input Parameter: 898 . ts - timestepping context 899 900 Output Parameter: 901 . rktype - type of RK-scheme 902 903 Level: intermediate 904 905 .seealso: TSRKGetType() 906 @*/ 907 PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype) 908 { 909 PetscErrorCode ierr; 910 911 PetscFunctionBegin; 912 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 913 ierr = PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));CHKERRQ(ierr); 914 PetscFunctionReturn(0); 915 } 916 917 static PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype) 918 { 919 TS_RK *rk = (TS_RK*)ts->data; 920 921 PetscFunctionBegin; 922 *rktype = rk->tableau->name; 923 PetscFunctionReturn(0); 924 } 925 static PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype) 926 { 927 TS_RK *rk = (TS_RK*)ts->data; 928 PetscErrorCode ierr; 929 PetscBool match; 930 RKTableauLink link; 931 932 PetscFunctionBegin; 933 if (rk->tableau) { 934 ierr = PetscStrcmp(rk->tableau->name,rktype,&match);CHKERRQ(ierr); 935 if (match) PetscFunctionReturn(0); 936 } 937 for (link = RKTableauList; link; link=link->next) { 938 ierr = PetscStrcmp(link->tab.name,rktype,&match);CHKERRQ(ierr); 939 if (match) { 940 if (ts->setupcalled) {ierr = TSRKTableauReset(ts);CHKERRQ(ierr);} 941 rk->tableau = &link->tab; 942 if (ts->setupcalled) {ierr = TSRKTableauSetUp(ts);CHKERRQ(ierr);} 943 PetscFunctionReturn(0); 944 } 945 } 946 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype); 947 PetscFunctionReturn(0); 948 } 949 950 static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y) 951 { 952 TS_RK *rk = (TS_RK*)ts->data; 953 954 PetscFunctionBegin; 955 *ns = rk->tableau->s; 956 if(Y) *Y = rk->Y; 957 PetscFunctionReturn(0); 958 } 959 960 961 /* ------------------------------------------------------------ */ 962 /*MC 963 TSRK - ODE and DAE solver using Runge-Kutta schemes 964 965 The user should provide the right hand side of the equation 966 using TSSetRHSFunction(). 967 968 Notes: 969 The default is TSRK3, it can be changed with TSRKSetType() or -ts_rk_type 970 971 Level: beginner 972 973 .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3, 974 TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister() 975 976 M*/ 977 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 978 { 979 TS_RK *rk; 980 PetscErrorCode ierr; 981 982 PetscFunctionBegin; 983 ierr = TSRKInitializePackage();CHKERRQ(ierr); 984 985 ts->ops->reset = TSReset_RK; 986 ts->ops->destroy = TSDestroy_RK; 987 ts->ops->view = TSView_RK; 988 ts->ops->load = TSLoad_RK; 989 ts->ops->setup = TSSetUp_RK; 990 ts->ops->adjointsetup = TSAdjointSetUp_RK; 991 ts->ops->step = TSStep_RK; 992 ts->ops->interpolate = TSInterpolate_RK; 993 ts->ops->evaluatestep = TSEvaluateStep_RK; 994 ts->ops->rollback = TSRollBack_RK; 995 ts->ops->setfromoptions = TSSetFromOptions_RK; 996 ts->ops->getstages = TSGetStages_RK; 997 ts->ops->adjointstep = TSAdjointStep_RK; 998 999 ts->ops->adjointintegral = TSAdjointCostIntegral_RK; 1000 ts->ops->forwardintegral = TSForwardCostIntegral_RK; 1001 1002 ierr = PetscNewLog(ts,&rk);CHKERRQ(ierr); 1003 ts->data = (void*)rk; 1004 1005 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);CHKERRQ(ierr); 1006 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);CHKERRQ(ierr); 1007 1008 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 1009 PetscFunctionReturn(0); 1010 } 1011