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