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; 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 495 PetscReal h = ts->time_step; // already obtained from checkpointing 496 ierr = TSPreStep(ts);CHKERRQ(ierr); 497 for (i=s-1; i>=0; i--) { 498 Mat J,Jp; 499 rk->stage_time = t + h*(1.0-c[i]); 500 ierr = VecCopy(ts->vec_sensi,VecSensiTemp);CHKERRQ(ierr); 501 ierr = VecScale(VecSensiTemp,-h*b[i]); 502 for (j=i+1; j<s; j++) { 503 ierr = VecAXPY(VecSensiTemp,-h*A[j*s+i],VecDeltaLam[j]); 504 } 505 ierr = TSGetRHSJacobian(ts,&J,&Jp,NULL,NULL);CHKERRQ(ierr); 506 ierr = TSComputeRHSJacobian(ts,rk->stage_time,Y[i],J,Jp);CHKERRQ(ierr); 507 ierr = MatMultTranspose(J,VecSensiTemp,VecDeltaLam[i]);CHKERRQ(ierr); 508 } 509 510 for (j=0; j<s; j++) w[j] = 1.0; 511 ierr = VecMAXPY(ts->vec_sensi,s,w,VecDeltaLam);CHKERRQ(ierr); 512 513 ts->ptime += ts->time_step; 514 ts->steps++; 515 rk->status = TS_STEP_COMPLETE; 516 ierr = PetscObjectComposedDataSetReal((PetscObject)ts->vec_sensi,explicit_stage_time_id,ts->ptime);CHKERRQ(ierr); 517 PetscFunctionReturn(0); 518 } 519 520 #undef __FUNCT__ 521 #define __FUNCT__ "TSInterpolate_RK" 522 static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X) 523 { 524 TS_RK *rk = (TS_RK*)ts->data; 525 PetscInt s = rk->tableau->s,pinterp = rk->tableau->pinterp,i,j; 526 PetscReal h; 527 PetscReal tt,t; 528 PetscScalar *b; 529 const PetscReal *B = rk->tableau->binterp; 530 PetscErrorCode ierr; 531 532 PetscFunctionBegin; 533 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 534 switch (rk->status) { 535 case TS_STEP_INCOMPLETE: 536 case TS_STEP_PENDING: 537 h = ts->time_step; 538 t = (itime - ts->ptime)/h; 539 break; 540 case TS_STEP_COMPLETE: 541 h = ts->time_step_prev; 542 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 543 break; 544 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 545 } 546 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 547 for (i=0; i<s; i++) b[i] = 0; 548 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 549 for (i=0; i<s; i++) { 550 b[i] += h * B[i*pinterp+j] * tt; 551 } 552 } 553 ierr = VecCopy(rk->Y[0],X);CHKERRQ(ierr); 554 ierr = VecMAXPY(X,s,b,rk->YdotRHS);CHKERRQ(ierr); 555 ierr = PetscFree(b);CHKERRQ(ierr); 556 PetscFunctionReturn(0); 557 } 558 559 /*------------------------------------------------------------*/ 560 #undef __FUNCT__ 561 #define __FUNCT__ "TSReset_RK" 562 static PetscErrorCode TSReset_RK(TS ts) 563 { 564 TS_RK *rk = (TS_RK*)ts->data; 565 PetscInt s; 566 PetscErrorCode ierr; 567 568 PetscFunctionBegin; 569 if (!rk->tableau) PetscFunctionReturn(0); 570 s = rk->tableau->s; 571 ierr = VecDestroyVecs(s,&rk->Y);CHKERRQ(ierr); 572 ierr = VecDestroyVecs(s,&rk->YdotRHS);CHKERRQ(ierr); 573 if(ts->reverse_mode) { 574 ierr = VecDestroyVecs(s,&rk->VecDeltaLam);CHKERRQ(ierr); 575 ierr = VecDestroy(&rk->VecSensiTemp);CHKERRQ(ierr); 576 } 577 ierr = PetscFree(rk->work);CHKERRQ(ierr); 578 PetscFunctionReturn(0); 579 } 580 581 #undef __FUNCT__ 582 #define __FUNCT__ "TSDestroy_RK" 583 static PetscErrorCode TSDestroy_RK(TS ts) 584 { 585 PetscErrorCode ierr; 586 587 PetscFunctionBegin; 588 ierr = TSReset_RK(ts);CHKERRQ(ierr); 589 ierr = PetscFree(ts->data);CHKERRQ(ierr); 590 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);CHKERRQ(ierr); 591 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);CHKERRQ(ierr); 592 PetscFunctionReturn(0); 593 } 594 595 596 #undef __FUNCT__ 597 #define __FUNCT__ "DMCoarsenHook_TSRK" 598 static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx) 599 { 600 PetscFunctionBegin; 601 PetscFunctionReturn(0); 602 } 603 604 #undef __FUNCT__ 605 #define __FUNCT__ "DMRestrictHook_TSRK" 606 static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 607 { 608 PetscFunctionBegin; 609 PetscFunctionReturn(0); 610 } 611 612 613 #undef __FUNCT__ 614 #define __FUNCT__ "DMSubDomainHook_TSRK" 615 static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx) 616 { 617 PetscFunctionBegin; 618 PetscFunctionReturn(0); 619 } 620 621 #undef __FUNCT__ 622 #define __FUNCT__ "DMSubDomainRestrictHook_TSRK" 623 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 624 { 625 626 PetscFunctionBegin; 627 PetscFunctionReturn(0); 628 } 629 /* 630 #undef __FUNCT__ 631 #define __FUNCT__ "RKSetAdjCoe" 632 static PetscErrorCode RKSetAdjCoe(RKTableau tab) 633 { 634 PetscReal *A,*b; 635 PetscInt s,i,j; 636 PetscErrorCode ierr; 637 638 PetscFunctionBegin; 639 s = tab->s; 640 ierr = PetscMalloc2(s*s,&A,s,&b);CHKERRQ(ierr); 641 642 for (i=0; i<s; i++) 643 for (j=0; j<s; j++) { 644 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]; 645 ierr = PetscPrintf(PETSC_COMM_WORLD,"Coefficients: A[%D][%D]=%.6f\n",i,j,A[i*s+j]);CHKERRQ(ierr); 646 } 647 648 for (i=0; i<s; i++) b[i] = (tab->b[s-1-i]==0)? 0: -tab->b[s-1-i]; 649 650 ierr = PetscMemcpy(tab->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 651 ierr = PetscMemcpy(tab->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 652 ierr = PetscFree2(A,b);CHKERRQ(ierr); 653 PetscFunctionReturn(0); 654 } 655 */ 656 #undef __FUNCT__ 657 #define __FUNCT__ "TSSetUp_RK" 658 static PetscErrorCode TSSetUp_RK(TS ts) 659 { 660 TS_RK *rk = (TS_RK*)ts->data; 661 RKTableau tab; 662 PetscInt s; 663 PetscErrorCode ierr; 664 DM dm; 665 666 PetscFunctionBegin; 667 if (!rk->tableau) { 668 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 669 } 670 tab = rk->tableau; 671 s = tab->s; 672 // old 673 //if (ts->reverse_mode) { 674 // ierr = RKSetAdjCoe(tab);CHKERRQ(ierr); 675 //} 676 ierr = VecDuplicateVecs(ts->vec_sol,s,&rk->Y);CHKERRQ(ierr); 677 ierr = VecDuplicateVecs(ts->vec_sol,s,&rk->YdotRHS);CHKERRQ(ierr); 678 if (ts->reverse_mode) { 679 ierr = VecDuplicateVecs(ts->vec_sensi,s,&rk->VecDeltaLam);CHKERRQ(ierr); 680 ierr = VecDuplicate(ts->vec_sensi,&rk->VecSensiTemp);CHKERRQ(ierr); 681 } 682 ierr = PetscMalloc1(s,&rk->work);CHKERRQ(ierr); 683 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 684 if (dm) { 685 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 686 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 687 } 688 PetscFunctionReturn(0); 689 } 690 691 /*------------------------------------------------------------*/ 692 693 #undef __FUNCT__ 694 #define __FUNCT__ "TSSetFromOptions_RK" 695 static PetscErrorCode TSSetFromOptions_RK(PetscOptions *PetscOptionsObject,TS ts) 696 { 697 PetscErrorCode ierr; 698 char rktype[256]; 699 700 PetscFunctionBegin; 701 ierr = PetscOptionsHead(PetscOptionsObject,"RK ODE solver options");CHKERRQ(ierr); 702 { 703 RKTableauLink link; 704 PetscInt count,choice; 705 PetscBool flg; 706 const char **namelist; 707 ierr = PetscStrncpy(rktype,TSRKDefault,sizeof(rktype));CHKERRQ(ierr); 708 for (link=RKTableauList,count=0; link; link=link->next,count++) ; 709 ierr = PetscMalloc1(count,&namelist);CHKERRQ(ierr); 710 for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 711 ierr = PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rktype,&choice,&flg);CHKERRQ(ierr); 712 ierr = TSRKSetType(ts,flg ? namelist[choice] : rktype);CHKERRQ(ierr); 713 ierr = PetscFree(namelist);CHKERRQ(ierr); 714 } 715 ierr = PetscOptionsTail();CHKERRQ(ierr); 716 PetscFunctionReturn(0); 717 } 718 719 #undef __FUNCT__ 720 #define __FUNCT__ "PetscFormatRealArray" 721 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 722 { 723 PetscErrorCode ierr; 724 PetscInt i; 725 size_t left,count; 726 char *p; 727 728 PetscFunctionBegin; 729 for (i=0,p=buf,left=len; i<n; i++) { 730 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 731 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 732 left -= count; 733 p += count; 734 *p++ = ' '; 735 } 736 p[i ? 0 : -1] = 0; 737 PetscFunctionReturn(0); 738 } 739 740 #undef __FUNCT__ 741 #define __FUNCT__ "TSView_RK" 742 static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer) 743 { 744 TS_RK *rk = (TS_RK*)ts->data; 745 RKTableau tab = rk->tableau; 746 PetscBool iascii; 747 PetscErrorCode ierr; 748 TSAdapt adapt; 749 750 PetscFunctionBegin; 751 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 752 if (iascii) { 753 TSRKType rktype; 754 char buf[512]; 755 ierr = TSRKGetType(ts,&rktype);CHKERRQ(ierr); 756 ierr = PetscViewerASCIIPrintf(viewer," RK %s\n",rktype);CHKERRQ(ierr); 757 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 758 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 759 ierr = PetscViewerASCIIPrintf(viewer,"FSAL: %s\n",tab->FSAL ? "yes" : "no");CHKERRQ(ierr); 760 } 761 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 762 ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 763 PetscFunctionReturn(0); 764 } 765 766 #undef __FUNCT__ 767 #define __FUNCT__ "TSLoad_RK" 768 static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer) 769 { 770 PetscErrorCode ierr; 771 TSAdapt tsadapt; 772 773 PetscFunctionBegin; 774 ierr = TSGetAdapt(ts,&tsadapt);CHKERRQ(ierr); 775 ierr = TSAdaptLoad(tsadapt,viewer);CHKERRQ(ierr); 776 PetscFunctionReturn(0); 777 } 778 779 #undef __FUNCT__ 780 #define __FUNCT__ "TSRKSetType" 781 /*@C 782 TSRKSetType - Set the type of RK scheme 783 784 Logically collective 785 786 Input Parameter: 787 + ts - timestepping context 788 - rktype - type of RK-scheme 789 790 Level: intermediate 791 792 .seealso: TSRKGetType(), TSRK, TSRK2, TSRK3, TSRKPRSSP2, TSRK4, TSRK5 793 @*/ 794 PetscErrorCode TSRKSetType(TS ts,TSRKType rktype) 795 { 796 PetscErrorCode ierr; 797 798 PetscFunctionBegin; 799 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 800 ierr = PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));CHKERRQ(ierr); 801 PetscFunctionReturn(0); 802 } 803 804 #undef __FUNCT__ 805 #define __FUNCT__ "TSRKGetType" 806 /*@C 807 TSRKGetType - Get the type of RK scheme 808 809 Logically collective 810 811 Input Parameter: 812 . ts - timestepping context 813 814 Output Parameter: 815 . rktype - type of RK-scheme 816 817 Level: intermediate 818 819 .seealso: TSRKGetType() 820 @*/ 821 PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype) 822 { 823 PetscErrorCode ierr; 824 825 PetscFunctionBegin; 826 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 827 ierr = PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));CHKERRQ(ierr); 828 PetscFunctionReturn(0); 829 } 830 831 #undef __FUNCT__ 832 #define __FUNCT__ "TSRKGetType_RK" 833 PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype) 834 { 835 TS_RK *rk = (TS_RK*)ts->data; 836 PetscErrorCode ierr; 837 838 PetscFunctionBegin; 839 if (!rk->tableau) { 840 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 841 } 842 *rktype = rk->tableau->name; 843 PetscFunctionReturn(0); 844 } 845 #undef __FUNCT__ 846 #define __FUNCT__ "TSRKSetType_RK" 847 PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype) 848 { 849 TS_RK *rk = (TS_RK*)ts->data; 850 PetscErrorCode ierr; 851 PetscBool match; 852 RKTableauLink link; 853 854 PetscFunctionBegin; 855 if (rk->tableau) { 856 ierr = PetscStrcmp(rk->tableau->name,rktype,&match);CHKERRQ(ierr); 857 if (match) PetscFunctionReturn(0); 858 } 859 for (link = RKTableauList; link; link=link->next) { 860 ierr = PetscStrcmp(link->tab.name,rktype,&match);CHKERRQ(ierr); 861 if (match) { 862 ierr = TSReset_RK(ts);CHKERRQ(ierr); 863 rk->tableau = &link->tab; 864 PetscFunctionReturn(0); 865 } 866 } 867 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype); 868 PetscFunctionReturn(0); 869 } 870 871 #undef __FUNCT__ 872 #define __FUNCT__ "TSGetStages_RK" 873 static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y) 874 { 875 TS_RK *rk = (TS_RK*)ts->data; 876 877 PetscFunctionBegin; 878 *ns = rk->tableau->s; 879 if(Y) *Y = rk->Y; 880 PetscFunctionReturn(0); 881 } 882 883 884 /* ------------------------------------------------------------ */ 885 /*MC 886 TSRK - ODE and DAE solver using Runge-Kutta schemes 887 888 The user should provide the right hand side of the equation 889 using TSSetRHSFunction(). 890 891 Notes: 892 The default is TSRK3, it can be changed with TSRKSetType() or -ts_rk_type 893 894 Level: beginner 895 896 .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3, 897 TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister() 898 899 M*/ 900 #undef __FUNCT__ 901 #define __FUNCT__ "TSCreate_RK" 902 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 903 { 904 TS_RK *th; 905 PetscErrorCode ierr; 906 907 PetscFunctionBegin; 908 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 909 ierr = TSRKInitializePackage();CHKERRQ(ierr); 910 #endif 911 912 ts->ops->reset = TSReset_RK; 913 ts->ops->destroy = TSDestroy_RK; 914 ts->ops->view = TSView_RK; 915 ts->ops->load = TSLoad_RK; 916 ts->ops->setup = TSSetUp_RK; 917 ts->ops->step = TSStep_RK; 918 ts->ops->interpolate = TSInterpolate_RK; 919 ts->ops->evaluatestep = TSEvaluateStep_RK; 920 ts->ops->setfromoptions = TSSetFromOptions_RK; 921 ts->ops->getstages = TSGetStages_RK; 922 ts->ops->stepadj = TSStepAdj_RK; 923 924 ierr = PetscNewLog(ts,&th);CHKERRQ(ierr); 925 ts->data = (void*)th; 926 927 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);CHKERRQ(ierr); 928 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);CHKERRQ(ierr); 929 PetscFunctionReturn(0); 930 } 931