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