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 ierr = TSPreStep(ts);CHKERRQ(ierr); 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],&accept);CHKERRQ(ierr); 494 if (!accept) break; 495 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 496 } 497 if(!accept) continue; 498 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 499 rk->status = TS_STEP_PENDING; 500 501 /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 502 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 503 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 504 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 505 ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 506 if (accept) { 507 if (ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation*/ 508 rk->ptime = ts->ptime; 509 rk->time_step = ts->time_step; 510 } 511 /* ignore next_scheme for now */ 512 ts->ptime += ts->time_step; 513 ts->time_step = next_time_step; 514 ts->steps++; 515 rk->status = TS_STEP_COMPLETE; 516 ierr = PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);CHKERRQ(ierr); 517 break; 518 } else { /* Roll back the current step */ 519 for (j=0; j<s; j++) w[j] = -h*b[j]; 520 ierr = VecMAXPY(ts->vec_sol,s,w,rk->YdotRHS);CHKERRQ(ierr); 521 ts->time_step = next_time_step; 522 rk->status = TS_STEP_INCOMPLETE; 523 } 524 } 525 if (rk->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 526 PetscFunctionReturn(0); 527 } 528 529 #undef __FUNCT__ 530 #define __FUNCT__ "TSAdjointSetUp_RK" 531 static PetscErrorCode TSAdjointSetUp_RK(TS ts) 532 { 533 TS_RK *rk = (TS_RK*)ts->data; 534 RKTableau tab; 535 PetscInt s; 536 PetscErrorCode ierr; 537 538 PetscFunctionBegin; 539 if (ts->adjointsetupcalled++) PetscFunctionReturn(0); 540 tab = rk->tableau; 541 s = tab->s; 542 ierr = VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecDeltaLam);CHKERRQ(ierr); 543 if(ts->vecs_sensip) { 544 ierr = VecDuplicateVecs(ts->vecs_sensip[0],s*ts->numcost,&rk->VecDeltaMu);CHKERRQ(ierr); 545 } 546 ierr = VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&rk->VecSensiTemp);CHKERRQ(ierr); 547 PetscFunctionReturn(0); 548 } 549 550 #undef __FUNCT__ 551 #define __FUNCT__ "TSAdjointStep_RK" 552 static PetscErrorCode TSAdjointStep_RK(TS ts) 553 { 554 TS_RK *rk = (TS_RK*)ts->data; 555 RKTableau tab = rk->tableau; 556 const PetscInt s = tab->s; 557 const PetscReal *A = tab->A,*b = tab->b,*c = tab->c; 558 PetscScalar *w = rk->work; 559 Vec *Y = rk->Y,*VecDeltaLam = rk->VecDeltaLam,*VecDeltaMu = rk->VecDeltaMu,*VecSensiTemp = rk->VecSensiTemp; 560 PetscInt i,j,nadj; 561 PetscReal t; 562 PetscErrorCode ierr; 563 PetscReal h = ts->time_step; 564 Mat J,Jp; 565 566 PetscFunctionBegin; 567 t = ts->ptime; 568 rk->status = TS_STEP_INCOMPLETE; 569 h = ts->time_step; 570 ierr = TSPreStep(ts);CHKERRQ(ierr); 571 for (i=s-1; i>=0; i--) { 572 rk->stage_time = t + h*(1.0-c[i]); 573 for (nadj=0; nadj<ts->numcost; nadj++) { 574 ierr = VecCopy(ts->vecs_sensi[nadj],VecSensiTemp[nadj]);CHKERRQ(ierr); 575 ierr = VecScale(VecSensiTemp[nadj],-h*b[i]);CHKERRQ(ierr); 576 for (j=i+1; j<s; j++) { 577 ierr = VecAXPY(VecSensiTemp[nadj],-h*A[j*s+i],VecDeltaLam[nadj*s+j]);CHKERRQ(ierr); 578 } 579 } 580 /* Stage values of lambda */ 581 ierr = TSGetRHSJacobian(ts,&J,&Jp,NULL,NULL);CHKERRQ(ierr); 582 ierr = TSComputeRHSJacobian(ts,rk->stage_time,Y[i],J,Jp);CHKERRQ(ierr); 583 if (ts->vec_costintegral) { 584 ierr = TSAdjointComputeDRDYFunction(ts,rk->stage_time,Y[i],ts->vecs_drdy);CHKERRQ(ierr); 585 } 586 for (nadj=0; nadj<ts->numcost; nadj++) { 587 ierr = MatMultTranspose(J,VecSensiTemp[nadj],VecDeltaLam[nadj*s+i]);CHKERRQ(ierr); 588 if (ts->vec_costintegral) { 589 ierr = VecAXPY(VecDeltaLam[nadj*s+i],-h*b[i],ts->vecs_drdy[nadj]);CHKERRQ(ierr); 590 } 591 } 592 593 /* Stage values of mu */ 594 if(ts->vecs_sensip) { 595 ierr = TSAdjointComputeRHSJacobian(ts,rk->stage_time,Y[i],ts->Jacp);CHKERRQ(ierr); 596 if (ts->vec_costintegral) { 597 ierr = TSAdjointComputeDRDPFunction(ts,rk->stage_time,Y[i],ts->vecs_drdp);CHKERRQ(ierr); 598 } 599 600 for (nadj=0; nadj<ts->numcost; nadj++) { 601 ierr = MatMultTranspose(ts->Jacp,VecSensiTemp[nadj],VecDeltaMu[nadj*s+i]);CHKERRQ(ierr); 602 if (ts->vec_costintegral) { 603 ierr = VecAXPY(VecDeltaMu[nadj*s+i],-h*b[i],ts->vecs_drdp[nadj]);CHKERRQ(ierr); 604 } 605 } 606 } 607 } 608 609 for (j=0; j<s; j++) w[j] = 1.0; 610 for (nadj=0; nadj<ts->numcost; nadj++) { 611 ierr = VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecDeltaLam[nadj*s]);CHKERRQ(ierr); 612 if(ts->vecs_sensip) { 613 ierr = VecMAXPY(ts->vecs_sensip[nadj],s,w,&VecDeltaMu[nadj*s]);CHKERRQ(ierr); 614 } 615 } 616 ts->steps++; 617 rk->status = TS_STEP_COMPLETE; 618 PetscFunctionReturn(0); 619 } 620 621 #undef __FUNCT__ 622 #define __FUNCT__ "TSInterpolate_RK" 623 static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X) 624 { 625 TS_RK *rk = (TS_RK*)ts->data; 626 PetscInt s = rk->tableau->s,pinterp = rk->tableau->pinterp,i,j; 627 PetscReal h; 628 PetscReal tt,t; 629 PetscScalar *b; 630 const PetscReal *B = rk->tableau->binterp; 631 PetscErrorCode ierr; 632 633 PetscFunctionBegin; 634 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 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->time_step_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<pinterp; j++,tt*=t) { 650 for (i=0; i<s; i++) { 651 b[i] += h * B[i*pinterp+j] * tt; 652 } 653 } 654 ierr = VecCopy(rk->Y[0],X);CHKERRQ(ierr); 655 ierr = VecMAXPY(X,s,b,rk->YdotRHS);CHKERRQ(ierr); 656 ierr = PetscFree(b);CHKERRQ(ierr); 657 PetscFunctionReturn(0); 658 } 659 660 /*------------------------------------------------------------*/ 661 #undef __FUNCT__ 662 #define __FUNCT__ "TSReset_RK" 663 static PetscErrorCode TSReset_RK(TS ts) 664 { 665 TS_RK *rk = (TS_RK*)ts->data; 666 PetscInt s; 667 PetscErrorCode ierr; 668 669 PetscFunctionBegin; 670 if (!rk->tableau) PetscFunctionReturn(0); 671 s = rk->tableau->s; 672 ierr = VecDestroyVecs(s,&rk->Y);CHKERRQ(ierr); 673 ierr = VecDestroyVecs(s,&rk->YdotRHS);CHKERRQ(ierr); 674 if (rk->VecCostIntegral0) { 675 ierr = VecDestroy(&rk->VecCostIntegral0);CHKERRQ(ierr); 676 } 677 ierr = VecDestroyVecs(s*ts->numcost,&rk->VecDeltaLam);CHKERRQ(ierr); 678 ierr = VecDestroyVecs(s*ts->numcost,&rk->VecDeltaMu);CHKERRQ(ierr); 679 ierr = VecDestroyVecs(ts->numcost,&rk->VecSensiTemp);CHKERRQ(ierr); 680 ierr = PetscFree(rk->work);CHKERRQ(ierr); 681 PetscFunctionReturn(0); 682 } 683 684 #undef __FUNCT__ 685 #define __FUNCT__ "TSDestroy_RK" 686 static PetscErrorCode TSDestroy_RK(TS ts) 687 { 688 PetscErrorCode ierr; 689 690 PetscFunctionBegin; 691 ierr = TSReset_RK(ts);CHKERRQ(ierr); 692 ierr = PetscFree(ts->data);CHKERRQ(ierr); 693 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);CHKERRQ(ierr); 694 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);CHKERRQ(ierr); 695 PetscFunctionReturn(0); 696 } 697 698 699 #undef __FUNCT__ 700 #define __FUNCT__ "DMCoarsenHook_TSRK" 701 static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx) 702 { 703 PetscFunctionBegin; 704 PetscFunctionReturn(0); 705 } 706 707 #undef __FUNCT__ 708 #define __FUNCT__ "DMRestrictHook_TSRK" 709 static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 710 { 711 PetscFunctionBegin; 712 PetscFunctionReturn(0); 713 } 714 715 716 #undef __FUNCT__ 717 #define __FUNCT__ "DMSubDomainHook_TSRK" 718 static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx) 719 { 720 PetscFunctionBegin; 721 PetscFunctionReturn(0); 722 } 723 724 #undef __FUNCT__ 725 #define __FUNCT__ "DMSubDomainRestrictHook_TSRK" 726 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 727 { 728 729 PetscFunctionBegin; 730 PetscFunctionReturn(0); 731 } 732 /* 733 #undef __FUNCT__ 734 #define __FUNCT__ "RKSetAdjCoe" 735 static PetscErrorCode RKSetAdjCoe(RKTableau tab) 736 { 737 PetscReal *A,*b; 738 PetscInt s,i,j; 739 PetscErrorCode ierr; 740 741 PetscFunctionBegin; 742 s = tab->s; 743 ierr = PetscMalloc2(s*s,&A,s,&b);CHKERRQ(ierr); 744 745 for (i=0; i<s; i++) 746 for (j=0; j<s; j++) { 747 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]; 748 ierr = PetscPrintf(PETSC_COMM_WORLD,"Coefficients: A[%D][%D]=%.6f\n",i,j,A[i*s+j]);CHKERRQ(ierr); 749 } 750 751 for (i=0; i<s; i++) b[i] = (tab->b[s-1-i]==0)? 0: -tab->b[s-1-i]; 752 753 ierr = PetscMemcpy(tab->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 754 ierr = PetscMemcpy(tab->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 755 ierr = PetscFree2(A,b);CHKERRQ(ierr); 756 PetscFunctionReturn(0); 757 } 758 */ 759 #undef __FUNCT__ 760 #define __FUNCT__ "TSSetUp_RK" 761 static PetscErrorCode TSSetUp_RK(TS ts) 762 { 763 TS_RK *rk = (TS_RK*)ts->data; 764 RKTableau tab; 765 PetscInt s; 766 PetscErrorCode ierr; 767 DM dm; 768 769 PetscFunctionBegin; 770 if (!rk->VecCostIntegral0 && ts->vec_costintegral && ts->costintegralfwd) { /* back up cost integral */ 771 ierr = VecDuplicate(ts->vec_costintegral,&rk->VecCostIntegral0);CHKERRQ(ierr); 772 } 773 if (!rk->tableau) { 774 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 775 } 776 tab = rk->tableau; 777 s = tab->s; 778 ierr = VecDuplicateVecs(ts->vec_sol,s,&rk->Y);CHKERRQ(ierr); 779 ierr = VecDuplicateVecs(ts->vec_sol,s,&rk->YdotRHS);CHKERRQ(ierr); 780 ierr = PetscMalloc1(s,&rk->work);CHKERRQ(ierr); 781 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 782 if (dm) { 783 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 784 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 785 } 786 PetscFunctionReturn(0); 787 } 788 789 790 /*------------------------------------------------------------*/ 791 792 #undef __FUNCT__ 793 #define __FUNCT__ "TSSetFromOptions_RK" 794 static PetscErrorCode TSSetFromOptions_RK(PetscOptionItems *PetscOptionsObject,TS ts) 795 { 796 PetscErrorCode ierr; 797 char rktype[256]; 798 799 PetscFunctionBegin; 800 ierr = PetscOptionsHead(PetscOptionsObject,"RK ODE solver options");CHKERRQ(ierr); 801 { 802 RKTableauLink link; 803 PetscInt count,choice; 804 PetscBool flg; 805 const char **namelist; 806 ierr = PetscStrncpy(rktype,TSRKDefault,sizeof(rktype));CHKERRQ(ierr); 807 for (link=RKTableauList,count=0; link; link=link->next,count++) ; 808 ierr = PetscMalloc1(count,&namelist);CHKERRQ(ierr); 809 for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 810 ierr = PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rktype,&choice,&flg);CHKERRQ(ierr); 811 ierr = TSRKSetType(ts,flg ? namelist[choice] : rktype);CHKERRQ(ierr); 812 ierr = PetscFree(namelist);CHKERRQ(ierr); 813 } 814 ierr = PetscOptionsTail();CHKERRQ(ierr); 815 PetscFunctionReturn(0); 816 } 817 818 #undef __FUNCT__ 819 #define __FUNCT__ "PetscFormatRealArray" 820 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 821 { 822 PetscErrorCode ierr; 823 PetscInt i; 824 size_t left,count; 825 char *p; 826 827 PetscFunctionBegin; 828 for (i=0,p=buf,left=len; i<n; i++) { 829 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 830 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 831 left -= count; 832 p += count; 833 *p++ = ' '; 834 } 835 p[i ? 0 : -1] = 0; 836 PetscFunctionReturn(0); 837 } 838 839 #undef __FUNCT__ 840 #define __FUNCT__ "TSView_RK" 841 static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer) 842 { 843 TS_RK *rk = (TS_RK*)ts->data; 844 RKTableau tab = rk->tableau; 845 PetscBool iascii; 846 PetscErrorCode ierr; 847 TSAdapt adapt; 848 849 PetscFunctionBegin; 850 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 851 if (iascii) { 852 TSRKType rktype; 853 char buf[512]; 854 ierr = TSRKGetType(ts,&rktype);CHKERRQ(ierr); 855 ierr = PetscViewerASCIIPrintf(viewer," RK %s\n",rktype);CHKERRQ(ierr); 856 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 857 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 858 ierr = PetscViewerASCIIPrintf(viewer,"FSAL: %s\n",tab->FSAL ? "yes" : "no");CHKERRQ(ierr); 859 } 860 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 861 ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 862 PetscFunctionReturn(0); 863 } 864 865 #undef __FUNCT__ 866 #define __FUNCT__ "TSLoad_RK" 867 static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer) 868 { 869 PetscErrorCode ierr; 870 TSAdapt tsadapt; 871 872 PetscFunctionBegin; 873 ierr = TSGetAdapt(ts,&tsadapt);CHKERRQ(ierr); 874 ierr = TSAdaptLoad(tsadapt,viewer);CHKERRQ(ierr); 875 PetscFunctionReturn(0); 876 } 877 878 #undef __FUNCT__ 879 #define __FUNCT__ "TSRKSetType" 880 /*@C 881 TSRKSetType - Set the type of RK scheme 882 883 Logically collective 884 885 Input Parameter: 886 + ts - timestepping context 887 - rktype - type of RK-scheme 888 889 Level: intermediate 890 891 .seealso: TSRKGetType(), TSRK, TSRK2, TSRK3, TSRKPRSSP2, TSRK4, TSRK5 892 @*/ 893 PetscErrorCode TSRKSetType(TS ts,TSRKType rktype) 894 { 895 PetscErrorCode ierr; 896 897 PetscFunctionBegin; 898 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 899 ierr = PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));CHKERRQ(ierr); 900 PetscFunctionReturn(0); 901 } 902 903 #undef __FUNCT__ 904 #define __FUNCT__ "TSRKGetType" 905 /*@C 906 TSRKGetType - Get the type of RK scheme 907 908 Logically collective 909 910 Input Parameter: 911 . ts - timestepping context 912 913 Output Parameter: 914 . rktype - type of RK-scheme 915 916 Level: intermediate 917 918 .seealso: TSRKGetType() 919 @*/ 920 PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype) 921 { 922 PetscErrorCode ierr; 923 924 PetscFunctionBegin; 925 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 926 ierr = PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));CHKERRQ(ierr); 927 PetscFunctionReturn(0); 928 } 929 930 #undef __FUNCT__ 931 #define __FUNCT__ "TSRKGetType_RK" 932 PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype) 933 { 934 TS_RK *rk = (TS_RK*)ts->data; 935 PetscErrorCode ierr; 936 937 PetscFunctionBegin; 938 if (!rk->tableau) { 939 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 940 } 941 *rktype = rk->tableau->name; 942 PetscFunctionReturn(0); 943 } 944 #undef __FUNCT__ 945 #define __FUNCT__ "TSRKSetType_RK" 946 PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype) 947 { 948 TS_RK *rk = (TS_RK*)ts->data; 949 PetscErrorCode ierr; 950 PetscBool match; 951 RKTableauLink link; 952 953 PetscFunctionBegin; 954 if (rk->tableau) { 955 ierr = PetscStrcmp(rk->tableau->name,rktype,&match);CHKERRQ(ierr); 956 if (match) PetscFunctionReturn(0); 957 } 958 for (link = RKTableauList; link; link=link->next) { 959 ierr = PetscStrcmp(link->tab.name,rktype,&match);CHKERRQ(ierr); 960 if (match) { 961 ierr = TSReset_RK(ts);CHKERRQ(ierr); 962 rk->tableau = &link->tab; 963 PetscFunctionReturn(0); 964 } 965 } 966 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype); 967 PetscFunctionReturn(0); 968 } 969 970 #undef __FUNCT__ 971 #define __FUNCT__ "TSGetStages_RK" 972 static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y) 973 { 974 TS_RK *rk = (TS_RK*)ts->data; 975 976 PetscFunctionBegin; 977 *ns = rk->tableau->s; 978 if(Y) *Y = rk->Y; 979 PetscFunctionReturn(0); 980 } 981 982 983 /* ------------------------------------------------------------ */ 984 /*MC 985 TSRK - ODE and DAE solver using Runge-Kutta schemes 986 987 The user should provide the right hand side of the equation 988 using TSSetRHSFunction(). 989 990 Notes: 991 The default is TSRK3, it can be changed with TSRKSetType() or -ts_rk_type 992 993 Level: beginner 994 995 .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3, 996 TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister() 997 998 M*/ 999 #undef __FUNCT__ 1000 #define __FUNCT__ "TSCreate_RK" 1001 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 1002 { 1003 TS_RK *th; 1004 PetscErrorCode ierr; 1005 1006 PetscFunctionBegin; 1007 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1008 ierr = TSRKInitializePackage();CHKERRQ(ierr); 1009 #endif 1010 1011 ts->ops->reset = TSReset_RK; 1012 ts->ops->destroy = TSDestroy_RK; 1013 ts->ops->view = TSView_RK; 1014 ts->ops->load = TSLoad_RK; 1015 ts->ops->setup = TSSetUp_RK; 1016 ts->ops->adjointsetup = TSAdjointSetUp_RK; 1017 ts->ops->step = TSStep_RK; 1018 ts->ops->interpolate = TSInterpolate_RK; 1019 ts->ops->evaluatestep = TSEvaluateStep_RK; 1020 ts->ops->setfromoptions = TSSetFromOptions_RK; 1021 ts->ops->getstages = TSGetStages_RK; 1022 ts->ops->adjointstep = TSAdjointStep_RK; 1023 1024 ts->ops->adjointintegral = TSAdjointCostIntegral_RK; 1025 ts->ops->forwardintegral = TSForwardCostIntegral_RK; 1026 1027 ierr = PetscNewLog(ts,&th);CHKERRQ(ierr); 1028 ts->data = (void*)th; 1029 1030 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);CHKERRQ(ierr); 1031 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);CHKERRQ(ierr); 1032 PetscFunctionReturn(0); 1033 } 1034