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