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