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