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